decompose_LUP_Cmplx.h

Go to the documentation of this file.

#ifndef DECOMPOSE_LUP_CMPLX_INCLUDED
#define DECOMPOSE_LUP_CMPLX_INCLUDED

#include <valarray>
#include <cfloat>
#include <cmath>

#include "bridge_complex.h"

#include "IO/bridgeIO.h"

class Decompose_LUP_Cmplx
{
 public:
  Decompose_LUP_Cmplx(size_t N) : N(N), N2(2 * N), size(N * N2),
    m_lu(size), m_pivot(N) {}

  void set_matrix(const double *mat);

  // solve Ax = b: vec -> A^{-1} vec
  void solve(double *vec);

  // M -> A^{-1}
  void get_inverse(double *mat_inv);

  // M -> A^{-1} * M
  void mult_inverse(double *mat);

  // return det(A)
  dcomplex determinant();

 private:
  int N;
  int N2;
  int size;
  std::valarray<double> m_lu;
  // pivot index
  std::valarray<int> m_pivot;
  // # of pivot
  // even ->  1;
  // odd  -> -1;
  int m_sign;

  inline size_t re(int i, int j)
  {
    return N2 * i + 2 * j;
  }

  inline size_t im(int i, int j)
  {
    return N2 * i + 2 * j + 1;
  }

  inline size_t re(int i)
  {
    return 2 * i;
  }

  inline size_t im(int i)
  {
    return 2 * i + 1;
  }
};
#endif