COPASI API  4.16.103
Functions
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().

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 }
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)
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