COPASI API  4.16.103
Public Member Functions | Private Member Functions | Private Attributes | List of all members
CLinkMatrix Class Reference

#include <CLinkMatrix.h>

Inheritance diagram for CLinkMatrix:
Inheritance graph
[legend]
Collaboration diagram for CLinkMatrix:
Collaboration graph
[legend]

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_FLOAT64array ()
 
virtual const C_FLOAT64array () 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 elementTypeoperator() (const size_t &row, const size_t &col)
 
virtual const elementTypeoperator() (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_FLOAT64operator[] (size_t row)
 
virtual const C_FLOAT64operator[] (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_INTmSwapVector
 

Additional Inherited Members

- Public Types inherited from CMatrix< C_FLOAT64 >
typedef C_FLOAT64 elementType
 
- Protected Attributes inherited from CMatrix< C_FLOAT64 >
C_FLOAT64mArray
 
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.

20  :
22  mRowPivots(),
23  mPivotInverse(),
24  mSwapVector(),
25  mIndependent(0)
26 {}
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
size_t mIndependent
Definition: CLinkMatrix.h:200
CVector< C_INT > mSwapVector
Definition: CLinkMatrix.h:195
CVector< size_t > mPivotInverse
Definition: CLinkMatrix.h:190
CLinkMatrix::CLinkMatrix ( const CLinkMatrix src)

Copy constructor

Definition at line 28 of file CLinkMatrix.cpp.

28  :
34 {}
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
size_t mIndependent
Definition: CLinkMatrix.h:200
CVector< C_INT > mSwapVector
Definition: CLinkMatrix.h:195
CVector< size_t > mPivotInverse
Definition: CLinkMatrix.h:190
CLinkMatrix::~CLinkMatrix ( )
virtual

Destructor

Definition at line 36 of file CLinkMatrix.cpp.

37 {}

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
constC_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().

564 {
565  if (matrix.numCols() < mRowPivots.size())
566  {
567  return false;
568  }
569 
570  C_INT N = matrix.numRows();
571  C_INT LDA = matrix.numCols();
572  C_INT K1 = 1;
573  C_INT K2 = mRowPivots.size();
574 
575  dlaswp_(&N, matrix.array(), &LDA, &K1, &K2,
576  const_cast< C_INT * >(mSwapVector.array()),
577  const_cast< C_INT * >(&incr));
578 
579  return true;
580 }
#define C_INT
Definition: copasi.h:115
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
virtual size_t numRows() const
Definition: CMatrix.h:138
CVector< C_INT > mSwapVector
Definition: CLinkMatrix.h:195
int dlaswp_(integer *n, doublereal *a, integer *lda, integer *k1, integer *k2, integer *ipiv, integer *incx)
size_t size() const
Definition: CVector.h:100
CType * array()
Definition: CVector.h:139
virtual size_t numCols() const
Definition: CMatrix.h:144
virtual CType * array()
Definition: CMatrix.h:337
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().

116  {
117  if (vector.size() < mRowPivots.size())
118  {
119  return false;
120  }
121 
122  CVector< bool > Applied(mRowPivots.size());
123  Applied = false;
124 
125  CType Tmp;
126 
127  size_t i, imax = mRowPivots.size();
128  size_t to;
129  size_t from;
130 
131  for (i = 0; i < imax; i++)
132  if (!Applied[i])
133  {
134  to = i;
135  from = mRowPivots[to];
136 
137  if (from != i)
138  {
139  Tmp = vector[to];
140 
141  while (from != i)
142  {
143  vector[to] = vector[from];
144  Applied[to] = true;
145 
146  to = from;
147  from = mRowPivots[to];
148  }
149 
150  vector[to] = Tmp;
151  }
152 
153  Applied[to] = true;
154  }
155 
156  return true;
157  }
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
size_t size() const
Definition: CVector.h:100
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
constCVector< 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().

508 {
509  if (matrix.numRows() < pivots.size())
510  {
511  return false;
512  }
513 
514  CVector< bool > Applied(pivots.size());
515  Applied = false;
516 
517  CVector< C_FLOAT64 > Tmp(matrix.numCols());
518 
519  size_t i, imax = pivots.size();
520  size_t to;
521  size_t from;
522  size_t numCols = matrix.numCols();
523 
524  for (i = 0; i < imax; i++)
525  if (!Applied[i])
526  {
527  to = i;
528  from = pivots[to];
529 
530  if (from != i)
531  {
532  memcpy(Tmp.array(), matrix[to], sizeof(C_FLOAT64) * numCols);
533 
534  while (from != i)
535  {
536  memcpy(matrix[to], matrix[from], sizeof(C_FLOAT64) * numCols);
537  Applied[to] = true;
538 
539  to = from;
540  from = pivots[to];
541  }
542 
543  memcpy(matrix[to], Tmp.array(), sizeof(C_FLOAT64) * numCols);
544  }
545 
546  Applied[to] = true;
547  }
548 
549  return true;
550 }
virtual size_t numRows() const
Definition: CMatrix.h:138
size_t size() const
Definition: CVector.h:100
#define C_FLOAT64
Definition: copasi.h:92
virtual size_t numCols() const
Definition: CMatrix.h:144
virtual CType * array()
Definition: CMatrix.h:337
bool CLinkMatrix::build ( const CMatrix< C_FLOAT64 > &  matrix,
size_t  maxRank = C_INVALID_INDEX 
)

Build the link matrix for the given matrix

Parameters
constCMatrix< C_FLOAT64 > & matrix
size_tmaxRank (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().

45 {
46  bool success = true;
47 
48  CMatrix< C_FLOAT64 > M(matrix);
49 
50  C_INT NumCols = (C_INT) M.numCols();
51  C_INT NumRows = (C_INT) M.numRows();
52  C_INT LDA = std::max<C_INT>(1, NumCols);
53 
54  CVector< C_INT > JPVT(NumRows);
55  JPVT = 0;
56 
57  C_INT32 Dim = std::min(NumCols, NumRows);
58 
59  if (Dim == 0)
60  {
61  resize(NumRows, 0);
62 
63  mIndependent = 0;
64 
65  C_INT32 i;
66  mRowPivots.resize(NumRows);
67 
68  for (i = 0; i < NumRows; i++)
69  mRowPivots[i] = i;
70  }
71  else
72  {
74 
75  CVector< C_FLOAT64 > WORK(1);
76  C_INT LWORK = -1;
77  C_INT INFO;
78 
79  // QR factorization of the stoichiometry matrix
80  /*
81  * -- LAPACK routine (version 3.0) --
82  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
83  * Courant Institute, Argonne National Lab, and Rice University
84  * June 30, 1999
85  *
86  * Purpose
87  * =======
88  *
89  * DGEQP3 computes a QR factorization with column pivoting of a
90  * matrix A: A*P = Q*R using Level 3 BLAS.
91  *
92  * Arguments
93  * =========
94  *
95  * M (input) INTEGER
96  * The number of rows of the matrix A. M >= 0.
97  *
98  * N (input) INTEGER
99  * The number of columns of the matrix A. N >= 0.
100  *
101  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
102  * On entry, the M-by-N matrix A.
103  * On exit, the upper triangle of the array contains the
104  * min(M,N)-by-N upper trapezoidal matrix R; the elements below
105  * the diagonal, together with the array TAU, represent the
106  * orthogonal matrix Q as a product of min(M,N) elementary
107  * reflectors.
108  *
109  * LDA (input) INTEGER
110  * The leading dimension of the array A. LDA >= max(1,M).
111  *
112  * JPVT (input/output) INTEGER array, dimension (N)
113  * On entry, if JPVT(J).ne.0, the J-th column of A is permuted
114  * to the front of A*P (a leading column); if JPVT(J)=0,
115  * the J-th column of A is a free column.
116  * On exit, if JPVT(J)=K, then the J-th column of A*P was the
117  * the K-th column of A.
118  *
119  * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
120  * The scalar factors of the elementary reflectors.
121  *
122  * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
123  * On exit, if INFO=0, WORK(1) returns the optimal LWORK.
124  *
125  * LWORK (input) INTEGER
126  * The dimension of the array WORK. LWORK >= 3*N+1.
127  * For optimal performance LWORK >= 2*N+(N+1)*NB, where NB
128  * is the optimal blocksize.
129  *
130  * If LWORK = -1, then a workspace query is assumed; the routine
131  * only calculates the optimal size of the WORK array, returns
132  * this value as the first entry of the WORK array, and no error
133  * message related to LWORK is issued by XERBLA.
134  *
135  * INFO (output) INTEGER
136  * = 0: successful exit.
137  * < 0: if INFO = -i, the i-th argument had an illegal value.
138  *
139  * Further Details
140  * ===============
141  *
142  * The matrix Q is represented as a product of elementary reflectors
143  *
144  * Q = H(1) H(2) . . . H(k), where k = min(m,n).
145  *
146  * Each H(i) has the form
147  *
148  * H(i) = I - tau * v * v'
149  *
150  * where tau is a real/complex scalar, and v is a real/complex vector
151  * with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
152  * A(i+1:m,i), and tau in TAU(i).
153  *
154  * Based on contributions by
155  * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
156  * X. Sun, Computer Science Dept., Duke University, USA
157  *
158  */
159 
160  dgeqp3_(&NumCols, &NumRows, M.array(), &LDA,
161  JPVT.array(), TAU.array(), WORK.array(), &LWORK, &INFO);
162 
163  if (INFO < 0) fatalError();
164 
165  LWORK = (C_INT) WORK[0];
166  WORK.resize(LWORK);
167 
168  dgeqp3_(&NumCols, &NumRows, M.array(), &LDA,
169  JPVT.array(), TAU.array(), WORK.array(), &LWORK, &INFO);
170 
171  if (INFO < 0) fatalError();
172 
173  C_INT32 i;
174  mRowPivots.resize(NumRows);
175 
176  for (i = 0; i < NumRows; i++)
177  mRowPivots[i] = JPVT[i] - 1;
178 
179  /**
180  * Determine the rank of the matrix
181  * This code is copied from dgelsy.f
182  */
183 
184  C_INT independent = 0;
185 
186  if (fabs(M(0, 0)) != 0.0)
187  {
188  ++independent;
189 
190  C_INT imax = 1;
191  C_INT imin = 2;
192  C_INT mn = Dim;
193 
194  C_FLOAT64 * pIsmin = WORK.array();
195  C_FLOAT64 * pIsmax = pIsmin + Dim;
196 
197  *pIsmin = 1.0;
198  *pIsmax = 1.0;
199 
200  C_FLOAT64 smax = fabs(M(0, 0));
201  C_FLOAT64 smin = smax;
202  C_FLOAT64 RCOND = std::max(NumRows, NumCols) * smax * std::numeric_limits< C_FLOAT64 >::epsilon();
203 
204  C_FLOAT64 sminpr, s1, c1, smaxpr, s2, c2;
205 
206  while (independent < mn && independent < maxRank)
207  {
208  dlaic1_(&imin, &independent, pIsmin, &smin, &M(independent, 0), &M(independent, independent), &sminpr, &s1, &c1);
209  dlaic1_(&imax, &independent, pIsmax, &smax, &M(independent, 0), &M(independent, independent), &smaxpr, &s2, &c2);
210 
211  if (smaxpr * RCOND > sminpr)
212  {
213  break;
214  }
215 
216  for (size_t i = 0; i < independent - 1; ++i)
217  {
218  *(pIsmin + i) *= s1;
219  *(pIsmax + i) *= s2;
220  }
221 
222  *(pIsmin + independent) = c1;
223  *(pIsmax + independent) = c2;
224 
225  smin = sminpr;
226  smax = smaxpr;
227 
228  ++independent;
229  }
230  }
231 
232  mIndependent = independent;
233 
234  // Resize mL
235  resize(NumRows - independent, independent);
236 
237  if (NumRows != independent && independent != 0)
238  {
239  /* to take care of differences between fortran's and c's memory access,
240  we need to take the transpose, i.e.,the upper triangular */
241  char cL = 'U';
242  char cU = 'N'; /* values in the diagonal of R */
243 
244  // Calculate Row Echelon form of R.
245  // First invert R_1,1
246  /* int dtrtri_(char *uplo,
247  * char *diag,
248  * integer *n,
249  * doublereal * A,
250  * integer *lda,
251  * integer *info);
252  * -- LAPACK routine (version 3.0) --
253  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
254  * Courant Institute, Argonne National Lab, and Rice University
255  * March 31, 1993
256  *
257  * Purpose
258  * =======
259  *
260  * DTRTRI computes the inverse of a real upper or lower triangular
261  * matrix A.
262  *
263  * This is the Level 3 BLAS version of the algorithm.
264  *
265  * Arguments
266  * =========
267  *
268  * uplo (input) CHARACTER*1
269  * = 'U': A is upper triangular;
270  * = 'L': A is lower triangular.
271  *
272  * diag (input) CHARACTER*1
273  * = 'N': A is non-unit triangular;
274  * = 'U': A is unit triangular.
275  *
276  * n (input) INTEGER
277  * The order of the matrix A. n >= 0.
278  *
279  * A (input/output) DOUBLE PRECISION array, dimension (lda,n)
280  * On entry, the triangular matrix A. If uplo = 'U', the
281  * leading n-by-n upper triangular part of the array A contains
282  * the upper triangular matrix, and the strictly lower
283  * triangular part of A is not referenced. If uplo = 'L', the
284  * leading n-by-n lower triangular part of the array A contains
285  * the lower triangular matrix, and the strictly upper
286  * triangular part of A is not referenced. If diag = 'U', the
287  * diagonal elements of A are also not referenced and are
288  * assumed to be 1.
289  * On exit, the (triangular) inverse of the original matrix, in
290  * the same storage format.
291  *
292  * lda (input) INTEGER
293  * The leading dimension of the array A. lda >= max(1,n).
294  *
295  * info (output) INTEGER
296  * = 0: successful exit
297  * < 0: if info = -i, the i-th argument had an illegal value
298  * > 0: if info = i, A(i,i) is exactly zero. The triangular
299  * matrix is singular and its inverse can not be computed.
300  */
301  dtrtri_(&cL, &cU, &independent, M.array(), &LDA, &INFO);
302 
303  if (INFO < 0) fatalError();
304 
305  C_INT32 j, k;
306 
307  // Compute Link_0 = inverse(R_1,1) * R_1,2
308  // :TODO: Use dgemm
309  C_FLOAT64 * pTmp1 = array();
310  C_FLOAT64 * pTmp2;
311  C_FLOAT64 * pTmp3;
312 
313  for (j = 0; j < NumRows - independent; j++)
314  for (i = 0; i < independent; i++, pTmp1++)
315  {
316  pTmp2 = &M(j + independent, i);
317  pTmp3 = &M(i, i);
318 
319  // assert(&mL(j, i) == pTmp3);
320  *pTmp1 = 0.0;
321 
322  for (k = i; k < independent; k++, pTmp2++, pTmp3 += NumCols)
323  {
324  // assert(&M(j + independent, k) == pTmp2);
325  // assert(&M(k, i) == pTmp3);
326 
327  *pTmp1 += *pTmp3 * *pTmp2;
328  }
329 
330  if (fabs(*pTmp1) < 100.0 * std::numeric_limits< C_FLOAT64 >::epsilon()) *pTmp1 = 0.0;
331  }
332  }
333  }
334 
336 
337  return success;
338 }
#define C_INT
Definition: copasi.h:115
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
#define TAU
#define fatalError()
size_t mIndependent
Definition: CLinkMatrix.h:200
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define C_INT32
Definition: copasi.h:90
int dtrtri_(char *uplo, char *diag, integer *n, doublereal *a, integer *lda, integer *info)
int dgeqp3_(integer *m, integer *n, doublereal *a, integer *lda, integer *jpvt, doublereal *tau, doublereal *work, integer *lwork, integer *info)
int dlaic1_(integer *job, integer *j, doublereal *x, doublereal *sest, doublereal *w, doublereal *gamma, doublereal *sestpr, doublereal *s, doublereal *c__)
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
#define C_FLOAT64
Definition: copasi.h:92
void completePivotInformation()
virtual C_FLOAT64 * array()
Definition: CMatrix.h:337
#define min(a, b)
Definition: f2c.h:175
#define max(a, b)
Definition: f2c.h:176
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().

385 {
386  size_t * pPivot = mRowPivots.array();
387  size_t * pPivotEnd = pPivot + mRowPivots.size();
388  size_t Index = 0;
389 
390  for (; pPivot != pPivotEnd; ++pPivot, ++Index)
391  {
392  *pPivot = Index;
393  }
394 
396 }
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
size_t size() const
Definition: CVector.h:100
CType * array()
Definition: CVector.h:139
void 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().

341 {
342  // We need to convert the pivot vector into a swap vector.
344 
345  size_t * pCurrentIndex = mPivotInverse.array();
346  size_t * pEnd = pCurrentIndex + mRowPivots.size();
347 
348  for (size_t i = 0; pCurrentIndex != pEnd; ++pCurrentIndex, ++i)
349  {
350  *pCurrentIndex = i;
351  }
352 
353  CVector< size_t > CurrentColumn(mPivotInverse);
354  size_t * pCurrentColumn = CurrentColumn.array();
355 
357  C_INT * pJPVT = mSwapVector.array();
358  pCurrentIndex = mPivotInverse.array();
359  const size_t * pPivot = mRowPivots.array();
360 
361  for (; pCurrentIndex != pEnd; ++pPivot, ++pJPVT, ++pCurrentIndex, ++pCurrentColumn)
362  {
363  // Swap Index
364  size_t * pToIndex = & mPivotInverse[*pPivot];
365  size_t * pFromIndex = & mPivotInverse[*pCurrentColumn];
366 
367  // Swap Column
368  size_t * pToColumn = & CurrentColumn[*pToIndex];
369  size_t * pFromColumn = pCurrentColumn;
370 
371  // Swap
372  *pJPVT = *pToIndex + 1;
373 
374  size_t tmp = *pFromIndex;
375  *pFromIndex = *pToIndex;
376  *pToIndex = tmp;
377 
378  tmp = *pFromColumn;
379  *pFromColumn = *pToColumn;
380  *pToColumn = tmp;
381  }
382 }
#define C_INT
Definition: copasi.h:115
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
CVector< C_INT > mSwapVector
Definition: CLinkMatrix.h:195
size_t size() const
Definition: CVector.h:100
CType * array()
Definition: CVector.h:139
CVector< size_t > mPivotInverse
Definition: CLinkMatrix.h:190
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().

553 {
554  return applyColumnPivot(matrix, 1);
555 }
bool applyColumnPivot(CMatrix< C_FLOAT64 > &matrix, const C_INT &incr) const
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().

497 {
498  return applyRowPivot(matrix, mRowPivots);
499 }
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
bool applyRowPivot(CVectorCore< CType > &vector) const
Definition: CLinkMatrix.h:115
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().

492 {
493  return numRows();
494 }
virtual size_t numRows() const
Definition: CMatrix.h:138
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().

487 {
488  return mIndependent;
489 }
size_t mIndependent
Definition: CLinkMatrix.h:200
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.

40 {
41  return mRowPivots;
42 }
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
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
constCMatrix< 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().

447 {
448  bool success = true;
449 
450  if (m.numRows() != numCols())
451  {
452  return false;
453  }
454 
455  // p := L * m = (I, L0) * m = (m, L0 * m) = (p1, p2)
456  p.resize(mRowPivots.size(), m.numCols());
457  p = 0.0;
458 
459  // p1 := m
460  memcpy(p.array(), m.array(), sizeof(C_FLOAT64) * m.size());
461 
462  // p2 := L0 * m
463  // p2 (numDependent x m.numCols())
464  char T = 'N';
465  C_INT M = (C_INT) m.numCols(); /* LDA, LDC */
466  C_INT N = (C_INT) getNumDependent();
467  C_INT K = (C_INT) numCols();
468 
469  C_INT LDA = (C_INT) m.numCols();
470  C_INT LDB = (C_INT) numCols();
471  C_INT LDC = (C_INT) p.numCols();
472 
473  C_FLOAT64 Alpha = 1.0;
474  C_FLOAT64 Zero = 0.0;
475 
476  // DGEMM (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
477  // C := alpha A B + beta C
478  dgemm_(&T, &T, &M, &N, &K, &Alpha,
479  const_cast< C_FLOAT64 * >(m.array()), &LDA,
480  const_cast< C_FLOAT64 * >(array()), &LDB,
481  &Zero, p.array() + m.size(), &LDC);
482 
483  return success;
484 }
#define C_INT
Definition: copasi.h:115
size_t getNumDependent() const
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
#define K
int dgemm_(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
size_t size() const
Definition: CVector.h:100
#define C_FLOAT64
Definition: copasi.h:92
virtual size_t numCols() const
Definition: CMatrix.h:144
virtual C_FLOAT64 * array()
Definition: CMatrix.h:337
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
constC_FLOAT64 & alpha
constCMatrix< 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().

401 {
402  bool success = true;
403 
404  if (m.numCols() != mRowPivots.size())
405  {
406  return false;
407  }
408 
409  // p := alpha * (m1, m2) * (I, L0')' = alpha * m1 + alpha * m2 * L0
410 
411  p.resize(m.numRows(), getNumIndependent());
412 
413  char T = 'N';
414  C_INT M = (C_INT) p.numCols();
415  C_INT N = (C_INT) p.numRows();
416  C_INT K = (C_INT) numRows();
417 
418  C_INT LDA = (C_INT) numCols();
419  C_INT LDB = (C_INT) m.numCols();
420  C_INT LDC = (C_INT) p.numCols();
421 
422  // p := m1
423 
424  C_FLOAT64 *pRowP = p.array();
425  const C_FLOAT64 *pRowM = m.array();
426  C_FLOAT64 *pEnd = p.array() + p.size();
427 
428  for (; pRowP < pEnd; pRowP += LDC, pRowM += LDB)
429  {
430  memcpy(pRowP, pRowM, sizeof(C_FLOAT64) * LDC);
431  }
432 
433  // DGEMM (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
434  // C := alpha A B + beta C
435 
436  dgemm_(&T, &T, &M, &N, &K,
437  const_cast< C_FLOAT64 * >(&alpha),
438  const_cast< C_FLOAT64 * >(array()), &LDA,
439  const_cast< C_FLOAT64 * >(m.array()) + LDA, &LDB,
440  const_cast< C_FLOAT64 * >(&alpha), p.array(), &LDC);
441 
442  return success;
443 }
#define C_INT
Definition: copasi.h:115
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
#define K
virtual size_t numRows() const
Definition: CMatrix.h:138
const size_t & getNumIndependent() const
int dgemm_(char *transa, char *transb, integer *m, integer *n, integer *k, doublereal *alpha, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *beta, doublereal *c, integer *ldc)
size_t size() const
Definition: CVector.h:100
#define C_FLOAT64
Definition: copasi.h:92
virtual size_t numCols() const
Definition: CMatrix.h:144
virtual C_FLOAT64 * array()
Definition: CMatrix.h:337
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().

558 {
559  return applyColumnPivot(matrix, -1);
560 }
bool applyColumnPivot(CMatrix< C_FLOAT64 > &matrix, const C_INT &incr) const
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().

502 {
503  return applyRowPivot(matrix, mPivotInverse);
504 }
bool applyRowPivot(CVectorCore< CType > &vector) const
Definition: CLinkMatrix.h:115
CVector< size_t > mPivotInverse
Definition: CLinkMatrix.h:190

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: