COPASI API  4.16.103
CEigen.h
Go to the documentation of this file.
1 // Copyright (C) 2010 - 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 // 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) 2001 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /**
16  * File name: CEigen.h
17  *
18  * Programmer: Yongqun He
19  * Contact email: yohe@vt.edu
20  * Purpose: This is the .h file for the class CEigen.
21  * It is to calculate eigenvalues and eigenvectors of a matrix.
22  * It mainly uses the dgees_() subroutine of CLAPACK
23  *
24  */
25 
26 #ifndef COPASI_CEigen
27 #define COPASI_CEigen
28 
29 #include <iostream>
30 #include <iomanip>
31 
32 #include "utilities/CMatrix.h"
33 #include "utilities/CVector.h"
35 
36 class CEigen: public CCopasiContainer
37 {
38 private:
39  // variables for stability analysis
40 
41  /**
42  * the real part of the maximum eigenvalue
43  */
45 
46  /**
47  * the imaginary part of the maximum eigenvalue
48  */
50 
51  /**
52  * the number of eigenvalues with positive real part
53  */
54  size_t mNposreal;
55 
56  /**
57  * the number of eigenvalues with negative real part
58  */
59  size_t mNnegreal;
60 
61  /**
62  * the number of real eigenvalues
63  */
64  size_t mNreal;
65 
66  /**
67  * the number of imaginary eigenvalue numbers
68  */
69  size_t mNimag;
70 
71  /**
72  *
73  */
74  size_t mNcplxconj;
75 
76  /**
77  * the number of eigenvalues with value of zero
78  */
79  size_t mNzero;
80 
81  /**
82  * the stiffness of eigenvalues
83  */
85 
86  /**
87  * the hierary of the eigenvalues
88  */
90 
91  /**
92  * largest real part of a complex eigenvalue.
93  * If there is no complex eigenvalue it will be the smallest ev.
94  */
96 
97  /**
98  * imaginary part of complex with largest real part
99  */
101 
102  /**
103  * frequency of complex with largest real part
104  */
106 
107  /**
108  * An index that is supposed to indicate the presence of oscillations
109  */
111 
112  /**
113  * A heuristical index that is supposed to indicate the presence of oscillations, based on Eigenvalues
114  */
116 
117  /**
118  * A bifurcation test function, it is 0 for a Hopf-bifurcation. Definition according to Kuznetsov.
119  * This is only a necessary, not a sufficient condition for a H.-bifurcation
120  */
122 
123  /**
124  * A bifurcation test function, it is 0 for a Fold-type bifurcation.
125  */
127 
128  /**
129  * A bifurcation test function, defined according to the "Bifurcation Discovery Tool", Chickarmane et al 2005
130  */
132 
133  /**
134  * A bifurcation test function, defined according to the "Bifurcation Discovery Tool", Chickarmane et al 2005
135  */
137 
138  /**
139  * The resolution of needed for the stability analysis
140  */
142 
143  //there are 15 parameters in the dgees() subroutine
144 
145  /**
146  * #1: (input) characer*1
147  * = 'N': Schur vectors are not computed
148  * = 'V': Schur vectors are computed
149  */
150  char mJobvs;
151 
152  /**
153  * #2: (input) characer*1
154  * = 'N': Eigenvalues are not ordered
155  * = 'S': Eigenvalues are ordered
156  */
157  char mSort;
158 
159  /**
160  * #3: (input) Logical function of two double precision arguments
161  * It must be declared external
162  * = 'S': Select is used to select eigenvalues to sort to the top left
163  * of the Schur form
164  * = 'N': Eigenvalues are ordered Select is not refereced.
165  */
167 
168  /**
169  * #4: (input) The order of the matrix A
170  */
172 
173  /**
174  * #5: (input/output) The double precision array, dimension (LDA,N)
175  * On entry, the N-by-N matrix A
176  * On exit, A has been overwritten by its real Schur form T
177  * Use a pointer variable here instead of array since we don't know
178  * the dimension yet.
179  */
181 
182  /**
183  * #6: (input) The leading dimension of the array A. LDA >= max(1,N)
184  */
186 
187  /**
188  * #7: (output) an integer
189  * if Sort = 'N', its value is 0
190  * if Sort = 'S', its value = number of eigenvalues (after sorting)
191  * for which mSelect is true.
192  */
194 
195  /**
196  * #8: array with dimension (mN)
197  */
199 
200  /**
201  * #9: array with dimension (mN)
202  */
204 
205  /**
206  * #10: (output) array with dimension (mLdvs, mN)
207  * If mJobvs='V', mVS contains the orthogoanl matrix Z of Schur vectors
208  * If mJobvs='N', mVS is not referenced
209  */
211 
212  /**
213  * #11: an integer, the leading dimension of the array VS. mLdvs >= 1;
214  * if mJobvs='V', mLdvs >= N.
215  */
217 
218  /**
219  * #12: (workspace/output) double precision array, dimension (mLWork)
220  * On exit, if mInfo=0; mWork(1) contains the optimal mLWork
221  */
223 
224  /**
225  * #13: (input) Dimension of array Work, its value >= max(1,3*mN).
226  * For good performance, it must generally be larger
227  */
229 
230  /**
231  * #14: (workspace) Logical array, dimension (N)
232  * Not referenced if mSort = 'N'
233  */
235 
236  /**
237  * #15: (output) an integer
238  * =0: successful exit
239  * <0: if mInfo=-i, the ith argument had an illegal value
240  * >0: if mInfo=i, and i is
241  * <=N: the QR algorithm failed to compute all the eigenvalues;
242  * =N+1: the eigenvalues could not be reordered because some
243  * eigenvalues were too close to separate (ill-conditioned)
244  * =N+2: after reordering, roundoff changed values of some
245  * complex eigenvalues so that leading eigenvalues in the
246  * Schur form no longer satisfy mSelect=.True. This could
247  * caused by underflow due to scaling
248  */
250 
251 public:
252  /**
253  * Default constructor
254  * @param const std::string & name (default: "NoName")
255  * @param const CCopasiContainer * pParent (default: NULL)
256  */
257  CEigen(const std::string & name = "NoName",
258  const CCopasiContainer * pParent = NULL);
259 
260  /**
261  * Copy constructor
262  * @param const CMetab & src
263  * @param const CCopasiContainer * pParent (default: NULL)
264  */
265  CEigen(const CEigen & src,
266  const CCopasiContainer * pParent = NULL);
267  /**
268  * Destructor
269  */
270  virtual ~CEigen();
271 
272  /**
273  * This is the output method for any object. The default implementation
274  * provided with CCopasiObject uses the ostream operator<< of the object
275  * to print the object.To overide this default behaviour one needs to
276  * reimplement the virtual print function.
277  * @param std::ostream * ostream
278  */
279  virtual void print(std::ostream * ostream) const;
280 
281  /**
282  * initialize variables for eigenvalue calculations
283  */
284  void initialize();
285 
286  /**
287  * cleanup()
288  */
289  void cleanup();
290 
291  /**
292  * Eigenvalue calculations
293  * @param const C_FLOAT64 * matrix
294  * @param const unsigned C_INT32 & dim
295  */
296  void calcEigenValues(const CMatrix< C_FLOAT64 > & matrix);
297 
298  /**
299  * Calculate various eigenvalue statistics
300  */
301  void stabilityAnalysis(const C_FLOAT64 & resolution);
302 
303  /**
304  * Get the max eigenvalue real part
305  */
306  const C_FLOAT64 & getMaxrealpart() const;
307 
308  /**
309  * Get the max eigenvalue imaginary part
310  */
311  const C_FLOAT64 & getMaximagpart() const;
312 
313  /**
314  * Get the number of zero eigenvalues
315  */
316  const size_t & getNzero() const;
317 
318  /**
319  * Get the eigenvalue stiffness
320  */
321  const C_FLOAT64 & getStiffness() const;
322 
323  /**
324  * Get the eigenvalue hierarchy
325  */
326  const C_FLOAT64 & getHierarchy() const;
327 
328  /**
329  * Return number of real eigenvalues WeiSun 3/28/02
330  */
331  const size_t & getNreal() const;
332 
333  /**
334  * Return the number of imaginary eigenvalue numbers
335  */
336  const size_t & getNimag() const;
337 
338  const size_t & getNcplxconj() const;
339 
340  /**
341  * Return the number of eigenvalues with positive real part
342  */
343  const size_t & getNposreal() const;
344 
345  /**
346  * Return the number of eigenvalues with negative real part
347  */
348  const size_t & getNnegreal() const;
349 
350  const CVector< C_FLOAT64 > & getI() const;
351 
352  const CVector< C_FLOAT64 > & getR() const;
353 
354  // Friend functions
355  friend std::ostream &operator<<(std::ostream &os, const CEigen &A);
356 
357 private:
358  /**
359  * Initialize the contained CCopasiObjects
360  */
361  void initObjects();
362 };
363 
364 #endif
C_FLOAT64 mMaximagpart
Definition: CEigen.h:49
#define C_INT
Definition: copasi.h:115
C_FLOAT64 mBifurcationIndicator_Fold
Definition: CEigen.h:126
CVector< C_FLOAT64 > mI
Definition: CEigen.h:203
C_FLOAT64 mMaxRealOfComplex
Definition: CEigen.h:95
C_INT32 * mpSelect
Definition: CEigen.h:166
C_INT mLDA
Definition: CEigen.h:185
C_FLOAT64 mMaxrealpart
Definition: CEigen.h:44
C_INT mInfo
Definition: CEigen.h:249
void calcEigenValues(const CMatrix< C_FLOAT64 > &matrix)
Definition: CEigen.cpp:243
C_INT mLWork
Definition: CEigen.h:228
C_FLOAT64 mBifurcationIndicator_Hopf_BDT
Definition: CEigen.h:131
C_FLOAT64 mFreqOfMaxComplex
Definition: CEigen.h:105
C_FLOAT64 * mpVS
Definition: CEigen.h:210
size_t mNcplxconj
Definition: CEigen.h:74
Definition: CEigen.h:36
virtual ~CEigen()
Definition: CEigen.cpp:144
C_INT mLdvs
Definition: CEigen.h:216
C_FLOAT64 mHierarchy
Definition: CEigen.h:89
C_FLOAT64 mOscillationIndicator_EV
Definition: CEigen.h:115
CVector< C_FLOAT64 > mWork
Definition: CEigen.h:222
CMatrix< C_FLOAT64 > mA
Definition: CEigen.h:180
C_FLOAT64 mBifurcationIndicator_Hopf
Definition: CEigen.h:121
#define C_INT32
Definition: copasi.h:90
CVector< C_FLOAT64 > mR
Definition: CEigen.h:198
char mJobvs
Definition: CEigen.h:150
const CVector< C_FLOAT64 > & getI() const
Definition: CEigen.cpp:687
size_t mNnegreal
Definition: CEigen.h:59
C_FLOAT64 mOscillationIndicator
Definition: CEigen.h:110
void initObjects()
Definition: CEigen.cpp:152
size_t mNreal
Definition: CEigen.h:64
C_FLOAT64 mBifurcationIndicator_Fold_BDT
Definition: CEigen.h:136
#define C_LOGICAL
Definition: copasi.h:122
void cleanup()
Definition: CEigen.cpp:240
const size_t & getNimag() const
Definition: CEigen.cpp:661
C_FLOAT64 mResolution
Definition: CEigen.h:141
C_FLOAT64 mImagOfMaxComplex
Definition: CEigen.h:100
void initialize()
Definition: CEigen.cpp:227
void stabilityAnalysis(const C_FLOAT64 &resolution)
Definition: CEigen.cpp:454
size_t mNimag
Definition: CEigen.h:69
const size_t & getNposreal() const
Definition: CEigen.cpp:674
C_FLOAT64 mStiffness
Definition: CEigen.h:84
const CVector< C_FLOAT64 > & getR() const
Definition: CEigen.cpp:690
Header file of class CCopasiContainer.
const C_FLOAT64 & getHierarchy() const
Definition: CEigen.cpp:220
virtual void print(std::ostream *ostream) const
Definition: CEigen.cpp:150
C_INT mSdim
Definition: CEigen.h:193
const C_FLOAT64 & getMaxrealpart() const
Definition: CEigen.cpp:196
C_INT mN
Definition: CEigen.h:171
const size_t & getNcplxconj() const
Definition: CEigen.cpp:666
C_LOGICAL * mpBWork
Definition: CEigen.h:234
#define C_FLOAT64
Definition: copasi.h:92
const C_FLOAT64 & getMaximagpart() const
Definition: CEigen.cpp:202
const C_FLOAT64 & getStiffness() const
Definition: CEigen.cpp:214
CEigen(const std::string &name="NoName", const CCopasiContainer *pParent=NULL)
Definition: CEigen.cpp:50
const size_t & getNnegreal() const
Definition: CEigen.cpp:682
size_t mNzero
Definition: CEigen.h:79
const size_t & getNzero() const
Definition: CEigen.cpp:208
const size_t & getNreal() const
Definition: CEigen.cpp:653
size_t mNposreal
Definition: CEigen.h:54
char mSort
Definition: CEigen.h:157
friend std::ostream & operator<<(std::ostream &os, const CEigen &A)
Definition: CEigen.cpp:693