MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
umfpackwrap.h
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/umfpackwrap.h,v 1.51 2017/01/12 14:44:25 masarati Exp $ */
2 /*
3  * HmFe (C) is a FEM analysis code.
4  *
5  * Copyright (C) 1996-2017
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-2017
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-2017
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 UmfpackSparseSolutionManager_hh
68 #define UmfpackSparseSolutionManager_hh
69 
70 #ifdef USE_UMFPACK
71 
72 #include <iostream>
73 #include <vector>
74 
75 extern "C" {
76 #include <umfpack.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 
87 /* UmfpackSolver - begin */
88 
89 class UmfpackSolver: public LinearSolver {
90 public:
91  enum Scale {
92  SCALE_NONE = UMFPACK_SCALE_NONE,
93  SCALE_MAX = UMFPACK_SCALE_MAX,
94  SCALE_SUM = UMFPACK_SCALE_SUM,
95 
96  SCALE_UNDEF
97  };
98 
99 private:
100  integer iSize;
101  mutable doublereal *Axp;
102  mutable integer *Aip;
103  mutable integer *App;
104 
105  void *Symbolic;
106  mutable doublereal Control[UMFPACK_CONTROL];
107  mutable doublereal Info[UMFPACK_INFO];
108  mutable void *Numeric;
109  bool bHaveCond;
110 
111  bool bPrepareSymbolic(void);
112 
113  void Factor(void);
114  void Solve(bool bTranspose) const;
115 
116 public:
117  UmfpackSolver(const integer &size,
118  const doublereal &dPivot,
119  const doublereal &dDropTolerance,
120  const unsigned blockSize,
121  Scale scale = SCALE_UNDEF,
122  integer iMaxIter=-1);
123  ~UmfpackSolver(void);
124 
125  void Reset(void);
126  void Solve(void) const;
127  void SolveT(void) const;
128 
130  std::vector<doublereal>& Ax,
131  std::vector<integer>& Ar,
132  std::vector<integer>& Ac,
133  std::vector<integer>& Ap) const;
134 
135  virtual bool bGetConditionNumber(doublereal& dCond);
136 };
137 
138 /* UmfpackSolver - end */
139 
140 /* UmfpackSparseSolutionManager - begin */
141 
142 class UmfpackSparseSolutionManager: public SolutionManager {
143 protected:
144  mutable SpMapMatrixHandler A;
145 
146  std::vector<doublereal> x;
147  std::vector<doublereal> b;
148 
149  mutable MyVectorHandler xVH, bVH;
150 
151  ScaleOpt scale;
152  MatrixScaleBase* pMatScale;
153 
154  std::vector<doublereal> Ax;
155  std::vector<integer> Ai;
156  std::vector<integer> Adummy;
157  std::vector<integer> Ap;
158 
159  /* Passa in forma di Compressed Column (callback per solve,
160  * richiesto da SpMap e CC Matrix Handler) */
161  virtual void MakeCompressedColumnForm(void);
162 
163  template <typename MH>
164  void ScaleMatrixAndRightHandSide(MH &mh);
165 
166  template <typename MH>
167  MatrixScale<MH>& GetMatrixScale();
168 
169  void ScaleSolution(void);
170 
171  /* Backward Substitution */
172  void BackSub(doublereal t_iniz = 0.);
173 
174 public:
175  UmfpackSparseSolutionManager(integer Dim,
176  doublereal dPivot = -1.,
177  doublereal dDropTolerance = 0.,
178  const unsigned blockSize = 0,
179  const ScaleOpt& scale = ScaleOpt(),
180  integer iMaxIter=-1);
181  virtual ~UmfpackSparseSolutionManager(void);
182 #ifdef DEBUG
183  virtual void IsValid(void) const {
184  NO_OP;
185  };
186 #endif /* DEBUG */
187 
188  /* Inizializzatore generico */
189  virtual void MatrReset(void);
190 
191  /* Risolve il sistema Backward Substitution; fattorizza se necessario */
192  virtual void Solve(void);
193  virtual void SolveT(void);
194 
195  /* Rende disponibile l'handler per la matrice */
196  virtual MatrixHandler* pMatHdl(void) const;
197 
198  /* Rende disponibile l'handler per il termine noto */
199  virtual MyVectorHandler* pResHdl(void) const;
200 
201  /* Rende disponibile l'handler per la soluzione */
202  virtual MyVectorHandler* pSolHdl(void) const;
203 };
204 
205 /* UmfpackSparseSolutionManager - end */
206 
207 /* UmfpackSparseCCSolutionManager - begin */
208 
209 template <class CC>
210 class UmfpackSparseCCSolutionManager: public UmfpackSparseSolutionManager {
211 protected:
212  bool CCReady;
213  CC *Ac;
214 
215  virtual void MatrReset(void);
216  virtual void MakeCompressedColumnForm(void);
217 
218 public:
219  UmfpackSparseCCSolutionManager(integer Dim,
220  doublereal dPivot = -1.,
221  doublereal dDropTolerance = 0.,
222  const unsigned& blockSize = 0,
223  const ScaleOpt& scale = ScaleOpt(),
224  integer iMaxIter=-1);
225  virtual ~UmfpackSparseCCSolutionManager(void);
226 
227  /* Inizializzatore "speciale" */
228  virtual void MatrInitialize(void);
229 
230  /* Rende disponibile l'handler per la matrice */
231  virtual MatrixHandler* pMatHdl(void) const;
232 };
233 
234 /* UmfpackSparseCCSolutionManager - end */
235 
236 #endif /* USE_UMFPACK */
237 
238 #endif /* UmfpackSparseSolutionManager_hh */
239 
virtual void Reset(void)
Definition: ls.cc:68
virtual VectorHandler * pResHdl(void) const =0
virtual void SolveT(void) const
Definition: ls.cc:74
virtual bool bGetConditionNumber(doublereal &dCond)
Definition: ls.cc:132
#define NO_OP
Definition: myassert.h:74
virtual void SolveT(void)
Definition: solman.cc:78
virtual MatrixHandler * pMatHdl(void) const =0
virtual void MatrReset(void)=0
virtual void Solve(void)=0
virtual void Solve(void) const =0
virtual VectorHandler * pSolHdl(void) const =0
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
virtual void MakeCompactForm(SparseMatrixHandler &mh, std::vector< doublereal > &Ax, std::vector< integer > &Ar, std::vector< integer > &Ac, std::vector< integer > &Ap) const
Definition: ls.cc:123