CLinkMatrix Class Reference

#include <CLinkMatrix.h>

Inheritance diagram for CLinkMatrix:

Collaboration diagram for CLinkMatrix:

Public Member Functions
template<class CType >
bool	applyRowPivot (CVectorCore< CType > &vector) const
bool	build (const CMatrix< C_FLOAT64 > &matrix, size_t maxRank=C_INVALID_INDEX)
void	clearPivoting ()
	CLinkMatrix ()
	CLinkMatrix (const CLinkMatrix &src)
bool	doColumnPivot (CMatrix< C_FLOAT64 > &matrix) const
bool	doRowPivot (CMatrix< C_FLOAT64 > &matrix) const
size_t	getNumDependent () const
const size_t &	getNumIndependent () const
const CVector< size_t > &	getRowPivots () const
bool	leftMultiply (const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
bool	rightMultiply (const C_FLOAT64 &alpha, const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
bool	undoColumnPivot (CMatrix< C_FLOAT64 > &matrix) const
bool	undoRowPivot (CMatrix< C_FLOAT64 > &matrix) const
virtual	~CLinkMatrix ()
Public Member Functions inherited from CMatrix< C_FLOAT64 >
bool	applyPivot (const CVector< size_t > &pivot)
virtual C_FLOAT64 *	array ()
virtual const C_FLOAT64 *	array () const
	CMatrix (size_t rows=0, size_t cols=0)
	CMatrix (const CMatrix< C_FLOAT64 > &src)
virtual size_t	numCols () const
virtual size_t	numRows () const
virtual elementType &	operator() (const size_t &row, const size_t &col)
virtual const elementType &	operator() (const size_t &row, const size_t &col) const
virtual CMatrix< C_FLOAT64 > &	operator= (const CMatrix< C_FLOAT64 > &rhs)
virtual CMatrix< C_FLOAT64 > &	operator= (const C_FLOAT64 &value)
virtual C_FLOAT64 *	operator[] (size_t row)
virtual const C_FLOAT64 *	operator[] (size_t row) const
virtual void	resize (size_t rows, size_t cols, const bool &copy=false)
virtual size_t	size () const
virtual	~CMatrix ()

Private Member Functions
bool	applyColumnPivot (CMatrix< C_FLOAT64 > &matrix, const C_INT &incr) const
bool	applyRowPivot (CMatrix< C_FLOAT64 > &matrix, const CVector< size_t > &pivots) const
void	completePivotInformation ()

Private Attributes
size_t	mIndependent
CVector< size_t >	mPivotInverse
CVector< size_t >	mRowPivots
CVector< C_INT >	mSwapVector

Additional Inherited Members
Public Types inherited from CMatrix< C_FLOAT64 >
typedef C_FLOAT64	elementType
Protected Attributes inherited from CMatrix< C_FLOAT64 >
C_FLOAT64 *	mArray
size_t	mCols
size_t	mRows

Detailed Description

Definition at line 15 of file CLinkMatrix.h.

Constructor & Destructor Documentation

CLinkMatrix::CLinkMatrix

(

)

Default constructor

Definition at line 20 of file CLinkMatrix.cpp.

                        :
  CMatrix< C_FLOAT64 >(),
  mRowPivots(),
  mPivotInverse(),
  mSwapVector(),
  mIndependent(0)
{}

CLinkMatrix::CLinkMatrix

(

const CLinkMatrix &

src

)

Copy constructor

Definition at line 28 of file CLinkMatrix.cpp.

                                               :
  CMatrix< C_FLOAT64 >(src),
  mRowPivots(src.mRowPivots),
  mPivotInverse(src.mPivotInverse),
  mSwapVector(src.mSwapVector),
  mIndependent(src.mIndependent)
{}

CLinkMatrix::~CLinkMatrix ( )

virtual

Destructor

Definition at line 36 of file CLinkMatrix.cpp.

{}

Member Function Documentation

bool CLinkMatrix::applyColumnPivot ( CMatrix< C_FLOAT64 > &  matrix, const C_INT &  incr  ) const

private

Internal method performing apply and undo of column pivoting.

Parameters

CMatrix<	C_FLOAT64 > & matrix
const	C_INT & incr

Returns

bool success

Definition at line 562 of file CLinkMatrix.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), C_INT, dlaswp_(), mRowPivots, mSwapVector, CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), and CVectorCore< CType >::size().

