itkFEMLinearSystemWrapperItpack.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMLinearSystemWrapperItpack.h,v $
  Language:  C++
  Date:      $Date: 2009-01-30 21:10:11 $
  Version:   $Revision: 1.15 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

#ifndef __itkFEMLinearSystemWrapperItpack_h
#define __itkFEMLinearSystemWrapperItpack_h

#include "itkFEMSolution.h"
#include "itkFEMLinearSystemWrapper.h"
#include "itkFEMItpackSparseMatrix.h"
#include <vector>

typedef long      integer;
typedef double    doublereal;

extern "C" {
  typedef  
  int (*ItkItpackSolverFunction)(integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *, integer *, doublereal *, 
                                 integer *, doublereal *, integer *);
}


namespace itk {
namespace fem {

class LinearSystemWrapperItpack : public LinearSystemWrapper
{
public:

  typedef LinearSystemWrapperItpack Self;

  typedef LinearSystemWrapper Superclass;

  typedef ItpackSparseMatrix MatrixRepresentation;

  typedef std::vector<MatrixRepresentation> MatrixHolder;

  /* auto pointer to vector of matrices typedef */
  /* typedef std::auto_ptr<MatrixHolder> MatrixArrayPtr; */

  /* typedef std::auto_ptr<double> VectorRepresentation; */
  typedef double * VectorRepresentation;

  typedef std::vector<VectorRepresentation> VectorHolder;

  /* auto pointer to vector of vectors typedef */
  /* typedef std::auto_ptr<VectorHolder> VectorArrayPtr; */

  /* pointer to array of unsigned int typedef */
  /* typedef std::auto_ptr<unsigned int> UnsignedIntegerArrayPtr; */

  /* -----------------------------------------------------------------
   * 
   * Routines for setting/reporting itpack parameters 
   *
   * ----------------------------------------------------------------- 
   */

  void SetMaximumNumberIterations(int i) {   m_IPARM[0] = i; }

  int  GetMaximumNumberIterations() { return m_IPARM[0]; }

  //void SetErrorReportingLevel(int i) {   m_IPARM[1] = i; }

  int  GetErrorReportingLevel() { return m_IPARM[1]; }

  void SetCommunicationSwitch(int i) {   m_IPARM[2] = i; }

  int  GetCommunicationSwitch() { return m_IPARM[2]; }

  //void SetOutputNumber(int i) {   m_IPARM[3] = i; }

  int  GetOutputNumber() { return m_IPARM[3]; }

  void SetSymmetricMatrixFlag(int i) {   m_IPARM[4] = i; }

  int  GetSymmetricMatrixFlag() { return m_IPARM[4]; }

  void SetAdaptiveSwitch(int i) {   m_IPARM[5] = i; }

  int  GetAdaptiveSwitch() { return m_IPARM[5]; }

  void SetAdaptiveCaseSwitch(int i) {   m_IPARM[6] = i; }

  int  GetAdaptiveCaseSwitch() { return m_IPARM[6]; }

  void SetWorkspaceUsed(int i) {   m_IPARM[7] = i; }

  int  GetWorkspaceUsed() { return m_IPARM[7]; }

  void SetRedBlackOrderingSwitch(int i) {   m_IPARM[8] = i; }

  int  GetRedBlackOrderingSwitch() { return m_IPARM[8]; }

  void SetRemoveSwitch(int i) {   m_IPARM[9] = i; }

  int  GetRemoveSwitch() { return m_IPARM[9]; }

  void SetTimingSwitch(int i) {   m_IPARM[10] = i; }

  int  GetTimingSwitch() { return m_IPARM[10]; }

  void SetErrorAnalysisSwitch(int i) {   m_IPARM[11] = i; }

  int  GetErrorAnalysisSwitch() { return m_IPARM[11]; }

  void   SetAccuracy(double i) {   m_RPARM[0] = i; }

  double GetAccuracy()    { return m_RPARM[0]; }

  void   SetLargestJacobiEigenvalueEstimate(double i) {   m_RPARM[1] = i; }

  double GetLargestJacobiEigenvalueEstimate()    { return m_RPARM[1]; }

  void   SetSmallestJacobiEigenvalueEstimate(double i) {   m_RPARM[2] = i; }

  double GetSmallestJacobiEigenvalueEstimate()    { return m_RPARM[2]; }

  void   SetDampingFactor(double i) {   m_RPARM[3] = i; }

  double GetDampingFactor()    { return m_RPARM[3]; }

  void   SetOverrelaxationParameter(double i) {   m_RPARM[4] = i; }

