MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
umfpackwrap.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/umfpackwrap.cc,v 1.65 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 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
68 
69 #ifdef USE_UMFPACK
70 #include "solman.h"
71 #include "spmapmh.h"
72 #include "ccmh.h"
73 #include "dirccmh.h"
74 #include "umfpackwrap.h"
75 #include "dgeequ.h"
76 #include <cstring>
77 
78 /* in some cases, int and int32_t differ */
79 
80 #ifdef USE_UMFPACK_LONG
81 
82 #define UMFPACKWRAP_defaults umfpack_dl_defaults
83 #define UMFPACKWRAP_free_symbolic umfpack_dl_free_symbolic
84 #define UMFPACKWRAP_free_numeric umfpack_dl_free_numeric
85 
86 #define UMFPACKWRAP_symbolic(size, app, aip, axp, sym, ctrl, info) \
87  umfpack_dl_symbolic(size, size, app, aip, axp, sym, ctrl, info)
88 
89 #define UMFPACKWRAP_report_info umfpack_dl_report_info
90 #define UMFPACKWRAP_report_status umfpack_dl_report_status
91 #define UMFPACKWRAP_numeric umfpack_dl_numeric
92 #define UMFPACKWRAP_solve umfpack_dl_solve
93 
94 #else // ! USE_UMFPACK_LONG
95 
96 #define UMFPACKWRAP_defaults umfpack_di_defaults
97 #define UMFPACKWRAP_free_symbolic umfpack_di_free_symbolic
98 #define UMFPACKWRAP_free_numeric umfpack_di_free_numeric
99 
100 #define UMFPACKWRAP_symbolic(size, app, aip, axp, sym, ctrl, info) \
101  umfpack_di_symbolic(size, size, app, aip, axp, sym, ctrl, info)
102 
103 #define UMFPACKWRAP_report_info umfpack_di_report_info
104 #define UMFPACKWRAP_report_status umfpack_di_report_status
105 #define UMFPACKWRAP_numeric umfpack_di_numeric
106 #define UMFPACKWRAP_solve umfpack_di_solve
107 
108 #endif // ! USE_UMFPACK_LONG
109 
110 /* required factorization type (A * x = b) */
111 #define SYS_VALUE UMFPACK_A
112 #define SYS_VALUET UMFPACK_Aat
113 
114 /* UmfpackSolver - begin */
115 
116 UmfpackSolver::UmfpackSolver(const integer &size,
117  const doublereal &dPivot,
118  const doublereal &dDropTolerance,
119  const unsigned blockSize,
120  Scale scale,
121  integer iMaxIter)
122 : LinearSolver(0),
123 iSize(size),
124 Axp(0),
125 Aip(0),
126 App(0),
127 Symbolic(0),
128 Numeric(0),
129 bHaveCond(false)
130 {
131  // silence static analyzers
132  memset(&Info[0], 0, sizeof(Info));
133  UMFPACKWRAP_defaults(Control);
134 
135  if (dPivot != -1. && (dPivot >= 0. && dPivot <= 1.)) {
136  /*
137  * 1.0: true partial pivoting
138  * 0.0: treated as 1.0
139  *
140  * default: 0.1
141  */
142  Control[UMFPACK_PIVOT_TOLERANCE] = dPivot;
143  }
144 
145  if (dDropTolerance != 0.) {
146 #ifdef UMFPACK_DROPTOL
147  ASSERT(dDropTolerance > 0.);
148  Control[UMFPACK_DROPTOL] = dDropTolerance;
149 #endif
150  }
151 
152  if (blockSize > 0) {
153  Control[UMFPACK_BLOCK_SIZE] = blockSize;
154  }
155 
156  if (scale != SCALE_UNDEF) {
157  Control[UMFPACK_SCALE] = scale;
158  }
159 
160  if (iMaxIter >= 0) {
161  Control[UMFPACK_IRSTEP] = iMaxIter;
162  }
163 }
164 
165 UmfpackSolver::~UmfpackSolver(void)
166 {
167  UMFPACKWRAP_free_symbolic(&Symbolic);
168  ASSERT(Symbolic == 0);
169 
170  UMFPACKWRAP_free_numeric(&Numeric);
171  ASSERT(Numeric == 0);
172 }
173 
174 void
176 {
177  bHaveCond = false;
178 
179  if (Numeric) {
180  UMFPACKWRAP_free_numeric(&Numeric);
181  ASSERT(Numeric == 0);
182  }
183 
184  bHasBeenReset = true;
185 }
186 
187 void
188 UmfpackSolver::Solve(void) const
189 {
190  Solve(false);
191 }
192 
193 void
194 UmfpackSolver::SolveT(void) const
195 {
196  Solve(true);
197 }
198 
199 void
200 UmfpackSolver::Solve(bool bTranspose) const
201 {
202  if (bHasBeenReset) {
203  const_cast<UmfpackSolver *>(this)->Factor();
204  bHasBeenReset = false;
205  }
206 
207  int status;
208 
209  /*
210  * NOTE: Axp, Aip, App should have been set by * MakeCompactForm()
211  */
212 
213 #ifdef UMFPACK_REPORT
214  doublereal t = t_iniz;
215 #endif /* UMFPACK_REPORT */
216 
217  status = UMFPACKWRAP_solve((bTranspose ? SYS_VALUET : SYS_VALUE),
218  App, Aip, Axp,
220  Numeric, Control, Info);
221 
222  if (status != UMFPACK_OK) {
223  UMFPACKWRAP_report_info(Control, Info);
224  UMFPACKWRAP_report_status(Control, status) ;
225  silent_cerr("UMFPACKWRAP_solve failed" << std::endl);
226 
227  /* de-allocate memory */
228  UMFPACKWRAP_free_numeric(&Numeric);
229  ASSERT(Numeric == 0);
230 
232  }
233 
234  if (Control[UMFPACK_IRSTEP] > 0 && Info[UMFPACK_IR_TAKEN] >= Control[UMFPACK_IRSTEP]) {
235  silent_cerr("warning: UMFPACK_IR_TAKEN = " << Info[UMFPACK_IR_TAKEN]
236  << " >= UMFPACK_IRSTEP = " << Control[UMFPACK_IRSTEP] << std::endl
237  << "\tThe flag \"max iterations\" of the linear solver is too small or the condition number of the Jacobian matrix is too high" << std::endl);
238  }
239 
240 #ifdef UMFPACK_REPORT
241  UMFPACKWRAP_report_info(Control, Info);
242  doublereal t1 = umfpack_timer() - t; /* ?!? not used! */
243 #endif /* UMFPACK_REPORT */
244 }
245 
246 void
247 UmfpackSolver::Factor(void)
248 {
249  int status;
250 
251  /*
252  * NOTE: Axp, Aip, App should have been set by * MakeCompactForm()
253  */
254 
255  if (Symbolic == 0 && !bPrepareSymbolic()) {
257  }
258 
259 #ifdef UMFPACK_REPORT
260  UMFPACKWRAP_report_symbolic ("Symbolic factorization of A",
261  Symbolic, Control) ;
262  UMFPACKWRAP_report_info(Control, Info);
263  doublereal t1 = umfpack_timer() - t; /* ?!? not used! */
264 #endif /* UMFPACK_REPORT */
265 
266  status = UMFPACKWRAP_numeric(App, Aip, Axp, Symbolic,
267  &Numeric, Control, Info);
268  if (status == UMFPACK_ERROR_different_pattern) {
269  UMFPACKWRAP_free_symbolic(&Symbolic);
270  if (!bPrepareSymbolic()) {
272  }
273  status = UMFPACKWRAP_numeric(App, Aip, Axp, Symbolic,
274  &Numeric, Control, Info);
275  }
276 
277  if (status != UMFPACK_OK) {
278  UMFPACKWRAP_report_info(Control, Info);
279  UMFPACKWRAP_report_status(Control, status);
280  silent_cerr("UMFPACKWRAP_numeric failed" << std::endl);
281 
282  /* de-allocate memory */
283  UMFPACKWRAP_free_symbolic(&Symbolic);
284  UMFPACKWRAP_free_numeric(&Numeric);
285  ASSERT(Numeric == 0);
286 
288  } else {
289  bHaveCond = true;
290  }
291 
292 #ifdef UMFPACK_REPORT
293  UMFPACKWRAP_report_numeric ("Numeric factorization of A",
294  Numeric, Control);
295  UMFPACKWRAP_report_info(Control, Info);
296  t1 = umfpack_timer() - t; /* ?!? not used! */
297 #endif /* UMFPACK_REPORT */
298 }
299 
300 void
301 UmfpackSolver::MakeCompactForm(SparseMatrixHandler& mh,
302  std::vector<doublereal>& Ax,
303  std::vector<integer>& Ai,
304  std::vector<integer>& Ac,
305  std::vector<integer>& Ap) const
306 {
307  if (!bHasBeenReset) {
308  return;
309  }
310 
311  mh.MakeCompressedColumnForm(Ax, Ai, Ap, 0);
312 
313  Axp = &(Ax[0]);
314  Aip = &(Ai[0]);
315  App = &(Ap[0]);
316 }
317 
318 bool
319 UmfpackSolver::bPrepareSymbolic(void)
320 {
321  int status;
322 
323  status = UMFPACKWRAP_symbolic(iSize, App, Aip, Axp,
324  &Symbolic, Control, Info);
325  if (status != UMFPACK_OK) {
326  UMFPACKWRAP_report_info(Control, Info) ;
327  UMFPACKWRAP_report_status(Control, status);
328  silent_cerr("UMFPACKWRAP_symbolic failed" << std::endl);
329 
330  /* de-allocate memory */
331  UMFPACKWRAP_free_symbolic(&Symbolic);
332  ASSERT(Symbolic == 0);
333 
334  return false;
335  }
336 
337  return true;
338 }
339 
340 bool UmfpackSolver::bGetConditionNumber(doublereal& dCond)
341 {
342  if (!bHaveCond) {
343  return false;
344  }
345 
346  dCond = 1. / Info[UMFPACK_RCOND];
347 
348  return true;
349 }
350 
351 /* UmfpackSolver - end */
352 
353 /* UmfpackSparseSolutionManager - begin */
354 
355 UmfpackSparseSolutionManager::UmfpackSparseSolutionManager(integer Dim,
356  doublereal dPivot,
357  doublereal dDropTolerance,
358  const unsigned blockSize,
359  const ScaleOpt& s,
360  integer iMaxIter)
361 : A(Dim),
362 x(Dim),
363 b(Dim),
364 xVH(Dim, &x[0]),
365 bVH(Dim, &b[0]),
366 scale(s),
367 pMatScale(0)
368 {
369  UmfpackSolver::Scale uscale = UmfpackSolver::SCALE_UNDEF;
370 
371  switch (scale.algorithm) {
372  case SCALEA_UNDEF:
373  scale.when = SCALEW_NEVER; // Do not scale twice
374  break;
375 
376  case SCALEA_NONE:
377  uscale = UmfpackSolver::SCALE_NONE;
378  scale.when = SCALEW_NEVER;
379  break;
380 
381  case SCALEA_ROW_MAX:
382  uscale = UmfpackSolver::SCALE_MAX;
383  scale.when = SCALEW_NEVER; // Do not scale twice! Use built in scaling from Umfpack
384  break;
385 
386  case SCALEA_ROW_SUM:
387  uscale = UmfpackSolver::SCALE_SUM;
388  scale.when = SCALEW_NEVER; // Do not scale twice! Use built in scaling from Umfpack
389  break;
390 
391  default:
392  // Allocate MatrixScale<T> on demand
393  uscale = UmfpackSolver::SCALE_NONE; // Do not scale twice!
394  }
395 
396  SAFENEWWITHCONSTRUCTOR(pLS, UmfpackSolver,
397  UmfpackSolver(Dim, dPivot, dDropTolerance, blockSize, uscale, iMaxIter));
398 
399  (void)pLS->pdSetResVec(&b[0]);
400  (void)pLS->pdSetSolVec(&x[0]);
401  pLS->SetSolutionManager(this);
402 }
403 
404 UmfpackSparseSolutionManager::~UmfpackSparseSolutionManager(void)
405 {
406  if (pMatScale) {
407  SAFEDELETE(pMatScale);
408  pMatScale = 0;
409  }
410 }
411 
412 void
413 UmfpackSparseSolutionManager::MatrReset(void)
414 {
415  pLS->Reset();
416 }
417 
418 void
419 UmfpackSparseSolutionManager::MakeCompressedColumnForm(void)
420 {
421  ScaleMatrixAndRightHandSide<SpMapMatrixHandler>(A);
422 
423  pLS->MakeCompactForm(A, Ax, Ai, Adummy, Ap);
424 }
425 
426 template <typename MH>
427 void UmfpackSparseSolutionManager::ScaleMatrixAndRightHandSide(MH& mh)
428 {
429  if (scale.when != SCALEW_NEVER) {
430  MatrixScale<MH>& rMatScale = GetMatrixScale<MH>();
431 
432  if (pLS->bReset()) {
433  if (!rMatScale.bGetInitialized()
434  || scale.when == SolutionManager::SCALEW_ALWAYS) {
435  // (re)compute
436  rMatScale.ComputeScaleFactors(mh);
437  }
438  // in any case scale matrix and right-hand-side
439  rMatScale.ScaleMatrix(mh);
440 
441  if (silent_err) {
442  rMatScale.Report(std::cerr);
443  }
444  }
445 
446  rMatScale.ScaleRightHandSide(bVH);
447  }
448 }
449 
450 template <typename MH>
451 MatrixScale<MH>& UmfpackSparseSolutionManager::GetMatrixScale()
452 {
453  if (pMatScale == 0) {
454  pMatScale = MatrixScale<MH>::Allocate(scale);
455  }
456 
457  // Will throw std::bad_cast if the type does not match
458  return dynamic_cast<MatrixScale<MH>&>(*pMatScale);
459 }
460 
461 void UmfpackSparseSolutionManager::ScaleSolution(void)
462 {
463  if (scale.when != SCALEW_NEVER) {
464  ASSERT(pMatScale != 0);
465  // scale solution
466  pMatScale->ScaleSolution(xVH);
467  }
468 }
469 
470 /* Risolve il sistema Fattorizzazione + Backward Substitution */
471 void
472 UmfpackSparseSolutionManager::Solve(void)
473 {
474  MakeCompressedColumnForm();
475 
476  pLS->Solve();
477 
478  ScaleSolution();
479 }
480 
481 /* Risolve il sistema (trasposto) Fattorizzazione + Backward Substitution */
482 void
483 UmfpackSparseSolutionManager::SolveT(void)
484 {
485  MakeCompressedColumnForm();
486 
487  pLS->SolveT();
488 
489  ScaleSolution();
490 }
491 
492 /* Rende disponibile l'handler per la matrice */
494 UmfpackSparseSolutionManager::pMatHdl(void) const
495 {
496  return &A;
497 }
498 
499 /* Rende disponibile l'handler per il termine noto */
501 UmfpackSparseSolutionManager::pResHdl(void) const
502 {
503  return &bVH;
504 }
505 
506 /* Rende disponibile l'handler per la soluzione */
508 UmfpackSparseSolutionManager::pSolHdl(void) const
509 {
510  return &xVH;
511 }
512 
513 /* UmfpackSparseSolutionManager - end */
514 
515 template <class CC>
516 UmfpackSparseCCSolutionManager<CC>::UmfpackSparseCCSolutionManager(integer Dim,
517  doublereal dPivot,
518  doublereal dDropTolerance,
519  const unsigned& blockSize,
520  const ScaleOpt& scale,
521  integer iMaxIter)
522 : UmfpackSparseSolutionManager(Dim, dPivot, dDropTolerance, blockSize, scale, iMaxIter),
523 CCReady(false),
524 Ac(0)
525 {
526  NO_OP;
527 }
528 
529 template <class CC>
530 UmfpackSparseCCSolutionManager<CC>::~UmfpackSparseCCSolutionManager(void)
531 {
532  if (Ac) {
533  SAFEDELETE(Ac);
534  }
535 }
536 
537 template <class CC>
538 void
539 UmfpackSparseCCSolutionManager<CC>::MatrReset(void)
540 {
541  pLS->Reset();
542 }
543 
544 template <class CC>
545 void
546 UmfpackSparseCCSolutionManager<CC>::MakeCompressedColumnForm(void)
547 {
548  if (!CCReady) {
549  pLS->MakeCompactForm(A, Ax, Ai, Adummy, Ap);
550 
551  if (Ac == 0) {
552  SAFENEWWITHCONSTRUCTOR(Ac, CC, CC(Ax, Ai, Ap));
553  }
554 
555  CCReady = true;
556  }
557 
558  ScaleMatrixAndRightHandSide(*Ac);
559 }
560 
561 /* Inizializzatore "speciale" */
562 template <class CC>
563 void
564 UmfpackSparseCCSolutionManager<CC>::MatrInitialize()
565 {
566  CCReady = false;
567 
568  if (Ac) {
569  // If a DirCColMatrixHandler is in use and matrix scaling is enabled
570  // an uncaught exception (MatrixHandler::ErrRebuildMatrix) will be thrown
571  // if zero entries in the matrix become nonzero.
572  // For that reason we have to reinitialize Ac!
573  SAFEDELETE(Ac);
574  Ac = 0;
575  }
576 
577  MatrReset();
578 }
579 
580 /* Rende disponibile l'handler per la matrice */
581 template <class CC>
583 UmfpackSparseCCSolutionManager<CC>::pMatHdl(void) const
584 {
585  if (!CCReady) {
586  return &A;
587  }
588 
589  ASSERT(Ac != 0);
590  return Ac;
591 }
592 
593 template class UmfpackSparseCCSolutionManager<CColMatrixHandler<0> >;
594 template class UmfpackSparseCCSolutionManager<DirCColMatrixHandler<0> >;
595 
596 /* UmfpackSparseCCSolutionManager - end */
597 
598 #endif /* USE_UMFPACK */
599 
virtual integer MakeCompressedColumnForm(doublereal *const Ax, integer *const Ai, integer *const Ap, int offset=0) const =0
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
bool bGetInitialized() const
Definition: dgeequ.h:59
bool ComputeScaleFactors(const T &mh)
Definition: dgeequ.h:305
#define NO_OP
Definition: myassert.h:74
std::ostream & Report(std::ostream &os) const
Definition: dgeequ.cc:52
void Reset(scalar_func_type &d)
Definition: gradient.h:2836
Definition: mbdyn.h:76
VectorHandler & ScaleRightHandSide(VectorHandler &bVH) const
Definition: dgeequ.h:216
#define ASSERT(expression)
Definition: colamd.c:977
doublereal * pdSol
Definition: ls.h:75
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
static MatrixScale< T > * Allocate(const SolutionManager::ScaleOpt &scale)
Definition: dgeequ.h:313
doublereal * pdRhs
Definition: ls.h:74
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710
T & ScaleMatrix(T &mh) const
Definition: dgeequ.h:274