COPASI API  4.16.103
CLinkMatrix.h
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 #ifndef COPASI_CLinkMatrix
7 #define COPASI_CLinkMatrix
8 
11 
12 template <class CType> class CCopasiVector;
13 class CCopasiObject;
14 
15 class CLinkMatrix: public CMatrix< C_FLOAT64 >
16 {
17 public:
18  /**
19  * Default constructor
20  */
21  CLinkMatrix();
22 
23  /**
24  * Copy constructor
25  */
26  CLinkMatrix(const CLinkMatrix & src);
27 
28  /**
29  * Destructor
30  */
31  virtual ~CLinkMatrix();
32 
33  /**
34  * Build the link matrix for the given matrix
35  * @param const CMatrix< C_FLOAT64 > & matrix
36  * @param size_t maxRank (default: C_INVALID_INDEX)
37  * @return bool success
38  */
39  bool build(const CMatrix< C_FLOAT64 > & matrix, size_t maxRank = C_INVALID_INDEX);
40 
41  /**
42  * Clear the pivot and swap vectors
43  */
44  void clearPivoting();
45 
46  /**
47  * Right multiply the given matrix M with L, i.e., P = alpha M * L.
48  * Note the columns of M must be in the same order as L.
49  * @param const C_FLOAT64 & alpha
50  * @param const CMatrix< C_FLOAT64> & M
51  * @param CMatrix< C_FLOAT64> & P
52  * @result bool success;
53  */
54  bool rightMultiply(const C_FLOAT64 & alpha, const CMatrix< C_FLOAT64> & M, CMatrix< C_FLOAT64> & P) const;
55 
56  /**
57  * Left multiply the given matrix M with L, i.e., P = L * M
58  * @param const CMatrix< C_FLOAT64> & M
59  * @param CMatrix< C_FLOAT64> & P
60  * @result bool success;
61  */
62  bool leftMultiply(const CMatrix< C_FLOAT64> & M, CMatrix< C_FLOAT64> & P) const;
63 
64  /**
65  * Retrieve the pivot vector used to create the link matrix
66  * @return const CVector< size_t > & rowPivots
67  */
68  const CVector< size_t > & getRowPivots() const;
69 
70  /**
71  * Retrieve the number of linear independent rows of the input matrix
72  * @return const size_t & numIndependent
73  */
74  const size_t & getNumIndependent() const;
75 
76  /**
77  * Retrieve the number of linear dependent rows of the input matrix
78  * @return size_t numDependent
79  */
80  size_t getNumDependent() const;
81 
82  /**
83  * Do the row pivot
84  * @param CMatrix< C_FLOAT64 > & matrix
85  * @return bool success
86  */
87  bool doRowPivot(CMatrix< C_FLOAT64 > & matrix) const;
88 
89  /**
90  * Undo the row pivot
91  * @param CMatrix< C_FLOAT64 > & matrix
92  * @return bool success
93  */
94  bool undoRowPivot(CMatrix< C_FLOAT64 > & matrix) const;
95 
96  /**
97  * Do the column pivot
98  * @param CMatrix< C_FLOAT64 > & matrix
99  * @return bool success
100  */
101  bool doColumnPivot(CMatrix< C_FLOAT64 > & matrix) const;
102 
103  /**
104  * Undo the column pivot
105  * @param CMatrix< C_FLOAT64 > & matrix
106  * @return bool success
107  */
108  bool undoColumnPivot(CMatrix< C_FLOAT64 > & matrix) const;
109 
110  /**
111  * Apply the row pivot
112  * @param CVectorCore< class CType > & vector
113  * @return bool success
114  */
115  template <class CType> bool applyRowPivot(CVectorCore< CType > & vector) const
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  }
158 
159 private:
160  /**
161  * Internal method performing apply and undo of column pivoting.
162  * @param CMatrix< C_FLOAT64 > & matrix
163  * @param const C_INT & incr
164  * @return bool success
165  */
166  bool applyColumnPivot(CMatrix< C_FLOAT64 > & matrix, const C_INT & incr) const;
167 
168  /**
169  * Internal method performing apply and undo of row pivoting.
170  * @param CMatrix< C_FLOAT64 > & matrix
171  * @param const CVector< size_t > & pivots
172  * @return bool success
173  */
174  bool applyRowPivot(CMatrix< C_FLOAT64 > & matrix,
175  const CVector< size_t > & pivots) const;
176 
177  /**
178  * Complete the pivot information.
179  */
181 
182  /**
183  * The row pivoting performed to create the link matrix
184  */
186 
187  /**
188  * The pivot vector used for undoing row swapping
189  */
191 
192  /**
193  * The swap vector used for column swapping
194  */
196 
197  /**
198  * The number of linear independent rows.
199  */
200  size_t mIndependent;
201 };
202 
204 {
205 public:
207 
208 private:
209  const CLinkMatrix * mpA;
210  const size_t * mpNumIndependent;
211  static const C_FLOAT64 mZero;
212  static const C_FLOAT64 mUnit;
213 
214 public:
215  /**
216  * Default constructor
217  * @param const const CLinkMatrix & A
218  * @param const size_t & mNumIndependent
219  */
220  CLinkMatrixView(const CLinkMatrix & A);
221 
222  /**
223  * Destructor.
224  */
226 
227  /**
228  * Assignment operator
229  * @param const CLinkMatrixView & rhs
230  * @return CLinkMatrixView & lhs
231  */
233 
234  /**
235  * The number of rows of the matrix.
236  * @return size_t rows
237  */
238  size_t numRows() const;
239 
240  /**
241  * The number of columns of the matrix
242  * @return size_t cols
243  */
244  size_t numCols() const;
245 
246  /**
247  * Retrieve a matrix element using the c-style indexing.
248  * @param const size_t & row
249  * @param const size_t & col
250  * @return elementType & element
251  */
252  inline C_FLOAT64 & operator()(const size_t & row,
253  const size_t & col) const
254  {
255  if (row >= *mpNumIndependent)
256  return *const_cast< C_FLOAT64 * >(&(*mpA)(row - *mpNumIndependent, col));
257  else if (row != col)
258  return *const_cast< C_FLOAT64 * >(&mZero);
259  else
260  return *const_cast< C_FLOAT64 * >(&mUnit);
261  }
262 
263  /**
264  * Output stream operator
265  * @param ostream & os
266  * @param const CLinkMatrixView & A
267  * @return ostream & os
268  */
269  friend std::ostream &operator<<(std::ostream &os,
270  const CLinkMatrixView & A);
271 };
272 
273 #endif // COPASI_CLinkMatrix
#define C_INT
Definition: copasi.h:115
C_FLOAT64 & operator()(const size_t &row, const size_t &col) const
Definition: CLinkMatrix.h:252
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
virtual ~CLinkMatrix()
Definition: CLinkMatrix.cpp:36
friend std::ostream & operator<<(std::ostream &os, const CLinkMatrixView &A)
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
#define C_INVALID_INDEX
Definition: copasi.h:222
void clearPivoting()
const size_t & getNumIndependent() const
size_t numCols() const
C_FLOAT64 elementType
Definition: CLinkMatrix.h:206
bool rightMultiply(const C_FLOAT64 &alpha, const CMatrix< C_FLOAT64 > &M, CMatrix< C_FLOAT64 > &P) const
CVector< C_INT > mSwapVector
Definition: CLinkMatrix.h:195
CLinkMatrixView(const CLinkMatrix &A)
static const C_FLOAT64 mUnit
Definition: CLinkMatrix.h:212
const CVector< size_t > & getRowPivots() const
Definition: CLinkMatrix.cpp:39
bool doRowPivot(CMatrix< C_FLOAT64 > &matrix) const
const CLinkMatrix * mpA
Definition: CLinkMatrix.h:209
bool applyRowPivot(CVectorCore< CType > &vector) const
Definition: CLinkMatrix.h:115
size_t numRows() const
bool doColumnPivot(CMatrix< C_FLOAT64 > &matrix) const
size_t size() const
Definition: CVector.h:100
#define C_FLOAT64
Definition: copasi.h:92
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
CVector< size_t > mPivotInverse
Definition: CLinkMatrix.h:190
CLinkMatrixView & operator=(const CLinkMatrixView &rhs)
bool undoRowPivot(CMatrix< C_FLOAT64 > &matrix) const