Calculate the accuracy of BLR LU using multiplication of L and U factors.

blr-lu-mm.cpp

#include "hicma/hicma.h"

#include <cassert>
#include <cstdint>
#include <tuple>
#include <vector>
#include <iostream>
#include <fstream>

using namespace hicma;
using namespace std;

int main(int argc, char** argv) {
  timing::start("Overall");
  hicma::initialize();
  
  int64_t nleaf = atoi(argv[1]);
  int64_t rank = atoi(argv[2]);
  int64_t N = atoi(argv[3]);
  int64_t admis = atoi(argv[4]);
  int64_t basis = 0;
  int64_t nblocks = N / nleaf;
  assert(basis == NORMAL_BASIS || basis == SHARED_BASIS);

  /* Default parameters for statistics */
  double beta = 0.1;
  double nu = 0.5;//in matern, nu=0.5 exp (half smooth), nu=inf sqexp (inifinetly smooth)
  double noise = 1.e-1;
  double sigma = 1.0;
  
  starsh::exp_kernel_prepare(N, beta, nu, noise, sigma, 3);  

  std::vector<std::vector<double>> randx{get_sorted_random_vector(N)};
  print("Being compression");
  timing::start("Hierarchical compression");
  Hierarchical A(starsh::exp_kernel_fill, randx, N, N, rank, nleaf, admis,
               nblocks, nblocks, basis);
  double comp_time = timing::stop("Hierarchical compression");
  
  Hierarchical A_c = A;

  timing::start("LU decomposition");
  Hierarchical L, U;
  std::tie(L, U) = getrf(A);
  double fact_time = timing::stop("LU decomposition");

  print("Begin multiplication of L and U.");
  timing::start("BLR LU gemm");
  Hierarchical L1(identity, randx, N, N, rank, nleaf, admis,
                 nblocks, nblocks, basis);
  Hierarchical U1(identity, randx, N, N, rank, nleaf, admis,
                 nblocks, nblocks, basis);
  for (int i = 0; i < nblocks; ++i) {
    for (int j = 0; j < nblocks; ++j) {
      if (i == j ) {
        L1(i, j) = L(i,j);
        U1(i, j) = U(i,j);
      }
      else if (j < i) {             // lower triangle
        L1(i, j) = L(i,j);
      }
      else {                    // upper triangle
        U1(i, j) = U(i,j);
      }
      
    }
  }
  auto C = gemm(L1, U1, 1, 0);
  
  timing::start("BLR LU gemm");
  
  print("Calculate l2 error.");
  double error = l2_error(A_c, C);

  std::ofstream file;
  file.open("blr-lu-exp3d-mkl-threads.csv", std::ios::app | std::ios::out);
  file << N << "," << nleaf << "," << rank << "," << admis << ","  << error
       << "," << fact_time << "," << comp_time << ","
       << std::getenv("MKL_NUM_THREADS")
       <<  std::endl;
  file.close();
  
  return 0;
}