Referenced by doColumnPivot(), and undoColumnPivot().

{
  if (matrix.numCols() < mRowPivots.size())
    {
      return false;
    }

  C_INT N = matrix.numRows();
  C_INT LDA = matrix.numCols();
  C_INT K1 = 1;
  C_INT K2 = mRowPivots.size();

  dlaswp_(&N, matrix.array(), &LDA, &K1, &K2,
          const_cast< C_INT * >(mSwapVector.array()),
          const_cast< C_INT * >(&incr));

  return true;
}

template<class CType >

bool CLinkMatrix::applyRowPivot ( CVectorCore< CType > &  vector ) const

inline

Apply the row pivot

Parameters

CVectorCore<

class CType > & vector

Returns

bool success

Definition at line 115 of file CLinkMatrix.h.

References mRowPivots, and CVectorCore< CType >::size().

Referenced by CModel::buildRedStoi(), doRowPivot(), and undoRowPivot().

  {
    if (vector.size() < mRowPivots.size())
      {
        return false;
      }

    CVector< bool > Applied(mRowPivots.size());
    Applied = false;

    CType Tmp;

    size_t i, imax = mRowPivots.size();
    size_t to;
    size_t from;

    for (i = 0; i < imax; i++)
      if (!Applied[i])
        {
          to = i;
          from = mRowPivots[to];

          if (from != i)
            {
              Tmp = vector[to];

              while (from != i)
                {
                  vector[to] = vector[from];
                  Applied[to] = true;

                  to = from;
                  from = mRowPivots[to];
                }

              vector[to] = Tmp;
            }

          Applied[to] = true;
        }

    return true;
  }

bool CLinkMatrix::applyRowPivot ( CMatrix< C_FLOAT64 > &  matrix, const CVector< size_t > &  pivots  ) const

private

Internal method performing apply and undo of row pivoting.

Parameters

CMatrix<	C_FLOAT64 > & matrix
const	CVector< size_t > & pivots

Returns

bool success

Definition at line 506 of file CLinkMatrix.cpp.

References CMatrix< CType >::array(), C_FLOAT64, CMatrix< C_FLOAT64 >::numCols(), CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), and CVectorCore< CType >::size().

{
  if (matrix.numRows() < pivots.size())
    {
      return false;
    }

  CVector< bool > Applied(pivots.size());
  Applied = false;

  CVector< C_FLOAT64 > Tmp(matrix.numCols());

  size_t i, imax = pivots.size();
  size_t to;
  size_t from;
  size_t numCols = matrix.numCols();

  for (i = 0; i < imax; i++)
    if (!Applied[i])
      {
        to = i;
        from = pivots[to];

        if (from != i)
          {
            memcpy(Tmp.array(), matrix[to], sizeof(C_FLOAT64) * numCols);

            while (from != i)
              {
                memcpy(matrix[to], matrix[from], sizeof(C_FLOAT64) * numCols);
                Applied[to] = true;

                to = from;
                from = pivots[to];
              }

            memcpy(matrix[to], Tmp.array(), sizeof(C_FLOAT64) * numCols);
          }

        Applied[to] = true;
      }

  return true;
}

bool CLinkMatrix::build	(	const CMatrix< C_FLOAT64 > &	matrix,
		size_t	maxRank = C_INVALID_INDEX
	)

Build the link matrix for the given matrix

Parameters

const	CMatrix< C_FLOAT64 > & matrix
size_t	maxRank (default: C_INVALID_INDEX)

Returns

bool success

Determine the rank of the matrix This code is copied from dgelsy.f

Definition at line 44 of file CLinkMatrix.cpp.

References CVectorCore< CType >::array(), CMatrix< C_FLOAT64 >::array(), CMatrix< CType >::array(), C_FLOAT64, C_INT, C_INT32, completePivotInformation(), dgeqp3_(), dlaic1_(), dtrtri_(), fatalError, max, min, mIndependent, mRowPivots, CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), CMatrix< C_FLOAT64 >::resize(), CVector< CType >::resize(), and TAU.

