itkfemlinearsystemwrapperitpack.cxx

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

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMLinearSystemWrapperItpack.cxx,v $
  Language:  C++
  Date:      $Date: 2006-11-07 23:23:16 $
  Version:   $Revision: 1.21 $

  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.

=========================================================================*/
#include "itkNumericTraits.h"
#include "itkFEMLinearSystemWrapperItpack.h"
#include "itkFEMException.h"
#include "itpack.h"


namespace itk {
namespace fem {



/**
 * constructor 
 */
LinearSystemWrapperItpack::LinearSystemWrapperItpack()
{

  /* fill m_Methods with pointers to solver functions */
  m_Methods[0] = jcg_;
  m_Methods[1] = jsi_;
  m_Methods[2] = sor_;
  m_Methods[3] = ssorcg_;
  m_Methods[4] = ssorsi_;
  m_Methods[5] = rscg_;
  m_Methods[6] = rssi_;
  m_Method = 0;                   /* set default method to jcg_ */

  /* Set solving parameters */
  dfault_( &(m_IPARM[0]) , &(m_RPARM[0]) );

  // We don't want the solver routines to
  // overwrite the parameters.
  m_IPARM[2]=1;


  /* m_IPARM[0] = 500; */  /* number of iterations */
  m_IPARM[1] = -1;   /* no error message output */
  m_IPARM[4] = 1;    /* non-symmetric matrix */

  /* itpack recommended (but not default) value */
#undef min
  m_RPARM[7] = 500.0 * NumericTraits<double>::min();

  m_MaximumNonZeroValues = 0;
  m_Matrices = 0;
  m_Vectors = 0;
  m_Solutions = 0;

}


void LinearSystemWrapperItpack::InitializeMatrix(unsigned int matrixIndex)
{
  /* error checking */
  if (!m_Order) 
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeMatrix", "System order not set");
  }
  if (matrixIndex >= m_NumberOfMatrices)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeMatrix", "m_Matrices", matrixIndex);
  }
  if (!m_MaximumNonZeroValues)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeMatrix", "Maximum number of non zeros not set");
  }

  // allocate if necessay
  if (m_Matrices == 0)
  {
    m_Matrices = new MatrixHolder(m_NumberOfMatrices);
  }

  /* Set required variables */
  (*m_Matrices)[matrixIndex].Clear();
  (*m_Matrices)[matrixIndex].SetOrder(m_Order);
  (*m_Matrices)[matrixIndex].SetMaxNonZeroValues( m_MaximumNonZeroValues );

  return;

}


bool LinearSystemWrapperItpack::IsMatrixInitialized(unsigned int matrixIndex)
{
  if (!m_Matrices) return false;
  if ( !(*m_Matrices)[matrixIndex].GetOrder() ) return false;
  if ( !(*m_Matrices)[matrixIndex].GetMaxNonZeroValues() ) return false;


  return true;
}

void LinearSystemWrapperItpack::InitializeVector(unsigned int vectorIndex)
{

  /* error checking */
  if (!m_Order)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeVector", "System order not set");
  }
  if (vectorIndex >= m_NumberOfVectors)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeVector", "m_Vectors", vectorIndex);
  }

  /* allocate if necessay */
  if (m_Vectors == 0)
  {
    m_Vectors = new VectorHolder(m_NumberOfVectors);
  }
  
  /* delete old vector */
  if ( (*m_Vectors)[vectorIndex] != 0 )
  {
    delete [] (*m_Vectors)[vectorIndex];
  }

  /* insert new vector */
  (*m_Vectors)[vectorIndex] = new doublereal [m_Order];

  /* fill with zeros */
  for (unsigned int i=0; i<m_Order; i++)
  {
    (*m_Vectors)[vectorIndex][i] = 0.0;
  }

  return;

}


bool LinearSystemWrapperItpack::IsVectorInitialized(unsigned int vectorIndex)
{
  if (!m_Vectors) return false;
  if ( !(*m_Vectors)[vectorIndex] ) return false;

  return true;
}


