MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
y12wrap.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbwrap/y12wrap.cc,v 1.61 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  * Y12 C++ WRAPPER *
35  * *
36  *****************************************************************************/
37 
38 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
39 
40 #ifdef USE_Y12
41 
42 #include <cstring>
43 #include "spmh.h"
44 #include "spmapmh.h"
45 #include "dirccmh.h"
46 #include "ccmh.h"
47 #include "y12lib.h"
48 #include "y12wrap.h"
49 
50 /* Y12Solver - begin */
51 
52 /* Costruttore: si limita ad allocare la memoria */
53 Y12Solver::Y12Solver(integer iMatOrd, integer iWorkSpaceSize,
54  doublereal* pdTmpRhs,
55  integer iPivotParam,
56  bool bDupInd)
57 : LinearSolver(0),
58 iMaxSize(iWorkSpaceSize),
59 iCurSize(iWorkSpaceSize),
60 piRow(0),
61 piCol(0),
62 pdMat(0),
63 pir(0),
64 pic(0),
65 bDuplicateIndices(bDupInd),
66 iN(iMatOrd),
67 iNonZeroes(0),
68 piHA(NULL),
69 pdPIVOT(NULL)
70 {
71  (void)pdSetResVec(pdTmpRhs);
72  (void)pdSetSolVec(pdTmpRhs);
73 
74  ASSERT(pdTmpRhs != NULL);
75  ASSERT(iN > 0);
76 
77  if (bDuplicateIndices) {
78  /*
79  * NOTE: Y12 alters the index arrays :(
80  */
81  iRow.reserve(iCurSize);
82  iCol.reserve(iCurSize);
83  }
84 
85  SAFENEWARR(piHA, integer, 11*iN);
86  SAFENEWARR(pdPIVOT, doublereal, iN);
87 
88 #ifdef DEBUG
89  for (integer iCnt = 0; iCnt < 11*iN; iCnt++) {
90  piHA[iCnt] = 0;
91  }
92 
93  for (integer iCnt = 0; iCnt < iN; iCnt++) {
94  pdPIVOT[iCnt] = 0;
95  }
96 #endif /* DEBUG */
97 
98  iIFLAG[I_1] = 0;
99  iIFLAG[I_2] = 3; /* recommended row number for pivoting */
100  iIFLAG[I_3] = iPivotParam;
101  iIFLAG[I_4] = 0;
102  iIFLAG[I_5] = 2; /* store non-zero elements of L */
103 
104  dAFLAG[I_1] = 8.; /* Should be 4.<dAFLAG[0]<16. for stability */
105  dAFLAG[I_2] = 0.; /* Should be 0.<dAFLAG[1]<1.e-12 */
106  dAFLAG[I_3] = 1.e6; /* Should be dAFLAG[2]>1.e5 */
107  dAFLAG[I_4] = 0.; /* FIXME: Should be 0 < dAFLAG[3]<1.e-12 */
108 }
109 
110 /* Distruttore */
111 Y12Solver::~Y12Solver(void)
112 {
113  if (pdPIVOT != NULL) {
114  SAFEDELETEARR(pdPIVOT);
115  }
116  if (piHA != NULL) {
117  SAFEDELETEARR(piHA);
118  }
119 }
120 
121 #ifdef DEBUG
122 void
123 Y12Solver::IsValid(void) const
124 {
125  ASSERT(iCurSize > 0);
126  ASSERT(piRow != NULL);
127  ASSERT(piCol != NULL);
128  ASSERT(pdMat != NULL);
129  ASSERT(iN > 0);
130 
131  ASSERT(piHA != NULL);
132  ASSERT(pdPIVOT != NULL);
133 
134 #ifdef DEBUG_MEMMANAGER
135  ASSERT(defaultMemoryManager.fIsBlock(piHA, 11*iN*sizeof(integer)));
136  ASSERT(defaultMemoryManager.fIsBlock(pdPIVOT, 1*iN*sizeof(doublereal)));
137 #endif /* DEBUG_MEMMANAGER */
138 }
139 #endif /* DEBUG */
140 
141 /* Fattorizza la matrice */
142 void
143 Y12Solver::Factor(void)
144 {
145 #ifdef DEBUG
146  IsValid();
147 #endif /* DEBUG */
148 
149  ASSERT(iNonZeroes > 0);
150 
151  /* Sets parameters */
152  integer iIFAIL = 0;
153 
154  /*
155  * must be set to 0 before the first call to a routine
156  * of the y12m package
157  */
158  iIFLAG[I_1] = 0;
159  iIFLAG[I_5] = 2;
160 
161  if (bDuplicateIndices) {
162  /*
163  * NOTE: Y12 alters the index arrays :(
164  *
165  * FIXME: make it stl-ish
166  */
167  iRow.resize(iCurSize);
168  iCol.resize(iCurSize);
169 
170  pir = &iRow[0];
171  pic = &iCol[0];
172 
173 #ifdef HAVE_MEMMOVE
174  memmove(pir, piRow, sizeof(integer)*iNonZeroes);
175  memmove(pic, piCol, sizeof(integer)*iNonZeroes);
176 #else /* ! HAVE_MEMMOVE */
177  for (unsigned i = 0; i < iNonZeroes; i++) {
178  pir[i] = piRow[i];
179  pic[i] = piCol[i];
180  }
181 #endif /* ! HAVE_MEMMOVE */
182 
183  } else {
184  pir = piRow;
185  pic = piCol;
186  }
187 
188  y12prefactor(&iN, &iNonZeroes, pdMat,
189  pic, &iCurSize,
190  pir, &iCurSize,
191  piHA, &iN,
192  dAFLAG, iIFLAG, &iIFAIL);
193 
194  if (iIFAIL != 0) {
195  silent_cerr("Y12Solver (y12prefactor): "
196  "error during pre-factorization, code "
197  << iIFAIL << ":" << std::endl);
198  PutError(iIFAIL);
199  throw Y12Solver::ErrFactorization(iIFAIL, MBDYN_EXCEPT_ARGS);
200  }
201 
202  /* actual factorization */
203  y12factor(&iN, &iNonZeroes, pdMat,
204  pic, &iCurSize,
205  pir, &iCurSize,
206  pdPIVOT, LinearSolver::pdRhs,
207  piHA, &iN,
208  dAFLAG, iIFLAG, &iIFAIL);
209 
210  if (iIFAIL != 0) {
211  silent_cerr("Y12Solver (y12factor): "
212  "error during factorization, code "
213  << iIFAIL << ":" << std::endl);
214  PutError(iIFAIL);
215  throw Y12Solver::ErrFactorization(iIFAIL, MBDYN_EXCEPT_ARGS);
216  }
217 
218  if (dAFLAG[7] < 1.e-12) {
219  silent_cerr("Y12Solver (y12factor):"
220  " warning, possible bad conditioning of matrix"
221  << std::endl);
222  }
223 }
224 
225 /* Risolve */
226 void
227 Y12Solver::Solve(void) const
228 {
229 #ifdef DEBUG
230  IsValid();
231 #endif /* DEBUG */
232 
233  if (bHasBeenReset) {
234  const_cast<Y12Solver *>(this)->Factor();
235  }
236 
237  integer iIFAIL = 0;
238 
239  y12solve(&iN, pdMat, &iCurSize, LinearSolver::pdRhs,
240  pdPIVOT, pic,
241  piHA, &iN,
242  iIFLAG, &iIFAIL);
243 
244  if (iIFAIL != 0) {
245  silent_cerr("Y12Solver (y12solve): "
246  "error during back substitution, code "
247  << iIFAIL << ":" << std::endl);
248  PutError(iIFAIL);
249  throw Y12Solver::ErrFactorization(iIFAIL, MBDYN_EXCEPT_ARGS);
250  }
251 
252  if (bHasBeenReset) {
253  iIFLAG[I_5] = 3;
254  bHasBeenReset = false;
255  }
256 }
257 
258 /* Index Form */
259 void
260 Y12Solver::MakeCompactForm(SparseMatrixHandler& mh,
261  std::vector<doublereal>& Ax,
262  std::vector<integer>& Ar, std::vector<integer>& Ac,
263  std::vector<integer>& Ap) const
264 {
265  if (!bHasBeenReset) {
266  return;
267  }
268 
269  iNonZeroes = mh.MakeIndexForm(Ax, Ar, Ac, Ap, 1);
270  ASSERT(iNonZeroes > 0);
271 
272  pdMat = &Ax[0];
273  piRow = &Ar[0];
274  piCol = &Ac[0];
275 
276  /* iCurSize should be between 3 and 5 times iNonZeroes ... */
277  if (iCurSize > 5*iNonZeroes) {
278  iCurSize = 5*iNonZeroes;
279 
280  } else if (iCurSize < 3*iNonZeroes) {
281  if (iMaxSize < 5*iNonZeroes) {
282  iCurSize = iMaxSize;
283  } else {
284  iCurSize = 5*iNonZeroes;
285  }
286  }
287 }
288 
289 void
290 Y12Solver::PutError(integer rc) const
291 {
292  silent_cerr(std::endl);
293 
294  switch (rc) {
295  case 1:
296  silent_cerr("\tThe coefficient matrix A is not" << std::endl
297  << "\tfactorized, i.e. the call of subroutine" << std::endl
298  << "\tY12MD was not preceded by a call of Y12MC" << std::endl
299  << "\tduring the solution of Ax=b or during" << std::endl
300  << "\tthe solution of the first system in a" << std::endl
301  << "\tsequence ( Ax1 = b1 , Ax2 = b2,.....,Axp =" << std::endl
302  << "\tbp) of systems with the same coefficient" << std::endl
303  << "\tmatrix. This will work in all cases only" << std::endl
304  << "\tif the user sets IFLAG(1) .ge. 0 before" << std::endl
305  << "\tthe call of package Y12M (i.e. before the" << std::endl
306  << "\tfirst call of a subroutine of this" << std::endl
307  << "\tpackage)." << std::endl);
308  break;
309 
310  case 2:
311  silent_cerr("\tThe coefficient matrix A is not ordered," << std::endl
312  << "\ti.e. the call of subroutine Y12MC was not" << std::endl
313  << "\tpreceded by a call of Y12MB. This will" << std::endl
314  << "\twork in all cases only if the user sets" << std::endl
315  << "\tIFLAG(1) .ge. 0 before the call of package" << std::endl
316  << "\tY12M (i.e. before the first call of a" << std::endl
317  << "\tsubroutine of this package)." << std::endl);
318  break;
319 
320  case 3:
321  silent_cerr("\tA pivotal element abs(a(i,i;j)) < AFLAG(4)" << std::endl
322  << "\t* AFLAG(6) is selected. When AFLAG(4) is" << std::endl
323  << "\tsufficiently small this is an indication" << std::endl
324  << "\tthat the coefficient matrix is numerically" << std::endl
325  << "\tsingular." << std::endl);
326  break;
327 
328  case 4:
329  silent_cerr("\tAFLAG(5), the growth factor, is larger" << std::endl
330  << "\tthan AFLAG(3). When AFLAG(3) is" << std::endl
331  << "\tsufficiently large this indicates that the" << std::endl
332  << "\telements of the coefficient matrix A grow" << std::endl
333  << "\tso quickly during the factorization that" << std::endl
334  << "\tthe continuation of the computation is not" << std::endl
335  << "\tjustified. The choice of a smaller" << std::endl
336  << "\tstability factor, AFLAG(1), may give" << std::endl
337  << "\tbetter results in this case." << std::endl);
338  break;
339 
340  case 5:
341  silent_cerr("\tThe length NN of arrays A and SNR is not" << std::endl
342  << "\tsufficient. Larger values of NN (and" << std::endl
343  << "\tpossibly of NN1) should be used." << std::endl);
344  break;
345 
346  case 6:
347  silent_cerr("\tThe length NN1 of array RNR is not" << std::endl
348  << "\tsufficient. Larger values of NN1 (and" << std::endl
349  << "\tpossibly of NN) should be used." << std::endl);
350  break;
351 
352  case 7:
353  silent_cerr("\tA row without non-zero elements in its" << std::endl
354  << "\tactive part is found during the" << std::endl
355  << "\tdecomposition. If the drop-tolerance," << std::endl
356  << "\tAFLAG(2), is sufficiently small, then" << std::endl
357  << "\tIFAIL = 7 indicates that the matrix is" << std::endl
358  << "\tnumerically singular. If a large value of" << std::endl
359  << "\tthe drop-tolerance AFLAG(2) is used and if" << std::endl
360  << "\tIFAIL = 7 on exit, this is not certain. A" << std::endl
361  << "\trun with a smaller value of AFLAG(2)" << std::endl
362  << "\tand/or a careful check of the parameters" << std::endl
363  << "\tAFLAG(8) and AFLAG(5) is recommended in" << std::endl
364  << "\tthe latter case." << std::endl);
365  break;
366 
367  case 8:
368  silent_cerr("\tA column without non-zero elements in its" << std::endl
369  << "\tactive part is found during the" << std::endl
370  << "\tdecomposition. If the drop-tolerance," << std::endl
371  << "\tAFLAG(2), is sufficiently small, then" << std::endl
372  << "\tIFAIL = 8 indicates that the matrix is" << std::endl
373  << "\tnumerically singular. If a large value of" << std::endl
374  << "\tthe drop-tolerance AFLAG(2) is used and if" << std::endl
375  << "\tIFAIL = 8 on exit, this is not certain. A" << std::endl
376  << "\trun with a smaller value of AFLAG(2)" << std::endl
377  << "\tand/or a careful check of the parameters" << std::endl
378  << "\tAFLAG(8) and AFLAG(5) is recommended in" << std::endl
379  << "\tthe latter case." << std::endl);
380  break;
381 
382  case 9:
383  silent_cerr("\tA pivotal element is missing. This may" << std::endl
384  << "\toccur if AFLAG(2) > 0 and IFLAG(4) = 2" << std::endl
385  << "\t(i.e. some system after the first one in a" << std::endl
386  << "\tsequence of systems with the same" << std::endl
387  << "\tstructure is solved using a positive value" << std::endl
388  << "\tfor the drop-tolerance). The value of the" << std::endl
389  << "\tdrop-tolerance AFLAG(2), should be" << std::endl
390  << "\tdecreased and the coefficient matrix of" << std::endl
391  << "\tthe system refactorized. This error may" << std::endl
392  << "\talso occur when one of the special pivotal" << std::endl
393  << "\tstrategies (IFLAG(3)=0 or IFLAG(3)=2) is" << std::endl
394  << "\tused and the matrix is not suitable for" << std::endl
395  << "\tsuch a strategy." << std::endl);
396  break;
397 
398  case 10:
399  silent_cerr("\tSubroutine Y12MF is called with IFLAG(5) =" << std::endl
400  << "\t1 (i.e. with a request to remove the" << std::endl
401  << "\tnon-zero elements of the lower triangular" << std::endl
402  << "\tmatrix L). IFLAG(5)=2 must be" << std::endl
403  << "\tinitialized instead of IFLAG(5)=1." << std::endl);
404  break;
405 
406  case 11:
407  silent_cerr("\tThe coefficient matrix A contains at least" << std::endl
408  << "\ttwo elements in the same position (i,j)." << std::endl
409  << "\tThe input data should be examined" << std::endl
410  << "\tcarefully in this case." << std::endl);
411  break;
412 
413  case 12:
414  silent_cerr("\tThe number of equations in the system Ax=b" << std::endl
415  << "\tis smaller than 2 (i.e. N<2). The value" << std::endl
416  << "\tof N should be checked." << std::endl);
417  break;
418 
419  case 13:
420  silent_cerr("\tThe number of non-zero elements of the" << std::endl
421  << "\tcoefficient matrix is non-positive (i.e." << std::endl
422  << "\tZ.le.0 ). The value of the parameter Z" << std::endl
423  << "\t(renamed NZ in Y12MF) should be checked." << std::endl);
424  break;
425 
426  case 14:
427  silent_cerr("\tThe number of non-zero elements in the" << std::endl
428  << "\tcoefficient matrix is smaller than the" << std::endl
429  << "\tnumber of equations (i.e. Z < N ). If" << std::endl
430  << "\tthere is no mistake (i.e. if parameter Z," << std::endl
431  << "\trenamed NZ in Y12MF, is correctly assigned" << std::endl
432  << "\ton entry) then the coefficient matrix is" << std::endl
433  << "\tstructurally singular in this case." << std::endl);
434  break;
435 
436  case 15:
437  silent_cerr("\tThe length IHA of the first dimension of" << std::endl
438  << "\tarray HA is smaller than N. IHA.ge.N" << std::endl
439  << "\tshould be assigned." << std::endl);
440  break;
441 
442  case 16:
443  silent_cerr("\tThe value of parameter IFLAG(4) is not" << std::endl
444  << "\tassigned correctly. IFLAG(4) should be" << std::endl
445  << "\tequal to 0, 1 or 2. See more details in" << std::endl
446  << "\tthe description of this parameter." << std::endl);
447  break;
448 
449  case 17:
450  silent_cerr("\tA row without non-zero elements has been" << std::endl
451  << "\tfound in the coefficient matrix A of the" << std::endl
452  << "\tsystem before the Gaussian elimination is" << std::endl
453  << "\tinitiated. Matrix A is structurally" << std::endl
454  << "\tsingular." << std::endl);
455  break;
456 
457  case 18:
458  silent_cerr("\tA column without non-zero elements has" << std::endl
459  << "\tbeen found in the coefficient matrix A of" << std::endl
460  << "\tthe system before the Gaussian elimination" << std::endl
461  << "\tis initiated. Matrix A is structurally" << std::endl
462  << "\tsingular." << std::endl);
463  break;
464 
465  case 19:
466  silent_cerr("\tParameter IFLAG(2) is smaller than 1. The" << std::endl
467  << "\tvalue of IFLAG(2) should be a positive" << std::endl
468  << "\tinteger (IFLAG(2) = 3 is recommended)." << std::endl);
469  break;
470 
471  case 20:
472  silent_cerr("\tParameter IFLAG(3) is out of range." << std::endl
473  << "\tIFLAG(3) should be equal to 0, 1 or 2." << std::endl);
474  break;
475 
476  case 21:
477  silent_cerr("\tParameter IFLAG(5) is out of range." << std::endl
478  << "\tIFLAG(5) should be equal to 1, 2 or 3 (but" << std::endl
479  << "\twhen IFLAG(5) = 3 Y12MB and Y12MC should" << std::endl
480  << "\tnot be called; see also the message for" << std::endl
481  << "\tIFAIL = 22 below)." << std::endl);
482  break;
483 
484  case 22:
485  silent_cerr("\tEither subroutine Y12MB or subroutine" << std::endl
486  << "\tY12MC is called with IFLAG(5) = 3. Each of" << std::endl
487  << "\tthese subroutines should be called with" << std::endl
488  << "\tIFLAG(5) equal to 1 or 2." << std::endl);
489  break;
490 
491  case 23:
492  silent_cerr("\tThe number of allowed iterations" << std::endl
493  << "\t(parameter IFLAG(11) when Y12MF is used)" << std::endl
494  << "\tis smaller than 2. IFLAG(11) .ge. 2" << std::endl
495  << "\tshould be assigned." << std::endl);
496  break;
497 
498  case 24:
499  silent_cerr("\tAt least one element whose column number" << std::endl
500  << "\tis either larger than N or smaller than 1" << std::endl
501  << "\tis found." << std::endl);
502  break;
503 
504  case 25:
505  silent_cerr("\tAt least one element whose row number is" << std::endl
506  << "\teither larger than N or smaller than 1 is" << std::endl
507  << "\tfound." << std::endl);
508  break;
509 
510  default:
511  silent_cerr("\t Unhandled code." << std::endl);
512  break;
513  }
514 
515  silent_cerr(std::endl);
516 }
517 
518 /* Y12Solver - end */
519 
520 
521 /* Y12SparseSolutionManager - begin: code */
522 
523 /* Costruttore */
524 Y12SparseSolutionManager::Y12SparseSolutionManager(integer iSize,
525  integer iWorkSpaceSize,
526  const doublereal& dPivotFactor, bool bDupInd)
527 : iMatSize(iSize),
528 iColStart(iSize + 1),
529 dVec(iSize),
530 MH(iSize),
531 VH(iSize, &dVec[0])
532 {
533  ASSERT(iSize > 0);
534  ASSERT(((dPivotFactor >= 0.0) && (dPivotFactor <= 1.0)) || dPivotFactor == -1.0);
535 
536 
537  /* Valore di default */
538  if (iWorkSpaceSize == 0) {
539  /*
540  * y12 requires at least 3*numzeros to store factors
541  * for multiple backsubs
542  */
543  iWorkSpaceSize = 3*iSize*iSize;
544  }
545 
546  integer iPivot;
547  if (dPivotFactor == 0.) {
548  iPivot = 0;
549 
550  } else {
551  iPivot = 1;
552  }
553 
554  iRow.reserve(iWorkSpaceSize);
555  iCol.reserve(iWorkSpaceSize);
556  dMat.reserve(iWorkSpaceSize);
557 
558  SAFENEWWITHCONSTRUCTOR(SolutionManager::pLS,
559  Y12Solver,
560  Y12Solver(iMatSize, iWorkSpaceSize,
561  &dVec[0], iPivot, bDupInd));
562 
563  pLS->SetSolutionManager(this);
564 
565 #ifdef DEBUG
566  IsValid();
567 #endif /* DEBUG */
568 }
569 
570 
571 /* Distruttore; verificare la distruzione degli oggetti piu' complessi */
572 Y12SparseSolutionManager::~Y12SparseSolutionManager(void)
573 {
574 #ifdef DEBUG
575  IsValid();
576 #endif /* DEBUG */
577 
578  /* Dealloca arrays */
579 }
580 
581 #ifdef DEBUG
582 /* Test di validita' del manager */
583 void
584 Y12SparseSolutionManager::IsValid(void) const
585 {
586  ASSERT(iMatSize > 0);
587 
588 #ifdef DEBUG_MEMMANAGER
589  ASSERT(defaultMemoryManager.fIsPointerToBlock(pLS));
590 #endif /* DEBUG_MEMMANAGER */
591  pLS->IsValid();
592 }
593 #endif /* DEBUG */
594 
595 void
596 Y12SparseSolutionManager::MatrReset(void)
597 {
598  pLS->Reset();
599 }
600 
601 void
602 Y12SparseSolutionManager::MakeIndexForm(void)
603 {
604 #ifdef DEBUG
605  IsValid();
606 #endif /* DEBUG */
607 
608  /* FIXME: move this to the matrix handler! */
609  pLS->MakeCompactForm(MH, dMat, iRow, iCol, iColStart);
610 }
611 
612 /* Risolve il problema */
613 void
614 Y12SparseSolutionManager::Solve(void)
615 {
616 #ifdef DEBUG
617  IsValid();
618 #endif /* DEBUG */
619 
620  /* FIXME: move this to the matrix handler! */
621  MakeIndexForm();
622 
623 #if 0
624  std::cerr << "### after MakeIndexForm:" << std::endl
625  << "{col Ap[col]}={" << std::endl;
626  for (unsigned i = 0; i < iColStart.size(); i++) {
627  std::cerr << i << " " << iColStart[i] << std::endl;
628  }
629  std::cerr << "}" << std::endl;
630 
631  std::cerr << "{idx Ai[idx] col Ax[idx]}={" << std::endl;
632  unsigned c = 0;
633  for (unsigned i = 0; i < dMat.size(); i++) {
634  std::cerr << i << " " << iRow[i] << " " << c << " " << dMat[i] << std::endl;
635  if (i == iColStart[c]) {
636  c++;
637  }
638  }
639  std::cerr << "}" << std::endl;
640 #endif
641 
642  pLS->Solve();
643 
644 #if 0
645  std::cerr << "### after Solve:" << std::endl
646  << "{col Ap[col]}={" << std::endl;
647  for (unsigned i = 0; i < iColStart.size(); i++) {
648  std::cerr << i << " " << iColStart[i] << std::endl;
649  }
650  std::cerr << "}" << std::endl;
651 
652  std::cerr << "{idx Ai[idx] col Ax[idx]}={" << std::endl;
653  c = 0;
654  for (unsigned i = 0; i < dMat.size(); i++) {
655  std::cerr << i << " " << iRow[i] << " " << c << " " << dMat[i] << std::endl;
656  if (i == iColStart[c]) {
657  c++;
658  }
659  }
660  std::cerr << "}" << std::endl;
661 #endif
662 }
663 
664 /* Y12SparseSolutionManager - end */
665 
666 /* Y12SparseCCSolutionManager - begin */
667 
668 template <class CC>
669 Y12SparseCCSolutionManager<CC>::Y12SparseCCSolutionManager(integer Dim,
670  integer dummy, doublereal dPivot)
671 : Y12SparseSolutionManager(Dim, dummy, dPivot, true),
672 CCReady(false),
673 Ac(0)
674 {
675  NO_OP;
676 }
677 
678 template <class CC>
679 Y12SparseCCSolutionManager<CC>::~Y12SparseCCSolutionManager(void)
680 {
681  if (Ac) {
682  SAFEDELETE(Ac);
683  }
684 }
685 
686 template <class CC>
687 void
688 Y12SparseCCSolutionManager<CC>::MatrReset(void)
689 {
690  pLS->Reset();
691 }
692 
693 /* Risolve il sistema Fattorizzazione + Backward Substitution */
694 template <class CC>
695 void
696 Y12SparseCCSolutionManager<CC>::MakeIndexForm(void)
697 {
698  if (!CCReady) {
699  pLS->MakeCompactForm(MH, dMat, iRow, iCol, iColStart);
700 
701  if (Ac == 0) {
702  SAFENEWWITHCONSTRUCTOR(Ac, CC,
703  CC(dMat, iRow, iColStart));
704  }
705 
706  CCReady = true;
707  }
708 }
709 
710 /* Inizializzatore "speciale" */
711 template <class CC>
712 void
713 Y12SparseCCSolutionManager<CC>::MatrInitialize(void)
714 {
715  CCReady = false;
716 
717  MatrReset();
718 }
719 
720 /* Rende disponibile l'handler per la matrice */
721 template <class CC>
722 MatrixHandler*
723 Y12SparseCCSolutionManager<CC>::pMatHdl(void) const
724 {
725  if (!CCReady) {
726  return &MH;
727  }
728 
729  ASSERT(Ac != 0);
730  return Ac;
731 }
732 
733 template class Y12SparseCCSolutionManager<CColMatrixHandler<1> >;
734 template class Y12SparseCCSolutionManager<DirCColMatrixHandler<1> >;
735 
736 /* Y12SparseCCSolutionManager - end */
737 
738 #endif /* USE_Y12 */
739 
SparseMatrixHandler::MakeIndexForm
virtual integer MakeIndexForm(doublereal *const Ax, integer *const Arow, integer *const Acol, integer *const AcolSt, int offset=0) const =0
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
SAFEDELETEARR
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
y12solve
#define y12solve
Definition: y12lib.h:44
NO_OP
#define NO_OP
Definition: myassert.h:74
y12prefactor
#define y12prefactor
Definition: y12lib.h:42
false
Definition: mbdyn.h:76
true
Definition: mbdyn.h:77
y12factor
#define y12factor
Definition: y12lib.h:43
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
SAFENEWWITHCONSTRUCTOR
#define SAFENEWWITHCONSTRUCTOR(pnt, item, constructor)
Definition: mynewmem.h:698
defaultMemoryManager
#define defaultMemoryManager
Definition: mynewmem.h:691
c
static std::stack< cleanup * > c
Definition: cleanup.cc:59
LinearSolver::pdRhs
doublereal * pdRhs
Definition: ls.h:74
SAFENEWARR
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
SolutionManager::pLS
LinearSolver * pLS
Definition: solman.h:151
doublereal
double doublereal
Definition: colamd.c:52
integer
long int integer
Definition: colamd.c:51
SAFEDELETE
#define SAFEDELETE(pnt)
Definition: mynewmem.h:710