Referenced by CModel::buildLinkZero(), and CMCAMethod::createLinkMatrix().

{
  bool success = true;

  CMatrix< C_FLOAT64 > M(matrix);

  C_INT NumCols = (C_INT) M.numCols();
  C_INT NumRows = (C_INT) M.numRows();
  C_INT LDA = std::max<C_INT>(1, NumCols);

  CVector< C_INT > JPVT(NumRows);
  JPVT = 0;

  C_INT32 Dim = std::min(NumCols, NumRows);

  if (Dim == 0)
    {
      resize(NumRows, 0);

      mIndependent = 0;

      C_INT32 i;
      mRowPivots.resize(NumRows);

      for (i = 0; i < NumRows; i++)
        mRowPivots[i] = i;
    }
  else
    {
      CVector< C_FLOAT64 > TAU(Dim);

      CVector< C_FLOAT64 > WORK(1);
      C_INT LWORK = -1;
      C_INT INFO;

      // QR factorization of the stoichiometry matrix
      /*
       *  -- LAPACK routine (version 3.0) --
       *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
       *     Courant Institute, Argonne National Lab, and Rice University
       *     June 30, 1999
       *
       *  Purpose
       *  =======
       *
       *  DGEQP3 computes a QR factorization with column pivoting of a
       *  matrix A:  A*P = Q*R  using Level 3 BLAS.
       *
       *  Arguments
       *  =========
       *
       *  M       (input) INTEGER
       *          The number of rows of the matrix A. M >= 0.
       *
       *  N       (input) INTEGER
       *          The number of columns of the matrix A.  N >= 0.
       *
       *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
       *          On entry, the M-by-N matrix A.
       *          On exit, the upper triangle of the array contains the
       *          min(M,N)-by-N upper trapezoidal matrix R; the elements below
       *          the diagonal, together with the array TAU, represent the
       *          orthogonal matrix Q as a product of min(M,N) elementary
       *          reflectors.
       *
       *  LDA     (input) INTEGER
       *          The leading dimension of the array A. LDA >= max(1,M).
       *
       *  JPVT    (input/output) INTEGER array, dimension (N)
       *          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
       *          to the front of A*P (a leading column); if JPVT(J)=0,
       *          the J-th column of A is a free column.
       *          On exit, if JPVT(J)=K, then the J-th column of A*P was the
       *          the K-th column of A.
       *
       *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
       *          The scalar factors of the elementary reflectors.
       *
       *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
       *          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
       *
       *  LWORK   (input) INTEGER
       *          The dimension of the array WORK. LWORK >= 3*N+1.
       *          For optimal performance LWORK >= 2*N+(N+1)*NB, where NB
       *          is the optimal blocksize.
       *
       *          If LWORK = -1, then a workspace query is assumed; the routine
       *          only calculates the optimal size of the WORK array, returns
       *          this value as the first entry of the WORK array, and no error
       *          message related to LWORK is issued by XERBLA.
       *
       *  INFO    (output) INTEGER
       *          = 0: successful exit.
       *          < 0: if INFO = -i, the i-th argument had an illegal value.
       *
       *  Further Details
       *  ===============
       *
       *  The matrix Q is represented as a product of elementary reflectors
       *
       *     Q = H(1) H(2) . . . H(k), where k = min(m,n).
       *
       *  Each H(i) has the form
       *
       *     H(i) = I - tau * v * v'
       *
       *  where tau is a real/complex scalar, and v is a real/complex vector
       *  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
       *  A(i+1:m,i), and tau in TAU(i).
       *
       *  Based on contributions by
       *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
       *    X. Sun, Computer Science Dept., Duke University, USA
       *
       */

      dgeqp3_(&NumCols, &NumRows, M.array(), &LDA,
              JPVT.array(), TAU.array(), WORK.array(), &LWORK, &INFO);

      if (INFO < 0) fatalError();

      LWORK = (C_INT) WORK[0];
      WORK.resize(LWORK);

      dgeqp3_(&NumCols, &NumRows, M.array(), &LDA,
              JPVT.array(), TAU.array(), WORK.array(), &LWORK, &INFO);

      if (INFO < 0) fatalError();

      C_INT32 i;
      mRowPivots.resize(NumRows);

      for (i = 0; i < NumRows; i++)
        mRowPivots[i] = JPVT[i] - 1;

      /**
       * Determine the rank of the matrix
       * This code is copied from dgelsy.f
       */

      C_INT independent = 0;

      if (fabs(M(0, 0)) != 0.0)
        {
          ++independent;

          C_INT imax = 1;
          C_INT imin = 2;
          C_INT mn = Dim;

          C_FLOAT64 * pIsmin = WORK.array();
          C_FLOAT64 * pIsmax = pIsmin + Dim;

          *pIsmin = 1.0;
          *pIsmax = 1.0;

          C_FLOAT64 smax = fabs(M(0, 0));
          C_FLOAT64 smin = smax;
          C_FLOAT64 RCOND = std::max(NumRows, NumCols) * smax * std::numeric_limits< C_FLOAT64 >::epsilon();

          C_FLOAT64 sminpr, s1, c1, smaxpr, s2, c2;

          while (independent < mn && independent < maxRank)
            {
              dlaic1_(&imin, &independent, pIsmin, &smin, &M(independent, 0), &M(independent, independent), &sminpr, &s1, &c1);
              dlaic1_(&imax, &independent, pIsmax, &smax, &M(independent, 0), &M(independent, independent), &smaxpr, &s2, &c2);

              if (smaxpr * RCOND > sminpr)
                {
                  break;
                }

              for (size_t i = 0; i < independent - 1; ++i)
                {
                  *(pIsmin + i) *= s1;
                  *(pIsmax + i) *= s2;
                }

              *(pIsmin + independent) = c1;
              *(pIsmax + independent) = c2;

              smin = sminpr;
              smax = smaxpr;

              ++independent;
            }
        }

      mIndependent = independent;

      // Resize mL
      resize(NumRows - independent, independent);

      if (NumRows != independent && independent != 0)
        {
          /* to take care of differences between fortran's and c's memory  access,
               we need to take the transpose, i.e.,the upper triangular */
          char cL = 'U';
          char cU = 'N'; /* values in the diagonal of R */

          // Calculate Row Echelon form of R.
          // First invert R_1,1
          /* int dtrtri_(char *uplo,
           *             char *diag,
           *             integer *n,
           *             doublereal * A,
           *             integer *lda,
           *             integer *info);
           *  -- LAPACK routine (version 3.0) --
           *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
           *     Courant Institute, Argonne National Lab, and Rice University
           *     March 31, 1993
           *
           *  Purpose
           *  =======
           *
           *  DTRTRI computes the inverse of a real upper or lower triangular
           *  matrix A.
           *
           *  This is the Level 3 BLAS version of the algorithm.
           *
           *  Arguments
           *  =========
           *
           *  uplo    (input) CHARACTER*1
           *          = 'U':  A is upper triangular;
           *          = 'L':  A is lower triangular.
           *
           *  diag    (input) CHARACTER*1
           *          = 'N':  A is non-unit triangular;
           *          = 'U':  A is unit triangular.
           *
           *  n       (input) INTEGER
           *          The order of the matrix A.  n >= 0.
           *
           *  A       (input/output) DOUBLE PRECISION array, dimension (lda,n)
           *          On entry, the triangular matrix A.  If uplo = 'U', the
           *          leading n-by-n upper triangular part of the array A contains
           *          the upper triangular matrix, and the strictly lower
           *          triangular part of A is not referenced.  If uplo = 'L', the
           *          leading n-by-n lower triangular part of the array A contains
           *          the lower triangular matrix, and the strictly upper
           *          triangular part of A is not referenced.  If diag = 'U', the
           *          diagonal elements of A are also not referenced and are
           *          assumed to be 1.
           *          On exit, the (triangular) inverse of the original matrix, in
           *          the same storage format.
           *
           *  lda     (input) INTEGER
           *          The leading dimension of the array A.  lda >= max(1,n).
           *
           *  info    (output) INTEGER
           *          = 0: successful exit
           *          < 0: if info = -i, the i-th argument had an illegal value
           *          > 0: if info = i, A(i,i) is exactly zero.  The triangular
           *               matrix is singular and its inverse can not be computed.
           */
          dtrtri_(&cL, &cU, &independent, M.array(), &LDA, &INFO);

          if (INFO < 0) fatalError();

          C_INT32 j, k;

          // Compute Link_0 = inverse(R_1,1) * R_1,2
          // :TODO: Use dgemm
          C_FLOAT64 * pTmp1 = array();
          C_FLOAT64 * pTmp2;
          C_FLOAT64 * pTmp3;

          for (j = 0; j < NumRows - independent; j++)
            for (i = 0; i < independent; i++, pTmp1++)
              {
                pTmp2 = &M(j + independent, i);
                pTmp3 = &M(i, i);

                // assert(&mL(j, i) == pTmp3);
                *pTmp1 = 0.0;

                for (k = i; k < independent; k++, pTmp2++, pTmp3 += NumCols)
                  {
                    // assert(&M(j + independent, k) == pTmp2);
                    // assert(&M(k, i) == pTmp3);

                    *pTmp1 += *pTmp3 * *pTmp2;
                  }

                if (fabs(*pTmp1) < 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon()) *pTmp1 = 0.0;
              }
        }
    }

  completePivotInformation();

  return success;
}