  double GetOverrelaxationParameter()    { return m_RPARM[4]; }

  void   SetEstimatedSpectralRadiusSSOR(double i) {   m_RPARM[5] = i; }

  double GetEstimatedSpectralRadiusSSOR()    { return m_RPARM[5]; }

  void   SetEstimatedSpectralRadiusLU(double i) {   m_RPARM[6] = i; }

  double GetEstimatedSpectralRadiusLU()    { return m_RPARM[6]; }

  void   SetTolerance(double i) {   m_RPARM[7] = i; }

  double GetTolerance()    { return m_RPARM[7]; }

  void   SetTimeToConvergence(double i) {   m_RPARM[8] = i; }

  double GetTimeToConvergence()    { return m_RPARM[8]; }

  void   SetTimeForCall(double i) {   m_RPARM[9] = i; }

  double GetTimeForCall()    { return m_RPARM[9]; }

  void   SetDigitsInError(double i) {   m_RPARM[10] = i; }

  double GetDigitsInError()    { return m_RPARM[10]; }

  void   SetDigitsInResidual(double i) {   m_RPARM[11] = i; }

  double GetDigitsInResidual()    { return m_RPARM[11]; }

  void JacobianConjugateGradient() { m_Method = 0; }

  void JacobianSemiIterative() { m_Method = 1; }

  void SuccessiveOverrelaxation() { m_Method = 2; }

  void SymmetricSuccessiveOverrelaxationConjugateGradient() { m_Method = 3; }

  void SymmetricSuccessiveOverrelaxationSuccessiveOverrelaxation() { m_Method = 4; }

  void ReducedSystemConjugateGradient() { m_Method = 5; }

  void ReducedSystemSemiIteration() { m_Method = 6; }


  virtual void SetMaximumNonZeroValuesInMatrix(unsigned int maxNonZeroValues) {m_MaximumNonZeroValues = maxNonZeroValues;}


  void ScaleMatrix(Float scale, unsigned int matrixIndex);


  LinearSystemWrapperItpack();

  ~LinearSystemWrapperItpack();


  /* memory management routines */
  virtual void  InitializeMatrix(unsigned int matrixIndex);  

  virtual bool  IsMatrixInitialized(unsigned int matrixIndex);

  virtual void  DestroyMatrix(unsigned int matrixIndex);

  virtual void  InitializeVector(unsigned int vectorIndex);

  virtual bool  IsVectorInitialized(unsigned int vectorIndex);

  virtual void  DestroyVector(unsigned int vectorIndex);

  virtual void  InitializeSolution(unsigned int solutionIndex);

  virtual bool  IsSolutionInitialized(unsigned int solutionIndex);

  virtual void  DestroySolution(unsigned int solutionIndex);

  /* assembly & solving routines */
  virtual Float GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex) const;

  virtual void  SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex);

  virtual void  AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex);

  virtual void GetColumnsOfNonZeroMatrixElementsInRow( unsigned int row, ColumnArray& cols, unsigned int matrixIndex );

  virtual Float GetVectorValue(unsigned int i, unsigned int vectorIndex) const;

  virtual void  SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex);

  virtual void  AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex);

  virtual Float GetSolutionValue(unsigned int i, unsigned int solutionIndex) const;

  virtual void  SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex);

  virtual void  AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex);

  virtual void  Solve(void);


  /* matrix & vector manipulation routines */
  virtual void  SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2);

  virtual void  SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2);

  virtual void  SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2);

  virtual void  CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex);

  virtual void  CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex);

  virtual void  MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex, unsigned int rightMatrixIndex);

  virtual void  MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);

private:

  MatrixHolder *m_Matrices;

  VectorHolder *m_Vectors;

  VectorHolder *m_Solutions;

  //UnsignedIntegerArrayPtr m_MaximumNonZeroValues;
  unsigned int m_MaximumNonZeroValues;

  ItkItpackSolverFunction    m_Methods[7];

  integer m_Method;

  integer m_IPARM[12];

  doublereal m_RPARM[12];

};

class FEMExceptionItpackSolver : public FEMException
{
public:

  typedef long      integer;

  FEMExceptionItpackSolver(const char *file, unsigned int lineNumber, std::string location, integer errorCode);

  virtual ~FEMExceptionItpackSolver() throw() {}

  itkTypeMacro(FEMExceptionItpackSolver,FEMException);

};
}} // end namespace itk::fem

#endif // #ifndef __itkFEMLinearSystemWrapperItpack_h