00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef __itkFEMLinearSystemWrapperItpack_h
00019 #define __itkFEMLinearSystemWrapperItpack_h
00020
00021 #include "itkFEMSolution.h"
00022 #include "itkFEMLinearSystemWrapper.h"
00023 #include "itkFEMItpackSparseMatrix.h"
00024 #include <vector>
00025
00026
00027 namespace itk {
00028 namespace fem {
00029
00036 class LinearSystemWrapperItpack : public LinearSystemWrapper
00037 {
00038 public:
00039
00041 typedef LinearSystemWrapperItpack Self;
00042
00044 typedef LinearSystemWrapper Superclass;
00045
00047 typedef ItpackSparseMatrix MatrixRepresentation;
00048
00050 typedef std::vector<MatrixRepresentation> MatrixHolder;
00051
00052
00053
00054
00056
00057 typedef double * VectorRepresentation;
00058
00060 typedef std::vector<VectorRepresentation> VectorHolder;
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00079 void SetMaximumNumberIterations(int i) { m_IPARM[0] = i; }
00080
00084 int GetMaximumNumberIterations() { return m_IPARM[0]; }
00085
00086
00087
00091 int GetErrorReportingLevel() { return m_IPARM[1]; }
00092
00097 void SetCommunicationSwitch(int i) { m_IPARM[2] = i; }
00098
00102 int GetCommunicationSwitch() { return m_IPARM[2]; }
00103
00104
00105
00109 int GetOutputNumber() { return m_IPARM[3]; }
00110
00115 void SetSymmetricMatrixFlag(int i) { m_IPARM[4] = i; }
00116
00120 int GetSymmetricMatrixFlag() { return m_IPARM[4]; }
00121
00126 void SetAdaptiveSwitch(int i) { m_IPARM[5] = i; }
00127
00131 int GetAdaptiveSwitch() { return m_IPARM[5]; }
00132
00137 void SetAdaptiveCaseSwitch(int i) { m_IPARM[6] = i; }
00138
00142 int GetAdaptiveCaseSwitch() { return m_IPARM[6]; }
00143
00149 void SetWorkspaceUsed(int i) { m_IPARM[7] = i; }
00150
00155 int GetWorkspaceUsed() { return m_IPARM[7]; }
00156
00161 void SetRedBlackOrderingSwitch(int i) { m_IPARM[8] = i; }
00162
00166 int GetRedBlackOrderingSwitch() { return m_IPARM[8]; }
00167
00172 void SetRemoveSwitch(int i) { m_IPARM[9] = i; }
00173
00177 int GetRemoveSwitch() { return m_IPARM[9]; }
00178
00183 void SetTimingSwitch(int i) { m_IPARM[10] = i; }
00184
00188 int GetTimingSwitch() { return m_IPARM[10]; }
00189
00194 void SetErrorAnalysisSwitch(int i) { m_IPARM[11] = i; }
00195
00199 int GetErrorAnalysisSwitch() { return m_IPARM[11]; }
00200
00205 void SetAccuracy(double i) { m_RPARM[0] = i; }
00206
00210 double GetAccuracy() { return m_RPARM[0]; }
00211
00216 void SetLargestJacobiEigenvalueEstimate(double i) { m_RPARM[1] = i; }
00217
00221 double GetLargestJacobiEigenvalueEstimate() { return m_RPARM[1]; }
00222
00227 void SetSmallestJacobiEigenvalueEstimate(double i) { m_RPARM[2] = i; }
00228
00232 double GetSmallestJacobiEigenvalueEstimate() { return m_RPARM[2]; }
00233
00238 void SetDampingFactor(double i) { m_RPARM[3] = i; }
00239
00243 double GetDampingFactor() { return m_RPARM[3]; }
00244
00249 void SetOverrelaxationParameter(double i) { m_RPARM[4] = i; }
00250
00254 double GetOverrelaxationParameter() { return m_RPARM[4]; }
00255
00260 void SetEstimatedSpectralRadiusSSOR(double i) { m_RPARM[5] = i; }
00261
00265 double GetEstimatedSpectralRadiusSSOR() { return m_RPARM[5]; }
00266
00271 void SetEstimatedSpectralRadiusLU(double i) { m_RPARM[6] = i; }
00272
00276 double GetEstimatedSpectralRadiusLU() { return m_RPARM[6]; }
00277
00282 void SetTolerance(double i) { m_RPARM[7] = i; }
00283
00287 double GetTolerance() { return m_RPARM[7]; }
00288
00293 void SetTimeToConvergence(double i) { m_RPARM[8] = i; }
00294
00298 double GetTimeToConvergence() { return m_RPARM[8]; }
00299
00304 void SetTimeForCall(double i) { m_RPARM[9] = i; }
00305
00309 double GetTimeForCall() { return m_RPARM[9]; }
00310
00315 void SetDigitsInError(double i) { m_RPARM[10] = i; }
00316
00320 double GetDigitsInError() { return m_RPARM[10]; }
00321
00326 void SetDigitsInResidual(double i) { m_RPARM[11] = i; }
00327
00331 double GetDigitsInResidual() { return m_RPARM[11]; }
00332
00336 void JacobianConjugateGradient() { m_Method = 0; }
00337
00341 void JacobianSemiIterative() { m_Method = 1; }
00342
00346 void SuccessiveOverrelaxation() { m_Method = 2; }
00347
00352 void SymmetricSuccessiveOverrelaxationConjugateGradient() { m_Method = 3; }
00353
00358 void SymmetricSuccessiveOverrelaxationSuccessiveOverrelaxation() { m_Method = 4; }
00359
00363 void ReducedSystemConjugateGradient() { m_Method = 5; }
00364
00367 void ReducedSystemSemiIteration() { m_Method = 6; }
00368
00369
00370
00371
00372
00373
00374
00375
00376
00382 virtual void SetMaximumNonZeroValuesInMatrix(unsigned int maxNonZeroValues) {m_MaximumNonZeroValues = maxNonZeroValues;}
00383
00384
00385 void ScaleMatrix(Float scale, unsigned int matrixIndex);
00386
00387
00388
00389
00390
00391
00392
00393
00394
00398 LinearSystemWrapperItpack();
00399
00403 ~LinearSystemWrapperItpack();
00404
00405
00406
00407 virtual void InitializeMatrix(unsigned int matrixIndex);
00408
00409 virtual bool IsMatrixInitialized(unsigned int matrixIndex);
00410
00411 virtual void DestroyMatrix(unsigned int matrixIndex);
00412
00413 virtual void InitializeVector(unsigned int vectorIndex);
00414
00415 virtual bool IsVectorInitialized(unsigned int vectorIndex);
00416
00417 virtual void DestroyVector(unsigned int vectorIndex);
00418
00419 virtual void InitializeSolution(unsigned int solutionIndex);
00420
00421 virtual bool IsSolutionInitialized(unsigned int solutionIndex);
00422
00423 virtual void DestroySolution(unsigned int solutionIndex);
00424
00425
00426 virtual Float GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex) const;
00427
00428 virtual void SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex);
00429
00430 virtual void AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex);
00431
00432 virtual void GetColumnsOfNonZeroMatrixElementsInRow( unsigned int row, ColumnArray& cols, unsigned int matrixIndex );
00433
00434 virtual Float GetVectorValue(unsigned int i, unsigned int vectorIndex) const;
00435
00436 virtual void SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex);
00437
00438 virtual void AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex);
00439
00440 virtual Float GetSolutionValue(unsigned int i, unsigned int solutionIndex) const;
00441
00442 virtual void SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex);
00443
00444 virtual void AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex);
00445
00446 virtual void Solve(void);
00447
00448
00449
00450 virtual void SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2);
00451
00452 virtual void SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2);
00453
00454 virtual void SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2);
00455
00456 virtual void CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex);
00457
00458 virtual void CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex);
00459
00460 virtual void MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex, unsigned int rightMatrixIndex);
00461
00462 virtual void MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);
00463
00464 private:
00465
00466
00468 typedef long integer;
00469 typedef double doublereal;
00470
00472 MatrixHolder *m_Matrices;
00473
00475 VectorHolder *m_Vectors;
00476
00478 VectorHolder *m_Solutions;
00479
00481
00482 unsigned int m_MaximumNonZeroValues;
00483
00485 int (*m_Methods[7])(integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *, integer *, doublereal *,
00486 integer *, doublereal *, integer *);
00487
00489 integer m_Method;
00490
00492 integer m_IPARM[12];
00493
00495 doublereal m_RPARM[12];
00496
00497 };
00498
00505 class FEMExceptionItpackSolver : public FEMException
00506 {
00507 public:
00508
00510 typedef long integer;
00511
00517 FEMExceptionItpackSolver(const char *file, unsigned int lineNumber, std::string location, integer errorCode);
00518
00520 virtual ~FEMExceptionItpackSolver() throw() {}
00521
00523 itkTypeMacro(FEMExceptionItpackSolver,FEMException);
00524
00525 };
00526
00527
00528
00529 }}
00530
00531 #endif // #ifndef __itkFEMLinearSystemWrapperItpack_h
00532
00533
00534