MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
superluwrap.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/superluwrap.cc,v 1.20 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  * SuperLU C++ WRAPPER *
35  * *
36  *****************************************************************************/
37 
38 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
39 
40 
41 #ifdef USE_SUPERLU
42 #ifndef USE_SUPERLU_MT /* SUPERLU and SUPERLU_MT are incompatible */
43 
44 #include <sys/types.h>
45 #include <unistd.h>
46 #include <signal.h>
47 
48 #include "spmh.h"
49 #include "spmapmh.h"
50 #include "dirccmh.h"
51 #include "ccmh.h"
52 
53 
54 #include "superluwrap.h"
55 
56 extern "C" {
57 #include <dsp_defs.h>
58 #include <util.h>
59 }
60 
61 struct SuperLUSolverData {
62  SuperMatrix A,
63  AC,
64  L,
65  U,
66  B;
67 
68  SuperLUStat_t Gstat;
69  std::vector<int> perm_c, /* column permutation vector */
70  perm_r, /* row permutations from partial pivoting */
71  etree ; /* elimination tree of A'*A */
72 
73  superlu_options_t options;
74 };
75 
76 
77 /* SuperLUSolver - begin */
78 
79 /* Costruttore: si limita ad allocare la memoria */
80 SuperLUSolver::SuperLUSolver(integer iMatOrd,
81  const doublereal &dPivot,
82  unsigned ptype)
83 : Aip(0),
84 App(0),
85 Axp(0),
86 iN(iMatOrd),
87 iNonZeroes(0),
88 dPivotFactor(dPivot),
89 permutation(ptype),
90 bFirstSol(true),
91 bRegenerateMatrix(true)
92 {
93  ASSERT(iN > 0);
94  ASSERTMSGBREAK(permutation & (SUPERLU_COLAMD|SUPERLU_MMDATA),
95  "SuperLU scalar solver unknown permutation strategy");
96 
97  SAFENEW(sld, SuperLUSolverData);
98 
99  /*
100  * This means it's the first run
101  * FIXME: create a dependence on the library's internals
102  */
103  sld->A.Store = NULL;
104 
105 }
106 
107 /* Distruttore */
108 SuperLUSolver::~SuperLUSolver(void)
109 {
110 }
111 
112 
113 #ifdef DEBUG
114 void
115 SuperLUSolver::IsValid(void) const
116 {
117  ASSERT(Aip != NULL);
118  ASSERT(App != NULL);
119  ASSERT(Axp != NULL);
120  ASSERT(iN > 0);
121 }
122 #endif /* DEBUG */
123 
124 /* Fattorizza la matrice */
125 void
126 SuperLUSolver::Factor(void)
127 {
128 #ifdef DEBUG
129  IsValid();
130 #endif /* DEBUG */
131 
132  ASSERT(iNonZeroes > 0);
133 
134  doublereal drop_tol = 0.0;
135  void *work = NULL;
136  int info = 0, lwork = 0;
137  int panel_size = sp_ienv(1),
138  relax = sp_ienv(2);
139 
140  if (bRegenerateMatrix) {
141  /* NOTE: we could use sld->A.Store == NULL */
142  if (bFirstSol) {
143  set_default_options(&sld->options);
144  sld->options.DiagPivotThresh = dPivotFactor;
145  //colperm_t permc_spec = MMD_AT_PLUS_A;
146  switch (permutation) {
147  case SUPERLU_MMDATA:
148  sld->options.ColPerm = MMD_ATA;
149  break;
150  case SUPERLU_COLAMD:
151  default:
152  sld->options.ColPerm = COLAMD;
153  break;
154  }
155  sld->options.Fact = DOFACT;
156  sld->options.PrintStat = NO;
157  ASSERT(Astore == NULL);
158 
159  /* ---------------------------------------------------
160  * Allocate storage and initialize statistics variables.
161  * ---------------------------------------------------*/
162  /* Set up the dense matrix data structure for B. */
163  dCreate_Dense_Matrix(&sld->B, iN, 1,
165  iN, SLU_DN, SLU_D, SLU_GE);
166 
167  /* Set up the sparse matrix data structure for A. */
168  dCreate_CompCol_Matrix(&sld->A, iN, iN, iNonZeroes,
169  Axp, Aip, App, SLU_NC, SLU_D, SLU_GE);
170 
171  StatInit(&sld->Gstat);
172 
173  sld->perm_c.resize(iN);
174  sld->perm_r.resize(iN);
175  sld->etree.resize(iN);
176 
177  bFirstSol = false; /* never change this again */
178 
179  } else {
180  sld->options.Fact = DOFACT;
181  NCformat *Astore = (NCformat *) sld->A.Store;
182 
183  ASSERT(Astore);
184 
185  Astore->nnz = iNonZeroes;
186  Astore->nzval = Axp;
187  Astore->rowind = Aip;
188  Astore->colptr = App;
189 
190  Destroy_CompCol_Permuted(&sld->AC);
191  Destroy_SuperNode_Matrix(&sld->L);
192  Destroy_CompCol_Matrix(&sld->U);
193  }
194 
195  int *pc = &(sld->perm_c[0]);
196  get_perm_c(sld->options.ColPerm, &sld->A, pc);
197 
198  bRegenerateMatrix = false;
199  } else {
200  sld->options.Fact = SamePattern;
201  Destroy_CompCol_Permuted(&sld->AC);
202  Destroy_SuperNode_Matrix(&sld->L);
203  Destroy_CompCol_Matrix(&sld->U);
204  }
205 
206 
207  int *pr = &(sld->perm_r[0]),
208  *pc = &(sld->perm_c[0]),
209  *et = &(sld->etree[0]);
210 
211 
212  sp_preorder(&sld->options, &sld->A, pc, et, &sld->AC);
213  dgstrf(&sld->options, &sld->AC, drop_tol, relax, panel_size, et, work, lwork, pc, pr,
214  &sld->L, &sld->U, &sld->Gstat, &info);
215 
216 }
217 
218 /* Risolve */
219 void
220 SuperLUSolver::Solve(void) const
221 {
222 #ifdef DEBUG
223  IsValid();
224 #endif /* DEBUG */
225 
226  if (bHasBeenReset) {
227  const_cast<SuperLUSolver *>(this)->Factor();
228  bHasBeenReset = false;
229  }
230 
231  /* ------------------------------------------------------------
232  * Solve the system A*X=B, overwriting B with X.
233  * ------------------------------------------------------------*/
234  trans_t trans = NOTRANS;
235  int info = 0;
236 
237  int *pr = &(sld->perm_r[0]),
238  *pc = &(sld->perm_c[0]);
239 
240  dgstrs(trans, &sld->L, &sld->U, pc, pr,
241  &sld->B, &sld->Gstat, &info);
242 }
243 
244 /* Index Form */
245 void
246 SuperLUSolver::MakeCompactForm(SparseMatrixHandler& mh,
247  std::vector<doublereal>& Ax,
248  std::vector<integer>& Ar, std::vector<integer>& Ac,
249  std::vector<integer>& Ap) const
250 {
251  /* no need to rebuild matrix */
252  if (!bHasBeenReset) {
253  return;
254  }
255 
256  iNonZeroes = mh.MakeCompressedColumnForm(Ax, Ar, Ap, 0);
257  ASSERT(iNonZeroes > 0);
258 
259  Axp = &Ax[0];
260  Aip = &Ar[0];
261  App = &Ap[0];
262 
263  /* rebuild matrix ... (CC is broken) */
264  bRegenerateMatrix = true;
265 
266 #if 0
267  Destroy_CompCol_Matrix(&sld->A);
268 #endif
269 }
270 
271 /* SuperLUSolver - end */
272 
273 
274 /* SuperLUSparseSolutionManager - begin: code */
275 
276 /* Costruttore */
277 SuperLUSparseSolutionManager::SuperLUSparseSolutionManager(integer iSize,
278  const doublereal& dPivotFactor, unsigned ptype)
279 : iMatSize(iSize),
280 Ap(iSize + 1),
281 xb(iSize),
282 MH(iSize),
283 VH(iSize, &xb[0])
284 {
285  ASSERT(iSize > 0);
286  ASSERT((dPivotFactor >= 0.0) && (dPivotFactor <= 1.0));
287 
289  SuperLUSolver,
290  SuperLUSolver(iMatSize, dPivotFactor, ptype));
291 
292  pLS->pdSetResVec(&(xb[0]));
293  pLS->pdSetSolVec(&(xb[0]));
294  pLS->SetSolutionManager(this);
295 
296 #ifdef DEBUG
297  IsValid();
298 #endif /* DEBUG */
299 }
300 
301 
302 /* Distruttore; verificare la distruzione degli oggetti piu' complessi */
303 SuperLUSparseSolutionManager::~SuperLUSparseSolutionManager(void)
304 {
305 #ifdef DEBUG
306  IsValid();
307 #endif /* DEBUG */
308 
309  /* Dealloca roba, tra cui i thread */
310 }
311 
312 #ifdef DEBUG
313 /* Test di validita' del manager */
314 void
315 SuperLUSparseSolutionManager::IsValid(void) const
316 {
317  ASSERT(iMatSize > 0);
318 
319 #ifdef DEBUG_MEMMANAGER
320  ASSERT(defaultMemoryManager.fIsPointerToBlock(VH.pdGetVec()));
321  ASSERT(defaultMemoryManager.fIsPointerToBlock(pLS));
322 #endif /* DEBUG_MEMMANAGER */
323 
324  ASSERT((VH.IsValid(), 1));
325  ASSERT((pLS->IsValid(), 1));
326 }
327 #endif /* DEBUG */
328 
329 /* Inizializza il gestore delle matrici */
330 void
331 SuperLUSparseSolutionManager::MatrReset(void)
332 {
333 #ifdef DEBUG
334  IsValid();
335 #endif /* DEBUG */
336 
337  pLS->Reset();
338 }
339 
340 void
341 SuperLUSparseSolutionManager::MakeCompressedColumnForm(void)
342 {
343 #ifdef DEBUG
344  IsValid();
345 #endif /* DEBUG */
346 
347  /* FIXME: move this to the matrix handler! */
348  pLS->MakeCompactForm(MH, Ax, Ai, Adummy, Ap);
349 }
350 
351 /* Risolve il problema */
352 void
353 SuperLUSparseSolutionManager::Solve(void)
354 {
355 #ifdef DEBUG
356  IsValid();
357 #endif /* DEBUG */
358 
359  /* FIXME: move this to the matrix handler! */
360  MakeCompressedColumnForm();
361 
362 #if 0
363  std::cerr << "### after MakeIndexForm:" << std::endl
364  << "{col Ap[col]}={" << std::endl;
365  for (unsigned i = 0; i < Ap.size(); i++) {
366  std::cerr << i << " " << Ap[i] << std::endl;
367  }
368  std::cerr << "}" << std::endl;
369 
370  std::cerr << "{idx Ai[idx] col Ax[idx]}={" << std::endl;
371  unsigned c = 0;
372  for (unsigned i = 0; i < Ax.size(); i++) {
373  std::cerr << i << " " << Ai[i] << " " << c << " " << Ax[i] << std::endl;
374  if (i == Ap[c]) {
375  c++;
376  }
377  }
378  std::cerr << "}" << std::endl;
379 #endif
380 
381  pLS->Solve();
382 
383 #if 0
384  std::cerr << "### after Solve:" << std::endl
385  << "{col Ap[col]}={" << std::endl;
386  for (unsigned i = 0; i < Ap.size(); i++) {
387  std::cerr << i << " " << Ap[i] << std::endl;
388  }
389  std::cerr << "}" << std::endl;
390 
391  std::cerr << "{idx Ai[idx] col Ax[idx]}={" << std::endl;
392  c = 0;
393  for (unsigned i = 0; i < Ax.size(); i++) {
394  std::cerr << i << " " << Ai[i] << " " << c << " " << Ax[i] << std::endl;
395  if (i == Ap[c]) {
396  c++;
397  }
398  }
399  std::cerr << "}" << std::endl;
400 #endif
401 }
402 
403 /* SuperLUSparseSolutionManager - end */
404 
405 /* SuperLUSparseCCSolutionManager - begin */
406 
407 template <class CC>
408 SuperLUSparseCCSolutionManager<CC>::SuperLUSparseCCSolutionManager(integer Dim,
409  const doublereal &dPivot, unsigned ptype)
410 : SuperLUSparseSolutionManager(Dim, dPivot, ptype),
411 CCReady(false),
412 Ac(0)
413 {
414  NO_OP;
415 }
416 
417 template <class CC>
418 SuperLUSparseCCSolutionManager<CC>::~SuperLUSparseCCSolutionManager(void)
419 {
420  if (Ac) {
421  SAFEDELETE(Ac);
422  }
423 }
424 
425 template <class CC>
426 void
427 SuperLUSparseCCSolutionManager<CC>::MatrReset(void)
428 {
429  pLS->Reset();
430 }
431 
432 /* Risolve il sistema Fattorizzazione + Backward Substitution */
433 template <class CC>
434 void
435 SuperLUSparseCCSolutionManager<CC>::MakeCompressedColumnForm(void)
436 {
437  if (!CCReady) {
438  pLS->MakeCompactForm(MH, Ax, Ai, Adummy, Ap);
439 
440  if (Ac == 0) {
441  SAFENEWWITHCONSTRUCTOR(Ac, CC, CC(Ax, Ai, Ap));
442  }
443 
444  CCReady = true;
445  }
446 }
447 
448 /* Inizializzatore "speciale" */
449 template <class CC>
450 void
451 SuperLUSparseCCSolutionManager<CC>::MatrInitialize()
452 {
453  CCReady = false;
454 
455  MatrReset();
456 }
457 
458 /* Rende disponibile l'handler per la matrice */
459 template <class CC>
461 SuperLUSparseCCSolutionManager<CC>::pMatHdl(void) const
462 {
463  if (!CCReady) {
464  return &MH;
465  }
466 
467  ASSERT(Ac != 0);
468  return Ac;
469 }
470 
471 template class SuperLUSparseCCSolutionManager<CColMatrixHandler<0> >;
472 template class SuperLUSparseCCSolutionManager<DirCColMatrixHandler<0> >;
473 
474 /* SuperLUSparseCCSolutionManager - end */
475 
476 #endif /* !USE_SUPERLU_MT */
477 #endif /* USE_SUPERLU */
#define ASSERTMSGBREAK(expr, msg)
Definition: myassert.h:222
virtual integer MakeCompressedColumnForm(doublereal *const Ax, integer *const Ai, integer *const Ap, int offset=0) const =0
#define NO_OP
Definition: myassert.h:74
#define SAFENEW(pnt, item)
Definition: mynewmem.h:695
Definition: mbdyn.h:76
Definition: mbdyn.h:77
#define ASSERT(expression)
Definition: colamd.c:977
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
#define defaultMemoryManager
Definition: mynewmem.h:691
struct option options[]
Definition: ann_in.c:46
static std::stack< cleanup * > c
Definition: cleanup.cc:59
doublereal * pdRhs
Definition: ls.h:74
LinearSolver * pLS
Definition: solman.h:151
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710