CluX.cpp File Reference

#include <algorithm>
#include <cmath>
#include <limits>
#include "copasi.h"
#include "CluX.h"
#include "CMatrix.h"
#include "CVector.h"
#include "CProcessReport.h"
#include "lapack/lapackwrap.h"

Include dependency graph for CluX.cpp:

Go to the source code of this file.

Functions
bool	LUfactor (CMatrix< C_FLOAT64 > &A, CVector< size_t > &row, CVector< size_t > &col, CProcessReport *cb)

Function Documentation

bool LUfactor	(	CMatrix< C_FLOAT64 > &	A,
		CVector< size_t > &	row,
		CVector< size_t > &	col,
		CProcessReport *	cb
	)

Definition at line 30 of file CluX.cpp.

References CProcessReport::addItem(), C_FLOAT64, C_INVALID_INDEX, CProcessReport::finishItem(), min, CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), CProcessReport::progressItem(), CVector< CType >::resize(), and CVectorCore< CType >::size().

{
  size_t Rows = A.numRows();
  size_t Cols = A.numCols();
  size_t Dim = std::min(Rows, Cols);

  if (row.size() != Rows) row.resize(Rows);

  if (col.size() != Cols) col.resize(Cols);

  size_t i = 0, j = 0, k = 0;
  size_t rowP = 0, colP = Cols - 1;

  for (i = 0; i < Rows; i++)
    row[i] = i;

  for (i = 0; i < Cols; i++)
    col[i] = i;

  C_FLOAT64 Work = .5 * Dim * (Dim - 1);
  C_FLOAT64 Completion = 0;

  size_t hProcess;

  if (cb)
    hProcess = cb->addItem("LU decomposition...",
                           Completion,
                           &Work);

  C_FLOAT64 tmp;
  C_FLOAT64 * pTmp1;
  C_FLOAT64 * pTmp2;
  C_FLOAT64 * pTmp3;

  for (i = 0; i < Dim; i++)
    {
      Completion += Dim - i;

      if (cb && !cb->progressItem(hProcess)) return false;

      // find pivot in column j and  test for singularity.
      while (true)
        {
          rowP = C_INVALID_INDEX;
          tmp = 0.0;

          for (j = i, pTmp1 = &A(i, i); j < Rows; j++, pTmp1 += Cols)
            if (fabs(*pTmp1) > tmp) // (fabs(A(j, i)) > t)
              {
                rowP = j;
                tmp = fabs(*pTmp1); // fabs(A(j, i));
              }

          // rowP now has the index of maximum element
          // of column j, below the diagonal

          if (tmp < std::numeric_limits< C_FLOAT64 >::epsilon()) // now we have to swap colums to find a pivot
            {
              // Instead of blindly swapping with the last available
              // column we first check whether it is suitable. If not
              // we proceed with the lat but one ...
              while (i < colP)
                {
                  rowP = C_INVALID_INDEX;
                  tmp = 0.0;

                  for (j = i, pTmp1 = &A(i, colP); j < Rows; j++, pTmp1 += Cols)
                    if (fabs(*pTmp1) > tmp) // (fabs(A(j, colP)) > t)
                      {
                        rowP = j;
                        tmp = fabs(*pTmp1); // fabs(A(j, colP));
                      }

                  if (tmp >= std::numeric_limits< C_FLOAT64 >::epsilon()) break; // Found a suitable column

                  // and row

                  colP--; // proceed with previous column
                }

              if (i >= colP)
                {
                  if (cb) cb->finishItem(hProcess);

                  return true;
                }

              // Swap columns i and colP
              pTmp1 = &A(0, i);
              pTmp2 = &A(0, colP);

              for (k = 0; k < Rows; k++, pTmp1 += Cols, pTmp2 += Cols)
                {
                  tmp = *pTmp1;
                  *pTmp1 = *pTmp2;
                  *pTmp2 = tmp;
                }

              std::swap(col[i], col[colP]);
              colP--;
            }

          // Swap rows i and rowP
          pTmp1 = &A(i, 0);
          pTmp2 = &A(rowP, 0);

          for (k = 0; k < Cols; k++, pTmp1++, pTmp2++)
            {
              tmp = *pTmp1;
              *pTmp1 = *pTmp2;
              *pTmp2 = tmp;
            }

          std::swap(row[i], row[rowP]);
          break;
        }

      if (i < Rows - 1)    // compute elements i+1 <= k < M of i-th column
        {
          // note A(i,i) is guarranteed not to be zero
          tmp = 1.0 / A(i, i);

          pTmp1 = &A(i + 1, i);

          for (k = i + 1; k < Rows; k++, pTmp1 += Cols)
            * pTmp1 *= tmp;
        }

      if (i < Dim - 1)
        {
          // rank-1 update to trailing submatrix:   E = E - x*y;
          //
          // E is the region A(i+1:M, i+1:N)
          // x is the column vector A(i+1:M,i)
          // y is row vector A(i,i+1:N)

          size_t ii, jj;

          pTmp1 = &A(i + 1, i);
          pTmp3 = &A(i + 1, i + 1);

          for (ii = i + 1; ii < Rows; ii++, pTmp1 += Cols, pTmp3 += i + 1)
            {
              pTmp2 = &A(i, i + 1);

              for (jj = i + 1; jj < Cols; jj++, pTmp2++, pTmp3++)
                {
                  *pTmp3 -= *pTmp1 * *pTmp2;
                  //                  if (fabs(*pTmp3) < std::numeric_limits< C_FLOAT64 >::epsilon()) *pTmp3 = 0.0;
                }
            }
        }

#ifdef DEBUG_MATRIX
      DebugFile << rowP << "\tx\t" << colP + 1 << std::endl;
      DebugFile << "Row\t" << row << std::endl;
      DebugFile << "Col\t" << col << std::endl;
      DebugFile << A << std::endl;
#endif
    }

  if (cb) cb->finishItem(hProcess);

  return true;
}