20671555286

Committed 01 Jan 2026 09:54PM UTC coverage: 63.137% (-2.8%) from 65.933%

Build # 20671555286

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #311 from libKriging/perf-improvements

Perf improvements

Coverage Stats

232 of 890 new or added lines in 10 files covered. (26.07%)

3 existing lines in 3 files now uncovered.

8610 of 13637 relevant lines covered (63.14%)

78049.07 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/bench/bench-linearalgebra.cpp

// clang-format off
// MUST BE at the beginning before any other <cmath> include (e.g. in armadillo's headers)
#define _USE_MATH_DEFINES
#include <cmath>
// clang-format on

#include "libKriging/LinearAlgebra.hpp"
#include "libKriging/Covariance.hpp"
#include "libKriging/utils/lk_armadillo.hpp"

#include <chrono>
#include <iostream>
#include <iomanip>
#include <vector>
#include <numeric>
#include <algorithm>
#include <functional>

#ifdef _OPENMP
#include <omp.h>
#endif

// Statistics computation
struct Stats {
  double mean;
  double std;
  double min;
  double max;
  double median;
};

Stats compute_stats(const std::vector<double>& values) {
  if (values.empty()) {
    return {0.0, 0.0, 0.0, 0.0, 0.0};
  }

  double sum = std::accumulate(values.begin(), values.end(), 0.0);
  double mean = sum / values.size();

  double sq_sum = 0.0;
  for (double v : values) {
    sq_sum += (v - mean) * (v - mean);
  }
  double std = std::sqrt(sq_sum / values.size());

  std::vector<double> sorted = values;
  std::sort(sorted.begin(), sorted.end());

  double min = sorted.front();
  double max = sorted.back();
  double median = sorted[sorted.size() / 2];

  return {mean, std, min, max, median};
}

void print_stats(const std::string& operation, const Stats& stats) {
  std::cout << std::setw(25) << std::left << operation << " | ";
  std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.mean << " | ";
  std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.std << " | ";
  std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.min << " | ";
  std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.max << " | ";
  std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.median << std::endl;
}

void print_header() {
  std::cout << std::setw(25) << std::left << "Operation" << " | ";
  std::cout << std::setw(10) << std::right << "Mean (ms)" << " | ";
  std::cout << std::setw(10) << std::right << "Std (ms)" << " | ";
  std::cout << std::setw(10) << std::right << "Min (ms)" << " | ";
  std::cout << std::setw(10) << std::right << "Max (ms)" << " | ";
  std::cout << std::setw(10) << std::right << "Median (ms)" << std::endl;
}

void benchmark_configuration(arma::uword n, arma::uword d, int n_iterations) {
  std::cout << "\n";
  std::cout << "n=" << n << " d=" << d << " iterations=" << n_iterations << std::endl;

  // Prepare test data similar to Kriging usage
  arma::mat X = arma::randu(d, n);  // Column-major: d × n (transposed from observations)
  arma::vec theta = arma::ones(d);

  // Gaussian covariance function (matching Kriging usage)
  auto gauss_cov = [](const arma::vec& diff, const arma::vec& theta) -> double {
    double sum_sq = 0.0;
    for (arma::uword i = 0; i < diff.n_elem; i++) {
      double val = diff[i] / theta[i];
      sum_sq += val * val;
    }
    return std::exp(-0.5 * sum_sq);
  };

  std::vector<double> times_covMat_sym_X;
  std::vector<double> times_covMat_rect;
  std::vector<double> times_solve;
  std::vector<double> times_rsolve;

  for (int iter = 0; iter < n_iterations; iter++) {
    // Benchmark covMat_sym_X (symmetric covariance matrix - used in fit())
    {
      arma::mat R(n, n);
      arma::vec diag_vec;  // empty = use factor

      auto t1 = std::chrono::high_resolution_clock::now();
      LinearAlgebra::covMat_sym_X(&R, X, theta, gauss_cov, 1.0, diag_vec);
      auto t2 = std::chrono::high_resolution_clock::now();

      double elapsed = std::chrono::duration<double, std::milli>(t2 - t1).count();
      times_covMat_sym_X.push_back(elapsed);
    }

    // Benchmark covMat_rect (rectangular covariance - used in predict())
    {
      arma::uword n_pred = std::min(n, arma::uword(100));  // Typical prediction size
      arma::mat X_pred = arma::randu(d, n_pred);
      arma::mat R_rect(n, n_pred);

      auto t1 = std::chrono::high_resolution_clock::now();
      LinearAlgebra::covMat_rect(&R_rect, X, X_pred, theta, gauss_cov, 1.0);
      auto t2 = std::chrono::high_resolution_clock::now();

      double elapsed = std::chrono::duration<double, std::milli>(t2 - t1).count();
      times_covMat_rect.push_back(elapsed);
    }

    // Benchmark solve (used extensively in Kriging)
    {
      arma::mat A = arma::randu(n, n);
      A = A * A.t() + arma::eye(n, n) * 0.1;  // Make positive definite
      arma::mat B = arma::randu(n, 10);

      auto t1 = std::chrono::high_resolution_clock::now();
      arma::mat X_sol = LinearAlgebra::solve(A, B);
      auto t2 = std::chrono::high_resolution_clock::now();

      double elapsed = std::chrono::duration<double, std::milli>(t2 - t1).count();
      times_solve.push_back(elapsed);
    }

    // Benchmark rsolve (right-solve, used in Kriging)
    {
      arma::mat A = arma::randu(n, n);
      A = A * A.t() + arma::eye(n, n) * 0.1;  // Make positive definite
      arma::mat B = arma::randu(10, n);

      auto t1 = std::chrono::high_resolution_clock::now();
      arma::mat X_sol = LinearAlgebra::rsolve(A, B);
      auto t2 = std::chrono::high_resolution_clock::now();

      double elapsed = std::chrono::duration<double, std::milli>(t2 - t1).count();
      times_rsolve.push_back(elapsed);
    }
  }

  // Print results
  print_header();
  print_stats("covMat_sym_X", compute_stats(times_covMat_sym_X));
  print_stats("covMat_rect", compute_stats(times_covMat_rect));
  print_stats("solve", compute_stats(times_solve));
  print_stats("rsolve", compute_stats(times_rsolve));
  std::cout << std::endl;
}

