itkfemlinearsystemwrapperitpack.cxx

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

    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SetSolutionValue", "Indexed solution not yet allocated");
  }

  /* set value */
  (*m_Solutions)[solutionIndex][i] = value;

}


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

  (*m_Solutions)[solutionIndex][i] += value;
}


void LinearSystemWrapperItpack::Solve(void)
{

  /* error checking */
  if (!m_Order || !m_Matrices || !m_Vectors || !m_Solutions)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::Solve", "Not all necessary data members have been allocated");
  }
  if ( !(*m_Matrices)[0].GetOrder() ) {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::AddSolutionValue", "Primary matrix never filled");
  }

  /* itpack variables */
  integer N;
  integer NW;
  integer NCG;
  integer *IWKSP;
  doublereal *WKSP;
  integer IERR = 0;

  /* *******************************************************************
   * FIX ME: itpack does not allow for any non-zero diagonal elements
   * so "very small" numbers are inserted to allow for a solution
   *
  int i;
  doublereal fakeZero = 1.0e-16;

  //insert "fake" zeros 
  for (i=0; i<static_cast<int>(m_Order); i++)
  {
    if ( (*m_Matrices)[0].Get(i,i) == 0.0)
    {
      (*m_Matrices)[0].Set(i,i,fakeZero);
    }
  }
  // END FIXME 
   *********************************************************************/


  /* Initialize solution values (set to zero) */
  if (!this->IsSolutionInitialized(0)) 
  {
    this->InitializeSolution(0);
  }


 /* Set up itpack workspace variables
  *
  * N -> Order of system
  *
  * NCG -> temp var
  *
  * NW -> on input: length of wksp, on output: actual amount used
  *  jcg_()    - 4*N + NCG
  *  jsi_()    - 2*N
  *  sor_()    - N
  *  ssorcg_() - 6*N + NCG
  *  ssorsi_() - 5*N
  *  rscg_()   - N + 3*NB + NCG 
  *  rssi_()   - N + NB
  *
  * IWKSP -> temp variable used by itpack (integer[3*N])
  * WKSP -> temp variable used by itpack (doublereal[NW])
  */
  N = m_Order;
  if (m_IPARM[4] == 1)
  {
    NCG = 4 * m_IPARM[0];
  }
  else 
  {
    NCG = 2 * m_IPARM[0];
  }
  switch ( m_Method ) {
  case 0:
    NW = 4*N + NCG;
    break;
  case 1:
    NW = 2*N;
    break;
  case 2:
    NW = N;
    break;
  case 3:
    NW = 6*N + NCG;
    break;
  case 4:
    NW = 5*N;
    break;
  case 5:
    NW = N + m_IPARM[8] + NCG;
    break;
  case 6:
    NW = N + m_IPARM[8];
    break;
  }
  m_IPARM[7] = NW;
  IWKSP = new integer [ 3*N ];
  WKSP = new doublereal [ NW+2 ];

  integer i;
  for (i=0; i<NW; i++) 
  {
    WKSP[i] = 0.0;
  }
  for (i=0; i<(3*N); i++) 
  {
    IWKSP[i] = 0;
  }


  // Save maximum number of iteration, since it will
  // be overwritten.
  int max_num_iterations=m_IPARM[0];

  /* call to itpack solver routine */
  (*m_Methods[m_Method])( &N, (*m_Matrices)[0].GetIA(), (*m_Matrices)[0].GetJA(), (*m_Matrices)[0].GetA(), 
    (*m_Vectors)[0], (*m_Solutions)[0], &(IWKSP[0]), &NW, &(WKSP[0]), &(m_IPARM[0]), &(m_RPARM[0]), &IERR);

  m_IPARM[0]=max_num_iterations;

  /* remove exception throw on convergence failure */
  if (IERR < 100) 
  {
    if ( (IERR % 10) == 3 )
    {
      IERR = 0;
    }
  }

  /* clean up */
  delete [] IWKSP;
  delete [] WKSP;

  /* check for itpack error code */
  if (IERR > 0)
  {
    throw FEMExceptionItpackSolver(__FILE__, __LINE__, "LinearSystemWrapperItpack::Solve", IERR);
  }
  
}


