matrix_structure.hpp

Go to the documentation of this file.

#pragma once

#include "./mpi_structure.hpp"
#include <limits>
#include <iostream>
#include <cmath>
#include <cstdlib>

#include "config_structure.hpp"
#include "geometry_structure.hpp"
#include "vector_structure.hpp"

using namespace std;

const su2double eps = numeric_limits<su2double>::epsilon(); 
class CSysMatrix {
private:
    unsigned long nPoint,   
    nPointDomain,           
    nVar,                   
    nEqn;                   
    su2double *matrix;            
    su2double *ILU_matrix;         
    unsigned long *row_ptr;    
    unsigned long *col_ind;    
    unsigned long nnz;         
    su2double *block;             
    su2double *block_inverse;             
    su2double *block_weight;             
  su2double *prod_block_vector; 
    su2double *prod_row_vector;   
    su2double *aux_vector;         
  su2double *sum_vector;         
    su2double *invM;              
    bool *LineletBool;                          
    vector<unsigned long> *LineletPoint;        
    unsigned long nLinelet;                     
  su2double **UBlock, **invUBlock, **LBlock,
  **yVector, **zVector, **rVector, *LFBlock,
  *LyVector, *FzVector;           
  unsigned long max_nElem;
  
public:
  
    CSysMatrix(void);
  
    ~CSysMatrix(void);
  
  void Initialize(unsigned long nPoint, unsigned long nPointDomain, unsigned short nVar, unsigned short nEqn,
                  bool EdgeConnect, CGeometry *geometry, CConfig *config);
  
    void SetIndexes(unsigned long val_nPoint, unsigned long val_nPointDomain, unsigned short val_nVar, unsigned short val_nEq, unsigned long* val_row_ptr, unsigned long* val_col_ind, unsigned long val_nnz, CConfig *config);
  
    void SetValZero(void);
  
    su2double *GetBlock(unsigned long block_i, unsigned long block_j);
  
    su2double GetBlock(unsigned long block_i, unsigned long block_j, unsigned short iVar, unsigned short jVar);
  
    void SetBlock(unsigned long block_i, unsigned long block_j, su2double **val_block);
  
    void SetBlock(unsigned long block_i, unsigned long block_j, su2double *val_block);
  
    void AddBlock(unsigned long block_i, unsigned long block_j, su2double **val_block);
  
    void SubtractBlock(unsigned long block_i, unsigned long block_j, su2double **val_block);
  
    su2double *GetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j);
  
    void SetBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, su2double *val_block);
  

  void SetBlockTransposed_ILUMatrix(unsigned long block_i, unsigned long block_j, su2double *val_block);

    void SubtractBlock_ILUMatrix(unsigned long block_i, unsigned long block_j, su2double *val_block);
  
    void AddVal2Diag(unsigned long block_i, su2double val_matrix);
  
  void SetVal2Diag(unsigned long block_i, su2double val_matrix);

  void MatrixVectorProduct(su2double *matrix, su2double *vector, su2double *product);
  
  void MatrixMatrixProduct(su2double *matrix_a, su2double *matrix_b, su2double *product);
  
    void DeleteValsRowi(unsigned long i);
  
  su2double MatrixDeterminant(su2double **a, unsigned long n);

  void MatrixCoFactor(su2double **a, unsigned long n, su2double **b) ;

  void MatrixTranspose(su2double **a, unsigned long n) ;
  
  void Gauss_Elimination(unsigned long block_i, su2double* rhs, bool transposed = false);
  
    void Gauss_Elimination(su2double* Block, su2double* rhs);

    void Gauss_Elimination_ILUMatrix(unsigned long block_i, su2double* rhs);
  
    void ProdBlockVector(unsigned long block_i, unsigned long block_j, const CSysVector & vec);
  
    void UpperProduct(CSysVector & vec, unsigned long row_i);
  
    void LowerProduct(CSysVector & vec, unsigned long row_i);
  
    void DiagonalProduct(CSysVector & vec, unsigned long row_i);

    void SendReceive_Solution(CSysVector & x, CGeometry *geometry, CConfig *config);

  void SendReceive_SolutionTransposed(CSysVector & x, CGeometry *geometry, CConfig *config);

    void RowProduct(const CSysVector & vec, unsigned long row_i);
  
    void MatrixVectorProduct(const CSysVector & vec, CSysVector & prod);
  
    void MatrixVectorProduct(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config);
    
  void MatrixVectorProductTransposed(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config);

    void GetMultBlockBlock(su2double *c, su2double *a, su2double *b);
    
    void GetMultBlockVector(su2double *c, su2double *a, su2double *b);
    
    void GetSubsBlock(su2double *c, su2double *a, su2double *b);
    
    void GetSubsVector(su2double *c, su2double *a, su2double *b);
  
  void InverseDiagonalBlock(unsigned long block_i, su2double *invBlock, bool transpose = false);
  
    void InverseDiagonalBlock_ILUMatrix(unsigned long block_i, su2double *invBlock);
 
    void InverseBlock(su2double *Block, su2double *invBlock);
  
  void BuildJacobiPreconditioner(bool transpose = false);
    
    void ComputeJacobiPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config);
  
  unsigned long Jacobi_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, su2double tol, unsigned long m, su2double *residual, bool monitoring, CGeometry *geometry, CConfig *config);
  
  void BuildILUPreconditioner(bool transposed = false);
  
    void ComputeILUPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config);
  
  unsigned long ILU0_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, su2double tol, unsigned long m, su2double *residual, bool monitoring, CGeometry *geometry, CConfig *config);

    void ComputeLU_SGSPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config);
  
  unsigned long LU_SGS_Smoother(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec, su2double tol, unsigned long m, su2double *residual, bool monitoring, CGeometry *geometry, CConfig *config);
  
  unsigned short BuildLineletPreconditioner(CGeometry *geometry, CConfig *config);

    void ComputeLineletPreconditioner(const CSysVector & vec, CSysVector & prod, CGeometry *geometry, CConfig *config);

  void ComputeResidual(const CSysVector & sol, const CSysVector & f, CSysVector & res);

};