int main(int argc, char* argv[]) {
  int n_iterations = 10;
  if (argc > 1) {
    n_iterations = std::atoi(argv[1]);
  }

  std::cout << "=============================================================================" << std::endl;
  std::cout << "              LinearAlgebra Performance Benchmark (OpenMP)                   " << std::endl;
  std::cout << "=============================================================================" << std::endl;

#ifdef _OPENMP
  int max_threads = omp_get_max_threads();
  int optimal_threads = (max_threads > 2) ? 2 : max_threads;
  std::cout << "OpenMP enabled: max_threads=" << max_threads
            << ", using optimal=" << optimal_threads << std::endl;
#else
  std::cout << "OpenMP disabled: serial execution" << std::endl;
#endif

  std::cout << "Testing covariance matrix operations with Kriging-realistic matrices" << std::endl;
  std::cout << "=============================================================================" << std::endl;

  // Test various sizes relevant to Kriging
  benchmark_configuration(100, 2, n_iterations);
  benchmark_configuration(100, 4, n_iterations);
  benchmark_configuration(200, 2, n_iterations);
  benchmark_configuration(200, 4, n_iterations);
  benchmark_configuration(500, 2, n_iterations);
  benchmark_configuration(500, 4, n_iterations);
  benchmark_configuration(1000, 2, n_iterations);
  benchmark_configuration(1000, 4, n_iterations);

  std::cout << "=============================================================================" << std::endl;
  std::cout << "Benchmark complete" << std::endl;
  std::cout << "=============================================================================" << std::endl;

  return 0;
}

