1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #include "itkNumericTraits.h"
19 #include "itpack.h"
20 #include "itkFEMLinearSystemWrapperItpack.h"
21
22 namespace itk
23 {
24 namespace fem
25 {
26 /**
27 * constructor
28 */
LinearSystemWrapperItpack()29 LinearSystemWrapperItpack::LinearSystemWrapperItpack()
30 {
31 /* fill m_Methods with pointers to solver functions */
32 m_Methods[0] = jcg_;
33 m_Methods[1] = jsi_;
34 m_Methods[2] = sor_;
35 m_Methods[3] = ssorcg_;
36 m_Methods[4] = ssorsi_;
37 m_Methods[5] = rscg_;
38 m_Methods[6] = rssi_;
39 m_Method = 0; /* set default method to jcg_ */
40
41 /* Set solving parameters */
42 dfault_( &( m_IPARM[0] ), &( m_RPARM[0] ) );
43
44 // We don't want the solver routines to
45 // overwrite the parameters.
46 m_IPARM[2] = 1;
47
48 /* m_IPARM[0] = 500; */ /* number of iterations */
49 m_IPARM[1] = -1; /* no error message output */
50 m_IPARM[4] = 1; /* non-symmetric matrix */
51
52 /* itpack recommended (but not default) value */
53 #undef min
54 m_RPARM[7] = 500.0 * NumericTraits<double>::min();
55
56 m_MaximumNonZeroValues = 0;
57 m_Matrices = nullptr;
58 m_Vectors = nullptr;
59 m_Solutions = nullptr;
60 }
61
InitializeMatrix(unsigned int matrixIndex)62 void LinearSystemWrapperItpack::InitializeMatrix(unsigned int matrixIndex)
63 {
64 /* error checking */
65 if( !m_Order )
66 {
67 throw FEMExceptionLinearSystem(__FILE__,
68 __LINE__,
69 "LinearSystemWrapperItpack::InitializeMatrix",
70 "System order not set");
71 }
72 if( matrixIndex >= m_NumberOfMatrices )
73 {
74 throw FEMExceptionLinearSystemBounds(__FILE__,
75 __LINE__,
76 "LinearSystemWrapperItpack::InitializeMatrix",
77 "m_Matrices",
78 matrixIndex);
79 }
80 if( !m_MaximumNonZeroValues )
81 {
82 throw FEMExceptionLinearSystem(__FILE__,
83 __LINE__,
84 "LinearSystemWrapperItpack::InitializeMatrix",
85 "Maximum number of non zeros not set");
86 }
87
88 // allocate if necessay
89 if( m_Matrices == nullptr )
90 {
91 m_Matrices = new MatrixHolder(m_NumberOfMatrices);
92 }
93
94 /* Set required variables */
95 ( *m_Matrices )[matrixIndex].Clear();
96 ( *m_Matrices )[matrixIndex].SetOrder(m_Order);
97 ( *m_Matrices )[matrixIndex].SetMaxNonZeroValues(m_MaximumNonZeroValues);
98 }
99
IsMatrixInitialized(unsigned int matrixIndex)100 bool LinearSystemWrapperItpack::IsMatrixInitialized(unsigned int matrixIndex)
101 {
102 if( !m_Matrices )
103 {
104 return false;
105 }
106 if( !( *m_Matrices )[matrixIndex].GetOrder() )
107 {
108 return false;
109 }
110 if( !( *m_Matrices )[matrixIndex].GetMaxNonZeroValues() )
111 {
112 return false;
113 }
114
115 return true;
116 }
117
InitializeVector(unsigned int vectorIndex)118 void LinearSystemWrapperItpack::InitializeVector(unsigned int vectorIndex)
119 {
120 /* error checking */
121 if( !m_Order )
122 {
123 throw FEMExceptionLinearSystem(__FILE__,
124 __LINE__,
125 "LinearSystemWrapperItpack::InitializeVector",
126 "System order not set");
127 }
128 if( vectorIndex >= m_NumberOfVectors )
129 {
130 throw FEMExceptionLinearSystemBounds(__FILE__,
131 __LINE__,
132 "LinearSystemWrapperItpack::InitializeVector",
133 "m_Vectors",
134 vectorIndex);
135 }
136
137 /* allocate if necessay */
138 if( m_Vectors == nullptr )
139 {
140 m_Vectors = new VectorHolder(m_NumberOfVectors);
141 }
142
143 /* delete old vector */
144 delete[] ( *m_Vectors )[vectorIndex];
145
146 /* insert new vector */
147 ( *m_Vectors )[vectorIndex] = new doublereal[m_Order];
148 /* fill with zeros */
149 for( unsigned int i = 0; i < m_Order; i++ )
150 {
151 ( *m_Vectors )[vectorIndex][i] = 0.0;
152 }
153 }
154
IsVectorInitialized(unsigned int vectorIndex)155 bool LinearSystemWrapperItpack::IsVectorInitialized(unsigned int vectorIndex)
156 {
157 if( !m_Vectors )
158 {
159 return false;
160 }
161 if( !( *m_Vectors )[vectorIndex] )
162 {
163 return false;
164 }
165
166 return true;
167 }
168
InitializeSolution(unsigned int solutionIndex)169 void LinearSystemWrapperItpack::InitializeSolution(unsigned int solutionIndex)
170 {
171 /* FIX ME: exceptions */
172 if( !m_Order )
173 {
174 throw FEMExceptionLinearSystem(__FILE__,
175 __LINE__,
176 "LinearSystemWrapperItpack::InitializeSolution",
177 "System order not set");
178 }
179 if( solutionIndex >= m_NumberOfSolutions )
180 {
181 throw FEMExceptionLinearSystemBounds(__FILE__,
182 __LINE__,
183 "LinearSystemWrapperItpack::InitializeSolution",
184 "m_Solutions",
185 solutionIndex);
186 }
187
188 // allocate if necessay
189 if( m_Solutions == nullptr )
190 {
191 m_Solutions = new VectorHolder(m_NumberOfSolutions);
192 }
193
194 /* delete old vector */
195 delete[] ( *m_Solutions )[solutionIndex];
196
197 /* insert new vector */
198 ( *m_Solutions )[solutionIndex] = new doublereal[m_Order];
199 /* fill with zeros */
200 for( unsigned int i = 0; i < m_Order; i++ )
201 {
202 ( *m_Solutions )[solutionIndex][i] = 0.0;
203 }
204 }
205
IsSolutionInitialized(unsigned int solutionIndex)206 bool LinearSystemWrapperItpack::IsSolutionInitialized(unsigned int solutionIndex)
207 {
208 if( !m_Solutions )
209 {
210 return false;
211 }
212 if( !( *m_Solutions )[solutionIndex] )
213 {
214 return false;
215 }
216
217 return true;
218 }
219
DestroyMatrix(unsigned int matrixIndex)220 void LinearSystemWrapperItpack::DestroyMatrix(unsigned int matrixIndex)
221 {
222 /* FIX ME: exceptions */
223 if( m_Matrices )
224 {
225 if( matrixIndex >= m_NumberOfMatrices )
226 {
227 throw FEMExceptionLinearSystemBounds(__FILE__,
228 __LINE__,
229 "LinearSystemWrapperItpack::DestroyMatrix",
230 "m_Matrices",
231 matrixIndex);
232 }
233
234 ( *m_Matrices )[matrixIndex].Clear();
235 }
236 }
237
DestroyVector(unsigned int vectorIndex)238 void LinearSystemWrapperItpack::DestroyVector(unsigned int vectorIndex)
239 {
240 /* FIXME: exceptions */
241 if( m_Vectors )
242 {
243 if( vectorIndex >= m_NumberOfVectors )
244 {
245 throw FEMExceptionLinearSystemBounds(__FILE__,
246 __LINE__,
247 "LinearSystemWrapperItpack::DestroyVector",
248 "m_Vectors",
249 vectorIndex);
250 }
251
252 /* delete vector */
253 delete[] ( *m_Vectors )[vectorIndex];
254 ( *m_Vectors )[vectorIndex] = nullptr;
255 }
256 }
257
DestroySolution(unsigned int solutionIndex)258 void LinearSystemWrapperItpack::DestroySolution(unsigned int solutionIndex)
259 {
260 // FIXME: exceptions
261 if( m_Solutions )
262 {
263 if( solutionIndex >= m_NumberOfSolutions )
264 {
265 throw FEMExceptionLinearSystemBounds(__FILE__,
266 __LINE__,
267 "LinearSystemWrapperItpack::DestroySolution",
268 "m_Solutions",
269 solutionIndex);
270 }
271
272 /* delete vector */
273 delete[] ( *m_Solutions )[solutionIndex];
274 ( *m_Solutions )[solutionIndex] = nullptr;
275 }
276 }
277
GetMatrixValue(unsigned int i,unsigned int j,unsigned int matrixIndex) const278 LinearSystemWrapperItpack::Float LinearSystemWrapperItpack::GetMatrixValue(unsigned int i,
279 unsigned int j,
280 unsigned int matrixIndex) const
281 {
282 /* error checking */
283 if( !m_Matrices )
284 {
285 throw FEMExceptionLinearSystem(__FILE__,
286 __LINE__,
287 "LinearSystemWrapperItpack::GetMatrixValue",
288 "No matrices have been allocated");
289 }
290 if( matrixIndex >= m_NumberOfMatrices )
291 {
292 throw FEMExceptionLinearSystemBounds(__FILE__,
293 __LINE__,
294 "LinearSystemWrapperItpack::GetMatrixValue",
295 "m_Matrices",
296 matrixIndex);
297 }
298 if( ( i >= m_Order ) || ( j >= m_Order ) )
299 {
300 throw FEMExceptionLinearSystemBounds(__FILE__,
301 __LINE__,
302 "LinearSystemWrapperItpack::GetMatrixValue",
303 "m_Matrices[]",
304 i,
305 j);
306 }
307
308 /* return value */
309 return ( *m_Matrices )[matrixIndex].Get(i, j);
310 }
311
SetMatrixValue(unsigned int i,unsigned int j,Float value,unsigned int matrixIndex)312 void LinearSystemWrapperItpack::SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex)
313 {
314 /* error checking */
315 if( !m_Matrices )
316 {
317 throw FEMExceptionLinearSystem(__FILE__,
318 __LINE__,
319 "LinearSystemWrapperItpack::SetMatrixValue",
320 "No matrices have been allocated");
321 }
322 if( ( i >= m_Order ) || ( j >= m_Order ) )
323 {
324 throw FEMExceptionLinearSystemBounds(__FILE__,
325 __LINE__,
326 "LinearSystemWrapperItpack::SetMatrixValue",
327 "m_Matrices[]",
328 i,
329 j);
330 }
331 if( matrixIndex >= m_NumberOfMatrices )
332 {
333 throw FEMExceptionLinearSystemBounds(__FILE__,
334 __LINE__,
335 "LinearSystemWrapperItpack::SetMatrixValue",
336 "m_Matrices",
337 matrixIndex);
338 }
339
340 /* set value */
341 ( ( *m_Matrices )[matrixIndex] ).Set(i, j, value);
342 }
343
AddMatrixValue(unsigned int i,unsigned int j,Float value,unsigned int matrixIndex)344 void LinearSystemWrapperItpack::AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex)
345 {
346 // FIXME: error checking
347 if( !m_Matrices )
348 {
349 throw FEMExceptionLinearSystem(__FILE__,
350 __LINE__,
351 "LinearSystemWrapperItpack::AddMatrixValue",
352 "No matrices have been allocated");
353 }
354 if( ( i >= m_Order ) || ( j >= m_Order ) )
355 {
356 throw FEMExceptionLinearSystemBounds(__FILE__,
357 __LINE__,
358 "LinearSystemWrapperItpack::AddMatrixValue",
359 "m_Matrices[]",
360 i,
361 j);
362 }
363 if( matrixIndex >= m_NumberOfMatrices )
364 {
365 throw FEMExceptionLinearSystemBounds(__FILE__,
366 __LINE__,
367 "LinearSystemWrapperItpack::AddMatrixValue",
368 "m_Matrices",
369 matrixIndex);
370 }
371
372 ( ( *m_Matrices )[matrixIndex] ).Add(i, j, value);
373 }
374
GetColumnsOfNonZeroMatrixElementsInRow(unsigned int row,ColumnArray & cols,unsigned int matrixIndex)375 void LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow(unsigned int row,
376 ColumnArray & cols,
377 unsigned int matrixIndex)
378 {
379 /* FIXME: error checking */
380 if( !m_Matrices )
381 {
382 throw FEMExceptionLinearSystem(__FILE__,
383 __LINE__,
384 "LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow",
385 "No matrices have been allocated");
386 }
387 if( row >= this->m_Order )
388 {
389 throw FEMExceptionLinearSystemBounds(__FILE__,
390 __LINE__,
391 "LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow",
392 "m_Matrices[]",
393 row);
394 }
395 if( matrixIndex >= m_NumberOfMatrices )
396 {
397 throw FEMExceptionLinearSystemBounds(__FILE__,
398 __LINE__,
399 "LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow",
400 "m_Matrices",
401 matrixIndex);
402 }
403
404 cols.clear();
405
406 ItpackSparseMatrix *mat = &( *m_Matrices )[matrixIndex];
407
408 /* Check if matrix is in readable form */
409 if( mat->m_MatrixFinalized )
410 {
411 /* get search bounds in appropriate row */
412 unsigned int lower = mat->m_IA[row] - 1;
413 unsigned int upper = mat->m_IA[row + 1] - 1;
414 for( unsigned int j = lower; j < upper; j++ )
415 {
416 cols.push_back(mat->m_JA[j] - 1);
417 }
418 }
419 else /* Scan the linked list to obtain the correct indices. */
420 {
421 int wrk = mat->m_IA[row] - 1;
422 while( wrk > 0 )
423 {
424 cols.push_back(mat->m_JA[wrk] - 1);
425 wrk = mat->m_IWORK[wrk] - 1;
426 }
427 }
428 }
429
ScaleMatrix(Float scale,unsigned int matrixIndex)430 void LinearSystemWrapperItpack::ScaleMatrix(Float scale, unsigned int matrixIndex)
431 {
432 /* error checking */
433 if( !m_Matrices )
434 {
435 throw FEMExceptionLinearSystem(__FILE__,
436 __LINE__,
437 "LinearSystemWrapperItpack::ScaleMatrix",
438 "No matrices have been allocated");
439 }
440 if( matrixIndex >= m_NumberOfMatrices )
441 {
442 throw FEMExceptionLinearSystemBounds(__FILE__,
443 __LINE__,
444 "LinearSystemWrapperItpack::ScaleMatrix",
445 "m_Matrices",
446 matrixIndex);
447 }
448
449 int i;
450 doublereal *values = ( *m_Matrices )[matrixIndex].GetA();
451 for( i = 0; i < ( *m_Matrices )[matrixIndex].GetIA()[this->m_Order] - 1; i++ )
452 {
453 values[i] = values[i] * scale;
454 }
455 }
456
GetVectorValue(unsigned int i,unsigned int vectorIndex) const457 LinearSystemWrapperItpack::Float LinearSystemWrapperItpack::GetVectorValue(unsigned int i,
458 unsigned int vectorIndex) const
459 {
460 /* error checking */
461 if( !m_Vectors )
462 {
463 throw FEMExceptionLinearSystem(__FILE__,
464 __LINE__,
465 "LinearSystemWrapperItpack::GetVectorValue",
466 "No vectors have been allocated");
467 }
468 if( i >= m_Order )
469 {
470 throw FEMExceptionLinearSystemBounds(__FILE__,
471 __LINE__,
472 "LinearSystemWrapperItpack::GetVectorValue",
473 "m_Vectors[]",
474 i);
475 }
476 if( vectorIndex >= m_NumberOfVectors )
477 {
478 throw FEMExceptionLinearSystemBounds(__FILE__,
479 __LINE__,
480 "LinearSystemWrapperItpack::GetVectorValue",
481 "m_Vectors",
482 vectorIndex);
483 }
484 if( !( *m_Vectors )[vectorIndex] )
485 {
486 throw FEMExceptionLinearSystem(__FILE__,
487 __LINE__,
488 "LinearSystemWrapperItpack::GetVectorValue",
489 "Indexed vector not yet allocated");
490 }
491
492 /* return value */
493 return ( *m_Vectors )[vectorIndex][i];
494 }
495
SetVectorValue(unsigned int i,Float value,unsigned int vectorIndex)496 void LinearSystemWrapperItpack::SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex)
497 {
498 /* error checking */
499 if( !m_Vectors )
500 {
501 throw FEMExceptionLinearSystem(__FILE__,
502 __LINE__,
503 "LinearSystemWrapperItpack::SetVectorValue",
504 "No vectors have been allocated");
505 }
506 if( i >= m_Order )
507 {
508 throw FEMExceptionLinearSystemBounds(__FILE__,
509 __LINE__,
510 "LinearSystemWrapperItpack::SetVectorValue",
511 "m_Vectors[]",
512 i);
513 }
514 if( vectorIndex >= m_NumberOfVectors )
515 {
516 throw FEMExceptionLinearSystemBounds(__FILE__,
517 __LINE__,
518 "LinearSystemWrapperItpack::SetVectorValue",
519 "m_Vectors",
520 vectorIndex);
521 }
522 if( !( *m_Vectors )[vectorIndex] )
523 {
524 throw FEMExceptionLinearSystem(__FILE__,
525 __LINE__,
526 "LinearSystemWrapperItpack::SetVectorValue",
527 "Indexed vector not yet allocated");
528 }
529
530 /* set value */
531 ( *m_Vectors )[vectorIndex][i] = value;
532 }
533
AddVectorValue(unsigned int i,Float value,unsigned int vectorIndex)534 void LinearSystemWrapperItpack::AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex)
535 {
536 /*( error checking */
537 if( !m_Vectors )
538 {
539 throw FEMExceptionLinearSystem(__FILE__,
540 __LINE__,
541 "LinearSystemWrapperItpack::AddVectorValue",
542 "No vectors have been allocated");
543 }
544 if( i >= m_Order )
545 {
546 throw FEMExceptionLinearSystemBounds(__FILE__,
547 __LINE__,
548 "LinearSystemWrapperItpack::AddVectorValue",
549 "m_Vectors[]",
550 i);
551 }
552 if( vectorIndex >= m_NumberOfVectors )
553 {
554 throw FEMExceptionLinearSystemBounds(__FILE__,
555 __LINE__,
556 "LinearSystemWrapperItpack::AddVectorValue",
557 "m_Vectors",
558 vectorIndex);
559 }
560 if( !( *m_Vectors )[vectorIndex] )
561 {
562 throw FEMExceptionLinearSystem(__FILE__,
563 __LINE__,
564 "LinearSystemWrapperItpack::AddVectorValue",
565 "Indexed vector not yet allocated");
566 }
567
568 /* add value */
569 ( *m_Vectors )[vectorIndex][i] += value;
570 }
571
GetSolutionValue(unsigned int i,unsigned int solutionIndex) const572 LinearSystemWrapperItpack::Float LinearSystemWrapperItpack::GetSolutionValue(unsigned int i,
573 unsigned int solutionIndex) const
574 {
575 // FIXME: error checking
576 if( ( i >= m_Order ) || !m_Solutions || ( solutionIndex >= m_NumberOfSolutions ) || !( *m_Solutions )[solutionIndex] )
577 {
578 return 0.0;
579 }
580 return ( *m_Solutions )[solutionIndex][i];
581 }
582
SetSolutionValue(unsigned int i,Float value,unsigned int solutionIndex)583 void LinearSystemWrapperItpack::SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex)
584 {
585 /* error checking */
586 if( !m_Solutions )
587 {
588 throw FEMExceptionLinearSystem(__FILE__,
589 __LINE__,
590 "LinearSystemWrapperItpack::SetSolutionValue",
591 "No solutions have been allocated");
592 }
593 if( i >= m_Order )
594 {
595 throw FEMExceptionLinearSystemBounds(__FILE__,
596 __LINE__,
597 "LinearSystemWrapperItpack::SetSolutionValue",
598 "m_Solutions[]",
599 i);
600 }
601 if( solutionIndex >= m_NumberOfSolutions )
602 {
603 throw FEMExceptionLinearSystemBounds(__FILE__,
604 __LINE__,
605 "LinearSystemWrapperItpack::SetSolutionValue",
606 "m_Solutions",
607 solutionIndex);
608 }
609 if( !( *m_Solutions )[solutionIndex] )
610 {
611 throw FEMExceptionLinearSystem(__FILE__,
612 __LINE__,
613 "LinearSystemWrapperItpack::SetSolutionValue",
614 "Indexed solution not yet allocated");
615 }
616
617 /* set value */
618 ( *m_Solutions )[solutionIndex][i] = value;
619 }
620
AddSolutionValue(unsigned int i,Float value,unsigned int solutionIndex)621 void LinearSystemWrapperItpack::AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex)
622 {
623 /* error checking */
624 if( !m_Solutions )
625 {
626 throw FEMExceptionLinearSystem(__FILE__,
627 __LINE__,
628 "LinearSystemWrapperItpack::AddSolutionValue",
629 "No solutions have been allocated");
630 }
631 if( i >= m_Order )
632 {
633 throw FEMExceptionLinearSystemBounds(__FILE__,
634 __LINE__,
635 "LinearSystemWrapperItpack::AddSolutionValue",
636 "m_Solutions[]",
637 i);
638 }
639 if( solutionIndex >= m_NumberOfSolutions )
640 {
641 throw FEMExceptionLinearSystemBounds(__FILE__,
642 __LINE__,
643 "LinearSystemWrapperItpack::AddSolutionValue",
644 "m_Solutions",
645 solutionIndex);
646 }
647 if( !( *m_Solutions )[solutionIndex] )
648 {
649 throw FEMExceptionLinearSystem(__FILE__,
650 __LINE__,
651 "LinearSystemWrapperItpack::AddSolutionValue",
652 "Indexed solution not yet allocated");
653 }
654
655 ( *m_Solutions )[solutionIndex][i] += value;
656 }
657
Solve()658 void LinearSystemWrapperItpack::Solve()
659 {
660 /* error checking */
661 if( !m_Order || !m_Matrices || !m_Vectors || !m_Solutions )
662 {
663 throw FEMExceptionLinearSystem(__FILE__,
664 __LINE__,
665 "LinearSystemWrapperItpack::Solve",
666 "Not all necessary data members have been allocated");
667 }
668 if( !( *m_Matrices )[0].GetOrder() )
669 {
670 throw FEMExceptionLinearSystem(__FILE__,
671 __LINE__,
672 "LinearSystemWrapperItpack::AddSolutionValue",
673 "Primary matrix never filled");
674 }
675
676 /* itpack variables */
677 integer N;
678 integer NB;
679 integer NW;
680 integer NCG;
681 integer * IWKSP;
682 doublereal *WKSP;
683 integer IERR = 0;
684
685 /*********************************************************************
686 * FIX ME: itpack does not allow for any non-zero diagonal elements
687 * so "very small" numbers are inserted to allow for a solution
688 *
689 * int i;
690 * doublereal fakeZero = 1.0e-16;
691
692 * //insert "fake" zeros
693 * for (i=0; i<static_cast<int>(m_Order); i++)
694 * {
695 * if( (*m_Matrices)[0].Get(i,i) == 0.0)
696 * {
697 * (*m_Matrices)[0].Set(i,i,fakeZero);
698 * }
699 * }
700 * // END FIXME
701 *********************************************************************/
702
703 /* Initialize solution values (set to zero) */
704 if( !this->IsSolutionInitialized(0) )
705 {
706 this->InitializeSolution(0);
707 }
708
709 /* Set up itpack workspace variables
710 *
711 * N -> Order of system
712 *
713 * NCG -> temp var
714 *
715 * NW -> on input: length of wksp, on output: actual amount used
716 * jcg_() - 4*N + NCG
717 * jsi_() - 2*N
718 * sor_() - N
719 * ssorcg_() - 6*N + NCG
720 * ssorsi_() - 5*N
721 * rscg_() - N + 3*NB + NCG
722 * rssi_() - N + NB
723 *
724 * IWKSP -> temp variable used by itpack (integer[3*N])
725 * WKSP -> temp variable used by itpack (doublereal[NW])
726 */
727 N = m_Order;
728 if( m_IPARM[4] == 1 )
729 {
730 NCG = 4 * m_IPARM[0];
731 }
732 else
733 {
734 NCG = 2 * m_IPARM[0];
735 }
736 NB = m_IPARM[8] + N; // upper bound of what can be computed in prbndx_
737
738 switch( m_Method )
739 {
740 case 0:
741 NW = 4 * N + NCG;
742 break;
743 case 1:
744 NW = 2 * N;
745 break;
746 case 2:
747 NW = N;
748 break;
749 case 3:
750 NW = 6 * N + NCG;
751 break;
752 case 4:
753 NW = 5 * N;
754 break;
755 case 5:
756 NW = N + 3 * NB + NCG;
757 break;
758 case 6:
759 NW = N + NB;
760 break;
761 default:
762 std::ostringstream msg;
763 msg << "m_Method is" << m_Method << " but must be >=0 and <= 6";
764 throw FEMExceptionLinearSystem(__FILE__,
765 __LINE__,
766 "LinearSystemWrapperItpack::Solve",
767 msg.str().c_str());
768 }
769 m_IPARM[7] = NW;
770 IWKSP = new integer[3 * N];
771 WKSP = new doublereal[NW + 2];
772
773 integer i;
774 for( i = 0; i < NW; i++ )
775 {
776 WKSP[i] = 0.0;
777 }
778 for( i = 0; i < ( 3 * N ); i++ )
779 {
780 IWKSP[i] = 0;
781 }
782
783 // Save maximum number of iteration, since it will
784 // be overwritten.
785 int max_num_iterations = m_IPARM[0];
786
787 /* call to itpack solver routine */
788 ( *m_Methods[m_Method] )(&N, ( *m_Matrices )[0].GetIA(), ( *m_Matrices )[0].GetJA(), ( *m_Matrices )[0].GetA(),
789 ( *m_Vectors )[0], ( *m_Solutions )[0], &( IWKSP[0] ), &NW, &( WKSP[0] ), &( m_IPARM[0] ),
790 &( m_RPARM[0] ), &IERR);
791
792 m_IPARM[0] = max_num_iterations;
793
794 /* remove exception throw on convergence failure */
795 if( IERR < 100 )
796 {
797 if( ( IERR % 10 ) == 3 )
798 {
799 IERR = 0;
800 }
801 }
802
803 /* clean up */
804 delete[] IWKSP;
805 delete[] WKSP;
806
807 /* check for itpack error code */
808 if( IERR > 0 )
809 {
810 throw FEMExceptionItpackSolver(__FILE__, __LINE__, "LinearSystemWrapperItpack::Solve", IERR);
811 }
812 }
813
SwapMatrices(unsigned int matrixIndex1,unsigned int matrixIndex2)814 void LinearSystemWrapperItpack::SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2)
815 {
816 /* error checking */
817 if( !m_Matrices )
818 {
819 throw FEMExceptionLinearSystem(__FILE__,
820 __LINE__,
821 "LinearSystemWrapperItpack::SwapMatrices",
822 "No matrices allocated");
823 }
824 if( matrixIndex1 >= m_NumberOfMatrices )
825 {
826 throw FEMExceptionLinearSystemBounds(__FILE__,
827 __LINE__,
828 "LinearSystemWrapperItpack::SwapMatrices",
829 "m_Matrices",
830 matrixIndex1);
831 }
832 if( matrixIndex2 >= m_NumberOfMatrices )
833 {
834 throw FEMExceptionLinearSystemBounds(__FILE__,
835 __LINE__,
836 "LinearSystemWrapperItpack::SwapMatrices",
837 "m_Matrices",
838 matrixIndex2);
839 }
840
841 int n = ( *m_Matrices )[matrixIndex2].GetOrder();
842 int nz = ( *m_Matrices )[matrixIndex2].GetMaxNonZeroValues();
843 integer * ia = ( *m_Matrices )[matrixIndex2].GetIA();
844 integer * ja = ( *m_Matrices )[matrixIndex2].GetJA();
845 doublereal *a = ( *m_Matrices )[matrixIndex2].GetA();
846
847 ( *m_Matrices )[matrixIndex2].SetOrder( ( *m_Matrices )[matrixIndex1].GetOrder() );
848 ( *m_Matrices )[matrixIndex2].SetMaxNonZeroValues( ( *m_Matrices )[matrixIndex1].GetMaxNonZeroValues() );
849 ( *m_Matrices )[matrixIndex2].SetCompressedRow( ( *m_Matrices )[matrixIndex1].GetIA(),
850 ( *m_Matrices )[matrixIndex1].GetJA(),
851 ( *m_Matrices )[matrixIndex1].GetA() );
852
853 ( *m_Matrices )[matrixIndex1].SetOrder(n);
854 ( *m_Matrices )[matrixIndex1].SetMaxNonZeroValues(nz);
855 ( *m_Matrices )[matrixIndex1].SetCompressedRow(ia, ja, a);
856 }
857
SwapVectors(unsigned int vectorIndex1,unsigned int vectorIndex2)858 void LinearSystemWrapperItpack::SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2)
859 {
860 /* error checking */
861 if( !m_Vectors )
862 {
863 throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapVectors", "No vectors allocated");
864 }
865 if( vectorIndex1 >= m_NumberOfVectors )
866 {
867 throw FEMExceptionLinearSystemBounds(__FILE__,
868 __LINE__,
869 "LinearSystemWrapperItpack::SwapVectors",
870 "m_Vectors",
871 vectorIndex1);
872 }
873 if( vectorIndex2 >= m_NumberOfVectors )
874 {
875 throw FEMExceptionLinearSystemBounds(__FILE__,
876 __LINE__,
877 "LinearSystemWrapperItpack::SwapVectors",
878 "m_Vectors",
879 vectorIndex2);
880 }
881
882 VectorRepresentation temp = ( *m_Vectors )[vectorIndex1];
883
884 ( *m_Vectors )[vectorIndex1] = ( *m_Vectors )[vectorIndex2];
885 ( *m_Vectors )[vectorIndex2] = temp;
886 }
887
SwapSolutions(unsigned int solutionIndex1,unsigned int solutionIndex2)888 void LinearSystemWrapperItpack::SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2)
889 {
890 /* error checking */
891 if( !m_Solutions )
892 {
893 throw FEMExceptionLinearSystem(__FILE__,
894 __LINE__,
895 "LinearSystemWrapperItpack::SwapSolutions",
896 "No solutions allocated");
897 }
898 if( solutionIndex1 >= m_NumberOfSolutions )
899 {
900 throw FEMExceptionLinearSystemBounds(__FILE__,
901 __LINE__,
902 "LinearSystemWrapperItpack::SwapSolutions",
903 "m_Solutions",
904 solutionIndex1);
905 }
906 if( solutionIndex2 >= m_NumberOfSolutions )
907 {
908 throw FEMExceptionLinearSystemBounds(__FILE__,
909 __LINE__,
910 "LinearSystemWrapperItpack::SwapSolutions",
911 "m_Solutions",
912 solutionIndex2);
913 }
914
915 VectorRepresentation temp = ( *m_Solutions )[solutionIndex1];
916
917 ( *m_Solutions )[solutionIndex1] = ( *m_Solutions )[solutionIndex2];
918 ( *m_Solutions )[solutionIndex2] = temp;
919 }
920
CopySolution2Vector(unsigned solutionIndex,unsigned int vectorIndex)921 void LinearSystemWrapperItpack::CopySolution2Vector(unsigned solutionIndex,
922 unsigned int vectorIndex)
923 {
924 /* error checking */
925 if( !m_Vectors )
926 {
927 throw FEMExceptionLinearSystem(__FILE__,
928 __LINE__,
929 "LinearSystemWrapperItpack::CopySolution2Vector",
930 "No vectors allocated");
931 }
932 if( !m_Solutions )
933 {
934 throw FEMExceptionLinearSystem(__FILE__,
935 __LINE__,
936 "LinearSystemWrapperItpack::CopySolution2Vector",
937 "No solutions allocated");
938 }
939 if( vectorIndex >= m_NumberOfVectors )
940 {
941 throw FEMExceptionLinearSystemBounds(__FILE__,
942 __LINE__,
943 "LinearSystemWrapperItpack::CopySolution2Vector",
944 "m_Vectors",
945 vectorIndex);
946 }
947 if( solutionIndex >= m_NumberOfSolutions )
948 {
949 throw FEMExceptionLinearSystemBounds(__FILE__,
950 __LINE__,
951 "LinearSystemWrapperItpack::CopySolution2Vector",
952 "m_Solutions",
953 solutionIndex);
954 }
955
956 this->InitializeVector(vectorIndex);
957 /* copy values */
958 for( unsigned int i = 0; i < m_Order; i++ )
959 {
960 ( *m_Vectors )[vectorIndex][i] = ( *m_Solutions )[solutionIndex][i];
961 }
962 }
963
CopyVector2Solution(unsigned vectorIndex,unsigned int solutionIndex)964 void LinearSystemWrapperItpack::CopyVector2Solution(unsigned vectorIndex, unsigned int solutionIndex)
965 {
966 /* error checking */
967 if( !m_Vectors )
968 {
969 throw FEMExceptionLinearSystem(__FILE__,
970 __LINE__,
971 "LinearSystemWrapperItpack::CopySolution2Vector",
972 "No vectors allocated");
973 }
974 if( !m_Solutions )
975 {
976 throw FEMExceptionLinearSystem(__FILE__,
977 __LINE__,
978 "LinearSystemWrapperItpack::CopySolution2Vector",
979 "No solutions allocated");
980 }
981 if( vectorIndex >= m_NumberOfVectors )
982 {
983 throw FEMExceptionLinearSystemBounds(__FILE__,
984 __LINE__,
985 "LinearSystemWrapperItpack::CopySolution2Vector",
986 "m_Vectors",
987 vectorIndex);
988 }
989 if( solutionIndex >= m_NumberOfSolutions )
990 {
991 throw FEMExceptionLinearSystemBounds(__FILE__,
992 __LINE__,
993 "LinearSystemWrapperItpack::CopySolution2Vector",
994 "m_Solutions",
995 solutionIndex);
996 }
997
998 this->InitializeSolution(solutionIndex);
999 /* copy values */
1000 for( unsigned int i = 0; i < m_Order; i++ )
1001 {
1002 ( *m_Solutions )[solutionIndex][i] = ( *m_Vectors )[vectorIndex][i];
1003 }
1004 }
1005
MultiplyMatrixMatrix(unsigned int resultMatrixIndex,unsigned int leftMatrixIndex,unsigned int rightMatrixIndex)1006 void LinearSystemWrapperItpack::MultiplyMatrixMatrix(unsigned int resultMatrixIndex,
1007 unsigned int leftMatrixIndex,
1008 unsigned int rightMatrixIndex)
1009 {
1010 /* error checking */
1011 if( !m_Matrices )
1012 {
1013 throw FEMExceptionLinearSystem(__FILE__,
1014 __LINE__,
1015 "LinearSystemWrapperItpack::MultiplyMatrixMatrix",
1016 "No matrices allocated");
1017 }
1018 if( resultMatrixIndex >= m_NumberOfMatrices )
1019 {
1020 throw FEMExceptionLinearSystemBounds(__FILE__,
1021 __LINE__,
1022 "LinearSystemWrapperItpack::MultiplyMatrixMatrix",
1023 "m_Matrices",
1024 resultMatrixIndex);
1025 }
1026 if( leftMatrixIndex >= m_NumberOfMatrices )
1027 {
1028 throw FEMExceptionLinearSystemBounds(__FILE__,
1029 __LINE__,
1030 "LinearSystemWrapperItpack::MultiplyMatrixMatrix",
1031 "m_Matrices",
1032 leftMatrixIndex);
1033 }
1034 if( rightMatrixIndex >= m_NumberOfMatrices )
1035 {
1036 throw FEMExceptionLinearSystemBounds(__FILE__,
1037 __LINE__,
1038 "LinearSystemWrapperItpack::MultiplyMatrixMatrix",
1039 "m_Matrices",
1040 rightMatrixIndex);
1041 }
1042
1043 ( *m_Matrices )[leftMatrixIndex].mult( &( ( *m_Matrices )[rightMatrixIndex] ), &( ( *m_Matrices )[resultMatrixIndex] ) );
1044 }
1045
MultiplyMatrixVector(unsigned int resultVectorIndex,unsigned int matrixIndex,unsigned int vectorIndex)1046 void LinearSystemWrapperItpack::MultiplyMatrixVector(unsigned int resultVectorIndex,
1047 unsigned int matrixIndex,
1048 unsigned int vectorIndex)
1049 {
1050 /* error checking */
1051 if( !m_Matrices )
1052 {
1053 throw FEMExceptionLinearSystem(__FILE__,
1054 __LINE__,
1055 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1056 "No matrices allocated");
1057 }
1058 if( !m_Vectors )
1059 {
1060 throw FEMExceptionLinearSystem(__FILE__,
1061 __LINE__,
1062 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1063 "No vectors allocated");
1064 }
1065 if( resultVectorIndex >= m_NumberOfVectors )
1066 {
1067 throw FEMExceptionLinearSystemBounds(__FILE__,
1068 __LINE__,
1069 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1070 "m_Vectors",
1071 resultVectorIndex);
1072 }
1073 if( matrixIndex >= m_NumberOfMatrices )
1074 {
1075 throw FEMExceptionLinearSystemBounds(__FILE__,
1076 __LINE__,
1077 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1078 "m_Matrices",
1079 matrixIndex);
1080 }
1081 if( vectorIndex >= m_NumberOfVectors )
1082 {
1083 throw FEMExceptionLinearSystemBounds(__FILE__,
1084 __LINE__,
1085 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1086 "m_Vectors",
1087 vectorIndex);
1088 }
1089
1090 /* perform mult */
1091 ( *m_Matrices )[matrixIndex].mult( ( *m_Vectors )[vectorIndex], ( *m_Vectors )[resultVectorIndex] );
1092 }
1093
1094
MultiplyMatrixSolution(unsigned int resultVectorIndex,unsigned int matrixIndex,unsigned int solutionIndex)1095 void LinearSystemWrapperItpack::MultiplyMatrixSolution(unsigned int resultVectorIndex,
1096 unsigned int matrixIndex,
1097 unsigned int solutionIndex)
1098 {
1099 /* error checking */
1100 if( !m_Matrices )
1101 {
1102 throw FEMExceptionLinearSystem(__FILE__,
1103 __LINE__,
1104 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1105 "No matrices allocated");
1106 }
1107
1108 if( !m_Vectors )
1109 {
1110 throw FEMExceptionLinearSystem(__FILE__,
1111 __LINE__,
1112 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1113 "No vectors allocated");
1114 }
1115
1116 if( !m_Solutions)
1117 {
1118 throw FEMExceptionLinearSystem(__FILE__,
1119 __LINE__,
1120 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1121 "No solutions allocated");
1122 }
1123
1124 if( resultVectorIndex >= m_NumberOfVectors )
1125 {
1126 throw FEMExceptionLinearSystemBounds(__FILE__,
1127 __LINE__,
1128 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1129 "m_Vectors",
1130 resultVectorIndex);
1131 }
1132
1133 if( matrixIndex >= m_NumberOfMatrices )
1134 {
1135 throw FEMExceptionLinearSystemBounds(__FILE__,
1136 __LINE__,
1137 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1138 "m_Matrices",
1139 matrixIndex);
1140 }
1141
1142 if( solutionIndex >= m_NumberOfSolutions )
1143 {
1144 throw FEMExceptionLinearSystemBounds(__FILE__,
1145 __LINE__,
1146 "LinearSystemWrapperItpack::MultiplyMatrixVector",
1147 "m_Solutions",
1148 solutionIndex);
1149 }
1150
1151 /* perform multiplication */
1152 ( *m_Matrices )[matrixIndex].mult( ( *m_Solutions )[solutionIndex], ( *m_Vectors )[resultVectorIndex] );
1153 }
1154
1155
~LinearSystemWrapperItpack()1156 LinearSystemWrapperItpack::~LinearSystemWrapperItpack()
1157 {
1158 delete m_Matrices;
1159
1160 unsigned int i;
1161 if( m_Vectors )
1162 {
1163 for( i = 0; i < m_NumberOfVectors; i++ )
1164 {
1165 delete[] ( *m_Vectors )[i];
1166 }
1167 delete m_Vectors;
1168 }
1169
1170 if( m_Solutions )
1171 {
1172 for( i = 0; i < m_NumberOfSolutions; i++ )
1173 {
1174 delete[] ( *m_Solutions )[i];
1175 }
1176 delete m_Solutions;
1177 }
1178 }
1179
FEMExceptionItpackSolver(const char * file,unsigned int lineNumber,std::string location,integer errorCode)1180 FEMExceptionItpackSolver::FEMExceptionItpackSolver(const char *file, unsigned int lineNumber, std::string location,
1181 integer errorCode) :
1182 FEMException(file, lineNumber)
1183 {
1184 std::string solverError;
1185
1186 if( errorCode < 100 )
1187 {
1188 errorCode = errorCode % 10;
1189 }
1190
1191 switch( errorCode )
1192 {
1193 case 1:
1194 solverError = "Invalid order of system";
1195 break;
1196 case 2:
1197 solverError = "Workspace is not large enough";
1198 break;
1199 case 3:
1200 solverError = "Failure to converge before reaching maximum number of iterations";
1201 break;
1202 case 4:
1203 solverError = "Invalid order of black subsystem";
1204 break;
1205 case 101:
1206 solverError = "A diagonal element is not positive";
1207 break;
1208 case 102:
1209 solverError = "No diagonal element in a row";
1210 break;
1211 case 201:
1212 solverError = "Red-black indexing is not possible";
1213 break;
1214 case 301:
1215 solverError = "No entry in a row of the original matrix";
1216 break;
1217 case 302:
1218 solverError = "No entry in a row of the permuted matrix";
1219 break;
1220 case 303:
1221 solverError = "Sorting error in a row of the permuted matrix";
1222 break;
1223 case 401:
1224 solverError = "A diagonal element is not positive";
1225 break;
1226 case 402:
1227 solverError = "No diagonal element in a row";
1228 break;
1229 case 501:
1230 solverError = "Failure to converge before reaching maximum number of iterations";
1231 break;
1232 case 502:
1233 solverError = "Function does not change sign at endpoints";
1234 break;
1235 case 601:
1236 solverError = "Successive iterations are not monotone increasing";
1237 break;
1238 default:
1239 solverError = "Unknown error code returned";
1240 }
1241
1242 std::ostringstream buf;
1243 buf << "Error: " << solverError;
1244
1245 SetDescription( buf.str().c_str() );
1246
1247 SetLocation(location);
1248 }
1249
1250 } // end namespace fem
1251 } // end namespace itk
1252