void LinearSystemWrapperItpack::SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2)
{

  /* error checking */
  if (!m_Matrices)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapMatrices", "No matrices allocated");
  }
  if (matrixIndex1 >= m_NumberOfMatrices)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapMatrices", "m_Matrices", matrixIndex1);
  }
  if (matrixIndex2 >= m_NumberOfMatrices)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapMatrices", "m_Matrices", matrixIndex2);
  }

  int n = (*m_Matrices)[matrixIndex2].GetOrder();
  int nz = (*m_Matrices)[matrixIndex2].GetMaxNonZeroValues();
  integer* ia = (*m_Matrices)[matrixIndex2].GetIA();
  integer *ja = (*m_Matrices)[matrixIndex2].GetJA();
  doublereal *a = (*m_Matrices)[matrixIndex2].GetA();

  (*m_Matrices)[matrixIndex2].SetOrder( (*m_Matrices)[matrixIndex1].GetOrder() );
  (*m_Matrices)[matrixIndex2].SetMaxNonZeroValues ( (*m_Matrices)[matrixIndex1].GetMaxNonZeroValues() );
  (*m_Matrices)[matrixIndex2].SetCompressedRow((*m_Matrices)[matrixIndex1].GetIA(), 
    (*m_Matrices)[matrixIndex1].GetJA(), (*m_Matrices)[matrixIndex1].GetA() );

  (*m_Matrices)[matrixIndex1].SetOrder( n );
  (*m_Matrices)[matrixIndex1].SetMaxNonZeroValues ( nz );
  (*m_Matrices)[matrixIndex1].SetCompressedRow(ia, ja, a);


}


void LinearSystemWrapperItpack::SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2)
{
  /* error checking */
  if (!m_Vectors)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapVectors", "No vectors allocated");
  }
  if (vectorIndex1 >= m_NumberOfVectors)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapVectors", "m_Vectors", vectorIndex1);
  }
  if (vectorIndex2 >= m_NumberOfVectors)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapVectors", "m_Vectors", vectorIndex2);
  }

  VectorRepresentation temp = (*m_Vectors)[vectorIndex1];

  (*m_Vectors)[vectorIndex1] = (*m_Vectors)[vectorIndex2];
  (*m_Vectors)[vectorIndex2] = temp;
}


void LinearSystemWrapperItpack::SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2)
{
  /* error checking */
  if (!m_Solutions)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapSolutions", "No solutions allocated");
  }
  if (solutionIndex1 >= m_NumberOfSolutions)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapSolutions", "m_Solutions", solutionIndex1);
  }
  if (solutionIndex2 >= m_NumberOfSolutions)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::SwapSolutions", "m_Solutions", solutionIndex2);
  }

  VectorRepresentation temp = (*m_Solutions)[solutionIndex1];

  (*m_Solutions)[solutionIndex1] = (*m_Solutions)[solutionIndex2];
  (*m_Solutions)[solutionIndex2] = temp;
}


void LinearSystemWrapperItpack::CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex)
{

  /* error checking */
  if (!m_Vectors)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "No vectors allocated");
  }
  if (!m_Solutions)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "No solutions allocated");
  }
  if (vectorIndex >= m_NumberOfVectors)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "m_Vectors", vectorIndex);
  }
  if (solutionIndex >= m_NumberOfSolutions)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "m_Solutions", solutionIndex);
  }


  this->InitializeVector(vectorIndex);

  /* copy values */
  for (unsigned int i=0; i<m_Order; i++)
  {
    (*m_Vectors)[vectorIndex][i] = (*m_Solutions)[solutionIndex][i];
  }
}


void LinearSystemWrapperItpack::CopyVector2Solution(unsigned vectorIndex, unsigned int solutionIndex)
{

  /* error checking */
  if (!m_Vectors)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "No vectors allocated");
  }
  if (!m_Solutions)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "No solutions allocated");
  }
  if (vectorIndex >= m_NumberOfVectors)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "m_Vectors", vectorIndex);
  }
  if (solutionIndex >= m_NumberOfSolutions)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::CopySolution2Vector", "m_Solutions", solutionIndex);
  }


  this->InitializeSolution(solutionIndex);

  /* copy values */
  for (unsigned int i=0; i<m_Order; i++)
  {
    (*m_Solutions)[solutionIndex][i] = (*m_Vectors)[vectorIndex][i];
  }
}