void LinearSystemWrapperItpack::InitializeSolution(unsigned int solutionIndex)
{

  /* FIX ME: exceptions */
  if (!m_Order)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeSolution", "System order not set");
  }
  if (solutionIndex >= m_NumberOfSolutions)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::InitializeSolution", "m_Solutions", solutionIndex);
  }
  

  // allocate if necessay
  if (m_Solutions == 0)
  {
    m_Solutions = new VectorHolder(m_NumberOfSolutions);
  }

  /* delete old vector */
  if ( (*m_Solutions)[solutionIndex] != 0 )
  {
    delete [] (*m_Solutions)[solutionIndex];
  }

  /* insert new vector */
  (*m_Solutions)[solutionIndex] = new doublereal [m_Order];

  /* fill with zeros */
  for (unsigned int i=0; i<m_Order; i++)
  {
    (*m_Solutions)[solutionIndex][i] = 0.0;
  }

  return;

}


bool LinearSystemWrapperItpack::IsSolutionInitialized(unsigned int solutionIndex)
{
  if (!m_Solutions) return false;
  if ( !(*m_Solutions)[solutionIndex] ) return false;

  return true;
}

void LinearSystemWrapperItpack::DestroyMatrix(unsigned int matrixIndex)
{

  /* FIX ME: exceptions */
  if ( !m_Matrices ) return;
  if ( matrixIndex >= m_NumberOfMatrices )
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::DestroyMatrix", "m_Matrices", matrixIndex);
  }

  (*m_Matrices)[matrixIndex].Clear();
}


void LinearSystemWrapperItpack::DestroyVector(unsigned int vectorIndex)
{
  /* FIXME: exceptions */
  if (!m_Vectors) return;
  if (vectorIndex >= m_NumberOfVectors)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::DestroyVector", "m_Vectors", vectorIndex);
  }

  if ( !(*m_Vectors)[vectorIndex] ) return;

  /* delete vector */
  delete [] (*m_Vectors)[vectorIndex];
  (*m_Vectors)[vectorIndex] = 0;

}


void LinearSystemWrapperItpack::DestroySolution(unsigned int solutionIndex)
{
  // FIXME: exceptions
  if (!m_Solutions) return;
  if (solutionIndex >= m_NumberOfSolutions)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::DestroySolution", "m_Solutions", solutionIndex);
  }
  if ( !(*m_Solutions)[solutionIndex] ) return;

  /* delete vector */
  delete [] (*m_Solutions)[solutionIndex];
  (*m_Solutions)[solutionIndex] = 0;

}


LinearSystemWrapperItpack::Float LinearSystemWrapperItpack::GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex) const
{
  /* error checking */
  if (!m_Matrices)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetMatrixValue", "No matrices have been allocated");
  }
  if (matrixIndex >= m_NumberOfMatrices)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetMatrixValue", "m_Matrices", matrixIndex);
  }
  if ( (i >= m_Order) || (j >= m_Order) )
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetMatrixValue", "m_Matrices[]", i,j);
  }

  /* return value */
  return (*m_Matrices)[matrixIndex].Get(i,j);
}


void LinearSystemWrapperItpack::SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex)
{
  /* error checking */
  if (!m_Matrices) 
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetMatrixValue", "No matrices have been allocated");
  }
  if ( (i >= m_Order) || (j >= m_Order) ) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetMatrixValue", "m_Matrices[]", i, j);
  }
  if (matrixIndex >= m_NumberOfMatrices) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetMatrixValue", "m_Matrices", matrixIndex);
  }

  /* set value */
  ((*m_Matrices)[matrixIndex]).Set(i,j,value);
}


void LinearSystemWrapperItpack::AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex)
{
  // FIXME: error checking
  if (!m_Matrices) 
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddMatrixValue", "No matrices have been allocated");
  }
  if ( (i >= m_Order) || (j >= m_Order) ) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddMatrixValue", "m_Matrices[]", i, j);
  }
  if (matrixIndex >= m_NumberOfMatrices) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddMatrixValue", "m_Matrices", matrixIndex);
  }

  ((*m_Matrices)[matrixIndex]).Add(i,j,value);
}


void LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow( unsigned int row, ColumnArray& cols, unsigned int matrixIndex )
{

  /* FIXME: error checking */
  if (!m_Matrices) 
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow", "No matrices have been allocated");
  }
  if (row >= this->m_Order) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow", "m_Matrices[]", row);
  }
  if (matrixIndex >= m_NumberOfMatrices) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetColumnsOfNonZeroMatrixElementsInRow", "m_Matrices", matrixIndex);
  }

  cols.clear();
  
  ItpackSparseMatrix* mat=&(*m_Matrices)[matrixIndex];

  /* Check if matrix is in readable form */
  if (mat->m_MatrixFinalized )
  {
    /* get search bounds in appropriate row */
    unsigned int lower = mat->m_IA[row]-1;
    unsigned int upper = mat->m_IA[row+1]-1;
    
    for(unsigned int j=lower; j<upper; j++)
    {
      cols.push_back(mat->m_JA[j]-1);
    }
  }
  else /* Scan the linked list to obtain the correct indices. */
  {
    int wrk=mat->m_IA[row]-1;
    while(wrk>0)
    {
      cols.push_back(mat->m_JA[wrk]-1);
      wrk=mat->m_IWORK[wrk]-1;
    }
  }

}


void LinearSystemWrapperItpack::ScaleMatrix(Float scale, unsigned int matrixIndex)
{
  /* error checking */
  if (!m_Matrices) 
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::ScaleMatrix", "No matrices have been allocated");
  }
  if (matrixIndex >= m_NumberOfMatrices) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::ScaleMatrix", "m_Matrices", matrixIndex);
  }

  int i;
  doublereal *values = (*m_Matrices)[matrixIndex].GetA();
  for (i=0; i<(*m_Matrices)[matrixIndex].GetIA()[this->m_Order] - 1; i++) 
  {
    values[i] = values[i]*scale;
  }
  
}


LinearSystemWrapperItpack::Float LinearSystemWrapperItpack::GetVectorValue(unsigned int i, unsigned int vectorIndex) const
{
  /* error checking */
  if (!m_Vectors) 
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetVectorValue", "No vectors have been allocated");
  }
  if (i >= m_Order) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetVectorValue", "m_Vectors[]", i);
  }
  if (vectorIndex >= m_NumberOfVectors) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetVectorValue", "m_Vectors", vectorIndex);
  }
  if ( !(*m_Vectors)[vectorIndex] )
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::GetVectorValue", "Indexed vector not yet allocated");
  }

  /* return value */
  return (*m_Vectors)[vectorIndex][i];
}


void LinearSystemWrapperItpack::SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex)
{
  /* error checking */
  if (!m_Vectors) 
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetVectorValue", "No vectors have been allocated");
  }
  if (i >= m_Order) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetVectorValue", "m_Vectors[]", i);
  }
  if (vectorIndex >= m_NumberOfVectors) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetVectorValue", "m_Vectors", vectorIndex);
  }
  if ( !(*m_Vectors)[vectorIndex] )
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetVectorValue", "Indexed vector not yet allocated");
  }

  /* set value */
  (*m_Vectors)[vectorIndex][i] = value;
}


void LinearSystemWrapperItpack::AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex)
{
  /*( error checking */
  if (!m_Vectors) 
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddVectorValue", "No vectors have been allocated");
  }
  if (i >= m_Order) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddVectorValue", "m_Vectors[]", i);
  }
  if (vectorIndex >= m_NumberOfVectors) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddVectorValue", "m_Vectors", vectorIndex);
  }
  if ( !(*m_Vectors)[vectorIndex] )
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddVectorValue", "Indexed vector not yet allocated");
  }

  /* add value */
  (*m_Vectors)[vectorIndex][i] += value;
}


LinearSystemWrapperItpack::Float LinearSystemWrapperItpack::GetSolutionValue(unsigned int i, unsigned int solutionIndex) const
{
  // FIXME: error checking
  if ( (i >= m_Order) || !m_Solutions || (solutionIndex >= m_NumberOfSolutions) || !(*m_Solutions)[solutionIndex] )
  {
    return 0.0;
  }
  return (*m_Solutions)[solutionIndex][i];
}


void LinearSystemWrapperItpack::SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex)
{
  /* error checking */
  if (!m_Solutions) 
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetSolutionValue", "No solutions have been allocated");
  }
  if (i >= m_Order) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetSolutionValue", "m_Solutions[]", i);
  }
  if (solutionIndex >= m_NumberOfSolutions) 
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetSolutionValue", "m_Solutions", solutionIndex);
  }
  if ( !(*m_Solutions)[solutionIndex] )
  {