PardisoSolver.h

#ifndef PARDIDSO_SOLVER_H
#define PARDIDSO_SOLVER_H
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iomanip>
 
#include <fstream>
#include <iostream>
 
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/SparseCore>
 
#include "Solver.h"
 
//#define PRINT_PAR
//#define PRINT_OMP
 
// get gflops manually from pardiso output
//#define MEAS_GFLOPS
 
typedef Eigen::VectorXd Vect;
typedef Eigen::SparseMatrix<double> SpMat;
 
using namespace Eigen;
 
 
/* PARDISO prototype. */
extern "C" void pardisoinit (void   *, int    *,   int *, int *, double *, int *);
extern "C" void pardiso     (void   *, int    *,   int *, int *,    int *, int *, 
                            double *, int    *,    int *, int *,   int *, int *,
                            int *, double *, double *, int *, double *);
extern "C" void pardiso_chkmatrix  (int *, int *, double *, int *, int *, int *);
extern "C" void pardiso_chkvec     (int *, int *, double *, int *);
extern "C" void pardiso_printstats (int *, int *, double *, int *, int *, int *,
                            double *, int *);
 
 
class PardisoSolver : public Solver{
 
private:
 
 
    int MPI_rank;           
    /* matrix size */
    int n;                  
    long unsigned int nnz;       
    SpMat Q;                
    int* ia;                
    int* ja;                
    double* a;              
    double* b;              
    double* x;              
    void *pt[64];           
    /* Pardiso control parameters. */
    int      iparm[64];
    double   dparm[64];
    int      maxfct, mnum, phase, error, msglvl, solver;
 
    int      num_procs;     
    /* Auxiliary variables. */
    int      i;
    int      k;
    long unsigned int l;
 
    double   ddum;              
    int      idum;              
    int     mtype;              
    int     init;               
public:
    PardisoSolver(int& MPI_rank);
 
 
    /* ======================================================================== */
    void symbolic_factorization(SpMat& Q, int& init);
 
 
    void factorize(SpMat& Q, double& log_det, double& t_priorLatChol);
 
    void factorize_w_constr(SpMat& Q, const MatrixXd& D, double& log_det, MatrixXd& V);
 
    void factorize_solve(SpMat& Q, Vect& rhs, Vect& sol, double &log_det, double& t_condLatChol, double& t_condLatSolve);
 
    void factorize_solve_w_constr(SpMat& Q, Vect& rhs, const MatrixXd& Dxy, double &log_det, Vect& sol, MatrixXd& V);
    
    void selected_inversion(SpMat& Q, Vect& inv_diag);
 
    void selected_inversion_w_constr(SpMat& Q, const MatrixXd& D, Vect& inv_diag, MatrixXd& V);
 
 
    void compute_inverse_pardiso(MatrixXd& H, MatrixXd& C);
 
 
    ~PardisoSolver();
 
}; // end class
 
 
#endif