COPASI API  4.16.103
CLinkMatrix.cpp
Go to the documentation of this file.
1 // Copyright (C) 2013 - 2015 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 #include <cmath>
7 
8 #include "copasi.h"
9 
10 #include "CLinkMatrix.h"
11 
12 #include "CCopasiVector.h"
13 
14 #include "lapack/blaswrap.h"
15 #include "lapack/lapackwrap.h"
16 
17 // Define to output debugging information
18 // #define DEBUG_MATRIX
19 
21  CMatrix< C_FLOAT64 >(),
22  mRowPivots(),
23  mPivotInverse(),
24  mSwapVector(),
25  mIndependent(0)
26 {}
27 
29  CMatrix< C_FLOAT64 >(src),
30  mRowPivots(src.mRowPivots),
31  mPivotInverse(src.mPivotInverse),
32  mSwapVector(src.mSwapVector),
33  mIndependent(src.mIndependent)
34 {}
35 
37 {}
38 
40 {
41  return mRowPivots;
42 }
43 
44 bool CLinkMatrix::build(const CMatrix< C_FLOAT64 > & matrix, size_t maxRank)
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 }
339 
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 }
383 
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 }
397 
399  const CMatrix< C_FLOAT64> & m,
400  CMatrix< C_FLOAT64> & p) const
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 }
444 
446  CMatrix< C_FLOAT64> & p) const
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 }
485 
486 const size_t & CLinkMatrix::getNumIndependent() const
487 {
488  return mIndependent;
489 }
490 
492 {
493  return numRows();
494 }
495 
497 {
498  return applyRowPivot(matrix, mRowPivots);
499 }
500 
502 {
503  return applyRowPivot(matrix, mPivotInverse);
504 }
505 
507  const CVector< size_t > & pivots) const
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 }
551 
553 {
554  return applyColumnPivot(matrix, 1);
555 }
556 
558 {
559  return applyColumnPivot(matrix, -1);
560 }
561 
563  const C_INT & incr) const
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 }
581 
582 //**********************************************************************
583 // CLinkMatrixView
584 //**********************************************************************
585 
588 
590  mpA(&A),
591  mpNumIndependent(&A.getNumIndependent())
593 
596 
599 {
600  mpA = rhs.mpA;
602 
603  return *this;
604 }
605 
607 {
608  return *mpNumIndependent + mpA->numRows();
609 }
610 
612 {
613  return mpA->numCols();
614 }
615 
616 std::ostream &operator<<(std::ostream &os,
617  const CLinkMatrixView & A)
618 {
619  size_t i, imax = A.numRows();
620  size_t j, jmax = A.numCols();
621  os << "Matrix(" << imax << "x" << jmax << ")" << std::endl;
622 
623  for (i = 0; i < imax; i++)
624  {
625  for (j = 0; j < jmax; j++)
626  os << "\t" << A(i, j);
627 
628  os << std::endl;
629  }
630 
631  return os;
632 }
#define C_INT
Definition: copasi.h:115
bool build(const CMatrix< C_FLOAT64 > &matrix, size_t maxRank=C_INVALID_INDEX)
Definition: CLinkMatrix.cpp:44
size_t getNumDependent() const
CVector< size_t > mRowPivots
Definition: CLinkMatrix.h:185
#define K
#define TAU
virtual ~CLinkMatrix()
Definition: CLinkMatrix.cpp:36
virtual size_t numRows() const
Definition: CMatrix.h:138
#define fatalError()
bool applyColumnPivot(CMatrix< C_FLOAT64 > &matrix, const C_INT &incr) const
bool leftMultiply(const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
size_t mIndependent
Definition: CLinkMatrix.h:200
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
void clearPivoting()
#define C_INT32
Definition: copasi.h:90
const size_t & getNumIndependent() const
size_t numCols() const
bool rightMultiply(const C_FLOAT64 &alpha, const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
std::ostream & operator<<(std::ostream &os, const CLinkMatrixView &A)
CVector< C_INT > mSwapVector
Definition: CLinkMatrix.h:195
int dtrtri_(char *uplo, char *diag, integer *n, doublereal *a, integer *lda, integer *info)
CLinkMatrixView(const CLinkMatrix &A)
static const C_FLOAT64 mUnit
Definition: CLinkMatrix.h:212
const CVector< size_t > & getRowPivots() const
Definition: CLinkMatrix.cpp:39
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
bool doRowPivot(CMatrix< C_FLOAT64 > &matrix) const
const CLinkMatrix * mpA
Definition: CLinkMatrix.h:209
bool applyRowPivot(CVectorCore< CType > &vector) const
Definition: CLinkMatrix.h:115
int dlaswp_(integer *n, doublereal *a, integer *lda, integer *k1, integer *k2, integer *ipiv, integer *incx)
int dgeqp3_(integer *m, integer *n, doublereal *a, integer *lda, integer *jpvt, doublereal *tau, doublereal *work, integer *lwork, integer *info)
size_t numRows() const
int dlaic1_(integer *job, integer *j, doublereal *x, doublereal *sest, doublereal *w, doublereal *gamma, doublereal *sestpr, doublereal *s, doublereal *c__)
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)
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
bool doColumnPivot(CMatrix< C_FLOAT64 > &matrix) const
size_t size() const
Definition: CVector.h:100
#define C_FLOAT64
Definition: copasi.h:92
CType * array()
Definition: CVector.h:139
virtual size_t size() const
Definition: CMatrix.h:132
bool undoColumnPivot(CMatrix< C_FLOAT64 > &matrix) const
void completePivotInformation()
static const C_FLOAT64 mZero
Definition: CLinkMatrix.h:211
const size_t * mpNumIndependent
Definition: CLinkMatrix.h:210
virtual size_t numCols() const
Definition: CMatrix.h:144
CVector< size_t > mPivotInverse
Definition: CLinkMatrix.h:190
CLinkMatrixView & operator=(const CLinkMatrixView &rhs)
bool undoRowPivot(CMatrix< C_FLOAT64 > &matrix) const
virtual CType * array()
Definition: CMatrix.h:337
#define min(a, b)
Definition: f2c.h:175
#define CONSTRUCTOR_TRACE
Definition: copasi.h:202
#define max(a, b)
Definition: f2c.h:176