void CLinkMatrix::clearPivoting

(

)

Clear the pivot and swap vectors

Definition at line 384 of file CLinkMatrix.cpp.

References CVectorCore< CType >::array(), completePivotInformation(), mRowPivots, and CVectorCore< CType >::size().

Referenced by CModel::compile().

{
  size_t * pPivot = mRowPivots.array();
  size_t * pPivotEnd = pPivot + mRowPivots.size();
  size_t Index = 0;

  for (; pPivot != pPivotEnd; ++pPivot, ++Index)
    {
      *pPivot = Index;
    }

  completePivotInformation();
}

void CLinkMatrix::completePivotInformation ( )

private

Complete the pivot information.

Definition at line 340 of file CLinkMatrix.cpp.

References CVectorCore< CType >::array(), C_INT, mPivotInverse, mRowPivots, mSwapVector, CVector< CType >::resize(), and CVectorCore< CType >::size().

Referenced by build(), and clearPivoting().

{
  // We need to convert the pivot vector into a swap vector.
  mPivotInverse.resize(mRowPivots.size());

  size_t * pCurrentIndex = mPivotInverse.array();
  size_t * pEnd  = pCurrentIndex + mRowPivots.size();

  for (size_t i = 0; pCurrentIndex != pEnd; ++pCurrentIndex, ++i)
    {
      *pCurrentIndex = i;
    }

  CVector< size_t > CurrentColumn(mPivotInverse);
  size_t * pCurrentColumn = CurrentColumn.array();

  mSwapVector.resize(mRowPivots.size());
  C_INT * pJPVT = mSwapVector.array();
  pCurrentIndex = mPivotInverse.array();
  const size_t * pPivot = mRowPivots.array();

  for (; pCurrentIndex != pEnd; ++pPivot, ++pJPVT, ++pCurrentIndex, ++pCurrentColumn)
    {
      // Swap Index
      size_t * pToIndex = & mPivotInverse[*pPivot];
      size_t * pFromIndex = & mPivotInverse[*pCurrentColumn];

      // Swap Column
      size_t * pToColumn = & CurrentColumn[*pToIndex];
      size_t * pFromColumn = pCurrentColumn;

      // Swap
      *pJPVT =  *pToIndex + 1;

      size_t tmp = *pFromIndex;
      *pFromIndex = *pToIndex;
      *pToIndex = tmp;

      tmp = *pFromColumn;
      *pFromColumn = *pToColumn;
      *pToColumn = tmp;
    }
}

bool CLinkMatrix::doColumnPivot

(

CMatrix< C_FLOAT64 > &

matrix

)

const

Do the column pivot

Parameters

CMatrix<

C_FLOAT64 > & matrix

Returns

bool success

Definition at line 552 of file CLinkMatrix.cpp.

References applyColumnPivot().

Referenced by CMCAMethod::calculateUnscaledConcentrationCC().

{
  return applyColumnPivot(matrix, 1);
}

bool CLinkMatrix::doRowPivot

(

CMatrix< C_FLOAT64 > &

matrix

)

const

Do the row pivot

Parameters

CMatrix<

C_FLOAT64 > & matrix

Returns

bool success

Definition at line 496 of file CLinkMatrix.cpp.

References applyRowPivot(), and mRowPivots.

Referenced by CModel::buildLinkZero(), and CMCAMethod::createLinkMatrix().

{
  return applyRowPivot(matrix, mRowPivots);
}

size_t CLinkMatrix::getNumDependent

(

)

const

Retrieve the number of linear dependent rows of the input matrix

Returns

size_t numDependent

Definition at line 491 of file CLinkMatrix.cpp.

References CMatrix< C_FLOAT64 >::numRows().

Referenced by leftMultiply().

{
  return numRows();
}

const size_t & CLinkMatrix::getNumIndependent

