Commit d7fa5caf authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Implemente GESV LAPACK inversion

parent 084e3c97

src/libnptm/lapack_calls.cpp

+43 −7
Original line number Diff line number Diff line
@@ -19,22 +19,21 @@ 
 * \brief Implementation of the interface with LAPACK libraries.

 */





#ifdef USE_LAPACK



#include <string>



#ifndef INCLUDE_TYPES_H_

#include "../include/types.h"

#endif



/*

#ifdef USE_LAPACK

#ifdef USE_MKL

#include <mkl_lapacke.h>

#else

#include <lapacke.h>

#endif

*/



#ifdef USE_LAPACK



#ifndef INCLUDE_LOGGING_H_

#include "../include/logging.h"
@@ -64,16 +63,18 @@ void zinvert(dcomplex **mat, np_int n, int &jer, const RuntimeSettings& rs) { 
  jer = 0;

  char buffer[128];

  string message;

  lapack_int info;

  lapack_int info, inc1 = 1;

  lcomplex *arr = &(mat[0][0]);

  lcomplex *arr_orig;

  lcomplex lapack_one = 1.0 + I * 0.0;

  np_int nn = n * n;

  if (rs.use_refinement) {

  if (rs.use_refinement && rs.invert_mode != RuntimeSettings::INV_MODE_RBT) {

    lapack_int inc1 = 1;

    arr_orig = new lcomplex[nn];

    zcopy_(&nn, arr, &inc1, arr_orig, &inc1);

  }

  if (rs.invert_mode == RuntimeSettings::INV_MODE_LU) {

    // >>> LU INVERSION SECTION <<<

    np_int* IPIV = new np_int[n];

    LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, arr, n, IPIV);

    LAPACKE_zgetri(LAPACK_ROW_MAJOR, n, arr, n, IPIV);
@@ -81,8 +82,43 @@ void zinvert(dcomplex **mat, np_int n, int &jer, const RuntimeSettings& rs) { 
    if (rs.use_refinement) {

      info = lapack_newton(rs, arr_orig, n, arr);

    }

  } // inversion mode switch

    // >>> END OF LU INVERSION SECTION <<<

  } else if (rs.invert_mode == RuntimeSettings::INV_MODE_GESV) {

    // >>> GESV INVERSION SECTION <<<

    lcomplex *id = new lcomplex[nn]();

    lapack_int *piv = new lapack_int[n];

    for (lapack_int i = 0; i < n; i++) {

      id[i * (n + 1)] = lapack_one;

    }

    LAPACKE_zgesv(LAPACK_ROW_MAJOR, n, n, arr, n, piv, id, n);

    if (info != LAPACK_SUCCESS) {

      message = "ERROR: call to zgesv_() returned info code " + to_string(info) + "!\n";

      rs.logger->err(message);

      exit(1);

    }

    delete[] piv; // free host memory

    if (rs.use_refinement) {

      info = lapack_newton(rs, arr_orig, n, id);

    }

    zcopy_(&nn, id, &inc1, arr, &inc1);

    delete[] id;

    // >>> END OF GESV INVERSION SECTION <<<

  } else if (rs.invert_mode == RuntimeSettings::INV_MODE_RBT) {

    // >>> RBT INVERSION SECTION <<<

    // RBT inversion not implemented

    message = "ERROR: not implemented!\n";

    rs.logger->err(message);

    exit(1);

    // >>> END OF RBT INVERSION SECTION <<<

  } else if (rs.invert_mode == RuntimeSettings::INV_MODE_SVD) {

    // >>> SVD INVERSION SECTION <<<

    // SVD inversion not implemented

    message = "ERROR: not implemented!\n";

    rs.logger->err(message);

    exit(1);

    // >>> END OF SVD INVERSION SECTION <<<

  } // inversion mode switch

  if (rs.use_refinement && rs.invert_mode != RuntimeSettings::INV_MODE_RBT) {

    delete[] arr_orig;

  }

}