itkFEMLinearSystemWrapper.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMLinearSystemWrapper.h,v $
  Language:  C++
  Date:      $Date: 2009-01-29 21:55:14 $
  Version:   $Revision: 1.23 $

  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 __itkFEMLinearSystemWrapper_h
#define __itkFEMLinearSystemWrapper_h 

#include "itkMacro.h"
#include "itkFEMSolution.h"
#include "itkFEMException.h"

#include <vector>
#include <typeinfo>
#include <string>

namespace itk {
namespace fem {


class LinearSystemWrapper : public Solution
{
public:
  typedef LinearSystemWrapper Self;
  typedef Solution            Superclass;
  typedef Self*               Pointer;
  typedef const Self*         ConstPointer;

  typedef std::vector<unsigned int> ColumnArray;

  LinearSystemWrapper() 
    : m_Order(0), m_NumberOfMatrices(1), m_NumberOfVectors(1), m_NumberOfSolutions(1) {}
      /* , m_PrimaryMatrixSetupFunction(0), m_PrimaryVectorSetupFunction(0), m_PrimarySolutionSetupFunction(0) {} */

  virtual ~LinearSystemWrapper() {};

  virtual void Clean( void );

  void SetSystemOrder(unsigned int N) { m_Order = N; }

  unsigned int GetSystemOrder() const { return m_Order; }

  void SetNumberOfMatrices(unsigned int nMatrices) { m_NumberOfMatrices = nMatrices; }

  //virtual void SetMaximumNonZeroValuesInMatrix(unsigned int maxNonZeroValues) = 0;

  unsigned int GetNumberOfMatrices() { return m_NumberOfMatrices; }

  void SetNumberOfVectors(unsigned int nVectors) { m_NumberOfVectors = nVectors; }

  unsigned int GetNumberOfVectors() { return m_NumberOfVectors; }

  void SetNumberOfSolutions(unsigned int nSolutions) { m_NumberOfSolutions = nSolutions; }

  unsigned int GetNumberOfSolutions() { return m_NumberOfSolutions; }

  virtual void InitializeMatrix(unsigned int matrixIndex = 0) = 0;


  virtual bool IsMatrixInitialized(unsigned int matrixIndex = 0) = 0;

  virtual void DestroyMatrix(unsigned int matrixIndex = 0) = 0;

  virtual void InitializeVector(unsigned int vectorIndex = 0) = 0;


  virtual bool IsVectorInitialized(unsigned int vectorIndex = 0) = 0;

  virtual void DestroyVector(unsigned int vectorIndex = 0) = 0;

  virtual void InitializeSolution(unsigned int solutionIndex = 0) = 0;

  virtual bool IsSolutionInitialized(unsigned int solutionIndex = 0) = 0;

  virtual void DestroySolution(unsigned int solutionIndex = 0) = 0;

  virtual Float GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex = 0) const = 0;

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

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

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

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

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

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

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

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

  virtual void Solve(void) = 0;

  virtual void SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2) = 0;

  virtual void CopyMatrix(unsigned int matrixIndex1, unsigned int matrixIndex2);

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

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


  virtual void ScaleMatrix(Float scale, unsigned int matrixIndex = 0);


  void ScaleVector(Float scale, unsigned int vectorIndex = 0);


  void ScaleSolution(Float scale, unsigned int solutionIndex = 0);

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

  virtual void AddMatrixMatrix(unsigned int matrixIndex1, unsigned int matrixIndex2);

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

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

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

  virtual void CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex) = 0;
  virtual void CopyVector(unsigned int vectorSource, unsigned int vectorDestination);

  virtual void OptimizeMatrixStorage(unsigned int matrixIndex, unsigned int tempMatrixIndex);

  virtual void ReverseCuthillMckeeOrdering(ColumnArray& newNumbering, unsigned int matrixIndex = 0);

protected:

  unsigned int m_Order;

  unsigned int m_NumberOfMatrices;

  unsigned int m_NumberOfVectors;

  unsigned int m_NumberOfSolutions;

  /*
   * Function used to prepare primary matrix for numerical solving 
   */
  //void (*m_PrimaryMatrixSetupFunction)(LinearSystemWrapper *lsw);

  /*
   * Function used to prepare primary vector for numerical solving 
   */
  /* void (*m_PrimaryVectorSetupFunction)(LinearSystemWrapper *lsw);*/

  /*
   * Function used to prepare primary matrix for numerical solving 
   */
  /* void (*m_PrimarySolutionSetupFunction)(LinearSystemWrapper *lsw); */

private:

  void CuthillMckeeOrdering(ColumnArray& newNumbering, int startingRow, unsigned int matrixIndex = 0);

  void FollowConnectionsCuthillMckeeOrdering(unsigned int rowNumber, ColumnArray& rowDegree, ColumnArray& newNumbering, unsigned int nextRowNumber, unsigned int matrixIndex = 0);

  LinearSystemWrapper(const LinearSystemWrapper&);

  const LinearSystemWrapper& operator= (const LinearSystemWrapper&);

};

class FEMExceptionLinearSystem : public FEMException
{
public:
  FEMExceptionLinearSystem(const char *file, unsigned int lineNumber, std::string location, std::string moreDescription);
 
  virtual ~FEMExceptionLinearSystem() throw() {}
  
  itkTypeMacro(FEMExceptionLinearSystem,FEMException);
  
};

class FEMExceptionLinearSystemBounds : public FEMException
{
public:
  FEMExceptionLinearSystemBounds(const char *file, unsigned int lineNumber, std::string location, std::string moreDescription, unsigned int index1);
 
  FEMExceptionLinearSystemBounds(const char *file, unsigned int lineNumber, std::string location, std::string moreDescription, unsigned int index1, unsigned int index2);

  virtual ~FEMExceptionLinearSystemBounds() throw() {}
  
  itkTypeMacro(FEMExceptionLinearSystem,FEMException);
  
};

}} // end namespace itk::fem

#endif // #ifndef __itkFEMLinearSystemWrapper_h