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
00029 typedef long integer;
00030 typedef double doublereal;
00031
00032 extern "C" {
00033 typedef
00034 int (*ItkItpackSolverFunction)(integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *, integer *, doublereal *,
00035 integer *, doublereal *, integer *);
00036 }
00037
00038
00039 namespace itk {
00040 namespace fem {
00041
00048 class LinearSystemWrapperItpack : public LinearSystemWrapper
00049 {
00050 public:
00051
00053 typedef LinearSystemWrapperItpack Self;
00054
00056 typedef LinearSystemWrapper Superclass;
00057
00059 typedef ItpackSparseMatrix MatrixRepresentation;
00060
00062 typedef std::vector<MatrixRepresentation> MatrixHolder;
00063
00064
00065
00066
00068
00069 typedef double * VectorRepresentation;
00070
00072 typedef std::vector<VectorRepresentation> VectorHolder;
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00091 void SetMaximumNumberIterations(int i) { m_IPARM[0] = i; }
00092
00096 int GetMaximumNumberIterations() { return m_IPARM[0]; }
00097
00098
00099
00103 int GetErrorReportingLevel() { return m_IPARM[1]; }
00104
00109 void SetCommunicationSwitch(int i) { m_IPARM[2] = i; }
00110
00114 int GetCommunicationSwitch() { return m_IPARM[2]; }
00115
00116
00117
00121 int GetOutputNumber() { return m_IPARM[3]; }
00122
00127 void SetSymmetricMatrixFlag(int i) { m_IPARM[4] = i; }
00128
00132 int GetSymmetricMatrixFlag() { return m_IPARM[4]; }
00133
00138 void SetAdaptiveSwitch(int i) { m_IPARM[5] = i; }
00139
00143 int GetAdaptiveSwitch() { return m_IPARM[5]; }
00144
00149 void SetAdaptiveCaseSwitch(int i) { m_IPARM[6] = i; }
00150
00154 int GetAdaptiveCaseSwitch() { return m_IPARM[6]; }
00155
00161 void SetWorkspaceUsed(int i) { m_IPARM[7] = i; }
00162
00167 int GetWorkspaceUsed() { return m_IPARM[7]; }
00168
00173 void SetRedBlackOrderingSwitch(int i) { m_IPARM[8] = i; }
00174
00178 int GetRedBlackOrderingSwitch() { return m_IPARM[8]; }
00179
00184 void SetRemoveSwitch(int i) { m_IPARM[9] = i; }
00185
00189 int GetRemoveSwitch() { return m_IPARM[9]; }
00190
00195 void SetTimingSwitch(int i) { m_IPARM[10] = i; }
00196
00200 int GetTimingSwitch() { return m_IPARM[10]; }
00201
00206 void SetErrorAnalysisSwitch(int i) { m_IPARM[11] = i; }
00207
00211 int GetErrorAnalysisSwitch() { return m_IPARM[11]; }
00212
00217 void SetAccuracy(double i) { m_RPARM[0] = i; }
00218
00222 double GetAccuracy() { return m_RPARM[0]; }
00223
00228 void SetLargestJacobiEigenvalueEstimate(double i) { m_RPARM[1] = i; }
00229
00233 double GetLargestJacobiEigenvalueEstimate() { return m_RPARM[1]; }
00234
00239 void SetSmallestJacobiEigenvalueEstimate(double i) { m_RPARM[2] = i; }
00240
00244 double GetSmallestJacobiEigenvalueEstimate() { return m_RPARM[2]; }
00245
00250 void SetDampingFactor(double i) { m_RPARM[3] = i; }
00251
00255 double GetDampingFactor() { return m_RPARM[3]; }
00256
00261 void SetOverrelaxationParameter(double i) { m_RPARM[4] = i; }
00262
00266 double GetOverrelaxationParameter() { return m_RPARM[4]; }
00267
00272 void SetEstimatedSpectralRadiusSSOR(double i) { m_RPARM[5] = i; }
00273
00277 double GetEstimatedSpectralRadiusSSOR() { return m_RPARM[5]; }
00278
00283 void SetEstimatedSpectralRadiusLU(double i) { m_RPARM[6] = i; }
00284
00288 double GetEstimatedSpectralRadiusLU() { return m_RPARM[6]; }
00289
00294 void SetTolerance(double i) { m_RPARM[7] = i; }
00295
00299 double GetTolerance() { return m_RPARM[7]; }
00300
00305 void SetTimeToConvergence(double i) { m_RPARM[8] = i; }
00306
00310 double GetTimeToConvergence() { return m_RPARM[8]; }
00311
00316 void SetTimeForCall(double i) { m_RPARM[9] = i; }
00317
00321 double GetTimeForCall() { return m_RPARM[9]; }
00322
00327 void SetDigitsInError(double i) { m_RPARM[10] = i; }
00328
00332 double GetDigitsInError() { return m_RPARM[10]; }
00333
00338 void SetDigitsInResidual(double i) { m_RPARM[11] = i; }
00339
00343 double GetDigitsInResidual() { return m_RPARM[11]; }
00344
00348 void JacobianConjugateGradient() { m_Method = 0; }
00349
00353 void JacobianSemiIterative() { m_Method = 1; }
00354
00358 void SuccessiveOverrelaxation() { m_Method = 2; }
00359
00364 void SymmetricSuccessiveOverrelaxationConjugateGradient() { m_Method = 3; }
00365
00370 void SymmetricSuccessiveOverrelaxationSuccessiveOverrelaxation() { m_Method = 4; }
00371
00375 void ReducedSystemConjugateGradient() { m_Method = 5; }
00376
00379 void ReducedSystemSemiIteration() { m_Method = 6; }
00380
00381
00394 virtual void SetMaximumNonZeroValuesInMatrix(unsigned int maxNonZeroValues) {m_MaximumNonZeroValues = maxNonZeroValues;}
00395
00396
00397 void ScaleMatrix(Float scale, unsigned int matrixIndex);
00398
00399
00410 LinearSystemWrapperItpack();
00411
00415 ~LinearSystemWrapperItpack();
00416
00417
00418
00419 virtual void InitializeMatrix(unsigned int matrixIndex);
00420
00421 virtual bool IsMatrixInitialized(unsigned int matrixIndex);
00422
00423 virtual void DestroyMatrix(unsigned int matrixIndex);
00424
00425 virtual void InitializeVector(unsigned int vectorIndex);
00426
00427 virtual bool IsVectorInitialized(unsigned int vectorIndex);
00428
00429 virtual void DestroyVector(unsigned int vectorIndex);
00430
00431 virtual void InitializeSolution(unsigned int solutionIndex);
00432
00433 virtual bool IsSolutionInitialized(unsigned int solutionIndex);
00434
00435 virtual void DestroySolution(unsigned int solutionIndex);
00436
00437
00438 virtual Float GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex) const;
00439
00440 virtual void SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex);
00441
00442 virtual void AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex);
00443
00444 virtual void GetColumnsOfNonZeroMatrixElementsInRow( unsigned int row, ColumnArray& cols, unsigned int matrixIndex );
00445
00446 virtual Float GetVectorValue(unsigned int i, unsigned int vectorIndex) const;
00447
00448 virtual void SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex);
00449
00450 virtual void AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex);
00451
00452 virtual Float GetSolutionValue(unsigned int i, unsigned int solutionIndex) const;
00453
00454 virtual void SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex);
00455
00456 virtual void AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex);
00457
00458 virtual void Solve(void);
00459
00460
00461
00462 virtual void SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2);
00463
00464 virtual void SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2);
00465
00466 virtual void SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2);
00467
00468 virtual void CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex);
00469
00470 virtual void CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex);
00471
00472 virtual void MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex, unsigned int rightMatrixIndex);
00473
00474 virtual void MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);
00475
00476 private:
00477
00479 MatrixHolder *m_Matrices;
00480
00482 VectorHolder *m_Vectors;
00483
00485 VectorHolder *m_Solutions;
00486
00488
00489 unsigned int m_MaximumNonZeroValues;
00490
00492 ItkItpackSolverFunction m_Methods[7];
00493
00495 integer m_Method;
00496
00498 integer m_IPARM[12];
00499
00501 doublereal m_RPARM[12];
00502
00503 };
00504
00511 class FEMExceptionItpackSolver : public FEMException
00512 {
00513 public:
00514
00516 typedef long integer;
00517
00523 FEMExceptionItpackSolver(const char *file, unsigned int lineNumber, std::string location, integer errorCode);
00524
00526 virtual ~FEMExceptionItpackSolver() throw() {}
00527
00529 itkTypeMacro(FEMExceptionItpackSolver,FEMException);
00530
00531 };
00532 }}
00533
00534 #endif // #ifndef __itkFEMLinearSystemWrapperItpack_h
00535