(

)

const

Retrieve the number of linear independent rows of the input matrix

Returns

const size_t & numIndependent

Definition at line 486 of file CLinkMatrix.cpp.

References mIndependent.

Referenced by CModel::buildLinkZero(), CMCAMethod::createLinkMatrix(), and rightMultiply().

{
  return mIndependent;
}

const CVector< size_t > & CLinkMatrix::getRowPivots

(

)

const

Retrieve the pivot vector used to create the link matrix

Returns

const CVector< size_t > & rowPivots

Definition at line 39 of file CLinkMatrix.cpp.

References mRowPivots.

{
  return mRowPivots;
}

bool CLinkMatrix::leftMultiply	(	const CMatrix< C_FLOAT64 > &	M,
		CMatrix< C_FLOAT64 > &	P
	)		const

Left multiply the given matrix M with L, i.e., P = L * M

Parameters

const	CMatrix< C_FLOAT64> & M
CMatrix<	C_FLOAT64> & P

Returns

bool success;

Definition at line 445 of file CLinkMatrix.cpp.

References CMatrix< CType >::array(), CMatrix< C_FLOAT64 >::array(), C_FLOAT64, C_INT, dgemm_(), getNumDependent(), K, mRowPivots, CMatrix< C_FLOAT64 >::numCols(), CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), CMatrix< CType >::resize(), CVectorCore< CType >::size(), and CMatrix< CType >::size().

Referenced by CMCAMethod::calculateUnscaledConcentrationCC().

{
  bool success = true;

  if (m.numRows() != numCols())
    {
      return false;
    }

  // p := L * m = (I, L0) * m = (m, L0 * m) = (p1, p2)
  p.resize(mRowPivots.size(), m.numCols());
  p = 0.0;

  // p1 := m
  memcpy(p.array(), m.array(), sizeof(C_FLOAT64) * m.size());

  // p2 := L0 * m
  // p2 (numDependent x m.numCols())
  char T = 'N';
  C_INT M = (C_INT) m.numCols(); /* LDA, LDC */
  C_INT N = (C_INT) getNumDependent();
  C_INT K = (C_INT) numCols();

  C_INT LDA = (C_INT) m.numCols();
  C_INT LDB = (C_INT) numCols();
  C_INT LDC = (C_INT) p.numCols();

  C_FLOAT64 Alpha = 1.0;
  C_FLOAT64 Zero = 0.0;

  // DGEMM (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
  // C := alpha A B + beta C
  dgemm_(&T, &T, &M, &N, &K, &Alpha,
         const_cast< C_FLOAT64 * >(m.array()), &LDA,
         const_cast< C_FLOAT64 * >(array()), &LDB,
         &Zero, p.array() + m.size(), &LDC);

  return success;
}

bool CLinkMatrix::rightMultiply	(	const C_FLOAT64 &	alpha,
		const CMatrix< C_FLOAT64 > &	M,
		CMatrix< C_FLOAT64 > &	P
	)		const

Right multiply the given matrix M with L, i.e., P = alpha M * L. Note the columns of M must be in the same order as L.

Parameters

const	C_FLOAT64 & alpha
const	CMatrix< C_FLOAT64> & M
CMatrix<	C_FLOAT64> & P