1	// clang-format off
2	// MUST BE at the beginning before any other <cmath> include (e.g. in armadillo's headers)
3	#define _USE_MATH_DEFINES
4	#include <cmath>
5	// clang-format on
6
7	#include "libKriging/LinearAlgebra.hpp"
8	#include "libKriging/Covariance.hpp"
9	#include "libKriging/utils/lk_armadillo.hpp"
10
11	#include <chrono>
12	#include <iostream>
13	#include <iomanip>
14	#include <vector>
15	#include <numeric>
16	#include <algorithm>
17	#include <functional>
18
19	#ifdef _OPENMP
20	#include <omp.h>
21	#endif
22
23	// Statistics computation
24	struct Stats {
25	double mean;
26	double std;
27	double min;
28	double max;
29	double median;
30	};
31
NEW 32	Stats compute_stats(const std::vector<double>& values) {	×
NEW 33	if (values.empty()) {	×
NEW 34	return {0.0, 0.0, 0.0, 0.0, 0.0};	×
35	}
36
NEW 37	double sum = std::accumulate(values.begin(), values.end(), 0.0);	×
NEW 38	double mean = sum / values.size();	×
39
NEW 40	double sq_sum = 0.0;	×
NEW 41	for (double v : values) {	×
NEW 42	sq_sum += (v - mean) * (v - mean);	×
43	}
NEW 44	double std = std::sqrt(sq_sum / values.size());	×
45
NEW 46	std::vector<double> sorted = values;	×
NEW 47	std::sort(sorted.begin(), sorted.end());	×
48
NEW 49	double min = sorted.front();	×
NEW 50	double max = sorted.back();	×
NEW 51	double median = sorted[sorted.size() / 2];	×
52
NEW 53	return {mean, std, min, max, median};	×
NEW 54	}	×
55
NEW 56	void print_stats(const std::string& operation, const Stats& stats) {	×
NEW 57	std::cout << std::setw(25) << std::left << operation << " \| ";	×
NEW 58	std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.mean << " \| ";	×
NEW 59	std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.std << " \| ";	×
NEW 60	std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.min << " \| ";	×
NEW 61	std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.max << " \| ";	×
NEW 62	std::cout << std::setw(10) << std::right << std::fixed << std::setprecision(3) << stats.median << std::endl;	×
NEW 63	}	×
64
NEW 65	void print_header() {	×
NEW 66	std::cout << std::setw(25) << std::left << "Operation" << " \| ";	×
NEW 67	std::cout << std::setw(10) << std::right << "Mean (ms)" << " \| ";	×
NEW 68	std::cout << std::setw(10) << std::right << "Std (ms)" << " \| ";	×
NEW 69	std::cout << std::setw(10) << std::right << "Min (ms)" << " \| ";	×
NEW 70	std::cout << std::setw(10) << std::right << "Max (ms)" << " \| ";	×
NEW 71	std::cout << std::setw(10) << std::right << "Median (ms)" << std::endl;	×
NEW 72	}	×
73
NEW 74	void benchmark_configuration(arma::uword n, arma::uword d, int n_iterations) {	×
NEW 75	std::cout << "\n";	×
NEW 76	std::cout << "n=" << n << " d=" << d << " iterations=" << n_iterations << std::endl;	×
77
78	// Prepare test data similar to Kriging usage
NEW 79	arma::mat X = arma::randu(d, n); // Column-major: d × n (transposed from observations)	×
NEW 80	arma::vec theta = arma::ones(d);	×
81
82	// Gaussian covariance function (matching Kriging usage)
NEW 83	auto gauss_cov = [](const arma::vec& diff, const arma::vec& theta) -> double {	×
NEW 84	double sum_sq = 0.0;	×
NEW 85	for (arma::uword i = 0; i < diff.n_elem; i++) {	×
NEW 86	double val = diff[i] / theta[i];	×
NEW 87	sum_sq += val * val;	×
88	}
NEW 89	return std::exp(-0.5 * sum_sq);	×
90	};
91
NEW 92	std::vector<double> times_covMat_sym_X;	×
NEW 93	std::vector<double> times_covMat_rect;	×
NEW 94	std::vector<double> times_solve;	×
NEW 95	std::vector<double> times_rsolve;	×
96
NEW 97	for (int iter = 0; iter < n_iterations; iter++) {	×
98	// Benchmark covMat_sym_X (symmetric covariance matrix - used in fit())
99	{
NEW 100	arma::mat R(n, n);	×
NEW 101	arma::vec diag_vec; // empty = use factor	×
102
NEW 103	auto t1 = std::chrono::high_resolution_clock::now();	×
NEW 104	LinearAlgebra::covMat_sym_X(&R, X, theta, gauss_cov, 1.0, diag_vec);	×
NEW 105	auto t2 = std::chrono::high_resolution_clock::now();	×
106
NEW 107	double elapsed = std::chrono::duration<double, std::milli>(t2 - t1).count();	×
NEW 108	times_covMat_sym_X.push_back(elapsed);	×
NEW 109	}	×
110
111	// Benchmark covMat_rect (rectangular covariance - used in predict())
112	{
NEW 113	arma::uword n_pred = std::min(n, arma::uword(100)); // Typical prediction size	×
NEW 114	arma::mat X_pred = arma::randu(d, n_pred);	×
NEW 115	arma::mat R_rect(n, n_pred);	×
116
NEW 117	auto t1 = std::chrono::high_resolution_clock::now();	×
NEW 118	LinearAlgebra::covMat_rect(&R_rect, X, X_pred, theta, gauss_cov, 1.0);	×
NEW 119	auto t2 = std::chrono::high_resolution_clock::now();	×
120
NEW 121	double elapsed = std::chrono::duration<double, std::milli>(t2 - t1).count();	×
NEW 122	times_covMat_rect.push_back(elapsed);	×
NEW 123	}	×
124
125	// Benchmark solve (used extensively in Kriging)
126	{
NEW 127	arma::mat A = arma::randu(n, n);	×
NEW 128	A = A * A.t() + arma::eye(n, n) * 0.1; // Make positive definite	×
NEW 129	arma::mat B = arma::randu(n, 10);	×
130
NEW 131	auto t1 = std::chrono::high_resolution_clock::now();	×
NEW 132	arma::mat X_sol = LinearAlgebra::solve(A, B);	×
NEW 133	auto t2 = std::chrono::high_resolution_clock::now();	×
134
NEW 135	double elapsed = std::chrono::duration<double, std::milli>(t2 - t1).count();	×
NEW 136	times_solve.push_back(elapsed);	×
NEW 137	}	×
138
139	// Benchmark rsolve (right-solve, used in Kriging)
140	{
NEW 141	arma::mat A = arma::randu(n, n);	×
NEW 142	A = A * A.t() + arma::eye(n, n) * 0.1; // Make positive definite	×
NEW 143	arma::mat B = arma::randu(10, n);	×
144
NEW 145	auto t1 = std::chrono::high_resolution_clock::now();	×
NEW 146	arma::mat X_sol = LinearAlgebra::rsolve(A, B);	×
NEW 147	auto t2 = std::chrono::high_resolution_clock::now();	×
148
NEW 149	double elapsed = std::chrono::duration<double, std::milli>(t2 - t1).count();	×
NEW 150	times_rsolve.push_back(elapsed);	×
NEW 151	}	×
152	}
153
154	// Print results
NEW 155	print_header();	×
NEW 156	print_stats("covMat_sym_X", compute_stats(times_covMat_sym_X));	×
NEW 157	print_stats("covMat_rect", compute_stats(times_covMat_rect));	×
NEW 158	print_stats("solve", compute_stats(times_solve));	×
NEW 159	print_stats("rsolve", compute_stats(times_rsolve));	×
NEW 160	std::cout << std::endl;	×
NEW 161	}	×
162
NEW 163	int main(int argc, char* argv[]) {	×
NEW 164	int n_iterations = 10;	×
NEW 165	if (argc > 1) {	×
NEW 166	n_iterations = std::atoi(argv[1]);	×
167	}
168
NEW 169	std::cout << "=============================================================================" << std::endl;	×
NEW 170	std::cout << " LinearAlgebra Performance Benchmark (OpenMP) " << std::endl;	×
NEW 171	std::cout << "=============================================================================" << std::endl;	×
172
173	#ifdef _OPENMP
NEW 174	int max_threads = omp_get_max_threads();	×
NEW 175	int optimal_threads = (max_threads > 2) ? 2 : max_threads;	×
NEW 176	std::cout << "OpenMP enabled: max_threads=" << max_threads	×
NEW 177	<< ", using optimal=" << optimal_threads << std::endl;	×
178	#else
179	std::cout << "OpenMP disabled: serial execution" << std::endl;
180	#endif
181
NEW 182	std::cout << "Testing covariance matrix operations with Kriging-realistic matrices" << std::endl;	×
NEW 183	std::cout << "=============================================================================" << std::endl;	×
184
185	// Test various sizes relevant to Kriging
NEW 186	benchmark_configuration(100, 2, n_iterations);	×
NEW 187	benchmark_configuration(100, 4, n_iterations);	×
NEW 188	benchmark_configuration(200, 2, n_iterations);	×
NEW 189	benchmark_configuration(200, 4, n_iterations);	×
NEW 190	benchmark_configuration(500, 2, n_iterations);	×
NEW 191	benchmark_configuration(500, 4, n_iterations);	×
NEW 192	benchmark_configuration(1000, 2, n_iterations);	×
NEW 193	benchmark_configuration(1000, 4, n_iterations);	×
194
NEW 195	std::cout << "=============================================================================" << std::endl;	×
NEW 196	std::cout << "Benchmark complete" << std::endl;	×
NEW 197	std::cout << "=============================================================================" << std::endl;	×
198
NEW 199	return 0;	×
200	}

libKriging / libKriging / 20671555286

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous