itkFEMLinearSystemWrapper.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#ifndef itkFEMLinearSystemWrapper_h
#define itkFEMLinearSystemWrapper_h

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

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

namespace itk
{
namespace fem
{
class ITKFEM_EXPORT 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 SetSystemOrder(unsigned int N)
  {
    m_Order = N;
  }

  unsigned int GetSystemOrder() const
  {
    return m_Order;
  }

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

  /*
   * Set the maximum number of entries permitted in a matrix
   * \param matrixIndex index of matrix to set value for
   * \param maxNonZeros maximum number of entries allowed in matrix
   * \note in general this function does nothing, however it may
   *       redefined by the derived wrapper if necessary
   */
  // virtual void SetMaximumNonZeroValuesInMatrix(unsigned int maxNonZeroValues)
  // = 0;

  unsigned int GetNumberOfMatrices() const
  {
    return m_NumberOfMatrices;
  }

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

  unsigned int GetNumberOfVectors() const
  {
    return m_NumberOfVectors;
  }

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

  unsigned int GetNumberOfSolutions() const
  {
    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 MultiplyMatrixSolution(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int solutionIndex);

  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 ITK_ABI_EXPORT FEMExceptionLinearSystem : public FEMException
{
public:
  FEMExceptionLinearSystem(const char *file, unsigned int lineNumber, std::string location, std::string moreDescription);

  virtual ~FEMExceptionLinearSystem()
  ITK_NOEXCEPT ITK_OVERRIDE;

  itkTypeMacro(FEMExceptionLinearSystem, FEMException);
};

class ITK_ABI_EXPORT 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()
  ITK_NOEXCEPT ITK_OVERRIDE;

  itkTypeMacro(FEMExceptionLinearSystem, FEMException);
};
}
}  // end namespace itk::fem

#endif // #ifndef itkFEMLinearSystemWrapper_h