itkfemlinearsystemwrapperitpack.h

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMLinearSystemWrapperItpack.h,v $
  Language:  C++
  Date:      $Date: 2008-05-23 20:25:28 $
  Version:   $Revision: 1.14 $

  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>

/** Array of pointers to available solver functions */
/** typedefs from f2c.h  */
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
 * \brief LinearSystemWrapper class that uses Itpack numeric library functions
 *        to define and solve a sparse linear system of equations
 * \sa LinearSystemWrapper
 */
class LinearSystemWrapperItpack : public LinearSystemWrapper
{
public:
 
  /** Standard "Self" typedef. */
  typedef LinearSystemWrapperItpack Self;

  /** Standard "Superclass" typedef. */
  typedef LinearSystemWrapper Superclass;

  /** matrix representatin typedef */
  typedef ItpackSparseMatrix MatrixRepresentation;

  /** vector of matrices typedef */
  typedef std::vector<MatrixRepresentation> MatrixHolder;

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

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

  /** vector of vector typedef */
  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 
   *
   * ----------------------------------------------------------------- 
   */

  /**
   * Set the maximum number of iterations
   * \param i maximum number of iterations that may be performed
   */
  void SetMaximumNumberIterations(int i) {   m_IPARM[0] = i; }

  /**
   * Get the maximum number iterations that may be performed
   */
  int  GetMaximumNumberIterations() { return m_IPARM[0]; }

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

  /**
   * Get a flag indicating the type of error reporting 
   */
  int  GetErrorReportingLevel() { return m_IPARM[1]; }

  /**
   * Set the communication switch - meaningless in this implementation
   * \param i flag value
   */
  void SetCommunicationSwitch(int i) {   m_IPARM[2] = i; }

  /**
   * Get the communication flag - meaningless in this implementation
   */
  int  GetCommunicationSwitch() { return m_IPARM[2]; }

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

  /**
   * Get the output number - meaningless in this implementation
   */
  int  GetOutputNumber() { return m_IPARM[3]; }

  /**
   * Set flag indicating symmetric matrix is being used
   * \param i 1=symmetric, 0=non-symmetric
   */
  void SetSymmetricMatrixFlag(int i) {   m_IPARM[4] = i; }

  /** 
   * Get flag indicating use of symmetric matrix (1=symmetric, 0=non-symmetric)
   */
  int  GetSymmetricMatrixFlag() { return m_IPARM[4]; }

  /**
   * Set flag for ???
   * \param i ??
   */
  void SetAdaptiveSwitch(int i) {   m_IPARM[5] = i; }

  /**
   * Get flag indicating ??
   */
  int  GetAdaptiveSwitch() { return m_IPARM[5]; }

  /**
   * Set flag for ??
   * \param i ??
   */
  void SetAdaptiveCaseSwitch(int i) {   m_IPARM[6] = i; }

  /**
   * Get flag indicating ??
   */
  int  GetAdaptiveCaseSwitch() { return m_IPARM[6]; }

  /**
   * Set size of workspace used by solver
   * \param i size of the workspace vector
   * \note this value is set by default
   */
  void SetWorkspaceUsed(int i) {   m_IPARM[7] = i; }

  /**
   * Get the size of the workspace used by solver
   * \note after solver is called this is the amount of workspace actually used
   */
  int  GetWorkspaceUsed() { return m_IPARM[7]; }

  /**
   * Set flag indicating use of red black ordering
   * \param i 1=red black ordering used, 0=not
   */
  void SetRedBlackOrderingSwitch(int i) {   m_IPARM[8] = i; }

  /**
   * Get the flag indicating use of red black ordering
   */
  int  GetRedBlackOrderingSwitch() { return m_IPARM[8]; }

  /**
   * Set flag indicating ??
   * \param i ??
   */
  void SetRemoveSwitch(int i) {   m_IPARM[9] = i; }

  /**
   * Get flag indicating ??
   */
  int  GetRemoveSwitch() { return m_IPARM[9]; }

  /**
   * Set the flag indicating use of timer routines - meaningless in this implementation
   * \param i flag
   */
  void SetTimingSwitch(int i) {   m_IPARM[10] = i; }

  /**
   * Get the flag indicating use of the timer routines - meaningless in this implementation
   */
  int  GetTimingSwitch() { return m_IPARM[10]; }

  /**
   * Set the flag for level of error reporting - meaningless in this implementation
   * \param i flag for level of error analysis
   */
  void SetErrorAnalysisSwitch(int i) {   m_IPARM[11] = i; }

  /**
   * Get the flag for level of error reporting - meaningless in this implementation
   */
  int  GetErrorAnalysisSwitch() { return m_IPARM[11]; }

  /**
   * Set the level of accuracy for an acceptable solution
   * \param i accuracy desired
   */
  void   SetAccuracy(double i) {   m_RPARM[0] = i; }

  /**
   * Get the level of accuracy
   */
  double GetAccuracy()    { return m_RPARM[0]; }

  /**
   * Set ??
   * \param i larges jacobian eigenvalue estimate
   */
  void   SetLargestJacobiEigenvalueEstimate(double i) {   m_RPARM[1] = i; }

  /**
   * Get ??
   */
  double GetLargestJacobiEigenvalueEstimate()    { return m_RPARM[1]; }

  /**
   * Set ??
   * \param i smalles jacobian eigenvalue estimate
   */
  void   SetSmallestJacobiEigenvalueEstimate(double i) {   m_RPARM[2] = i; }

  /**
   * Get ??
   */
  double GetSmallestJacobiEigenvalueEstimate()    { return m_RPARM[2]; }

  /**
   * Set the damping factor used by ??
   * \param i damping factor
   */
  void   SetDampingFactor(double i) {   m_RPARM[3] = i; }

  /**
   * Get the damping factor used by ??
   */
  double GetDampingFactor()    { return m_RPARM[3]; }

  /**
   * Set the over-relaxation parameter ??
   * \param i parameter
   */
  void   SetOverrelaxationParameter(double i) {   m_RPARM[4] = i; }

  /**
   * Get the over-relaxation parameter ??
   */
  double GetOverrelaxationParameter()    { return m_RPARM[4]; }

  /**
   * Set the ??
   * \param i ??