1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/kluwrap.h,v 1.5 2014/01/13 08:59:07 masarati Exp $ */ 2 /* 3 * HmFe (C) is a FEM analysis code. 4 * 5 * Copyright (C) 1996-2006 6 * 7 * Marco Morandini <morandini@aero.polimi.it> 8 * 9 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano 10 * via La Masa, 34 - 20156 Milano, Italy 11 * http://www.aero.polimi.it 12 * 13 * Changing this copyright notice is forbidden. 14 * This program is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 17 * 18 */ 19 /* December 2001 20 * Modified to add a Sparse matrix in row form and to implement methods 21 * to be used in the parallel MBDyn Solver. 22 * 23 * Copyright (C) 2001-2005 24 * 25 * Giuseppe Quaranta <quaranta@aero.polimi.it> 26 * 27 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano 28 * via La Masa, 34 - 20156 Milano, Italy 29 * http://www.aero.polimi.it 30 * 31 */ 32 /* 33 * MBDyn (C) is a multibody analysis code. 34 * http://www.mbdyn.org 35 * 36 * Copyright (C) 1996-2005 37 * 38 * Pierangelo Masarati <masarati@aero.polimi.it> 39 * Paolo Mantegazza <mantegazza@aero.polimi.it> 40 * 41 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano 42 * via La Masa, 34 - 20156 Milano, Italy 43 * http://www.aero.polimi.it 44 * 45 * Changing this copyright notice is forbidden. 46 * 47 * This program is free software; you can redistribute it and/or modify 48 * it under the terms of the GNU General Public License as published by 49 * the Free Software Foundation (version 2 of the License). 50 * 51 * 52 * This program is distributed in the hope that it will be useful, 53 * but WITHOUT ANY WARRANTY; without even the implied warranty of 54 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 55 * GNU General Public License for more details. 56 * 57 * You should have received a copy of the GNU General Public License 58 * along with this program; if not, write to the Free Software 59 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 60 */ 61 62 /* 63 * Umfpack is used by permission; please read its Copyright, 64 * License and Availability note. 65 */ 66 67 #ifndef KLUSparseSolutionManager_hh 68 #define KLUSparseSolutionManager_hh 69 70 #ifdef USE_KLU 71 72 #include <iostream> 73 #include <vector> 74 75 extern "C" { 76 #include <klu.h> 77 } 78 79 #include "myassert.h" 80 #include "mynewmem.h" 81 #include "ls.h" 82 #include "solman.h" 83 #include "spmapmh.h" 84 #include "ccmh.h" 85 #include "dgeequ.h" 86 #include "linsol.h" 87 /* KLUSolver - begin */ 88 89 class KLUSolver: public LinearSolver { 90 private: 91 integer iSize; 92 mutable doublereal *Axp; 93 mutable integer *Aip; 94 mutable integer *App; 95 96 klu_symbolic *Symbolic; 97 mutable klu_common Control; 98 mutable klu_numeric *Numeric; 99 100 bool bPrepareSymbolic(void); 101 102 void Factor(void); 103 104 public: 105 enum Scale { 106 SCALE_NONE = 0, 107 SCALE_SUM = 1, 108 SCALE_MAX = 2, 109 SCALE_UNDEF 110 }; 111 112 KLUSolver(const integer &size, const doublereal &dPivot, Scale scale = SCALE_UNDEF); 113 ~KLUSolver(void); 114 115 void Reset(void); 116 void Solve(void) const; 117 118 void MakeCompactForm(SparseMatrixHandler&, 119 std::vector<doublereal>& Ax, 120 std::vector<integer>& Ar, 121 std::vector<integer>& Ac, 122 std::vector<integer>& Ap) const; 123 124 virtual bool bGetConditionNumber(doublereal& dCond); 125 }; 126 127 /* KLUSolver - end */ 128 129 /* KLUSparseSolutionManager - begin */ 130 131 class KLUSparseSolutionManager: public SolutionManager { 132 protected: 133 mutable SpMapMatrixHandler A; 134 135 std::vector<doublereal> b; 136 137 mutable MyVectorHandler bVH; 138 139 ScaleOpt scale; 140 MatrixScaleBase* pMatScale; 141 142 std::vector<doublereal> Ax; 143 std::vector<integer> Ai; 144 std::vector<integer> Adummy; 145 std::vector<integer> Ap; 146 147 /* Passa in forma di Compressed Column (callback per solve, 148 * richiesto da SpMap e CC Matrix Handler) */ 149 virtual void MakeCompressedColumnForm(void); 150 151 template <typename MH> 152 void ScaleMatrixAndRightHandSide(MH& mh); 153 154 template <typename MH> 155 MatrixScale<MH>& GetMatrixScale(); 156 157 void ScaleSolution(void); 158 159 /* Backward Substitution */ 160 void BackSub(doublereal t_iniz = 0.); 161 162 public: 163 KLUSparseSolutionManager(integer Dim, 164 doublereal dPivot = -1., 165 const ScaleOpt& scale = ScaleOpt()); 166 virtual ~KLUSparseSolutionManager(void); 167 #ifdef DEBUG IsValid(void)168 virtual void IsValid(void) const { 169 NO_OP; 170 }; 171 #endif /* DEBUG */ 172 173 /* Inizializzatore generico */ 174 virtual void MatrReset(void); 175 176 /* Risolve il sistema Backward Substitution; fattorizza se necessario */ 177 virtual void Solve(void); 178 179 /* Rende disponibile l'handler per la matrice */ 180 virtual MatrixHandler* pMatHdl(void) const; 181 182 /* Rende disponibile l'handler per il termine noto */ 183 virtual MyVectorHandler* pResHdl(void) const; 184 185 /* Rende disponibile l'handler per la soluzione */ 186 virtual MyVectorHandler* pSolHdl(void) const; 187 }; 188 189 /* KLUSparseSolutionManager - end */ 190 191 /* KLUSparseCCSolutionManager - begin */ 192 193 template <class CC> 194 class KLUSparseCCSolutionManager: public KLUSparseSolutionManager { 195 protected: 196 bool CCReady; 197 CC *Ac; 198 199 virtual void MatrReset(void); 200 virtual void MakeCompressedColumnForm(void); 201 202 public: 203 KLUSparseCCSolutionManager(integer Dim, 204 doublereal dPivot = -1., 205 const ScaleOpt& scale = ScaleOpt()); 206 virtual ~KLUSparseCCSolutionManager(void); 207 208 /* Inizializzatore "speciale" */ 209 virtual void MatrInitialize(void); 210 211 /* Rende disponibile l'handler per la matrice */ 212 virtual MatrixHandler* pMatHdl(void) const; 213 }; 214 215 /* KLUSparseCCSolutionManager - end */ 216 217 #endif /* USE_UMFPACK */ 218 219 #endif /* KLUSparseSolutionManager_hh */ 220 221