COPASI API  4.16.103
CEigen.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2014 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.cpp
17  *
18  * Programmer: Yongqun He
19  * Contact email: yohe@vt.edu
20  * Purpose: This is the .cpp file for the class CEigen.
21  * It is to calculate eigenvalues and eigenvectors of a matrix.
22  *
23  */
24 
25 #ifdef SunOS
26 # include <ieeefp.h>
27 #endif
28 
29 #ifdef WIN32
30 # define _USE_MATH_DEFINES
31 #endif
32 
33 #include <cmath>
34 #include <complex>
35 
36 #include "copasi.h"
37 
38 #include "CEigen.h"
39 
41 #include "utilities/CReadConfig.h"
43 #include "utilities/CSort.h"
44 
45 #include "lapack/lapackwrap.h" //use CLAPACK
46 
47 /**
48  * Default constructor
49  */
50 CEigen::CEigen(const std::string & name,
51  const CCopasiContainer * pParent):
52  CCopasiContainer(name, pParent, "Eigen Values",
53  CCopasiObject::Container),
54  mMaxrealpart(0),
55  mMaximagpart(0),
56  mNposreal(0),
57  mNnegreal(0),
58  mNreal(0),
59  mNimag(0),
60  mNcplxconj(0),
61  mNzero(0),
62  mStiffness(0),
63  mHierarchy(0),
64 
65  mMaxRealOfComplex(0.0),
66  mImagOfMaxComplex(0.0),
67  mFreqOfMaxComplex(0.0),
68  mOscillationIndicator(0.0),
69  mOscillationIndicator_EV(0.0),
70  mBifurcationIndicator_Hopf(0.0),
71  mBifurcationIndicator_Fold(0.0),
72  mBifurcationIndicator_Hopf_BDT(0.0),
73  mBifurcationIndicator_Fold_BDT(0.0),
74 
75  mResolution(0),
76  mJobvs('N'),
77  mSort('N'),
78  mpSelect(NULL),
79  mN(0),
80  mA(),
81  mLDA(0),
82  mSdim(0),
83  mR(),
84  mI(),
85  mpVS(NULL),
86  mLdvs(1),
87  mWork(1),
88  mLWork(4096),
89  mpBWork(NULL),
90  mInfo(0)
91 {
93  initObjects();
94 }
95 
96 CEigen::CEigen(const CEigen & src,
97  const CCopasiContainer * pParent):
98  CCopasiContainer(src, pParent),
99  mMaxrealpart(src.mMaxrealpart),
100  mMaximagpart(src.mMaximagpart),
101  mNposreal(src.mNposreal),
102  mNnegreal(src.mNnegreal),
103  mNreal(src.mNreal),
104  mNimag(src.mNimag),
105  mNcplxconj(src.mNcplxconj),
106  mNzero(src.mNzero),
107  mStiffness(src.mStiffness),
108  mHierarchy(src.mHierarchy),
109 
110  mMaxRealOfComplex(src.mMaxRealOfComplex),
111  mImagOfMaxComplex(src.mImagOfMaxComplex),
112  mFreqOfMaxComplex(src.mFreqOfMaxComplex),
113  mOscillationIndicator(src.mOscillationIndicator),
114  mOscillationIndicator_EV(src.mOscillationIndicator_EV),
115  mBifurcationIndicator_Hopf(src.mBifurcationIndicator_Hopf),
116  mBifurcationIndicator_Fold(src.mBifurcationIndicator_Fold),
117  mBifurcationIndicator_Hopf_BDT(src.mBifurcationIndicator_Hopf_BDT),
118  mBifurcationIndicator_Fold_BDT(src.mBifurcationIndicator_Fold_BDT),
119 
120  mResolution(src.mResolution),
121  mJobvs(src.mJobvs),
122  mSort(src.mSort),
123  mpSelect(NULL),
124  mN(src.mN),
125  mA(src.mA),
126  mLDA(src.mLDA),
127  mSdim(src.mSdim),
128  mR(src.mR),
129  mI(src.mI),
130  mpVS(NULL),
131  mLdvs(src.mLdvs),
132  mWork(src.mWork),
133  mLWork(src.mLWork),
134  mpBWork(NULL),
135  mInfo(src.mInfo)
136 {
138  initObjects();
139 }
140 
141 /**
142  * Deconstructor
143  */
145 {
146  cleanup();
148 }
149 
150 void CEigen::print(std::ostream * ostream) const {(*ostream) << (*this);}
151 
153 {
155  addObjectReference("Maximum imaginary part", mMaximagpart, CCopasiObject::ValueDbl);
156  addObjectReference("# Positive eigenvalues", mNposreal, CCopasiObject::ValueInt);
157  addObjectReference("# Negative eigenvalues", mNnegreal, CCopasiObject::ValueInt);
158  addObjectReference("# Real eigenvalues", mNreal, CCopasiObject::ValueInt);
159  addObjectReference("# Imaginary eigenvalues", mNimag, CCopasiObject::ValueInt);
160  addObjectReference("# Complex conjugated eigenvalues", mNcplxconj, CCopasiObject::ValueInt);
161  addObjectReference("# Zero eigenvalues", mNzero, CCopasiObject::ValueInt);
165  addVectorReference("Vector of real part of eigenvalues", mR, CCopasiObject::ValueDbl);
166  addVectorReference("Vector of imaginary part of eigenvalues", mI, CCopasiObject::ValueDbl);
167 
168  addObjectReference("Maximum real part of complex eigenvalue", mMaxRealOfComplex, CCopasiObject::ValueDbl);
169  addObjectReference("Imaginary part of largest complex eigenvalue", mImagOfMaxComplex, CCopasiObject::ValueDbl);
170  addObjectReference("Linear Frequency of largest complex eigenvalue", mFreqOfMaxComplex, CCopasiObject::ValueDbl);
172  addObjectReference("EV-based oscillation indicator", mOscillationIndicator_EV, CCopasiObject::ValueDbl);
173  addObjectReference("Hopf bifurcation test function", mBifurcationIndicator_Hopf, CCopasiObject::ValueDbl);
174  addObjectReference("Fold bifurcation test function", mBifurcationIndicator_Fold, CCopasiObject::ValueDbl);
175  addObjectReference("Hopf bifurcation test function (BDT)", mBifurcationIndicator_Hopf_BDT, CCopasiObject::ValueDbl);
176  addObjectReference("Fold bifurcation test function (BDT)", mBifurcationIndicator_Fold_BDT, CCopasiObject::ValueDbl);
177 }
178 
179 /**
180  * return the matrix
181  */
182 //TNT::Matrix < C_FLOAT64 > CEigen::getMatrix()
183 //{
184 // return mMatrix;
185 //}
186 
187 /**
188  * Set the Matrix
189  */
190 //void CEigen::setMatrix(C_INT32 rows, C_INT32 cols)
191 //{
192 // mMatrix.newsize(rows, cols);
193 //}
194 
195 //Get the max eigenvalue real part
197 {
198  return mMaxrealpart;
199 }
200 
201 //Get the max eigenvalue imaginary part
203 {
204  return mMaximagpart;
205 }
206 
207 // Get the number of zero eigenvalues
208 const size_t & CEigen::getNzero() const
209 {
210  return mNzero;
211 }
212 
213 //Get the eigenvalue stiffness
215 {
216  return mStiffness;
217 }
218 
219 //Get the eigenvalue hierarchy
221 {
222  return mHierarchy;
223 }
224 
225 //initialize variables for eigenvalue calculations
226 //
228 {
229  cleanup();
230 
232  mNzero = mNcplxconj = 0;
233 
234  mLDA = mN > 1 ? mN : 1;
235 
236  mR.resize(mN);
237  mI.resize(mN);
238 }
239 
241 {}
242 
244 {
245  assert(matrix.numRows() == matrix.numCols());
246  mN = (C_INT) matrix.numRows();
247  initialize();
248 
249  if (!mN) return;
250 
251  // copy the jacobian into mA
252  mA.resize(matrix.numRows(), matrix.numCols());
253  C_FLOAT64 * pA = mA.array();
254  C_FLOAT64 * pAEnd = pA + mA.size();
255  const C_FLOAT64 * pMatrix = matrix.array();
256 
257  for (; pA != pAEnd; ++pA, ++pMatrix)
258  {
259  *pA = *pMatrix;
260 
261  if (!finite(*pA) && !isnan(*pA))
262  {
263  if (*pA > 0)
265  else
267  }
268  }
269 
270  // Querry for the work array size.
271  mLWork = -1;
272  dgees_(&mJobvs, // 'N'
273  &mSort, // 'N'
274  NULL, // NULL,
275  &mN, // n,
276  mA.array(),
277  & mLDA,
278  & mSdim, // output
279  mR.array(),
280  mI.array(),
281  mpVS,
282  & mLdvs,
283  mWork.array(),
284  & mLWork,
285  mpBWork, // NULL
286  &mInfo); // output
287 
288  if (mInfo != 0)
289  {
290  // Exception
292  }
293 
294  mLWork = (C_INT) mWork[0];
296 
297  // calculate the eigenvalues
298  /* int dgees_(char *jobvs,
299  * char *sort,
300  * L_fp select,
301  * integer *n,
302  * doublereal *a,
303  * integer *lda,
304  * integer *sdim,
305  * doublereal *wr,
306  * doublereal *wi,
307  * doublereal *vs,
308  * integer *ldvs,
309  * doublereal *work,
310  * integer *lwork,
311  * logical *bwork,
312  * integer *info);
313  * Arguments
314  * =========
315  *
316  * JOBVS (input) CHARACTER*1
317  * = 'N': Schur vectors are not computed;
318  * = 'V': Schur vectors are computed.
319  *
320  * SORT (input) CHARACTER*1
321  * Specifies whether or not to order the eigenvalues on the
322  * diagonal of the Schur form.
323  * = 'N': Eigenvalues are not ordered;
324  * = 'S': Eigenvalues are ordered (see SELECT).
325  *
326  * SELECT (input) LOGICAL FUNCTION of two DOUBLE PRECISION arguments
327  * SELECT must be declared EXTERNAL in the calling subroutine.
328  * If SORT = 'S', SELECT is used to select eigenvalues to sort
329  * to the top left of the Schur form.
330  * If SORT = 'N', SELECT is not referenced.
331  * An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
332  * SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
333  * conjugate pair of eigenvalues is selected, then both complex
334  * eigenvalues are selected.
335  * Note that a selected complex eigenvalue may no longer
336  * satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
337  * ordering may change the value of complex eigenvalues
338  * (especially if the eigenvalue is ill-conditioned); in this
339  * case INFO is set to N+2 (see INFO below).
340  *
341  * N (input) INTEGER
342  * The order of the matrix A. N >= 0.
343  *
344  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
345  * On entry, the N-by-N matrix A.
346  * On exit, A has been overwritten by its real Schur form T.
347  *
348  * LDA (input) INTEGER
349  * The leading dimension of the array A. LDA >= max(1,N).
350  *
351  * SDIM (output) INTEGER
352  * If SORT = 'N', SDIM = 0.
353  * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
354  * for which SELECT is true. (Complex conjugate
355  * pairs for which SELECT is true for either
356  * eigenvalue count as 2.)
357  *
358  * WR (output) DOUBLE PRECISION array, dimension (N)
359  * WI (output) DOUBLE PRECISION array, dimension (N)
360  * WR and WI contain the real and imaginary parts,
361  * respectively, of the computed eigenvalues in the same order
362  * that they appear on the diagonal of the output Schur form T.
363  * Complex conjugate pairs of eigenvalues will appear
364  * consecutively with the eigenvalue having the positive
365  * imaginary part first.
366  *
367  * VS (output) DOUBLE PRECISION array, dimension (LDVS,N)
368  * If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
369  * vectors.
370  * If JOBVS = 'N', VS is not referenced.
371  *
372  * LDVS (input) INTEGER
373  * The leading dimension of the array VS. LDVS >= 1; if
374  * JOBVS = 'V', LDVS >= N.
375  *
376  * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
377  * On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
378  *
379  * LWORK (input) INTEGER
380  * The dimension of the array WORK. LWORK >= max(1,3*N).
381  * For good performance, LWORK must generally be larger.
382  *
383  * If LWORK = -1, then a workspace query is assumed; the routine
384  * only calculates the optimal size of the WORK array, returns
385  * this value as the first entry of the WORK array, and no error
386  * message related to LWORK is issued by XERBLA.
387  *
388  * BWORK (workspace) LOGICAL array, dimension (N)
389  * Not referenced if SORT = 'N'.
390  *
391  * INFO (output) INTEGER
392  * = 0: successful exit
393  * < 0: if INFO = -i, the i-th argument had an illegal value.
394  * > 0: if INFO = i, and i is
395  * <= N: the QR algorithm failed to compute all the
396  * eigenvalues; elements 1:i-1 and i+1:N of WR and WI
397  * contain those eigenvalues which have converged; if
398  * JOBVS = 'V', VS contains the matrix which reduces A
399  * to its partially converged Schur form.
400  * = N+1: the eigenvalues could not be reordered because some
401  * eigenvalues were too close to separate (the problem
402  * is very ill-conditioned);
403  * = N+2: after reordering, roundoff changed values of some
404  * complex eigenvalues so that leading eigenvalues in
405  * the Schur form no longer satisfy SELECT=.TRUE. This
406  * could also be caused by underflow due to scaling.
407  *
408  */
409  dgees_(&mJobvs, // 'N'
410  &mSort, // 'N'
411  NULL, // NULL,
412  &mN, // n,
413  mA.array(),
414  & mLDA,
415  & mSdim, // output
416  mR.array(),
417  mI.array(),
418  mpVS,
419  & mLdvs,
420  mWork.array(),
421  & mLWork,
422  mpBWork, // NULL
423  &mInfo); // output
424 
425  if (mInfo != 0)
426  {
427  if (mInfo < 0)
428  {
429  // Exception
431  }
432  else if (mInfo <= mN)
433  {
434  // Warning
436  }
437  else if (mInfo == mN + 1)
438  {
439  // Warning
441  }
442  else if (mInfo == mN + 2)
443  {
444  // Warning
446  }
447  else // Catch all, should never happen.
448  fatalError();
449  }
450 
451  return;
452 }
453 
454 void CEigen::stabilityAnalysis(const C_FLOAT64 & resolution)
455 {
456  if (!mN) return;
457 
458  size_t mx, mn; // YH: n is the 4th parameter, not here
459  size_t i;
460  C_FLOAT64 distt, maxt, tott;
461  mResolution = resolution;
462 
463  // sort the eigenvalues
464  CVector<size_t> Pivot;
465 
467 
468  // The sort order is ascending however we need descending
469  size_t *pTo = Pivot.array();
470  size_t *pFrom = pTo + mN - 1;
471 
472  for (; pTo < pFrom; ++pTo, --pFrom)
473  {
474  size_t Tmp = *pFrom;
475  *pFrom = *pTo;
476  *pTo = Tmp;
477  }
478 
479  mR.applyPivot(Pivot);
480  mI.applyPivot(Pivot);
481 
482  // calculate various eigenvalue statistics
483  mMaxrealpart = mR[0];
484  mMaxRealOfComplex = mR[mN - 1] - 1;
485  mMaximagpart = fabs(mI[0]);
486  mImagOfMaxComplex = 0.0;
487 
488  for (i = 0; (C_INT) i < mN; i++)
489  {
490  // for the largest real part
491  if (mR[i] > mMaxrealpart)
492  mMaxrealpart = mR[i];
493 
494  // for the largest imaginary part
495  if (fabs(mI[i]) > mMaximagpart)
496  mMaximagpart = fabs(mI[i]);
497 
498  // for the largest complex eigenvalue
499  if ((fabs(mI[i]) > resolution) && (mR[i] > mMaxRealOfComplex))
500  {
501  mMaxRealOfComplex = mR[i];
502  mImagOfMaxComplex = fabs(mI[i]);
503  }
504 
505  if (fabs(mR[i]) > resolution)
506  {
507  // positive real part
508  if (mR[i] >= resolution)
509  mNposreal++;
510 
511  // negative real part
512  if (mR[i] <= -resolution)
513  mNnegreal++;
514 
515  if (fabs(mI[i]) > resolution)
516  {
517  // complex
518  mNcplxconj++;
519  }
520  else
521  {
522  mI[i] = 0.0;
523  // pure real
524  mNreal++;
525  }
526  }
527  else
528  {
529  mR[i] = 0.0;
530 
531  if (fabs(mI[i]) > resolution)
532  {
533  // pure imaginary
534  mNimag++;
535  }
536  else
537  {
538  mI[i] = 0.0;
539  // zero
540  mNzero++;
541  }
542  }
543  }
544 
545  if (mImagOfMaxComplex == 0)
546  mMaxRealOfComplex = mR[mN - 1]; //default value
547 
548  mFreqOfMaxComplex = mImagOfMaxComplex / (2 * M_PI);
549 
550  if (mNposreal > 0)
551  {
552  if (mR[0] > fabs(mR[mN - 1]))
553  mx = 0;
554  else
555  mx = mN - 1;
556 
557  if ((C_INT) mNposreal == mN)
558  mn = mNposreal - 1;
559  else if (mR[mNposreal - 1] < fabs(mR[mNposreal]))
560  mn = mNposreal - 1;
561  else
562  mn = mNposreal;
563  }
564  else
565  {
566  mx = mN - 1; // index of the largest absolute real part
567  mn = 0; // index of the smallest absolute real part
568  }
569 
570  mStiffness = fabs(mR[mx]) / fabs(mR[mn]);
571 
572  maxt = tott = fabs(1 / mR[mn]);
573  distt = 0.0;
574 
575  for (i = 1; (C_INT) i < mN; i++)
576  if (i != mn)
577  {
578  distt += maxt - fabs(1 / mR[i]);
579  tott += fabs(1 / mR[i]);
580  }
581 
582  mHierarchy = distt / tott / (mN - 1);
583 
584  //we calculate an oscillation indicator based on the eigenvalues. It is using a rather heuristical approach
585  if (mN < 2) // at least 2 Eigenvalues are required for oscillations
587  else
588  {
589  if (fabs(mI[0]) > resolution)
590  mOscillationIndicator_EV = mR[0] * (mR[0] > 0 ? fabs(mI[0]) / (0.01 + fabs(mI[0])) : 2 - fabs(mI[0]) / (0.01 + fabs(mI[0])));
591  else
592  mOscillationIndicator_EV = 2 * mR[1] - (mR[0] > 0 ? 2 * mR[0] : 0);
593  }
594 
595  //bifurcation test functions, according to Kuznetsov "Elements of applied bifurcation theory"
596  // this function can also be calculated without using the eigenvalues, but since we already
597  //have them we can use them
598 
599  std::complex<C_FLOAT64> tmpcpl = 1.0;
600  unsigned C_INT32 index_min = 0; // the index of the EV with smallest abs real part
601  C_FLOAT64 tmpmin = 1e300;
602 
603  for (i = 0; (C_INT) i < mN; i++)
604  {
605  tmpcpl *= std::complex<C_FLOAT64>(mR[i], mI[i]);
606 
607  if (fabs(mR[i]) < tmpmin)
608  {
609  tmpmin = fabs(mR[i]);
610  index_min = i;
611  }
612  }
613 
614  mBifurcationIndicator_Fold = std::abs<C_FLOAT64>(tmpcpl);
615 
616  tmpcpl = 1.0;
617 
618  for (i = 0; (C_INT) i < mN - 1; i++)
619  {
620  tmpcpl *= (std::complex<C_FLOAT64>(mR[i], mI[i]) + std::complex<C_FLOAT64>(mR[i + 1], mI[i + 1]));
621  }
622 
623  mBifurcationIndicator_Hopf = (tmpcpl.real());
624 
625  //now the test functions from the "Bifurcation discovery tool" (as described in the publication by Chickarmane et al.)
626 
627  //calculate the product of EVs excluding the minimal one
628  tmpcpl = 1.0;
629 
630  for (i = 0; (C_INT) i < mN; i++)
631  {
632  if (i != index_min)
633  tmpcpl *= std::complex<C_FLOAT64>(mR[i], mI[i]);
634  }
635 
636  mBifurcationIndicator_Fold_BDT = mBifurcationIndicator_Fold / (1 - 0.99 * exp(-std::abs<C_FLOAT64>(tmpcpl)));
637 
638  C_FLOAT64 tmp_product = 1.0;
639 
640  for (i = 0; (C_INT) i < mN; i++)
641  {
642  tmp_product *= mR[i] / (1 - 0.99 * exp(-mI[i]));
643  }
644 
645  mBifurcationIndicator_Hopf_BDT = tmp_product;
646 
647  mOscillationIndicator = 0.0; //not used at the moment.
648 }
649 
650 /**
651  * Return number of real eigenvalues WeiSun 3/28/02
652  */
653 const size_t & CEigen::getNreal() const
654 {
655  return mNreal;
656 }
657 
658 /**
659  * Return the number of imaginary eigenvalue numbers
660  */
661 const size_t & CEigen::getNimag() const
662 {
663  return mNimag;
664 }
665 
666 const size_t & CEigen::getNcplxconj() const
667 {
668  return mNcplxconj;
669 }
670 
671 /**
672  * Return the number of eigenvalues with positive real part
673  */
674 const size_t & CEigen::getNposreal() const
675 {
676  return mNposreal;
677 }
678 
679 /**
680  * Return the number of eigenvalues with negative real part
681  */
682 const size_t & CEigen::getNnegreal() const
683 {
684  return mNnegreal;
685 }
686 
688 {return mI;}
689 
691 {return mR;}
692 
693 std::ostream &operator<<(std::ostream &os, const CEigen &A)
694 {
695  os << std::endl;
696  os << "KINETIC STABILITY ANALYSIS";
697  os << std::endl;
698  os << "The linear stability analysis based on the eigenvalues" << std::endl;
699  os << "of the Jacobian matrix is only valid for steady states." << std::endl;
700  os << std::endl;
701  os << "Summary:" << std::endl;
702  os << "This state ";
703 
704  // Output statistics
705 
706  if (A.mMaxrealpart > A.mResolution)
707  os << "is unstable";
708  else if (A.mMaxrealpart < -A.mResolution)
709  os << "is asymptotically stable";
710  else
711  os << "'s stability is undetermined";
712 
713  if (A.mMaximagpart > A.mResolution)
714  {
715  os << "," << std::endl;
716  os << "transient states in its vicinity have oscillatory components";
717  }
718 
719  os << "." << std::endl;
720  os << std::endl;
721 
722  os << "Eigenvalue statistics:" << std::endl;
723  // Output Max Real Part
724  os << " Largest real part: ";
725  os << std::setprecision(6) << A.mMaxrealpart << std::endl;
726  // Output Max imaginary Part
727  os << " Largest absolute imaginary part: ";
728  os << std::setprecision(6) << A.mMaximagpart << std::endl;
729 
730  if (A.mImagOfMaxComplex > A.mResolution)
731  os << " The complex eigenvalues with the largest real part are: "
732  << A.mMaxRealOfComplex << " +|- " << A.mImagOfMaxComplex << "i" << std::endl;
733 
734  // Output Eigen-nreal
735  os.unsetf(std::ios_base::scientific);
736  os.unsetf(std::ios_base::showpoint);
737  os << " " << A.mNreal;
738  os << " are purely real" << std::endl;
739  // Output Eigen-nimage
740  os << " " << A.mNimag;
741  os << " are purely imaginary" << std::endl;
742  // Output Eigen-ncplxconj
743  os << " " << A.mNcplxconj;
744  os << " are complex" << std::endl;
745  // Output Eigen-nzero
746  os << " " << A.mNzero;
747  os << " are equal to zero" << std::endl;
748  // Output Eigen-nposreal
749  os << " " << A.mNposreal;
750  os << " have positive real part" << std::endl;
751  // Output Eigen-nnegreal
752  os << " " << A.mNnegreal;
753  os << " have negative real part" << std::endl;
754 
755  // Set point manipulators
756  os.setf(std::ios_base::showpoint);
757  // Output Eigne-stiffness
758  os << " stiffness = " << A.mStiffness << std::endl;
759  os << " time hierarchy = " << A.mHierarchy << std::endl;
760 
761  return os;
762 }
C_FLOAT64 mMaximagpart
Definition: CEigen.h:49
#define C_INT
Definition: copasi.h:115
C_FLOAT64 mBifurcationIndicator_Fold
Definition: CEigen.h:126
void sortWithPivot(RandomAccessIterator first, RandomAccessIterator last, CVector< size_t > &pivot)
Definition: CSort.h:77
CVector< C_FLOAT64 > mI
Definition: CEigen.h:203
C_FLOAT64 mMaxRealOfComplex
Definition: CEigen.h:95
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
virtual size_t numRows() const
Definition: CMatrix.h:138
#define fatalError()
C_FLOAT64 mOscillationIndicator_EV
Definition: CEigen.h:115
CVector< C_FLOAT64 > mWork
Definition: CEigen.h:222
CMatrix< C_FLOAT64 > mA
Definition: CEigen.h:180
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
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
void cleanup()
Definition: CEigen.cpp:240
#define DESTRUCTOR_TRACE
Definition: copasi.h:206
const size_t & getNimag() const
Definition: CEigen.cpp:661
C_FLOAT64 mResolution
Definition: CEigen.h:141
C_FLOAT64 mImagOfMaxComplex
Definition: CEigen.h:100
#define MCEigen
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
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
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
size_t size() const
Definition: CVector.h:100
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
CType * array()
Definition: CVector.h:139
const C_FLOAT64 & getStiffness() const
Definition: CEigen.cpp:214
virtual size_t size() const
Definition: CMatrix.h:132
int dgees_(char *jobvs, char *sort, L_fp select, integer *n, doublereal *a, integer *lda, integer *sdim, doublereal *wr, doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work, integer *lwork, logical *bwork, integer *info)
std::ostream & operator<<(std::ostream &os, const CEigen &A)
Definition: CEigen.cpp:693
CEigen(const std::string &name="NoName", const CCopasiContainer *pParent=NULL)
Definition: CEigen.cpp:50
CCopasiObject * addVectorReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
const size_t & getNnegreal() const
Definition: CEigen.cpp:682
size_t mNzero
Definition: CEigen.h:79
virtual size_t numCols() const
Definition: CMatrix.h:144
const size_t & getNzero() const
Definition: CEigen.cpp:208
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
const size_t & getNreal() const
Definition: CEigen.cpp:653
virtual CType * array()
Definition: CMatrix.h:337
size_t mNposreal
Definition: CEigen.h:54
char mSort
Definition: CEigen.h:157
#define CONSTRUCTOR_TRACE
Definition: copasi.h:202
bool applyPivot(const CVectorCore< size_t > &pivot)
Definition: CVector.h:172
#define max(a, b)
Definition: f2c.h:176