ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkFEMLinearSystemWrapperItpack.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 
00019 #ifndef __itkFEMLinearSystemWrapperItpack_h
00020 #define __itkFEMLinearSystemWrapperItpack_h
00021 
00022 #include "itkFEMSolution.h"
00023 #include "itkFEMLinearSystemWrapper.h"
00024 #include "itkFEMItpackSparseMatrix.h"
00025 #include <vector>
00026 
00030 typedef long   integer;
00031 typedef double doublereal;
00032 
00033 extern "C" {
00034 typedef
00035 int ( *ItkItpackSolverFunction )(integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *,
00036                                  integer *, doublereal *,
00037                                  integer *, doublereal *, integer *);
00038 }
00039 
00040 namespace itk
00041 {
00042 namespace fem
00043 {
00051 class LinearSystemWrapperItpack : public LinearSystemWrapper
00052 {
00053 public:
00054 
00056   typedef LinearSystemWrapperItpack Self;
00057 
00059   typedef LinearSystemWrapper Superclass;
00060 
00062   typedef ItpackSparseMatrix MatrixRepresentation;
00063 
00065   typedef std::vector<MatrixRepresentation> MatrixHolder;
00066 
00067   /* auto pointer to vector of matrices typedef */
00068   /* typedef std::auto_ptr<MatrixHolder> MatrixArrayPtr; */
00069 
00071   /* typedef std::auto_ptr<double> VectorRepresentation; */
00072   typedef double *VectorRepresentation;
00073 
00075   typedef std::vector<VectorRepresentation> VectorHolder;
00076 
00077   /* auto pointer to vector of vectors typedef */
00078   /* typedef std::auto_ptr<VectorHolder> VectorArrayPtr; */
00079 
00080   /* pointer to array of unsigned int typedef */
00081   /* typedef std::auto_ptr<unsigned int> UnsignedIntegerArrayPtr; */
00082 
00083   /* -----------------------------------------------------------------
00084    *
00085    * Routines for setting/reporting itpack parameters
00086    *
00087    * -----------------------------------------------------------------
00088    */
00089 
00094   void SetMaximumNumberIterations(int i)
00095   {
00096     m_IPARM[0] = i;
00097   }
00098 
00102   int  GetMaximumNumberIterations() const
00103   {
00104     return m_IPARM[0];
00105   }
00106 
00107   // void SetErrorReportingLevel(int i) {   m_IPARM[1] = i; }
00108 
00112   int  GetErrorReportingLevel() const
00113   {
00114     return m_IPARM[1];
00115   }
00116 
00121   void SetCommunicationSwitch(int i)
00122   {
00123     m_IPARM[2] = i;
00124   }
00125 
00129   int  GetCommunicationSwitch() const
00130   {
00131     return m_IPARM[2];
00132   }
00133 
00134   // void SetOutputNumber(int i) {   m_IPARM[3] = i; }
00135 
00139   int  GetOutputNumber() const
00140   {
00141     return m_IPARM[3];
00142   }
00143 
00148   void SetSymmetricMatrixFlag(int i)
00149   {
00150     m_IPARM[4] = i;
00151   }
00152 
00156   int  GetSymmetricMatrixFlag()
00157   {
00158     return m_IPARM[4];
00159   }
00160 
00165   void SetAdaptiveSwitch(int i)
00166   {
00167     m_IPARM[5] = i;
00168   }
00169 
00173   int  GetAdaptiveSwitch() const
00174   {
00175     return m_IPARM[5];
00176   }
00177 
00182   void SetAdaptiveCaseSwitch(int i)
00183   {
00184     m_IPARM[6] = i;
00185   }
00186 
00190   int  GetAdaptiveCaseSwitch() const
00191   {
00192     return m_IPARM[6];
00193   }
00194 
00200   void SetWorkspaceUsed(int i)
00201   {
00202     m_IPARM[7] = i;
00203   }
00204 
00209   int  GetWorkspaceUsed()
00210   {
00211     return m_IPARM[7];
00212   }
00213 
00218   void SetRedBlackOrderingSwitch(int i)
00219   {
00220     m_IPARM[8] = i;
00221   }
00222 
00226   int  GetRedBlackOrderingSwitch()
00227   {
00228     return m_IPARM[8];
00229   }
00230 
00235   void SetRemoveSwitch(int i)
00236   {
00237     m_IPARM[9] = i;
00238   }
00239 
00243   int  GetRemoveSwitch()
00244   {
00245     return m_IPARM[9];
00246   }
00247 
00252   void SetTimingSwitch(int i)
00253   {
00254     m_IPARM[10] = i;
00255   }
00256 
00260   int  GetTimingSwitch()
00261   {
00262     return m_IPARM[10];
00263   }
00264 
00269   void SetErrorAnalysisSwitch(int i)
00270   {
00271     m_IPARM[11] = i;
00272   }
00273 
00277   int  GetErrorAnalysisSwitch() const
00278   {
00279     return m_IPARM[11];
00280   }
00281 
00286   void   SetAccuracy(double i)
00287   {
00288     m_RPARM[0] = i;
00289   }
00290 
00294   double GetAccuracy() const
00295   {
00296     return m_RPARM[0];
00297   }
00298 
00303   void   SetLargestJacobiEigenvalueEstimate(double i)
00304   {
00305     m_RPARM[1] = i;
00306   }
00307 
00311   double GetLargestJacobiEigenvalueEstimate() const
00312   {
00313     return m_RPARM[1];
00314   }
00315 
00320   void   SetSmallestJacobiEigenvalueEstimate(double i)
00321   {
00322     m_RPARM[2] = i;
00323   }
00324 
00328   double GetSmallestJacobiEigenvalueEstimate()
00329   {
00330     return m_RPARM[2];
00331   }
00332 
00337   void   SetDampingFactor(double i)
00338   {
00339     m_RPARM[3] = i;
00340   }
00341 
00345   double GetDampingFactor() const
00346   {
00347     return m_RPARM[3];
00348   }
00349 
00354   void   SetOverrelaxationParameter(double i)
00355   {
00356     m_RPARM[4] = i;
00357   }
00358 
00362   double GetOverrelaxationParameter()
00363   {
00364     return m_RPARM[4];
00365   }
00366 
00371   void   SetEstimatedSpectralRadiusSSOR(double i)
00372   {
00373     m_RPARM[5] = i;
00374   }
00375 
00379   double GetEstimatedSpectralRadiusSSOR() const
00380   {
00381     return m_RPARM[5];
00382   }
00383 
00388   void   SetEstimatedSpectralRadiusLU(double i)
00389   {
00390     m_RPARM[6] = i;
00391   }
00392 
00396   double GetEstimatedSpectralRadiusLU() const
00397   {
00398     return m_RPARM[6];
00399   }
00400 
00405   void   SetTolerance(double i)
00406   {
00407     m_RPARM[7] = i;
00408   }
00409 
00413   double GetTolerance()
00414   {
00415     return m_RPARM[7];
00416   }
00417 
00422   void   SetTimeToConvergence(double i)
00423   {
00424     m_RPARM[8] = i;
00425   }
00426 
00430   double GetTimeToConvergence()
00431   {
00432     return m_RPARM[8];
00433   }
00434 
00439   void   SetTimeForCall(double i)
00440   {
00441     m_RPARM[9] = i;
00442   }
00443 
00447   double GetTimeForCall()
00448   {
00449     return m_RPARM[9];
00450   }
00451 
00456   void   SetDigitsInError(double i)
00457   {
00458     m_RPARM[10] = i;
00459   }
00460 
00464   double GetDigitsInError() const
00465   {
00466     return m_RPARM[10];
00467   }
00468 
00473   void   SetDigitsInResidual(double i)
00474   {
00475     m_RPARM[11] = i;
00476   }
00477 
00481   double GetDigitsInResidual() const
00482   {
00483     return m_RPARM[11];
00484   }
00485 
00489   void JacobianConjugateGradient()
00490   {
00491     m_Method = 0;
00492   }
00493 
00497   void JacobianSemiIterative()
00498   {
00499     m_Method = 1;
00500   }
00501 
00505   void SuccessiveOverrelaxation()
00506   {
00507     m_Method = 2;
00508   }
00509 
00514   void SymmetricSuccessiveOverrelaxationConjugateGradient()
00515   {
00516     m_Method = 3;
00517   }
00518 
00523   void SymmetricSuccessiveOverrelaxationSuccessiveOverrelaxation()
00524   {
00525     m_Method = 4;
00526   }
00527 
00531   void ReducedSystemConjugateGradient()
00532   {
00533     m_Method = 5;
00534   }
00535 
00538   void ReducedSystemSemiIteration()
00539   {
00540     m_Method = 6;
00541   }
00542 
00555   virtual void SetMaximumNonZeroValuesInMatrix(unsigned int maxNonZeroValues)
00556   {
00557     m_MaximumNonZeroValues =
00558       maxNonZeroValues;
00559   }
00560 
00561   void ScaleMatrix(Float scale, unsigned int matrixIndex);
00562 
00573   LinearSystemWrapperItpack();
00574 
00578   ~LinearSystemWrapperItpack();
00579 
00580   /* memory management routines */
00581   virtual void  InitializeMatrix(unsigned int matrixIndex);
00582 
00583   virtual bool  IsMatrixInitialized(unsigned int matrixIndex);
00584 
00585   virtual void  DestroyMatrix(unsigned int matrixIndex);
00586 
00587   virtual void  InitializeVector(unsigned int vectorIndex);
00588 
00589   virtual bool  IsVectorInitialized(unsigned int vectorIndex);
00590 
00591   virtual void  DestroyVector(unsigned int vectorIndex);
00592 
00593   virtual void  InitializeSolution(unsigned int solutionIndex);
00594 
00595   virtual bool  IsSolutionInitialized(unsigned int solutionIndex);
00596 
00597   virtual void  DestroySolution(unsigned int solutionIndex);
00598 
00599   /* assembly & solving routines */
00600   virtual Float GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex) const;
00601 
00602   virtual void  SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex);
00603 
00604   virtual void  AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex);
00605 
00606   virtual void GetColumnsOfNonZeroMatrixElementsInRow(unsigned int row, ColumnArray & cols, unsigned int matrixIndex);
00607 
00608   virtual Float GetVectorValue(unsigned int i, unsigned int vectorIndex) const;
00609 
00610   virtual void  SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex);
00611 
00612   virtual void  AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex);
00613 
00614   virtual Float GetSolutionValue(unsigned int i, unsigned int solutionIndex) const;
00615 
00616   virtual void  SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex);
00617 
00618   virtual void  AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex);
00619 
00620   virtual void  Solve(void);
00621 
00622   /* matrix & vector manipulation routines */
00623   virtual void  SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2);
00624 
00625   virtual void  SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2);
00626 
00627   virtual void  SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2);
00628 
00629   virtual void  CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex);
00630 
00631   virtual void  CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex);
00632 
00633   virtual void  MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex,
00634                                      unsigned int rightMatrixIndex);
00635 
00636   virtual void  MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);
00637 
00638 private:
00639 
00641   MatrixHolder *m_Matrices;
00642 
00644   VectorHolder *m_Vectors;
00645 
00647   VectorHolder *m_Solutions;
00648 
00651   // UnsignedIntegerArrayPtr m_MaximumNonZeroValues;
00652   unsigned int m_MaximumNonZeroValues;
00653 
00655   ItkItpackSolverFunction m_Methods[7];
00656 
00658   integer m_Method;
00659 
00661   integer m_IPARM[12];
00662 
00664   doublereal m_RPARM[12];
00665 };
00666 
00674 class ITK_ABI_EXPORT FEMExceptionItpackSolver : public FEMException
00675 {
00676 public:
00677 
00679   typedef long integer;
00680 
00686   FEMExceptionItpackSolver(const char *file, unsigned int lineNumber, std::string location, integer errorCode);
00687 
00689   virtual ~FEMExceptionItpackSolver()
00690   throw ( )
00691   {
00692   }
00693 
00695   itkTypeMacro(FEMExceptionItpackSolver, FEMException);
00696 };
00697 }
00698 }  // end namespace itk::fem
00699 
00700 #endif // #ifndef __itkFEMLinearSystemWrapperItpack_h
00701