Orfeo Toolbox  3.16
itkFEMLinearSystemWrapperItpack.cxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkFEMLinearSystemWrapperItpack.cxx,v $
5  Language: C++
6  Date: $Date: 2009-02-19 18:40:54 $
7  Version: $Revision: 1.23 $
8 
9  Copyright (c) Insight Software Consortium. All rights reserved.
10  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 #include "itkNumericTraits.h"
18 #include "itkFEMLinearSystemWrapperItpack.h"
19 #include "itkFEMException.h"
20 #include "itpack.h"
21 
22 
23 namespace itk {
24 namespace fem {
25 
29 LinearSystemWrapperItpack::LinearSystemWrapperItpack()
30 {
31 
32  /* fill m_Methods with pointers to solver functions */
33  m_Methods[0] = jcg_;
34  m_Methods[1] = jsi_;
35  m_Methods[2] = sor_;
36  m_Methods[3] = ssorcg_;
37  m_Methods[4] = ssorsi_;
38  m_Methods[5] = rscg_;
39  m_Methods[6] = rssi_;
40  m_Method = 0; /* set default method to jcg_ */
41 
42  /* Set solving parameters */
43  dfault_( &(m_IPARM[0]) , &(m_RPARM[0]) );
44 
45  // We don't want the solver routines to
46  // overwrite the parameters.
47  m_IPARM[2]=1;
48 
49 
50  /* m_IPARM[0] = 500; */ /* number of iterations */
51  m_IPARM[1] = -1; /* no error message output */
52  m_IPARM[4] = 1; /* non-symmetric matrix */
53 
54  /* itpack recommended (but not default) value */
55 #undef min
56  m_RPARM[7] = 500.0 * NumericTraits<double>::min();
57 
58  m_MaximumNonZeroValues = 0;
59  m_Matrices = 0;
60  m_Vectors = 0;
61  m_Solutions = 0;
62 
63 }
64 
65 
66 void LinearSystemWrapperItpack::InitializeMatrix(unsigned int matrixIndex)
67 {
68  /* error checking */
69  if (!m_Order)
70  {
71  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeMatrix", "System order not set");
72  }
73  if (matrixIndex >= m_NumberOfMatrices)
74  {
75  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeMatrix", "m_Matrices", matrixIndex);
76  }
77  if (!m_MaximumNonZeroValues)
78  {
79  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeMatrix", "Maximum number of non zeros not set");
80  }
81 
82  // allocate if necessay
83  if (m_Matrices == 0)
84  {
85  m_Matrices = new MatrixHolder(m_NumberOfMatrices);
86  }
87 
88  /* Set required variables */
89  (*m_Matrices)[matrixIndex].Clear();
90  (*m_Matrices)[matrixIndex].SetOrder(m_Order);
91  (*m_Matrices)[matrixIndex].SetMaxNonZeroValues( m_MaximumNonZeroValues );
92 
93  return;
94 
95 }
96 
97 
98 bool LinearSystemWrapperItpack::IsMatrixInitialized(unsigned int matrixIndex)
99 {
100  if (!m_Matrices) return false;
101  if ( !(*m_Matrices)[matrixIndex].GetOrder() ) return false;
102  if ( !(*m_Matrices)[matrixIndex].GetMaxNonZeroValues() ) return false;
103 
104 
105  return true;
106 }
107 
108 void LinearSystemWrapperItpack::InitializeVector(unsigned int vectorIndex)
109 {
110 
111  /* error checking */
112  if (!m_Order)
113  {
114  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeVector", "System order not set");
115  }
116  if (vectorIndex >= m_NumberOfVectors)
117  {
118  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeVector", "m_Vectors", vectorIndex);
119  }
120 
121  /* allocate if necessay */
122  if (m_Vectors == 0)
123  {
124  m_Vectors = new VectorHolder(m_NumberOfVectors);
125  }
126 
127  /* delete old vector */
128  if ( (*m_Vectors)[vectorIndex] != 0 )
129  {
130  delete [] (*m_Vectors)[vectorIndex];
131  }
132 
133  /* insert new vector */
134  (*m_Vectors)[vectorIndex] = new doublereal [m_Order];
135 
136  /* fill with zeros */
137  for (unsigned int i=0; i<m_Order; i++)
138  {
139  (*m_Vectors)[vectorIndex][i] = 0.0;
140  }
141 
142  return;
143 
144 }
145 
146 
147 bool LinearSystemWrapperItpack::IsVectorInitialized(unsigned int vectorIndex)
148 {
149  if (!m_Vectors) return false;
150  if ( !(*m_Vectors)[vectorIndex] ) return false;
151 
152  return true;
153 }
154 
155 
156 void LinearSystemWrapperItpack::InitializeSolution(unsigned int solutionIndex)
157 {
158 
159  /* FIX ME: exceptions */
160  if (!m_Order)
161  {
162  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeSolution", "System order not set");
163  }
164  if (solutionIndex >= m_NumberOfSolutions)
165  {
166  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeSolution", "m_Solutions", solutionIndex);
167  }
168 
169 
170  // allocate if necessay
171  if (m_Solutions == 0)
172  {
173  m_Solutions = new VectorHolder(m_NumberOfSolutions);
174  }
175 
176  /* delete old vector */
177  if ( (*m_Solutions)[solutionIndex] != 0 )
178  {
179  delete [] (*m_Solutions)[solutionIndex];
180  }
181 
182  /* insert new vector */
183  (*m_Solutions)[solutionIndex] = new doublereal [m_Order];
184 
185  /* fill with zeros */
186  for (unsigned int i=0; i<m_Order; i++)
187  {
188  (*m_Solutions)[solutionIndex][i] = 0.0;
189  }
190 
191  return;
192 
193 }
194 
195 
196 bool LinearSystemWrapperItpack::IsSolutionInitialized(unsigned int solutionIndex)
197 {
198  if (!m_Solutions) return false;
199  if ( !(*m_Solutions)[solutionIndex] ) return false;
200 
201  return true;
202 }
203 
204 void LinearSystemWrapperItpack::DestroyMatrix(unsigned int matrixIndex)
205 {
206 
207  /* FIX ME: exceptions */
208  if ( !m_Matrices ) return;
209  if ( matrixIndex >= m_NumberOfMatrices )
210  {
211  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::DestroyMatrix", "m_Matrices", matrixIndex);
212  }
213 
214  (*m_Matrices)[matrixIndex].Clear();
215 }
216 
217 
218 void LinearSystemWrapperItpack::DestroyVector(unsigned int vectorIndex)
219 {
220  /* FIXME: exceptions */
221  if (!m_Vectors) return;
222  if (vectorIndex >= m_NumberOfVectors)
223  {
224  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::DestroyVector", "m_Vectors", vectorIndex);
225  }
226 
227  if ( !(*m_Vectors)[vectorIndex] ) return;
228 
229  /* delete vector */
230  delete [] (*m_Vectors)[vectorIndex];
231  (*m_Vectors)[vectorIndex] = 0;
232 
233 }
234 
235 
236 void LinearSystemWrapperItpack::DestroySolution(unsigned int solutionIndex)
237 {
238  // FIXME: exceptions
239  if (!m_Solutions) return;
240  if (solutionIndex >= m_NumberOfSolutions)
241  {
242  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::DestroySolution", "m_Solutions", solutionIndex);
243  }
244  if ( !(*m_Solutions)[solutionIndex] ) return;
245 
246  /* delete vector */
247  delete [] (*m_Solutions)[solutionIndex];
248  (*m_Solutions)[solutionIndex] = 0;
249 
250 }
251 
252 
253 LinearSystemWrapperItpack::Float LinearSystemWrapperItpack::GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex) const
254 {
255  /* error checking */
256  if (!m_Matrices)
257  {
258  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetMatrixValue", "No matrices have been allocated");
259  }
260  if (matrixIndex >= m_NumberOfMatrices)
261  {
262  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetMatrixValue", "m_Matrices", matrixIndex);
263  }
264  if ( (i >= m_Order) || (j >= m_Order) )
265  {
266  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetMatrixValue", "m_Matrices[]", i,j);
267  }
268 
269  /* return value */
270  return (*m_Matrices)[matrixIndex].Get(i,j);
271 }
272 
273 
274 void LinearSystemWrapperItpack::SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex)
275 {
276  /* error checking */
277  if (!m_Matrices)
278  {
279  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetMatrixValue", "No matrices have been allocated");
280  }
281  if ( (i >= m_Order) || (j >= m_Order) )
282  {
283  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetMatrixValue", "m_Matrices[]", i, j);
284  }
285  if (matrixIndex >= m_NumberOfMatrices)
286  {
287  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetMatrixValue", "m_Matrices", matrixIndex);
288  }
289 
290  /* set value */
291  ((*m_Matrices)[matrixIndex]).Set(i,j,value);
292 }
293 
294 
295 void LinearSystemWrapperItpack::AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex)
296 {
297  // FIXME: error checking
298  if (!m_Matrices)
299  {
300  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddMatrixValue", "No matrices have been allocated");
301  }
302  if ( (i >= m_Order) || (j >= m_Order) )
303  {
304  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddMatrixValue", "m_Matrices[]", i, j);
305  }
306  if (matrixIndex >= m_NumberOfMatrices)
307  {
308  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddMatrixValue", "m_Matrices", matrixIndex);
309  }
310 
311  ((*m_Matrices)[matrixIndex]).Add(i,j,value);
312 }
313 
314 
315 void LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow( unsigned int row, ColumnArray& cols, unsigned int matrixIndex )
316 {
317 
318  /* FIXME: error checking */
319  if (!m_Matrices)
320  {
321  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow", "No matrices have been allocated");
322  }
323  if (row >= this->m_Order)
324  {
325  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow", "m_Matrices[]", row);
326  }
327  if (matrixIndex >= m_NumberOfMatrices)
328  {
329  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow", "m_Matrices", matrixIndex);
330  }
331 
332  cols.clear();
333 
334  ItpackSparseMatrix* mat=&(*m_Matrices)[matrixIndex];
335 
336  /* Check if matrix is in readable form */
337  if (mat->m_MatrixFinalized )
338  {
339  /* get search bounds in appropriate row */
340  unsigned int lower = mat->m_IA[row]-1;
341  unsigned int upper = mat->m_IA[row+1]-1;
342 
343  for(unsigned int j=lower; j<upper; j++)
344  {
345  cols.push_back(mat->m_JA[j]-1);
346  }
347  }
348  else /* Scan the linked list to obtain the correct indices. */
349  {
350  int wrk=mat->m_IA[row]-1;
351  while(wrk>0)
352  {
353  cols.push_back(mat->m_JA[wrk]-1);
354  wrk=mat->m_IWORK[wrk]-1;
355  }
356  }
357 
358 }
359 
360 
361 void LinearSystemWrapperItpack::ScaleMatrix(Float scale, unsigned int matrixIndex)
362 {
363  /* error checking */
364  if (!m_Matrices)
365  {
366  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::ScaleMatrix", "No matrices have been allocated");
367  }
368  if (matrixIndex >= m_NumberOfMatrices)
369  {
370  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::ScaleMatrix", "m_Matrices", matrixIndex);
371  }
372 
373  int i;
374  doublereal *values = (*m_Matrices)[matrixIndex].GetA();
375  for (i=0; i<(*m_Matrices)[matrixIndex].GetIA()[this->m_Order] - 1; i++)
376  {
377  values[i] = values[i]*scale;
378  }
379 
380 }
381 
382 
383 LinearSystemWrapperItpack::Float LinearSystemWrapperItpack::GetVectorValue(unsigned int i, unsigned int vectorIndex) const
384 {
385  /* error checking */
386  if (!m_Vectors)
387  {
388  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetVectorValue", "No vectors have been allocated");
389  }
390  if (i >= m_Order)
391  {
392  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetVectorValue", "m_Vectors[]", i);
393  }
394  if (vectorIndex >= m_NumberOfVectors)
395  {
396  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetVectorValue", "m_Vectors", vectorIndex);
397  }
398  if ( !(*m_Vectors)[vectorIndex] )
399  {
400  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetVectorValue", "Indexed vector not yet allocated");
401  }
402 
403  /* return value */
404  return (*m_Vectors)[vectorIndex][i];
405 }
406 
407 
408 void LinearSystemWrapperItpack::SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex)
409 {
410  /* error checking */
411  if (!m_Vectors)
412  {
413  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetVectorValue", "No vectors have been allocated");
414  }
415  if (i >= m_Order)
416  {
417  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetVectorValue", "m_Vectors[]", i);
418  }
419  if (vectorIndex >= m_NumberOfVectors)
420  {
421  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetVectorValue", "m_Vectors", vectorIndex);
422  }
423  if ( !(*m_Vectors)[vectorIndex] )
424  {
425  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetVectorValue", "Indexed vector not yet allocated");
426  }
427 
428  /* set value */
429  (*m_Vectors)[vectorIndex][i] = value;
430 }
431 
432 
433 void LinearSystemWrapperItpack::AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex)
434 {
435  /*( error checking */
436  if (!m_Vectors)
437  {
438  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddVectorValue", "No vectors have been allocated");
439  }
440  if (i >= m_Order)
441  {
442  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddVectorValue", "m_Vectors[]", i);
443  }
444  if (vectorIndex >= m_NumberOfVectors)
445  {
446  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddVectorValue", "m_Vectors", vectorIndex);
447  }
448  if ( !(*m_Vectors)[vectorIndex] )
449  {
450  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddVectorValue", "Indexed vector not yet allocated");
451  }
452 
453  /* add value */
454  (*m_Vectors)[vectorIndex][i] += value;
455 }
456 
457 
458 LinearSystemWrapperItpack::Float LinearSystemWrapperItpack::GetSolutionValue(unsigned int i, unsigned int solutionIndex) const
459 {
460  // FIXME: error checking
461  if ( (i >= m_Order) || !m_Solutions || (solutionIndex >= m_NumberOfSolutions) || !(*m_Solutions)[solutionIndex] )
462  {
463  return 0.0;
464  }
465  return (*m_Solutions)[solutionIndex][i];
466 }
467 
468 
469 void LinearSystemWrapperItpack::SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex)
470 {
471  /* error checking */
472  if (!m_Solutions)
473  {
474  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetSolutionValue", "No solutions have been allocated");
475  }
476  if (i >= m_Order)
477  {
478  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetSolutionValue", "m_Solutions[]", i);
479  }
480  if (solutionIndex >= m_NumberOfSolutions)
481  {
482  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetSolutionValue", "m_Solutions", solutionIndex);
483  }
484  if ( !(*m_Solutions)[solutionIndex] )
485  {
486  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetSolutionValue", "Indexed solution not yet allocated");
487  }
488 
489  /* set value */
490  (*m_Solutions)[solutionIndex][i] = value;
491 
492 }
493 
494 
495 void LinearSystemWrapperItpack::AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex)
496 {
497  /* error checking */
498  if (!m_Solutions)
499  {
500  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddSolutionValue", "No solutions have been allocated");
501  }
502  if (i >= m_Order)
503  {
504  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddSolutionValue", "m_Solutions[]", i);
505  }
506  if (solutionIndex >= m_NumberOfSolutions)
507  {
508  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddSolutionValue", "m_Solutions", solutionIndex);
509  }
510  if ( !(*m_Solutions)[solutionIndex] )
511  {
512  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddSolutionValue", "Indexed solution not yet allocated");
513  }
514 
515  (*m_Solutions)[solutionIndex][i] += value;
516 }
517 
518 
519 void LinearSystemWrapperItpack::Solve(void)
520 {
521 
522  /* error checking */
523  if (!m_Order || !m_Matrices || !m_Vectors || !m_Solutions)
524  {
525  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::Solve", "Not all necessary data members have been allocated");
526  }
527  if ( !(*m_Matrices)[0].GetOrder() ) {
528  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddSolutionValue", "Primary matrix never filled");
529  }
530 
531  /* itpack variables */
532  integer N;
533  integer NB;
534  integer NW;
535  integer NCG;
536  integer *IWKSP;
537  doublereal *WKSP;
538  integer IERR = 0;
539 
540  /* *******************************************************************
541  * FIX ME: itpack does not allow for any non-zero diagonal elements
542  * so "very small" numbers are inserted to allow for a solution
543  *
544  int i;
545  doublereal fakeZero = 1.0e-16;
546 
547  //insert "fake" zeros
548  for (i=0; i<static_cast<int>(m_Order); i++)
549  {
550  if ( (*m_Matrices)[0].Get(i,i) == 0.0)
551  {
552  (*m_Matrices)[0].Set(i,i,fakeZero);
553  }
554  }
555  // END FIXME
556  *********************************************************************/
557 
558 
559  /* Initialize solution values (set to zero) */
560  if (!this->IsSolutionInitialized(0))
561  {
562  this->InitializeSolution(0);
563  }
564 
565 
566  /* Set up itpack workspace variables
567  *
568  * N -> Order of system
569  *
570  * NCG -> temp var
571  *
572  * NW -> on input: length of wksp, on output: actual amount used
573  * jcg_() - 4*N + NCG
574  * jsi_() - 2*N
575  * sor_() - N
576  * ssorcg_() - 6*N + NCG
577  * ssorsi_() - 5*N
578  * rscg_() - N + 3*NB + NCG
579  * rssi_() - N + NB
580  *
581  * IWKSP -> temp variable used by itpack (integer[3*N])
582  * WKSP -> temp variable used by itpack (doublereal[NW])
583  */
584  N = m_Order;
585  if (m_IPARM[4] == 1)
586  {
587  NCG = 4 * m_IPARM[0];
588  }
589  else
590  {
591  NCG = 2 * m_IPARM[0];
592  }
593  NB = m_IPARM[8] + N; // upper bound of what can be computed in prbndx_
594  switch ( m_Method )
595  {
596  case 0:
597  NW = 4*N + NCG;
598  break;
599  case 1:
600  NW = 2*N;
601  break;
602  case 2:
603  NW = N;
604  break;
605  case 3:
606  NW = 6*N + NCG;
607  break;
608  case 4:
609  NW = 5*N;
610  break;
611  case 5:
612  NW = N + 3*NB + NCG;
613  break;
614  case 6:
615  NW = N + NB;
616  break;
617  }
618  m_IPARM[7] = NW;
619  IWKSP = new integer [ 3*N ];
620  WKSP = new doublereal [ NW+2 ];
621 
622  integer i;
623  for (i=0; i<NW; i++)
624  {
625  WKSP[i] = 0.0;
626  }
627  for (i=0; i<(3*N); i++)
628  {
629  IWKSP[i] = 0;
630  }
631 
632 
633  // Save maximum number of iteration, since it will
634  // be overwritten.
635  int max_num_iterations=m_IPARM[0];
636 
637  /* call to itpack solver routine */
638  (*m_Methods[m_Method])( &N, (*m_Matrices)[0].GetIA(), (*m_Matrices)[0].GetJA(), (*m_Matrices)[0].GetA(),
639  (*m_Vectors)[0], (*m_Solutions)[0], &(IWKSP[0]), &NW, &(WKSP[0]), &(m_IPARM[0]), &(m_RPARM[0]), &IERR);
640 
641  m_IPARM[0]=max_num_iterations;
642 
643  /* remove exception throw on convergence failure */
644  if (IERR < 100)
645  {
646  if ( (IERR % 10) == 3 )
647  {
648  IERR = 0;
649  }
650  }
651 
652  /* clean up */
653  delete [] IWKSP;
654  delete [] WKSP;
655 
656  /* check for itpack error code */
657  if (IERR > 0)
658  {
659  throw FEMExceptionItpackSolver(__FILE__, __LINE__, "LinearSystemWrapperItpack::Solve", IERR);
660  }
661 
662 }
663 
664 
665 void LinearSystemWrapperItpack::SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2)
666 {
667 
668  /* error checking */
669  if (!m_Matrices)
670  {
671  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapMatrices", "No matrices allocated");
672  }
673  if (matrixIndex1 >= m_NumberOfMatrices)
674  {
675  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapMatrices", "m_Matrices", matrixIndex1);
676  }
677  if (matrixIndex2 >= m_NumberOfMatrices)
678  {
679  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapMatrices", "m_Matrices", matrixIndex2);
680  }
681 
682  int n = (*m_Matrices)[matrixIndex2].GetOrder();
683  int nz = (*m_Matrices)[matrixIndex2].GetMaxNonZeroValues();
684  integer* ia = (*m_Matrices)[matrixIndex2].GetIA();
685  integer *ja = (*m_Matrices)[matrixIndex2].GetJA();
686  doublereal *a = (*m_Matrices)[matrixIndex2].GetA();
687 
688  (*m_Matrices)[matrixIndex2].SetOrder( (*m_Matrices)[matrixIndex1].GetOrder() );
689  (*m_Matrices)[matrixIndex2].SetMaxNonZeroValues ( (*m_Matrices)[matrixIndex1].GetMaxNonZeroValues() );
690  (*m_Matrices)[matrixIndex2].SetCompressedRow((*m_Matrices)[matrixIndex1].GetIA(),
691  (*m_Matrices)[matrixIndex1].GetJA(), (*m_Matrices)[matrixIndex1].GetA() );
692 
693  (*m_Matrices)[matrixIndex1].SetOrder( n );
694  (*m_Matrices)[matrixIndex1].SetMaxNonZeroValues ( nz );
695  (*m_Matrices)[matrixIndex1].SetCompressedRow(ia, ja, a);
696 
697 
698 }
699 
700 
701 void LinearSystemWrapperItpack::SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2)
702 {
703  /* error checking */
704  if (!m_Vectors)
705  {
706  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapVectors", "No vectors allocated");
707  }
708  if (vectorIndex1 >= m_NumberOfVectors)
709  {
710  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapVectors", "m_Vectors", vectorIndex1);
711  }
712  if (vectorIndex2 >= m_NumberOfVectors)
713  {
714  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapVectors", "m_Vectors", vectorIndex2);
715  }
716 
717  VectorRepresentation temp = (*m_Vectors)[vectorIndex1];
718 
719  (*m_Vectors)[vectorIndex1] = (*m_Vectors)[vectorIndex2];
720  (*m_Vectors)[vectorIndex2] = temp;
721 }
722 
723 
724 void LinearSystemWrapperItpack::SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2)
725 {
726  /* error checking */
727  if (!m_Solutions)
728  {
729  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapSolutions", "No solutions allocated");
730  }
731  if (solutionIndex1 >= m_NumberOfSolutions)
732  {
733  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapSolutions", "m_Solutions", solutionIndex1);
734  }
735  if (solutionIndex2 >= m_NumberOfSolutions)
736  {
737  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapSolutions", "m_Solutions", solutionIndex2);
738  }
739 
740  VectorRepresentation temp = (*m_Solutions)[solutionIndex1];
741 
742  (*m_Solutions)[solutionIndex1] = (*m_Solutions)[solutionIndex2];
743  (*m_Solutions)[solutionIndex2] = temp;
744 }
745 
746 
747 void LinearSystemWrapperItpack::CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex)
748 {
749 
750  /* error checking */
751  if (!m_Vectors)
752  {
753  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "No vectors allocated");
754  }
755  if (!m_Solutions)
756  {
757  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "No solutions allocated");
758  }
759  if (vectorIndex >= m_NumberOfVectors)
760  {
761  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "m_Vectors", vectorIndex);
762  }
763  if (solutionIndex >= m_NumberOfSolutions)
764  {
765  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "m_Solutions", solutionIndex);
766  }
767 
768 
769  this->InitializeVector(vectorIndex);
770 
771  /* copy values */
772  for (unsigned int i=0; i<m_Order; i++)
773  {
774  (*m_Vectors)[vectorIndex][i] = (*m_Solutions)[solutionIndex][i];
775  }
776 }
777 
778 
779 void LinearSystemWrapperItpack::CopyVector2Solution(unsigned vectorIndex, unsigned int solutionIndex)
780 {
781 
782  /* error checking */
783  if (!m_Vectors)
784  {
785  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "No vectors allocated");
786  }
787  if (!m_Solutions)
788  {
789  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "No solutions allocated");
790  }
791  if (vectorIndex >= m_NumberOfVectors)
792  {
793  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "m_Vectors", vectorIndex);
794  }
795  if (solutionIndex >= m_NumberOfSolutions)
796  {
797  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "m_Solutions", solutionIndex);
798  }
799 
800 
801  this->InitializeSolution(solutionIndex);
802 
803  /* copy values */
804  for (unsigned int i=0; i<m_Order; i++)
805  {
806  (*m_Solutions)[solutionIndex][i] = (*m_Vectors)[vectorIndex][i];
807  }
808 }
809 
810 
811 void LinearSystemWrapperItpack::MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex, unsigned int rightMatrixIndex)
812 {
813  /* error checking */
814  if (!m_Matrices)
815  {
816  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixMatrix", "No matrices allocated");
817  }
818  if (resultMatrixIndex >= m_NumberOfMatrices)
819  {
820  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixMatrix", "m_Matrices", resultMatrixIndex);
821  }
822  if (leftMatrixIndex >= m_NumberOfMatrices)
823  {
824  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixMatrix", "m_Matrices", leftMatrixIndex);
825  }
826  if (rightMatrixIndex >= m_NumberOfMatrices)
827  {
828  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixMatrix", "m_Matrices", rightMatrixIndex);
829  }
830 
831  (*m_Matrices)[leftMatrixIndex].mult( &((*m_Matrices)[rightMatrixIndex]), &((*m_Matrices)[resultMatrixIndex]) );
832 
833 }
834 
835 
836 void LinearSystemWrapperItpack::MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex)
837 {
838 
839  /* error checking */
840  if (!m_Matrices)
841  {
842  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "No matrices allocated");
843  }
844  if (!m_Vectors)
845  {
846  throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "No vectors allocated");
847  }
848  if (resultVectorIndex >= m_NumberOfVectors)
849  {
850  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "m_Vectors", resultVectorIndex);
851  }
852  if (matrixIndex >= m_NumberOfMatrices)
853  {
854  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "m_Matrices", matrixIndex);
855  }
856  if (vectorIndex >= m_NumberOfVectors)
857  {
858  throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "m_Vectors", vectorIndex);
859  }
860 
861  /* perform mult */
862  (*m_Matrices)[matrixIndex].mult( (*m_Vectors)[vectorIndex], (*m_Vectors)[resultVectorIndex] );
863 
864 }
865 
866 
867 LinearSystemWrapperItpack::~LinearSystemWrapperItpack(void)
868 {
869  delete m_Matrices;
870 
871  unsigned int i;
872  if (m_Vectors != 0)
873  {
874  for (i=0; i<m_NumberOfVectors; i++)
875  {
876  if ( (*m_Vectors)[i] != 0 )
877  {
878  delete [] (*m_Vectors)[i];
879  }
880  }
881  delete m_Vectors;
882  }
883 
884 
885  if (m_Solutions != 0)
886  {
887  for (i=0; i<m_NumberOfSolutions; i++)
888  {
889  if ( (*m_Solutions)[i] != 0 )
890  {
891  delete [] (*m_Solutions)[i];
892  }
893  }
894  delete m_Solutions;
895  }
896 
897 
898 
899 }
900 
901 
902 FEMExceptionItpackSolver::FEMExceptionItpackSolver(const char *file, unsigned int lineNumber, std::string location, integer errorCode) :
903  FEMException(file,lineNumber)
904 {
905  std::string solverError;
906 
907  if (errorCode < 100)
908  {
909  errorCode = errorCode % 10;
910  }
911 
912  switch (errorCode)
913  {
914  case 1 :
915  solverError = "Invalid order of system";
916  break;
917  case 2 :
918  solverError = "Workspace is not large enough";
919  break;
920  case 3 :
921  solverError = "Failure to converge before reaching maximum number of iterations";
922  break;
923  case 4 :
924  solverError = "Invalid order of black subsystem";
925  break;
926  case 101:
927  solverError = "A diagonal element is not positive";
928  break;
929  case 102 :
930  solverError = "No diagonal element in a row";
931  break;
932  case 201 :
933  solverError = "Red-black indexing is not possible";
934  break;
935  case 301 :
936  solverError = "No entry in a row of the original matrix";
937  break;
938  case 302 :
939  solverError = "No entry in a row of the permuted matrix";
940  break;
941  case 303 :
942  solverError = "Sorting error in a row of the permuted matrix";
943  break;
944  case 401 :
945  solverError = "A diagonal element is not positive";
946  break;
947  case 402 :
948  solverError = "No diagonal element in a row";
949  break;
950  case 501:
951  solverError = "Failure to converge before reaching maximum number of iterations";
952  break;
953  case 502 :
954  solverError = "Function does not change sign at endpoints";
955  break;
956  case 601 :
957  solverError = "Successive iterations are not monotone increasing";
958  break;
959  default :
960  solverError = "Unknown error code returned";
961  }
962 
963  OStringStream buf;
964  buf << "Error: " << solverError;
965 
966  SetDescription(buf.str().c_str());
967 
968  SetLocation(location);
969 
970 }
971 
972 }} // end namespace itk::fem
Generated at Sat Feb 2 2013 23:36:46 for Orfeo Toolbox with doxygen 1.8.1.1