Returns

bool success;

Definition at line 398 of file CLinkMatrix.cpp.

References CMatrix< CType >::array(), CMatrix< C_FLOAT64 >::array(), C_FLOAT64, C_INT, dgemm_(), getNumIndependent(), K, mRowPivots, CMatrix< CType >::numCols(), CMatrix< C_FLOAT64 >::numCols(), CMatrix< CType >::numRows(), CMatrix< C_FLOAT64 >::numRows(), CMatrix< CType >::resize(), CVectorCore< CType >::size(), and CMatrix< CType >::size().

Referenced by CMCAMethod::calculateUnscaledConcentrationCC().

{
  bool success = true;

  if (m.numCols() != mRowPivots.size())
    {
      return false;
    }

  // p := alpha * (m1, m2) * (I, L0')' = alpha * m1 + alpha * m2 * L0

  p.resize(m.numRows(), getNumIndependent());

  char T = 'N';
  C_INT M = (C_INT) p.numCols();
  C_INT N = (C_INT) p.numRows();
  C_INT K = (C_INT) numRows();

  C_INT LDA = (C_INT) numCols();
  C_INT LDB = (C_INT) m.numCols();
  C_INT LDC = (C_INT) p.numCols();

  // p := m1

  C_FLOAT64 *pRowP = p.array();
  const C_FLOAT64 *pRowM = m.array();
  C_FLOAT64 *pEnd = p.array() + p.size();

  for (; pRowP < pEnd; pRowP += LDC, pRowM += LDB)
    {
      memcpy(pRowP, pRowM, sizeof(C_FLOAT64) * LDC);
    }

  // DGEMM (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
  // C := alpha A B + beta C

  dgemm_(&T, &T, &M, &N, &K,
         const_cast< C_FLOAT64 * >(&alpha),
         const_cast< C_FLOAT64 * >(array()), &LDA,
         const_cast< C_FLOAT64 * >(m.array()) + LDA, &LDB,
         const_cast< C_FLOAT64 * >(&alpha), p.array(), &LDC);

  return success;
}

bool CLinkMatrix::undoColumnPivot

(

CMatrix< C_FLOAT64 > &

matrix

)

const

Undo the column pivot

Parameters

CMatrix<

C_FLOAT64 > & matrix

Returns

bool success

Definition at line 557 of file CLinkMatrix.cpp.

References applyColumnPivot().

Referenced by CMCAMethod::calculateUnscaledConcentrationCC().

{
  return applyColumnPivot(matrix, -1);
}

bool CLinkMatrix::undoRowPivot

(

CMatrix< C_FLOAT64 > &

matrix

)

const

Undo the row pivot

Parameters

CMatrix<

C_FLOAT64 > & matrix

Returns

bool success

Definition at line 501 of file CLinkMatrix.cpp.

References applyRowPivot(), and mPivotInverse.

Referenced by CMCAMethod::calculateUnscaledConcentrationCC().

{
  return applyRowPivot(matrix, mPivotInverse);
}

Member Data Documentation

size_t CLinkMatrix::mIndependent

private

The number of linear independent rows.

Definition at line 200 of file CLinkMatrix.h.

Referenced by build(), and getNumIndependent().

CVector< size_t > CLinkMatrix::mPivotInverse

private

The pivot vector used for undoing row swapping

Definition at line 190 of file CLinkMatrix.h.

Referenced by completePivotInformation(), and undoRowPivot().

CVector< size_t > CLinkMatrix::mRowPivots

private

The row pivoting performed to create the link matrix

Definition at line 185 of file CLinkMatrix.h.

Referenced by applyColumnPivot(), applyRowPivot(), build(), clearPivoting(), completePivotInformation(), doRowPivot(), getRowPivots(), leftMultiply(), and rightMultiply().

CVector< C_INT > CLinkMatrix::mSwapVector

private

The swap vector used for column swapping

Definition at line 195 of file CLinkMatrix.h.

Referenced by applyColumnPivot(), and completePivotInformation().

---

The documentation for this class was generated from the following files:

• CLinkMatrix.h
• CLinkMatrix.cpp