COPASI API  4.16.103
CluX.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2005 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <algorithm>
16 #include <cmath>
17 #include <limits>
18 
19 #include "copasi.h"
20 
21 #include "CluX.h"
22 #include "CMatrix.h"
23 #include "CVector.h"
24 #include "CProcessReport.h"
25 
26 #include "lapack/lapackwrap.h"
27 
28 // #define DEBUG_MATRIX
29 
31  CVector< size_t > & row,
32  CVector< size_t > & col,
33  CProcessReport * cb)
34 {
35  size_t Rows = A.numRows();
36  size_t Cols = A.numCols();
37  size_t Dim = std::min(Rows, Cols);
38 
39  if (row.size() != Rows) row.resize(Rows);
40 
41  if (col.size() != Cols) col.resize(Cols);
42 
43  size_t i = 0, j = 0, k = 0;
44  size_t rowP = 0, colP = Cols - 1;
45 
46  for (i = 0; i < Rows; i++)
47  row[i] = i;
48 
49  for (i = 0; i < Cols; i++)
50  col[i] = i;
51 
52  C_FLOAT64 Work = .5 * Dim * (Dim - 1);
53  C_FLOAT64 Completion = 0;
54 
55  size_t hProcess;
56 
57  if (cb)
58  hProcess = cb->addItem("LU decomposition...",
59  Completion,
60  &Work);
61 
62  C_FLOAT64 tmp;
63  C_FLOAT64 * pTmp1;
64  C_FLOAT64 * pTmp2;
65  C_FLOAT64 * pTmp3;
66 
67  for (i = 0; i < Dim; i++)
68  {
69  Completion += Dim - i;
70 
71  if (cb && !cb->progressItem(hProcess)) return false;
72 
73  // find pivot in column j and test for singularity.
74  while (true)
75  {
76  rowP = C_INVALID_INDEX;
77  tmp = 0.0;
78 
79  for (j = i, pTmp1 = &A(i, i); j < Rows; j++, pTmp1 += Cols)
80  if (fabs(*pTmp1) > tmp) // (fabs(A(j, i)) > t)
81  {
82  rowP = j;
83  tmp = fabs(*pTmp1); // fabs(A(j, i));
84  }
85 
86  // rowP now has the index of maximum element
87  // of column j, below the diagonal
88 
89  if (tmp < std::numeric_limits< C_FLOAT64 >::epsilon()) // now we have to swap colums to find a pivot
90  {
91  // Instead of blindly swapping with the last available
92  // column we first check whether it is suitable. If not
93  // we proceed with the lat but one ...
94  while (i < colP)
95  {
96  rowP = C_INVALID_INDEX;
97  tmp = 0.0;
98 
99  for (j = i, pTmp1 = &A(i, colP); j < Rows; j++, pTmp1 += Cols)
100  if (fabs(*pTmp1) > tmp) // (fabs(A(j, colP)) > t)
101  {
102  rowP = j;
103  tmp = fabs(*pTmp1); // fabs(A(j, colP));
104  }
105 
106  if (tmp >= std::numeric_limits< C_FLOAT64 >::epsilon()) break; // Found a suitable column
107 
108  // and row
109 
110  colP--; // proceed with previous column
111  }
112 
113  if (i >= colP)
114  {
115  if (cb) cb->finishItem(hProcess);
116 
117  return true;
118  }
119 
120  // Swap columns i and colP
121  pTmp1 = &A(0, i);
122  pTmp2 = &A(0, colP);
123 
124  for (k = 0; k < Rows; k++, pTmp1 += Cols, pTmp2 += Cols)
125  {
126  tmp = *pTmp1;
127  *pTmp1 = *pTmp2;
128  *pTmp2 = tmp;
129  }
130 
131  std::swap(col[i], col[colP]);
132  colP--;
133  }
134 
135  // Swap rows i and rowP
136  pTmp1 = &A(i, 0);
137  pTmp2 = &A(rowP, 0);
138 
139  for (k = 0; k < Cols; k++, pTmp1++, pTmp2++)
140  {
141  tmp = *pTmp1;
142  *pTmp1 = *pTmp2;
143  *pTmp2 = tmp;
144  }
145 
146  std::swap(row[i], row[rowP]);
147  break;
148  }
149 
150  if (i < Rows - 1) // compute elements i+1 <= k < M of i-th column
151  {
152  // note A(i,i) is guarranteed not to be zero
153  tmp = 1.0 / A(i, i);
154 
155  pTmp1 = &A(i + 1, i);
156 
157  for (k = i + 1; k < Rows; k++, pTmp1 += Cols)
158  * pTmp1 *= tmp;
159  }
160 
161  if (i < Dim - 1)
162  {
163  // rank-1 update to trailing submatrix: E = E - x*y;
164  //
165  // E is the region A(i+1:M, i+1:N)
166  // x is the column vector A(i+1:M,i)
167  // y is row vector A(i,i+1:N)
168 
169  size_t ii, jj;
170 
171  pTmp1 = &A(i + 1, i);
172  pTmp3 = &A(i + 1, i + 1);
173 
174  for (ii = i + 1; ii < Rows; ii++, pTmp1 += Cols, pTmp3 += i + 1)
175  {
176  pTmp2 = &A(i, i + 1);
177 
178  for (jj = i + 1; jj < Cols; jj++, pTmp2++, pTmp3++)
179  {
180  *pTmp3 -= *pTmp1 * *pTmp2;
181  // if (fabs(*pTmp3) < std::numeric_limits< C_FLOAT64 >::epsilon()) *pTmp3 = 0.0;
182  }
183  }
184  }
185 
186 #ifdef DEBUG_MATRIX
187  DebugFile << rowP << "\tx\t" << colP + 1 << std::endl;
188  DebugFile << "Row\t" << row << std::endl;
189  DebugFile << "Col\t" << col << std::endl;
190  DebugFile << A << std::endl;
191 #endif
192  }
193 
194  if (cb) cb->finishItem(hProcess);
195 
196  return true;
197 }
198 
199 #ifdef XXXX
201  CVector< size_t > & row,
202  CVector< size_t > & col,
203  CProcessReport * cb)
204 {
205  CVector< size_t > lCol;
206 
207  C_INT Rows = A.numRows();
208  C_INT Cols = A.numCols();
209 
210  if (row.size() != Rows) row.resize(Rows);
211 
212  if (col.size() != Cols) col.resize(Cols);
213 
214  if (lRow.size() != Rows) lRow.resize(Rows);
215 
216  if (lCol.size() != Rows) lCol.resize(Cols);
217 
218  size_t i = 0, j = 0, k = 0;
219  size_t ip, jp = Cols - 1;
220 
221  for (i = 0; i < Rows; i++)
222  {
223  row[i] = i;
224  lRow[i] = i;
225  }
226 
227  for (i = 0; i < Cols; i++)
228  {
229  col[i] = i;
230  lCol[i] = i;
231  }
232 
233  if (Rows == 0 || Cols == 0)
234  return false;
235 
236  CMatrix< C_FLOAT64 > T(Cols, Rows);
237 
238  for (i = 0; i < Cols; i++)
239  for (j = 0; j < Rows; j++)
240  T(i, j) = A(j, i);
241 
242  C_INT imax = std::min(Rows, Cols);
243  CVector<C_INT> RowPivot(imax);
244 
245  C_INT Info;
246 
247  bool found;
248  C_FLOAT64 tmp;
249  C_FLOAT64 * pTmp1;
250  C_FLOAT64 * pTmp2;
251  size_t Delta;
252 
253  C_INT subRows = Rows;
254  C_INT subCols = Cols;
255 
256  size_t iTop = 0;
257  size_t jTop = 0;
258  size_t LastCols = Cols;
259 
260  while (true)
261  {
262  /*
263  * -- LAPACK routine (version 3.0) --
264  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
265  * Courant Institute, Argonne National Lab, and Rice University
266  * March 31, 1993
267  *
268  * Purpose
269  * =======
270  *
271  * DGETRF computes an LU factorization of a general M-by-N matrix A
272  * using partial pivoting with row interchanges.
273  *
274  * The factorization has the form
275  * A = P * L * U
276  * where P is a permutation matrix, L is lower triangular with unit
277  * diagonal elements (lower trapezoidal if m > n), and U is upper
278  * triangular (upper trapezoidal if m < n).
279  *
280  * This is the right-looking Level 3 BLAS version of the algorithm.
281  *
282  * Arguments
283  * =========
284  *
285  * M (input) INTEGER
286  * The number of rows of the matrix A. M >= 0.
287  *
288  * N (input) INTEGER
289  * The number of columns of the matrix A. N >= 0.
290  *
291  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
292  * On entry, the M-by-N matrix to be factored.
293  * On exit, the factors L and U from the factorization
294  * A = P*L*U; the unit diagonal elements of L are not stored.
295  *
296  * LDA (input) INTEGER
297  * The leading dimension of the array A. LDA >= max(1,M).
298  *
299  * IPIV (output) INTEGER array, dimension (min(M,N))
300  * The pivot indices; for 1 <= i <= min(M,N), row i of the
301  * matrix was interchanged with row IPIV(i).
302  *
303  * INFO (output) INTEGER
304  * = 0: successful exit
305  * < 0: if INFO = -i, the i-th argument had an illegal value
306  * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
307  * has been completed, but the factor U is exactly
308  * singular, and division by zero will occur if it is used
309  * to solve a system of equations.
310  *
311  */
312  dgetrf_(&subRows, &subCols, &T(lCol[iTop], iTop), &Rows, RowPivot.array(), &Info);
313 
314  if (Info < 0) return false;
315 
316 #ifdef DEBUG_MATRIX
317  DebugFile << "After dgetrf_\n";
318  DebugFile << CTransposeView< CMatrix< C_FLOAT64 > >(T);
319 #endif
320 
321  for (i = 0; i < RowPivot.size(); i++)
322  if (i != RowPivot[i] - 1)
323  {
324  std::swap(row[i + iTop], row[RowPivot[i] - 1 + iTop]);
325 
326  pTmp1 = &T(0, i + iTop);
327  pTmp2 = &T(0, RowPivot[i] - 1 + iTop);
328 
329  for (k = 0; k < iTop; k++, pTmp1 += Rows, pTmp2 += Rows)
330  {
331  tmp = *pTmp1;
332  *pTmp1 = *pTmp2;
333  *pTmp2 = tmp;
334  }
335  }
336 
337 #ifdef DEBUG_MATRIX
338  DebugFile << row << std::endl;
339  DebugFile << "After row swap\n";
340  DebugFile << CTransposeView< CMatrix< C_FLOAT64 > >(T);
341 #endif
342 
343  if (Info == 0) break;
344 
345  iTop += Info - 1;
346 
347  // Search first non zero diagonal
348  for (i = iTop + 1, found = false; i < imax; i++)
349  if (T(lCol[i], i) != 0.0)
350  {
351  jTop = i;
352  found = true;
353  break;
354  }
355 
356  if (found)
357  {
358  for (i = Rows - 1; jTop <= i; i--)
359  {
360  for (j = Cols - 1; lCol[jTop] <= j; j--)
361  {
362  tmp = 0.0;
363 
364  pTmp1 = &T(iTop, i);
365  pTmp2 = &T(j, iTop);
366 
367  for (k = iTop; k < i && k < j; k++, pTmp1 += Rows, pTmp2++)
368  tmp += *pTmp1 * *pTmp2;
369 
370  if (i <= j) tmp += *pTmp2;
371 
372  T(j, i) = tmp;
373  }
374  }
375 
376  subRows = Rows - iTop;
377  subCols = Cols - lCol[iTop];
378 
379  pTmp1 = &T(lCol[iTop], 0);
380  pTmp2 = &T(lCol[jTop], 0);
381 
382  for (k = 0; k < Rows; k++, pTmp1++, pTmp2++)
383  {
384  tmp = *pTmp1;
385  *pTmp1 = *pTmp2;
386  *pTmp2 = tmp;
387  }
388 
389  std::swap(col[iTop], col[jTop]);
390 
391 #ifdef DEBUG_MATRIX
392  DebugFile << "After partial inversion: " << iTop << " , " << jTop << std::endl;
393  DebugFile << CTransposeView< CMatrix< C_FLOAT64 > >(T);
394  DebugFile << "After column swap\n";
395  DebugFile << CTransposeView< CMatrix< C_FLOAT64 > >(T);
396 #endif
397  }
398  else
399  {
400  for (j = iTop + 1, found = false; j < Cols; j++)
401  {
402  for (i = iTop; i < j && i < Rows; i++)
403  if (T(lCol[j], i) != 0)
404  {
405  found = true;
406  break;
407  }
408 
409  if (found) break;
410  }
411 
412  if (!found) break;
413 
414  jTop = j;
415 
416  subRows = Rows - iTop;
417  subCols = Cols - lCol[jTop];
418 
419  Delta = jTop - iTop;
420 
421  // Record the needed column swaps
422  for (j = jTop; j < LastCols; j += Delta)
423  for (i = j - Delta; i < j && i + Delta < LastCols; i++)
424  {
425  std::swap(col[i], col[i + Delta]);
426  std::swap(lCol[i], lCol[i + Delta]);
427  }
428 
429  LastCols -= Delta;
430  }
431 
432  RowPivot.resize(std::min(subRows, subCols));
433 
434 #ifdef DEBUG_MATRIX
435  DebugFile << lCol << std::endl;
436  DebugFile << col << std::endl;
437 #endif
438  }
439 
440  for (i = 0; i < Cols; i++)
441  for (j = 0; j < Rows; j++)
442  A(j, i) = T(lCol[i], j);
443 
444 #ifdef DEBUG_MATRIX
445  DebugFile << A;
446 #endif
447 
448  return true;
449 }
450 #endif // XXXX
int dgetrf_(integer *m, integer *n, doublereal *a, integer *lda, integer *ipiv, integer *info)
#define C_INT
Definition: copasi.h:115
virtual size_t numRows() const
Definition: CMatrix.h:138
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
#define C_INVALID_INDEX
Definition: copasi.h:222
virtual bool progressItem(const size_t &handle)
size_t addItem(const std::string &name, const std::string &value, const std::string *pEndValue=NULL)
bool LUfactor(CMatrix< C_FLOAT64 > &A, CVector< size_t > &row, CVector< size_t > &col, CProcessReport *cb)
Definition: CluX.cpp:30
virtual bool finishItem(const size_t &handle)
size_t size() const
Definition: CVector.h:100
#define C_FLOAT64
Definition: copasi.h:92
virtual size_t numCols() const
Definition: CMatrix.h:144
#define min(a, b)
Definition: f2c.h:175