1 // ============================================================================= 2 // PROJECT CHRONO - http://projectchrono.org 3 // 4 // Copyright (c) 2014 projectchrono.org 5 // All rights reserved. 6 // 7 // Use of this source code is governed by a BSD-style license that can be found 8 // in the LICENSE file at the top level of the distribution and at 9 // http://projectchrono.org/license-chrono.txt. 10 // 11 // ============================================================================= 12 // Authors: Dario Mangoni 13 // ============================================================================= 14 15 #ifndef CHPARDISOPROJECTENGINE_H 16 #define CHPARDISOPROJECTENGINE_H 17 18 #include "chrono/core/ChMatrix.h" 19 #include "chrono_pardisoproject/ChApiPardisoProject.h" 20 21 22 namespace chrono { 23 24 /// @addtogroup pardisoproject_module 25 /// @{ 26 27 /// Wrapper class for the PardisoProject direct linear solver. 28 /// This solver is not appropriate for VI and complementarity problems. 29 /// \warning {[WARNING for DEVELOPERS]: please consider that the solver does not handle C-like arrays so the matrices have to be set to one-indexed format manually. 30 /// Appropriate instructions have been put to avoid passing back the matrix in one-indexed state when not possible} 31 32 class ChApiPardisoProject ChPardisoProjectEngine { 33 public: 34 enum parproj_SYM { 35 STRUCTURAL_SYMMETRIC = 1, ///< real and structurally symmetric, supernode pivoting 36 SYMMETRIC_POSDEF = 2, ///< real and symmetric positive definite 37 SYMMETRIC_GENERAL = -2, ///< real and symmetric indefinite, diagonal or Bunch-Kaufman pivoting 38 UNSYMMETRIC = 11 ///< real and nonsymmetric, complete supernode pivoting 39 }; 40 41 enum parproj_PHASE { 42 END = -1, 43 ANALYZE = 11, 44 ANALYZE_FACTORIZE = 12, 45 FACTORIZE = 22, 46 SOLVE = 33, 47 FACTORIZE_SOLVE = 23, 48 COMPLETE = 13, 49 SELECTED_INVERSION = -22 50 }; 51 52 ChPardisoProjectEngine(parproj_SYM symmetry); 53 ~ChPardisoProjectEngine(); 54 55 /// Set the problem matrix and the right-hand side. 56 void SetProblem(const ChSparseMatrix& Z, ChVectorRef rhs, ChVectorRef sol); 57 58 /// Set the problem matrix. 59 void SetMatrix(const ChSparseMatrix& Z, bool isZeroIndexed = true); 60 void SetMatrix(int n, int *ia, int *ja, double *a, bool isZeroIndexed = true); 61 62 /// Set a new value for symmetry. \warning{This function triggers a Reinit()} 63 void SetMatrixSymmetry(parproj_SYM symmetry); 64 65 /// Set the right-hand side vector. 66 /// Note that it is the caller's responsibility to ensure that the size is appropriate. 67 void SetRhsVector(ChVectorRef b); 68 void SetRhsVector(double* b); 69 70 /// Set the solution vector. 71 /// Note that it is the caller's responsibility to ensure that the size is appropriate. 72 void SetSolutionVector(ChVectorRef x); 73 void SetSolutionVector(double* x); 74 75 /// Submit job to PardisoProject. 76 int PardisoProjectCall(parproj_PHASE job_call); 77 78 /// Check if the input matrix is in the appropriate form. 79 int CheckMatrix(bool print=true); 80 81 /// Check if the rhs vector is in the appropriate form. 82 int CheckRhsVectors(bool print=true); 83 84 /// Check if the input matrix has appropriate properties. 85 int CheckMatrixStats(bool print=true); 86 87 /// Return the value of the i-th IPARM coefficient; consider that is in zero-indexed format. GetIPARM(int id)88 int GetIPARM(int id) const { return iparm[id]; }; 89 /// Set the value of the i-th IPARM coefficient; consider that is in zero-indexed format. SetIPARM(int id,int val)90 void SetIPARM(int id, int val){ iparm[id] = val; }; 91 92 /// Return the value of the i-th DPARM coefficient; consider that is in zero-indexed format. GetDPARM(int id)93 double GetDPARM(int id) const { return dparm[id]; }; 94 /// Set the value of the i-th DPARM coefficient; consider that is in zero-indexed format. SetDPARM(int id,int val)95 void SetDPARM(int id, int val){ dparm[id] = val; }; 96 97 /// Return the value of the error flag. GetLastError()98 int GetLastError() {return error;}; 99 100 /// Set the index of the underlying arrays to zero-indexed. 101 void SetZeroIndexedFormat(); 102 /// Set the index of the underlying arrays to one-indexed. 103 void SetOneIndexedFormat(); 104 105 /// Set the solver type. \p directsparse = true for direct sparse solver, false for multi-recursive iterative solver. \warning{This function triggers a Reinit()} 106 void SetSolverType(bool directsparse = true); 107 108 /// Reinitialize the solver (e.g. when a new symmetry option is set) 109 void Reinit(); 110 111 /// Set the message level verbosity (0: no messages, 1:print stats). (Default: 0) 112 void SetMessageLevel(int msglvl); 113 114 /// Set the maximum number of numerical factorizations. SetMaxNumericalFactorization(int maxfct)115 void SetMaxNumericalFactorization(int maxfct){ maxfct = maxfct; }; 116 117 /// Set the maximum number of numerical factorizations. 118 void GetSchurComplement(ChSparseMatrix& Z, int nrows); 119 120 protected: 121 /// Shift the matrix indeces of the internal matrix arrays of a value of \p val 122 void shiftMatrixIndices(int* ext_ia, int* ext_ja, double* ext_a, int ext_n, int val, bool isOneIndexed); 123 124 /// Shift the matrix indeces of a value of \p val 125 void shiftInternalMatrixIndices(int val); 126 127 /// Get internal arrays of sparse matrix 128 bool getMatrixInternalArrays(ChSparseMatrix& sparseMat, int** ext_ia, int** ext_ja, double** ext_a); 129 130 private: 131 132 bool matOneIndexedFormat = false; 133 134 /* RHS and solution vectors. */ 135 double *b, *x; 136 int nrhs = 1; /* Number of right hand sides. */ 137 138 /* Internal solver memory pointer pt, */ 139 /* 32-bit: int pt[64]; 64-bit: long int pt[64] */ 140 /* or void *pt[64] should be OK on both architectures */ 141 void *pt[64]; 142 143 /* Pardiso control parameters. */ 144 int iparm[64]; 145 double dparm[64]; 146 int solver = 0; /* Solver type: 0: sparse direct solver; 1: multi-recursive iterative solver*/ 147 int maxfct = 1; /* Maximum number of numerical factorizations in memory.*/ 148 int mnum = 1; /* Which factorization to use. */ 149 int error = 0; /* Initialize error flag */ 150 int msglvl = 0; /* Print statistical information */ 151 152 /* Number of processors. */ 153 int num_procs; 154 155 /* Auxiliary variables. */ 156 char *var; 157 int i; 158 159 double ddum; /* Double dummy */ 160 int idum; /* Integer dummy. */ 161 162 // Matrix variables 163 int n = 0; 164 int *ia; 165 int *ja; 166 double *a; 167 168 parproj_SYM symmetry; 169 170 }; 171 172 /// @} pardisoproject_module 173 174 } // end of namespace chrono 175 176 #endif 177