MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
harwrap.h
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/harwrap.h,v 1.40 2017/01/12 14:44:25 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 /*****************************************************************************
33  * *
34  * HARWLIB C++ WRAPPER *
35  * *
36  *****************************************************************************/
37 
38 /*
39  * Wrapper for Harwell sparse LU solution
40  * http://www.netlib.org/harwell/
41  */
42 
43 /*
44  * Uso:
45  * l'HarwellSparseSolutionManager rende disponibile lo spazio per la matrice,
46  * la soluzione ed il termine noto. Inoltre rende disponibili
47  * metodi per l'handling della matrice e dei vettori
48  * (inserimento e lettura dei coefficienti). Infine consente la soluzione del
49  * problema lineare mediante fattorizzazione LU e successiva sostituzone.
50  *
51  * Uso consigliato:
52  * // creare un oggetto HarwellSparseSolutionManager:
53  * HarwellSparseSolutionManager SM(size, workspacesize, pivotfactor);
54  *
55  * // ottenere l'handling alla matrice ed al termine noto
56  * SparseMatrixHandler* pMH = SM.pMatHdl();
57  * VectorHandler* pRH = SM.pResHdl();
58  *
59  * // Ogni volta che si desidera riassemblare la matrice:
60  * // inizializzare lo SparseMatrixHandler
61  * SM.MatrReset();
62  *
63  * // Operare sulla matrice e sui vettori
64  * pMH->PutCoef(row, col, coef);
65  * coef = pMH->operator()(row, col);
66  * pRH->PutCoef(row, coef);
67  * coef = pRH->operator()(row);
68  *
69  * // Risolvere il problema; in questa fase si assume che:
70  * // - la matrice sia stata completata;
71  * // - il termine noto sia stato completato.
72  * // In tale caso viene eseguita la fattorizzazione e quindi la soluzione.
73  * // Se la fattorizzazione e' gia' stata compiuta, la matrice non e' piu'
74  * // stata modificata e si e' modificato il termine noto, la chiamata a
75  * // HarwellSparseSolutionManager::Solve()
76  * // semplicemente esegue una nuova sostituzione all'indietro.
77  * SM.Solve();
78  *
79  * // Se si desidera modificare la matrice per una nuova soluzione, occorre
80  * // inizializzare di nuovo lo SparseMatrixHandler, con:
81  * SM.MatrReset();
82  *
83  * // Per i parametri di inizializzazione e per eventuali modifiche
84  * // fare riferimento al codice sorgente ed alla libreria originaria
85  * // in FORTRAN od al suo porting in C, files <harwlib.f>, <harwlib.c>
86  * // Il pacchetto e' costituito dai files:
87  *
88  * <harwlib.f> libreria originaria in FORTRAN
89  * <harwlib.c> conversione in C mediante f2c
90  * <harwlib.h> header per <harwlib.c> con dichiarazione delle funzioni
91  * e dei common
92  * <harwrap.cc> wrapper in C++
93  * <harwrap.h> header per <harwrap.cc>.
94  * E' sufficiente includere questo per usare il wrapper.
95  */
96 
97 #ifndef HARWRAP_H
98 #define HARWRAP_H
99 
100 #ifdef USE_HARWELL
101 
102 #include <myassert.h>
103 #include <mynewmem.h>
104 #include <except.h>
105 #include <harwlib.h>
106 #include <solman.h>
107 #include <ls.h>
108 #include <submat.h>
109 #include <spmapmh.h>
110 
111 /* classi dichiarate in questo file */
112 class HarwellSolver; /* solutore */
113 class HarwellSparseSolutionManager; /* gestore della soluzione */
114 
115 
116 /* HarwellSolver - begin */
117 
118 /*
119  * Solutore LU per matrici sparse. usa spazio messo a disposizione da altri
120  * e si aspetta le matrici gia' bell'e preparate
121  */
122 
123 static char sLUClassName[] = "HarwellSolver";
124 
125 class HarwellSolver {
126  friend class HarwellSparseSolutionManager;
127 
128 
129 public:
130  class ErrFactorization : public MBDynErrBase {
131  private:
132  integer iErrCode;
133  public:
134  ErrFactorization(integer i, MBDYN_EXCEPT_ARGS_DECL) :
136  NO_OP;
137  };
138  integer iGetErrCode(void) const {
139  return iErrCode;
140  };
141  };
142 
143 private:
144  integer iMatSize;
145 
146  std::vector<integer>*const piRow;
147  std::vector<integer>*const piCol;
148  std::vector<doublereal>*const pdMat;
149 
150  integer iN; /* ordine della matrice */
151  integer iNonZeroes; /* coeff. non nulli */
152  doublereal* pdRhs; /* Soluzione e termine noto */
153 
154  doublereal dU; /* parametro di pivoting */
155  integer* piKeep; /* vettore di lavoro */
156  integer* piW; /* vettore di lavoro */
157  doublereal* pdW; /* vettore di lavoro */
158 
159 protected:
160  /* Costruttore: si limita ad allocare la memoria */
161  HarwellSolver(integer iMatOrd, integer iSize,
162  std::vector<integer>*const piTmpRow,
163  std::vector<integer>*const piTmpCol,
164  std::vector<doublereal>*const pdTmpMat,
165  doublereal* pdTmpRhs, doublereal dPivotFact)
166  : iMatSize(iSize),
167  piRow(piTmpRow),
168  piCol(piTmpCol),
169  pdMat(pdTmpMat),
170  iN(iMatOrd),
171  iNonZeroes(0),
172  pdRhs(pdTmpRhs),
173  dU(dPivotFact),
174  piKeep(NULL),
175  piW(NULL),
176  pdW(NULL) {
177  ASSERT(iMatSize > 0);
178  ASSERT(piRow != NULL);
179  ASSERT(piCol != NULL);
180  ASSERT(pdMat != NULL);
181  ASSERT(*ppiRow != NULL);
182  ASSERT(*ppiCol != NULL);
183  ASSERT(*ppdMat != NULL);
184  ASSERT(pdRhs != NULL);
185  ASSERT(iN > 0);
186 
187 #ifdef DEBUG_MEMMANAGER
188  ASSERT(piRow->Size()==iMatSize);
189  ASSERT(piCol->Size()==iMatSize);
190  ASSERT(pdMat->Size()==iMatSize);
191  ASSERT(defaultMemoryManager.fIsValid(pdRhs,
192  iN*sizeof(doublereal)));
193 #endif /* DEBUG_MEMMANAGER */
194 
195  SAFENEWARR(piKeep, integer, 5*iN);
196  SAFENEWARR(piW, integer, 8*iN);
197  SAFENEWARR(pdW, doublereal, iN);
198 
199 #ifdef DEBUG
200  for (integer iCnt = 0; iCnt < 5*iN; iCnt++) {
201  piKeep[iCnt] = 0;
202  }
203  for (integer iCnt = 0; iCnt < 8*iN; iCnt++) {
204  piW[iCnt] = 0;
205  }
206  for (integer iCnt = 0; iCnt < 1*iN; iCnt++) {
207  pdW[iCnt] = 0.;
208  }
209 #endif /* DEBUG */
210  };
211 
212  /* Distruttore */
213  ~HarwellSolver(void) {
214  if (pdW != NULL) {
215  SAFEDELETEARR(pdW);
216  }
217  if (piW != NULL) {
218  SAFEDELETEARR(piW);
219  }
220  if (piKeep != NULL) {
221  SAFEDELETEARR(piKeep);
222  }
223  };
224 
225 #ifdef DEBUG
226  void IsValid(void) const {
227  ASSERT(iMatSize > 0);
228  ASSERT(piRow != NULL);
229  ASSERT(piCol != NULL);
230  ASSERT(pdMat != NULL);
231  ASSERT(pdRhs != NULL);
232  ASSERT(iN > 0);
233 
234 #ifdef DEBUG_MEMMANAGER
235  ASSERT(piRow->Size()==iMatSize);
236  ASSERT(piCol->Size()==iMatSize);
237  ASSERT(pdMat->Size()==iMatSize);
238  ASSERT(defaultMemoryManager.fIsValid(pdRhs,
239  iN*sizeof(doublereal)));
240 #endif /* DEBUG_MEMMANAGER */
241 
242  ASSERT(piKeep != NULL);
243  ASSERT(piW != NULL);
244  ASSERT(pdW != NULL);
245 
246 #ifdef DEBUG_MEMMANAGER
247  ASSERT(defaultMemoryManager.fIsBlock(piKeep,
248  5*iN*sizeof(integer)));
249  ASSERT(defaultMemoryManager.fIsBlock(piW,
250  8*iN*sizeof(integer)));
251  ASSERT(defaultMemoryManager.fIsBlock(pdW,
252  1*iN*sizeof(doublereal)));
253 #endif /* DEBUG_MEMMANAGER */
254  };
255 #endif /* DEBUG */
256 
257  /* Fattorizza la matrice */
258  bool bLUFactor(void) {
259 #ifdef DEBUG
260  IsValid();
261 #endif /* DEBUG */
262 
263  ASSERT(iNonZeroes > 0);
264 
265  integer iLicn = iMatSize;
266  integer iLirn = iMatSize;
267  integer iFlag = 0;
268 
269  DEBUGCOUT("Calling ma28ad_()," << std::endl
270  << "iN = " << iN << std::endl
271  << "iNonZeroes = " << iNonZeroes << std::endl
272  << "pdMat = " << *ppdMat << std::endl
273  << "iLicn = " << iLicn << std::endl
274  << "piRow = " << *ppiRow << std::endl
275  << "iLirn = " << iLirn << std::endl
276  << "piCol = " << *ppiCol << std::endl
277  << "dU = " << dU << std::endl
278  << "piKeep = " << piKeep << std::endl
279  << "piW = " << piW << std::endl
280  << "iFlag = " << iFlag << std::endl);
281 
282  __FC_DECL__(ma28ad)(&iN, &iNonZeroes, &((*pdMat)[0]),
283  &iLicn, &((*piRow)[0]),
284  &iLirn, &((*piCol)[0]),
285  &dU, piKeep, piW, pdW,
286  &iFlag);
287 
288  if (iFlag < 0) {
289  silent_cerr(sLUClassName
290  << ": error during factorization, code "
291  << iFlag << std::endl);
292  throw ErrFactorization(iFlag);
293  }
294 
295  /* FIXME: handle iFlag > 0 ??? */
296  return true;
297  };
298 
299  /* Risolve */
300  void Solve(void) {
301 #ifdef DEBUG
302  IsValid();
303 #endif /* DEBUG */
304 
305  integer iLicn = iMatSize;
306  integer iMtype = 1;
307  __FC_DECL__(ma28cd)(&iN, &((*pdMat)[0]), &iLicn, &((*piCol)[0]),
308  piKeep, pdRhs, pdW, &iMtype);
309  };
310 };
311 
312 /* HarwellSolver - end */
313 
314 
315 /* HarwellSparseSolutionManager - begin */
316 
317 /*
318  * Gestisce la soluzione del problema; alloca le matrici occorrenti
319  * e gli oggetti dedicati alla gestione delle matrici ed alla soluzione
320  */
321 
322 class HarwellSparseSolutionManager : public SolutionManager {
323 protected:
324  integer iMatMaxSize; /* Dimensione max della matrice (per resize) */
325  integer iMatSize; /* ordine della matrice */
326  std::vector<integer> iRow; /* array di interi con:
327  * tabella di SparseData/indici di riga
328  * di HarwellSolver */
329  std::vector<integer> iCol; /* array di interi con:
330  * keys di SparseData/indici di colonna
331  * di HarwellSolver */
332  std::vector<integer> iColStart;
333  std::vector<doublereal> dMat; /* array di reali con la matrice */
334  std::vector<doublereal> dVec; /* array di reali con residuo/soluzione */
335 
336  mutable SpMapMatrixHandler MH; /* SparseMatrixHandler */
337 /* SparseMatrixHandler* pMH; puntatore a SparseMatrixHandler */
338  VectorHandler* pVH; /* puntatore a VectorHandler */
339  HarwellSolver* pLU; /* puntatore a HarwellSolver */
340 
341  flag fHasBeenReset; /* flag di matrice resettata */
342 
343  /* Prepara i vettori e la matrice per il solutore */
344  void PacVec(void);
345 
346  /* Fattorizza la matrice (non viene mai chiamato direttamente,
347  * ma da Solve se la matrice ancora non e' stata fattorizzata) */
348  void Factor(void);
349 
350 public:
351  /* Costruttore: usa e mette a disposizione matrici che gli sono date */
352  HarwellSparseSolutionManager(integer iSize,
353  integer iWorkSpaceSize = 0,
354  const doublereal& dPivotFactor = 1.0);
355 
356  /* Distruttore: dealloca le matrici e distrugge gli oggetti propri */
357  ~HarwellSparseSolutionManager(void);
358 
359  /* Usata per il debug */
360 #ifdef DEBUG
361  void IsValid(void) const;
362 #endif /* DEBUG */
363 
364  /* Inizializza il gestore delle matrici */
365  void MatrReset();
366 
367  /* Risolve il sistema */
368  void Solve(void);
369 
370  /* sposta il puntatore al vettore del residuo */
371  void pdSetResVec(doublereal* pRes){
372  pLU->pdRhs = pRes;
373  };
374 
375  /* sposta il puntatore al vettore del residuo */
376  void pdSetSolVec(doublereal* pSol){
377  pLU->pdRhs = pSol;
378  };
379 
380  /* Rende disponibile l'handler per la matrice */
381  MatrixHandler* pMatHdl(void) const {
382  return &MH;
383  };
384 
385  /* Rende disponibile l'handler per il termine noto */
386  VectorHandler* pResHdl(void) const {
387  ASSERT(pVH != NULL);
388  return pVH;
389  };
390 
391  /* Rende disponibile l'handler per la soluzione (e' lo stesso
392  * del termine noto, ma concettualmente sono separati) */
393  VectorHandler* pSolHdl(void) const {
394  ASSERT(pVH != NULL);
395  return pVH;
396  };
397 };
398 
399 /* HarwellSparseSolutionManager - end */
400 
401 #endif /* USE_HARWELL */
402 
403 #endif /* HARWRAP_H */
404 
virtual VectorHandler * pResHdl(void) const =0
long int flag
Definition: mbdyn.h:43
doublereal * pdSetSolVec(doublereal *pd)
Definition: solman.cc:101
#define MBDYN_EXCEPT_ARGS_PASSTHRU
Definition: except.h:55
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
#define MBDYN_EXCEPT_ARGS_DECL
Definition: except.h:43
int ma28cd(integer *n, doublereal *a, integer *licn, integer *icn, integer *ikeep, doublereal *rhs, doublereal *w, integer *mtype)
#define NO_OP
Definition: myassert.h:74
virtual MatrixHandler * pMatHdl(void) const =0
#define DEBUGCOUT(msg)
Definition: myassert.h:232
doublereal * pdSetResVec(doublereal *pd)
Definition: solman.cc:93
virtual void MatrReset(void)=0
virtual void Solve(void)=0
#define ASSERT(expression)
Definition: colamd.c:977
#define defaultMemoryManager
Definition: mynewmem.h:691
int ma28ad(integer *n, integer *nz, doublereal *a, integer *licn, integer *irn, integer *lirn, integer *icn, doublereal *u, integer *ikeep, integer *iw, doublereal *w, integer *iflag)
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
virtual VectorHandler * pSolHdl(void) const =0
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51