void LinearSystemWrapperItpack::MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex, unsigned int rightMatrixIndex)
{
  /* error checking */
  if (!m_Matrices)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixMatrix", "No matrices allocated");
  }
  if (resultMatrixIndex >= m_NumberOfMatrices)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixMatrix", "m_Matrices", resultMatrixIndex);
  }
  if (leftMatrixIndex >= m_NumberOfMatrices)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixMatrix", "m_Matrices", leftMatrixIndex);
  }
  if (rightMatrixIndex >= m_NumberOfMatrices)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixMatrix", "m_Matrices", rightMatrixIndex);
  }

  (*m_Matrices)[leftMatrixIndex].mult( &((*m_Matrices)[rightMatrixIndex]), &((*m_Matrices)[resultMatrixIndex]) );

}


void LinearSystemWrapperItpack::MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex)
{

  /* error checking */
  if (!m_Matrices)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "No matrices allocated");
  }
  if (!m_Vectors)
  {
    throw FEMExceptionLinearSystem(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "No vectors allocated");
  }
  if (resultVectorIndex >= m_NumberOfVectors)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "m_Vectors", resultVectorIndex);
  }
  if (matrixIndex >= m_NumberOfMatrices)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "m_Matrices", matrixIndex);
  }
  if (vectorIndex >= m_NumberOfVectors)
  {
    throw FEMExceptionLinearSystemBounds(__FILE__, __LINE__, "LinearSystemWrapperItpack::MultiplyMatrixVector", "m_Vectors", vectorIndex);
  }

  /* perform mult */
  (*m_Matrices)[matrixIndex].mult( (*m_Vectors)[vectorIndex], (*m_Vectors)[resultVectorIndex] );

}


LinearSystemWrapperItpack::~LinearSystemWrapperItpack(void)
{
  delete m_Matrices;

  unsigned int i;
  if (m_Vectors != 0)
  {
    for (i=0; i<m_NumberOfVectors; i++)
    {
      if ( (*m_Vectors)[i] != 0 )
      {
        delete [] (*m_Vectors)[i];
      }
    }
    delete m_Vectors;
  }

  
  if (m_Solutions != 0)
  {
    for (i=0; i<m_NumberOfSolutions; i++)
    {
      if ( (*m_Solutions)[i] != 0 )
      {
        delete [] (*m_Solutions)[i];
      }
    }
    delete m_Solutions;
  }
  


}


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

    if (errorCode < 100) 
    {
      errorCode = errorCode % 10;
    }

    switch (errorCode) 
    {
    case 1 :
      solverError = "Invalid order of system";
      break;
    case 2 :
      solverError = "Workspace is not large enough";
      break;
    case 3 :
      solverError = "Failure to converge before reaching maximum number of iterations";
      break;
    case 4 :
      solverError = "Invalid order of black subsystem";
      break;
    case 101:
      solverError = "A diagonal element is not positive";
      break;
    case 102 :
      solverError = "No diagonal element in a row";
      break;
    case 201 :
      solverError = "Red-black indexing is not possible";
      break;
    case 301 :
      solverError = "No entry in a row of the original matrix";
      break;
    case 302 :
      solverError = "No entry in a row of the permuted matrix";
      break;
    case 303 :
      solverError = "Sorting error in a row of the permuted matrix";
      break;
    case 401 :
      solverError = "A diagonal element is not positive";
      break;
    case 402 :
      solverError = "No diagonal element in a row";
      break;
    case 501:
      solverError = "Failure to converge before reaching maximum number of iterations";
      break;
    case 502 :
      solverError = "Function does not change sign at endpoints";
      break;
    case 601 :
      solverError = "Successive iterations are not monotone increasing";
      break;
    default :
      solverError = "Unknown error code returned";
    }

  OStringStream buf;
  buf << "Error: " << solverError;

  SetDescription(buf.str().c_str());

  SetLocation(location);

}


}}  // end namespace itk::fem