class CSysMatrixVectorProduct : public CMatrixVectorProduct {
private:
    CSysMatrix* sparse_matrix; 
    CGeometry* geometry; 
    CConfig* config; 
public:
  
    CSysMatrixVectorProduct(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref);
  
    ~CSysMatrixVectorProduct() {}
  
    void operator()(const CSysVector & u, CSysVector & v) const;
};

class CSysMatrixVectorProductTransposed : public CMatrixVectorProduct {
private:
  CSysMatrix* sparse_matrix; 
  CGeometry* geometry; 
  CConfig* config; 
public:

  CSysMatrixVectorProductTransposed(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref);

  ~CSysMatrixVectorProductTransposed() {}

  void operator()(const CSysVector & u, CSysVector & v) const;
};

class CJacobiPreconditioner : public CPreconditioner {
private:
    CSysMatrix* sparse_matrix; 
    CGeometry* geometry; 
    CConfig* config; 
public:
  
    CJacobiPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref);
  
    ~CJacobiPreconditioner() {}
  
    void operator()(const CSysVector & u, CSysVector & v) const;
};

class CJacobiTransposedPreconditioner : public CPreconditioner {
private:
  CSysMatrix* sparse_matrix; 
  CGeometry* geometry; 
  CConfig* config; 
public:

  CJacobiTransposedPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref);

  ~CJacobiTransposedPreconditioner() {}

  void operator()(const CSysVector & u, CSysVector & v) const;
};

class CILUPreconditioner : public CPreconditioner {
private:
    CSysMatrix* sparse_matrix; 
    CGeometry* geometry; 
    CConfig* config; 
public:
  
    CILUPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref);
  
    ~CILUPreconditioner() {}
  
    void operator()(const CSysVector & u, CSysVector & v) const;
};

class CLU_SGSPreconditioner : public CPreconditioner {
private:
    CSysMatrix* sparse_matrix; 
  CGeometry* geometry; 
    CConfig* config; 
public:
    
    CLU_SGSPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref);
    
    ~CLU_SGSPreconditioner() {}
    
    void operator()(const CSysVector & u, CSysVector & v) const;
};

class CLineletPreconditioner : public CPreconditioner {
private:
    CSysMatrix* sparse_matrix; 
  CGeometry* geometry; 
    CConfig* config; 
public:
    
    CLineletPreconditioner(CSysMatrix & matrix_ref, CGeometry *geometry_ref, CConfig *config_ref);
    
    ~CLineletPreconditioner() {}
    
    void operator()(const CSysVector & u, CSysVector & v) const;
};

#include "matrix_structure.inl"