#include <CEigen.h>

Inheritance diagram for CEigen:

Collaboration diagram for CEigen:

Public Member Functions
void	calcEigenValues (const CMatrix< C_FLOAT64 > &matrix)

	CEigen (const std::string &name="NoName", const CCopasiContainer *pParent=NULL)

	CEigen (const CEigen &src, const CCopasiContainer *pParent=NULL)

void	cleanup ()

const C_FLOAT64 &	getHierarchy () const

const CVector< C_FLOAT64 > &	getI () const

const C_FLOAT64 &	getMaximagpart () const

const C_FLOAT64 &	getMaxrealpart () const

const size_t &	getNcplxconj () const

const size_t &	getNimag () const

const size_t &	getNnegreal () const

const size_t &	getNposreal () const

const size_t &	getNreal () const

const size_t &	getNzero () const

const CVector< C_FLOAT64 > &	getR () const

const C_FLOAT64 &	getStiffness () const

void	initialize ()

virtual void	print (std::ostream *ostream) const

void	stabilityAnalysis (const C_FLOAT64 &resolution)

virtual	~CEigen ()

Public Member Functions inherited from CCopasiContainer
virtual bool	add (CCopasiObject *pObject, const bool &adopt=true)

	CCopasiContainer (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=CCopasiObject::Container)

	CCopasiContainer (const CCopasiContainer &src, const CCopasiContainer *pParent=NULL)

virtual std::string	getChildObjectUnits (const CCopasiObject *pObject) const

virtual const CObjectInterface *	getObject (const CCopasiObjectName &cn) const

virtual const objectMap &	getObjects () const

virtual std::string	getUnits () const

virtual const CCopasiObject *	getValueObject () const

virtual bool	remove (CCopasiObject *pObject)

virtual	~CCopasiContainer ()

Public Member Functions inherited from CCopasiObject
void	addDirectDependency (const CCopasiObject *pObject)

	CCopasiObject (const CCopasiObject &src, const CCopasiContainer *pParent=NULL)

void	clearDirectDependencies ()

void	clearRefresh ()

bool	dependsOn (DataObjectSet candidates, const DataObjectSet &context=DataObjectSet()) const

void	getAllDependencies (DataObjectSet &dependencies, const DataObjectSet &context) const

virtual CCopasiObjectName	getCN () const

virtual const DataObjectSet &	getDirectDependencies (const DataObjectSet &context=DataObjectSet()) const

virtual const std::string &	getKey () const

CCopasiContainer *	getObjectAncestor (const std::string &type) const

CCopasiDataModel *	getObjectDataModel ()

const CCopasiDataModel *	getObjectDataModel () const

virtual std::string	getObjectDisplayName (bool regular=true, bool richtext=false) const

const std::string &	getObjectName () const

CCopasiContainer *	getObjectParent () const

const std::string &	getObjectType () const

virtual const CObjectInterface::ObjectSet &	getPrerequisites () const

virtual Refresh *	getRefresh () const

UpdateMethod *	getUpdateMethod () const

virtual void *	getValuePointer () const

bool	hasCircularDependencies (DataObjectSet &candidates, DataObjectSet &verified, const DataObjectSet &context) const

bool	hasUpdateMethod () const

bool	isArray () const

bool	isContainer () const

bool	isDataModel () const

bool	isMatrix () const

bool	isNameVector () const

bool	isNonUniqueName () const

virtual bool	isPrerequisiteForContext (const CObjectInterface *pObject, const CMath::SimulationContextFlag &context, const CObjectInterface::ObjectSet &changedObjects) const

bool	isReference () const

bool	isRoot () const

bool	isSeparator () const

bool	isStaticString () const

bool	isValueBool () const

bool	isValueDbl () const

bool	isValueInt () const

bool	isValueInt64 () const

bool	isValueString () const

bool	isVector () const

virtual bool	mustBeDeleted (const DataObjectSet &deletedObjects) const

void	removeDirectDependency (const CCopasiObject *pObject)

void	setDirectDependencies (const DataObjectSet &directDependencies)

bool	setObjectName (const std::string &name)

virtual bool	setObjectParent (const CCopasiContainer *pParent)

void	setObjectValue (const C_FLOAT64 &value)

void	setObjectValue (const C_INT32 &value)

void	setObjectValue (const bool &value)

template<class CType >
void	setRefresh (CType pType, void(CType::method)(void))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const C_FLOAT64 &))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const C_INT32 &))

template<class CType >
void	setUpdateMethod (CType pType, void(CType::method)(const bool &))

virtual	~CCopasiObject ()

Public Member Functions inherited from CObjectInterface
	CObjectInterface ()

virtual	~CObjectInterface ()

Private Member Functions
void	initObjects ()

Private Attributes
CMatrix< C_FLOAT64 >	mA

C_FLOAT64	mBifurcationIndicator_Fold

C_FLOAT64	mBifurcationIndicator_Fold_BDT

C_FLOAT64	mBifurcationIndicator_Hopf

C_FLOAT64	mBifurcationIndicator_Hopf_BDT

C_FLOAT64	mFreqOfMaxComplex

C_FLOAT64	mHierarchy

CVector< C_FLOAT64 >	mI

C_FLOAT64	mImagOfMaxComplex

C_INT	mInfo

char	mJobvs

C_INT	mLDA

C_INT	mLdvs

C_INT	mLWork

C_FLOAT64	mMaximagpart

C_FLOAT64	mMaxRealOfComplex

C_FLOAT64	mMaxrealpart

C_INT	mN

size_t	mNcplxconj

size_t	mNimag

size_t	mNnegreal

size_t	mNposreal

size_t	mNreal

size_t	mNzero

C_FLOAT64	mOscillationIndicator

C_FLOAT64	mOscillationIndicator_EV

C_LOGICAL *	mpBWork

C_INT32 *	mpSelect

C_FLOAT64 *	mpVS

CVector< C_FLOAT64 >	mR

C_FLOAT64	mResolution

C_INT	mSdim

char	mSort

C_FLOAT64	mStiffness

CVector< C_FLOAT64 >	mWork

Friends
std::ostream &	operator<< (std::ostream &os, const CEigen &A)

Additional Inherited Members
Public Types inherited from CCopasiContainer
typedef std::multimap < std::string, CCopasiObject * >	objectMap

Public Types inherited from CCopasiObject
typedef std::set< const CCopasiObject * >	DataObjectSet

typedef std::vector< Refresh * >	DataUpdateSequence

Public Types inherited from CObjectInterface
typedef std::set< const CObjectInterface * >	ObjectSet

typedef std::vector < CObjectInterface * >	UpdateSequence

Static Public Member Functions inherited from CCopasiObject
static std::vector< Refresh * >	buildUpdateSequence (const DataObjectSet &objects, const DataObjectSet &uptoDateObjects, const DataObjectSet &context=DataObjectSet())

static void	setRenameHandler (CRenameHandler *rh)

Static Public Attributes inherited from CCopasiContainer
static const std::vector < CCopasiContainer * >	EmptyList

Protected Types inherited from CCopasiObject
enum	Flag { Container = 0x1, Vector = 0x2, Matrix = 0x4, NameVector = 0x8, Reference = 0x10, ValueBool = 0x20, ValueInt = 0x40, ValueInt64 = 0x80, ValueDbl = 0x100, NonUniqueName = 0x200, StaticString = 0x400, ValueString = 0x800, Separator = 0x1000, ModelEntity = 0x2000, Array = 0x4000, DataModel = 0x8000, Root = 0x10000, Gui = 0x20000 }

Protected Member Functions inherited from CCopasiContainer
template<class CType >
CCopasiObject *	addMatrixReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

template<class CType >
CCopasiObject *	addObjectReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

template<class CType >
CCopasiObject *	addVectorReference (const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)

void	initObjects ()

Protected Member Functions inherited from CCopasiObject
	CCopasiObject ()

	CCopasiObject (const std::string &name, const CCopasiContainer *pParent=NULL, const std::string &type="CN", const unsigned C_INT32 &flag=0)

Protected Attributes inherited from CCopasiContainer
objectMap	mObjects

Static Protected Attributes inherited from CCopasiObject
static CRenameHandler *	smpRenameHandler = NULL

Detailed Description

File name: CEigen.h

Programmer: Yongqun He Contact email: yohe@.nosp@m.vt.e.nosp@m.du Purpose: This is the .h file for the class CEigen. It is to calculate eigenvalues and eigenvectors of a matrix. It mainly uses the dgees_() subroutine of CLAPACK

Definition at line 36 of file CEigen.h.

Constructor & Destructor Documentation

CEigen::CEigen	(	const std::string &	name = `"NoName"`,
		const CCopasiContainer *	pParent = `NULL`
	)

Default constructor

Parameters

const	std::string & name (default: "NoName")
const	CCopasiContainer * pParent (default: NULL)

File name: CEigen.cpp

Programmer: Yongqun He Contact email: yohe@.nosp@m.vt.e.nosp@m.du Purpose: This is the .cpp file for the class CEigen. It is to calculate eigenvalues and eigenvectors of a matrix. Default constructor

Definition at line 50 of file CEigen.cpp.

References CONSTRUCTOR_TRACE, and initObjects().

                                                 :
   CCopasiContainer(name, pParent, "Eigen Values",
                    CCopasiObject::Container),
   mMaxrealpart(0),
   mMaximagpart(0),
   mNposreal(0),
   mNnegreal(0),
   mNreal(0),
   mNimag(0),
   mNcplxconj(0),
   mNzero(0),
   mStiffness(0),
   mHierarchy(0),
 
   mMaxRealOfComplex(0.0),
   mImagOfMaxComplex(0.0),
   mFreqOfMaxComplex(0.0),
   mOscillationIndicator(0.0),
   mOscillationIndicator_EV(0.0),
   mBifurcationIndicator_Hopf(0.0),
   mBifurcationIndicator_Fold(0.0),
   mBifurcationIndicator_Hopf_BDT(0.0),
   mBifurcationIndicator_Fold_BDT(0.0),
 
   mResolution(0),
   mJobvs('N'),
   mSort('N'),
   mpSelect(NULL),
   mN(0),
   mA(),
   mLDA(0),
   mSdim(0),
   mR(),
   mI(),
   mpVS(NULL),
   mLdvs(1),
   mWork(1),
   mLWork(4096),
   mpBWork(NULL),
   mInfo(0)
 {
   CONSTRUCTOR_TRACE;
   initObjects();
 }

CEigen::CEigen	(	const CEigen &	src,
		const CCopasiContainer *	pParent = `NULL`
	)

Copy constructor

Parameters

const	CMetab & src
const	CCopasiContainer * pParent (default: NULL)

Definition at line 96 of file CEigen.cpp.

References CONSTRUCTOR_TRACE, and initObjects().

                                                 :
   CCopasiContainer(src, pParent),
   mMaxrealpart(src.mMaxrealpart),
   mMaximagpart(src.mMaximagpart),
   mNposreal(src.mNposreal),
   mNnegreal(src.mNnegreal),
   mNreal(src.mNreal),
   mNimag(src.mNimag),
   mNcplxconj(src.mNcplxconj),
   mNzero(src.mNzero),
   mStiffness(src.mStiffness),
   mHierarchy(src.mHierarchy),
 
   mMaxRealOfComplex(src.mMaxRealOfComplex),
   mImagOfMaxComplex(src.mImagOfMaxComplex),
   mFreqOfMaxComplex(src.mFreqOfMaxComplex),
   mOscillationIndicator(src.mOscillationIndicator),
   mOscillationIndicator_EV(src.mOscillationIndicator_EV),
   mBifurcationIndicator_Hopf(src.mBifurcationIndicator_Hopf),
   mBifurcationIndicator_Fold(src.mBifurcationIndicator_Fold),
   mBifurcationIndicator_Hopf_BDT(src.mBifurcationIndicator_Hopf_BDT),
   mBifurcationIndicator_Fold_BDT(src.mBifurcationIndicator_Fold_BDT),
 
   mResolution(src.mResolution),
   mJobvs(src.mJobvs),
   mSort(src.mSort),
   mpSelect(NULL),
   mN(src.mN),
   mA(src.mA),
   mLDA(src.mLDA),
   mSdim(src.mSdim),
   mR(src.mR),
   mI(src.mI),
   mpVS(NULL),
   mLdvs(src.mLdvs),
   mWork(src.mWork),
   mLWork(src.mLWork),
   mpBWork(NULL),
   mInfo(src.mInfo)
 {
   CONSTRUCTOR_TRACE;
   initObjects();
 }

CEigen::~CEigen ( )

virtual

Destructor

Deconstructor

Definition at line 144 of file CEigen.cpp.

References cleanup(), and DESTRUCTOR_TRACE.

 {
   cleanup();
   DESTRUCTOR_TRACE;
 }

Member Function Documentation

void CEigen::calcEigenValues ( const CMatrix< C_FLOAT64 > & matrix )

Eigenvalue calculations

Parameters

const	C_FLOAT64 * matrix
const	unsigned C_INT32 & dim

Definition at line 243 of file CEigen.cpp.

References CVectorCore< CType >::array(), CMatrix< CType >::array(), C_FLOAT64, C_INT, dgees_(), CCopasiMessage::EXCEPTION, fatalError, initialize(), mA, max, MCEigen, mI, mInfo, mJobvs, mLDA, mLdvs, mLWork, mN, mpBWork, mpVS, mR, mSdim, mSort, mWork, CMatrix< CType >::numCols(), CMatrix< CType >::numRows(), CMatrix< CType >::resize(), CVector< CType >::resize(), CMatrix< CType >::size(), and CCopasiMessage::WARNING.

Referenced by CSteadyStateTask::process().

 {
   assert(matrix.numRows() == matrix.numCols());
   mN = (C_INT) matrix.numRows();
   initialize();
 
   if (!mN) return;
 
   // copy the jacobian into mA
   mA.resize(matrix.numRows(), matrix.numCols());
   C_FLOAT64 * pA = mA.array();
   C_FLOAT64 * pAEnd = pA + mA.size();
   const C_FLOAT64 * pMatrix = matrix.array();
 
   for (; pA != pAEnd; ++pA, ++pMatrix)
     {
       *pA = *pMatrix;
 
       if (!finite(*pA) && !isnan(*pA))
         {
           if (*pA > 0)
             *pA = std::numeric_limits< C_FLOAT64 >::max();
           else
             *pA = - std::numeric_limits< C_FLOAT64 >::max();
         }
     }
 
   // Querry for the work array size.
   mLWork = -1;
   dgees_(&mJobvs, // 'N'
          &mSort, // 'N'
          NULL, // NULL,
          &mN, // n,
          mA.array(),
          & mLDA,
          & mSdim, // output
          mR.array(),
          mI.array(),
          mpVS,
          & mLdvs,
          mWork.array(),
          & mLWork,
          mpBWork, // NULL
          &mInfo);        // output
 
   if (mInfo != 0)
     {
       // Exception
       CCopasiMessage(CCopasiMessage::EXCEPTION, MCEigen + 1, -mInfo);
     }
 
   mLWork = (C_INT) mWork[0];
   mWork.resize(mLWork);
 
   // calculate the eigenvalues
   /* 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);
    *  Arguments
    *  =========
    *
    *  JOBVS   (input) CHARACTER*1
    *          = 'N': Schur vectors are not computed;
    *          = 'V': Schur vectors are computed.
    *
    *  SORT    (input) CHARACTER*1
    *          Specifies whether or not to order the eigenvalues on the
    *          diagonal of the Schur form.
    *          = 'N': Eigenvalues are not ordered;
    *          = 'S': Eigenvalues are ordered (see SELECT).
    *
    *  SELECT  (input) LOGICAL FUNCTION of two DOUBLE PRECISION arguments
    *          SELECT must be declared EXTERNAL in the calling subroutine.
    *          If SORT = 'S', SELECT is used to select eigenvalues to sort
    *          to the top left of the Schur form.
    *          If SORT = 'N', SELECT is not referenced.
    *          An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
    *          SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
    *          conjugate pair of eigenvalues is selected, then both complex
    *          eigenvalues are selected.
    *          Note that a selected complex eigenvalue may no longer
    *          satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
    *          ordering may change the value of complex eigenvalues
    *          (especially if the eigenvalue is ill-conditioned); in this
    *          case INFO is set to N+2 (see INFO below).
    *
    *  N       (input) INTEGER
    *          The order of the matrix A. N >= 0.
    *
    *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
    *          On entry, the N-by-N matrix A.
    *          On exit, A has been overwritten by its real Schur form T.
    *
    *  LDA     (input) INTEGER
    *          The leading dimension of the array A.  LDA >= max(1,N).
    *
    *  SDIM    (output) INTEGER
    *          If SORT = 'N', SDIM = 0.
    *          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
    *                         for which SELECT is true. (Complex conjugate
    *                         pairs for which SELECT is true for either
    *                         eigenvalue count as 2.)
    *
    *  WR      (output) DOUBLE PRECISION array, dimension (N)
    *  WI      (output) DOUBLE PRECISION array, dimension (N)
    *          WR and WI contain the real and imaginary parts,
    *          respectively, of the computed eigenvalues in the same order
    *          that they appear on the diagonal of the output Schur form T.
    *          Complex conjugate pairs of eigenvalues will appear
    *          consecutively with the eigenvalue having the positive
    *          imaginary part first.
    *
    *  VS      (output) DOUBLE PRECISION array, dimension (LDVS,N)
    *          If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
    *          vectors.
    *          If JOBVS = 'N', VS is not referenced.
    *
    *  LDVS    (input) INTEGER
    *          The leading dimension of the array VS.  LDVS >= 1; if
    *          JOBVS = 'V', LDVS >= N.
    *
    *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
    *          On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
    *
    *  LWORK   (input) INTEGER
    *          The dimension of the array WORK.  LWORK >= max(1,3*N).
    *          For good performance, LWORK must generally be larger.
    *
    *          If LWORK = -1, then a workspace query is assumed; the routine
    *          only calculates the optimal size of the WORK array, returns
    *          this value as the first entry of the WORK array, and no error
    *          message related to LWORK is issued by XERBLA.
    *
    *  BWORK   (workspace) LOGICAL array, dimension (N)
    *          Not referenced if SORT = 'N'.
    *
    *  INFO    (output) INTEGER
    *          = 0: successful exit
    *          < 0: if INFO = -i, the i-th argument had an illegal value.
    *          > 0: if INFO = i, and i is
    *             <= N: the QR algorithm failed to compute all the
    *                   eigenvalues; elements 1:i-1 and i+1:N of WR and WI
    *                   contain those eigenvalues which have converged; if
    *                   JOBVS = 'V', VS contains the matrix which reduces A
    *                   to its partially converged Schur form.
    *             = N+1: the eigenvalues could not be reordered because some
    *                   eigenvalues were too close to separate (the problem
    *                   is very ill-conditioned);
    *             = N+2: after reordering, roundoff changed values of some
    *                   complex eigenvalues so that leading eigenvalues in
    *                   the Schur form no longer satisfy SELECT=.TRUE.  This
    *                   could also be caused by underflow due to scaling.
    *
    */
   dgees_(&mJobvs, // 'N'
          &mSort, // 'N'
          NULL, // NULL,
          &mN, // n,
          mA.array(),
          & mLDA,
          & mSdim, // output
          mR.array(),
          mI.array(),
          mpVS,
          & mLdvs,
          mWork.array(),
          & mLWork,
          mpBWork, // NULL
          &mInfo);        // output
 
   if (mInfo != 0)
     {
       if (mInfo < 0)
         {
           // Exception
           CCopasiMessage(CCopasiMessage::EXCEPTION, MCEigen + 1, -mInfo);
         }
       else if (mInfo <= mN)
         {
           // Warning
           CCopasiMessage(CCopasiMessage::WARNING, MCEigen + 2, mInfo);
         }
       else if (mInfo == mN + 1)
         {
           // Warning
           CCopasiMessage(CCopasiMessage::WARNING, MCEigen + 3, mInfo);
         }
       else if (mInfo == mN + 2)
         {
           // Warning
           CCopasiMessage(CCopasiMessage::EXCEPTION, MCEigen + 4, mInfo);
         }
       else // Catch all, should never happen.
         fatalError();
     }
 
   return;
 }

void CEigen::cleanup ( )

cleanup()

Definition at line 240 of file CEigen.cpp.

Referenced by initialize(), and ~CEigen().

241 {}

const C_FLOAT64 & CEigen::getHierarchy ( ) const

Get the eigenvalue hierarchy

Definition at line 220 of file CEigen.cpp.

References mHierarchy.

 {
   return mHierarchy;
 }

const CVector< C_FLOAT64 > & CEigen::getI ( ) const

Definition at line 687 of file CEigen.cpp.

References mI.

Referenced by StateSubwidget::loadJacobian(), operator<<(), and CSteadyStateTask::process().

688 {return mI;}

CEigen::mI

CVector< C_FLOAT64 > mI

Definition: CEigen.h:203

const C_FLOAT64 & CEigen::getMaximagpart ( ) const

Get the max eigenvalue imaginary part

Definition at line 202 of file CEigen.cpp.

References mMaximagpart.

 {
   return mMaximagpart;
 }

const C_FLOAT64 & CEigen::getMaxrealpart ( ) const

Get the max eigenvalue real part

return the matrix Set the Matrix

Definition at line 196 of file CEigen.cpp.

References mMaxrealpart.

 {
   return mMaxrealpart;
 }

const size_t & CEigen::getNcplxconj ( ) const

Definition at line 666 of file CEigen.cpp.

References mNcplxconj.

 {
   return mNcplxconj;
 }

const size_t & CEigen::getNimag ( ) const

Return the number of imaginary eigenvalue numbers

Definition at line 661 of file CEigen.cpp.

References mNimag.

Referenced by CLNATask::process().

 {
   return mNimag;
 }

const size_t & CEigen::getNnegreal ( ) const

Return the number of eigenvalues with negative real part

Definition at line 682 of file CEigen.cpp.

References mNnegreal.

 {
   return mNnegreal;
 }

const size_t & CEigen::getNposreal ( ) const

Return the number of eigenvalues with positive real part

Definition at line 674 of file CEigen.cpp.

References mNposreal.

Referenced by CLNATask::process().

 {
   return mNposreal;
 }

const size_t & CEigen::getNreal ( ) const

Return number of real eigenvalues WeiSun 3/28/02

Definition at line 653 of file CEigen.cpp.

References mNreal.

 {
   return mNreal;
 }

const size_t & CEigen::getNzero ( ) const

Get the number of zero eigenvalues

Definition at line 208 of file CEigen.cpp.

References mNzero.

Referenced by CLNATask::process().

 {
   return mNzero;
 }

const CVector< C_FLOAT64 > & CEigen::getR ( ) const

Definition at line 690 of file CEigen.cpp.

References mR.

Referenced by StateSubwidget::loadJacobian(), operator<<(), and CSteadyStateTask::process().

691 {return mR;}

CEigen::mR

CVector< C_FLOAT64 > mR

Definition: CEigen.h:198

const C_FLOAT64 & CEigen::getStiffness ( ) const

Get the eigenvalue stiffness

Definition at line 214 of file CEigen.cpp.

References mStiffness.

 {
   return mStiffness;
 }

void CEigen::initialize ( )

initialize variables for eigenvalue calculations

Definition at line 227 of file CEigen.cpp.

References cleanup(), mI, mLDA, mN, mNcplxconj, mNimag, mNnegreal, mNposreal, mNreal, mNzero, mR, and CVector< CType >::resize().

Referenced by calcEigenValues().

 {
   cleanup();
 
   mNreal = mNimag = mNposreal = mNnegreal =
                                   mNzero = mNcplxconj = 0;
 
   mLDA = mN > 1 ? mN : 1;
 
   mR.resize(mN);
   mI.resize(mN);
 }

void CEigen::initObjects ( )

private

Initialize the contained CCopasiObjects

Definition at line 152 of file CEigen.cpp.

References CCopasiContainer::addObjectReference(), CCopasiContainer::addVectorReference(), mBifurcationIndicator_Fold, mBifurcationIndicator_Fold_BDT, mBifurcationIndicator_Hopf, mBifurcationIndicator_Hopf_BDT, mFreqOfMaxComplex, mHierarchy, mI, mImagOfMaxComplex, mMaximagpart, mMaxRealOfComplex, mMaxrealpart, mNcplxconj, mNimag, mNnegreal, mNposreal, mNreal, mNzero, mOscillationIndicator, mOscillationIndicator_EV, mR, mResolution, mStiffness, CCopasiObject::ValueDbl, and CCopasiObject::ValueInt.

Referenced by CEigen().

 {
   addObjectReference("Maximum real part", mMaxrealpart, CCopasiObject::ValueDbl);
   addObjectReference("Maximum imaginary part", mMaximagpart, CCopasiObject::ValueDbl);
   addObjectReference("# Positive eigenvalues", mNposreal, CCopasiObject::ValueInt);
   addObjectReference("# Negative eigenvalues", mNnegreal, CCopasiObject::ValueInt);
   addObjectReference("# Real eigenvalues", mNreal, CCopasiObject::ValueInt);
   addObjectReference("# Imaginary eigenvalues", mNimag, CCopasiObject::ValueInt);
   addObjectReference("# Complex conjugated eigenvalues", mNcplxconj, CCopasiObject::ValueInt);
   addObjectReference("# Zero eigenvalues", mNzero, CCopasiObject::ValueInt);
   addObjectReference("Stiffness", mStiffness, CCopasiObject::ValueDbl);
   addObjectReference("Time hierachy", mHierarchy, CCopasiObject::ValueDbl);
   addObjectReference("Resolution", mResolution, CCopasiObject::ValueDbl);
   addVectorReference("Vector of real part of eigenvalues", mR, CCopasiObject::ValueDbl);
   addVectorReference("Vector of imaginary part of eigenvalues", mI, CCopasiObject::ValueDbl);
 
   addObjectReference("Maximum real part of complex eigenvalue", mMaxRealOfComplex, CCopasiObject::ValueDbl);
   addObjectReference("Imaginary part of largest complex eigenvalue", mImagOfMaxComplex, CCopasiObject::ValueDbl);
   addObjectReference("Linear Frequency of largest complex eigenvalue", mFreqOfMaxComplex, CCopasiObject::ValueDbl);
   addObjectReference("Oscillation indicator", mOscillationIndicator, CCopasiObject::ValueDbl);
   addObjectReference("EV-based oscillation indicator", mOscillationIndicator_EV, CCopasiObject::ValueDbl);
   addObjectReference("Hopf bifurcation test function", mBifurcationIndicator_Hopf, CCopasiObject::ValueDbl);
   addObjectReference("Fold bifurcation test function", mBifurcationIndicator_Fold, CCopasiObject::ValueDbl);
   addObjectReference("Hopf bifurcation test function (BDT)", mBifurcationIndicator_Hopf_BDT, CCopasiObject::ValueDbl);
   addObjectReference("Fold bifurcation test function (BDT)", mBifurcationIndicator_Fold_BDT, CCopasiObject::ValueDbl);
 }

void CEigen::print ( std::ostream * ostream ) const

virtual

This is the output method for any object. The default implementation provided with CCopasiObject uses the ostream operator<< of the object to print the object.To overide this default behaviour one needs to reimplement the virtual print function.

Parameters

std::ostream * ostream

Reimplemented from CCopasiObject.

Definition at line 150 of file CEigen.cpp.

150 {(*ostream) << (*this);}

void CEigen::stabilityAnalysis ( const C_FLOAT64 & resolution )

Calculate various eigenvalue statistics

Definition at line 454 of file CEigen.cpp.

References CVectorCore< CType >::applyPivot(), CVectorCore< CType >::array(), C_FLOAT64, C_INT, C_INT32, mBifurcationIndicator_Fold, mBifurcationIndicator_Fold_BDT, mBifurcationIndicator_Hopf, mBifurcationIndicator_Hopf_BDT, mFreqOfMaxComplex, mHierarchy, mI, mImagOfMaxComplex, mMaximagpart, mMaxRealOfComplex, mMaxrealpart, mN, mNcplxconj, mNimag, mNnegreal, mNposreal, mNreal, mNzero, mOscillationIndicator, mOscillationIndicator_EV, mR, mResolution, mStiffness, CVectorCore< CType >::size(), and sortWithPivot().

Referenced by CSteadyStateTask::process().

 {
   if (!mN) return;
 
   size_t mx, mn;            // YH: n is the 4th parameter, not here
   size_t i;
   C_FLOAT64 distt, maxt, tott;
   mResolution = resolution;
 
   // sort the eigenvalues
   CVector<size_t> Pivot;
 
   sortWithPivot(mR.array(), mR.array() + mR.size(), CompareDoubleWithNaN(), Pivot);
 
   // The sort order is ascending however we need descending
   size_t *pTo = Pivot.array();
   size_t *pFrom = pTo + mN - 1;
 
   for (; pTo < pFrom; ++pTo, --pFrom)
     {
       size_t Tmp = *pFrom;
       *pFrom = *pTo;
       *pTo = Tmp;
     }
 
   mR.applyPivot(Pivot);
   mI.applyPivot(Pivot);
 
   // calculate various eigenvalue statistics
   mMaxrealpart = mR[0];
   mMaxRealOfComplex = mR[mN - 1] - 1;
   mMaximagpart = fabs(mI[0]);
   mImagOfMaxComplex = 0.0;
 
   for (i = 0; (C_INT) i < mN; i++)
     {
       // for the largest real part
       if (mR[i] > mMaxrealpart)
         mMaxrealpart = mR[i];
 
       // for the largest imaginary part
       if (fabs(mI[i]) > mMaximagpart)
         mMaximagpart = fabs(mI[i]);
 
       // for the largest complex eigenvalue
       if ((fabs(mI[i]) > resolution) && (mR[i] > mMaxRealOfComplex))
         {
           mMaxRealOfComplex = mR[i];
           mImagOfMaxComplex = fabs(mI[i]);
         }
 
       if (fabs(mR[i]) > resolution)
         {
           // positive real part
           if (mR[i] >= resolution)
             mNposreal++;
 
           // negative real part
           if (mR[i] <= -resolution)
             mNnegreal++;
 
           if (fabs(mI[i]) > resolution)
             {
               // complex
               mNcplxconj++;
             }
           else
             {
               mI[i] = 0.0;
               // pure real
               mNreal++;
             }
         }
       else
         {
           mR[i] = 0.0;
 
           if (fabs(mI[i]) > resolution)
             {
               // pure imaginary
               mNimag++;
             }
           else
             {
               mI[i] = 0.0;
               // zero
               mNzero++;
             }
         }
     }
 
   if (mImagOfMaxComplex == 0)
     mMaxRealOfComplex = mR[mN - 1]; //default value
 
   mFreqOfMaxComplex = mImagOfMaxComplex / (2 * M_PI);
 
   if (mNposreal > 0)
     {
       if (mR[0] > fabs(mR[mN - 1]))
         mx = 0;
       else
         mx = mN - 1;
 
       if ((C_INT) mNposreal == mN)
         mn = mNposreal - 1;
       else if (mR[mNposreal - 1] < fabs(mR[mNposreal]))
         mn = mNposreal - 1;
       else
         mn = mNposreal;
     }
   else
     {
       mx = mN - 1; // index of the largest absolute real part
       mn = 0; // index of the smallest absolute real part
     }
 
   mStiffness = fabs(mR[mx]) / fabs(mR[mn]);
 
   maxt = tott = fabs(1 / mR[mn]);
   distt = 0.0;
 
   for (i = 1; (C_INT) i < mN; i++)
     if (i != mn)
       {
         distt += maxt - fabs(1 / mR[i]);
         tott += fabs(1 / mR[i]);
       }
 
   mHierarchy = distt / tott / (mN - 1);
 
   //we calculate an oscillation indicator based on the eigenvalues. It is using a rather heuristical approach
   if (mN < 2) // at least 2 Eigenvalues are required for oscillations
     mOscillationIndicator_EV = 0;
   else
     {
       if (fabs(mI[0]) > resolution)
         mOscillationIndicator_EV = mR[0] * (mR[0] > 0 ? fabs(mI[0]) / (0.01 + fabs(mI[0])) : 2 - fabs(mI[0]) / (0.01 + fabs(mI[0])));
       else
         mOscillationIndicator_EV = 2 * mR[1] - (mR[0] > 0 ? 2 * mR[0] : 0);
     }
 
   //bifurcation test functions, according to Kuznetsov "Elements of applied bifurcation theory"
   // this function can also be calculated without using the eigenvalues, but since we already
   //have them we can use them
 
   std::complex<C_FLOAT64> tmpcpl = 1.0;
   unsigned C_INT32 index_min = 0; // the index of the EV with smallest abs real part
   C_FLOAT64 tmpmin = 1e300;
 
   for (i = 0; (C_INT) i < mN; i++)
     {
       tmpcpl *= std::complex<C_FLOAT64>(mR[i], mI[i]);
 
       if (fabs(mR[i]) < tmpmin)
         {
           tmpmin = fabs(mR[i]);
           index_min = i;
         }
     }
 
   mBifurcationIndicator_Fold = std::abs<C_FLOAT64>(tmpcpl);
 
   tmpcpl = 1.0;
 
   for (i = 0; (C_INT) i < mN - 1; i++)
     {
       tmpcpl *= (std::complex<C_FLOAT64>(mR[i], mI[i]) + std::complex<C_FLOAT64>(mR[i + 1], mI[i + 1]));
     }
 
   mBifurcationIndicator_Hopf = (tmpcpl.real());
 
   //now the test functions from the "Bifurcation discovery tool" (as described in the publication by Chickarmane et al.)
 
   //calculate the product of EVs excluding the minimal one
   tmpcpl = 1.0;
 
   for (i = 0; (C_INT) i < mN; i++)
     {
       if (i != index_min)
         tmpcpl *= std::complex<C_FLOAT64>(mR[i], mI[i]);
     }
 
   mBifurcationIndicator_Fold_BDT = mBifurcationIndicator_Fold / (1 - 0.99 * exp(-std::abs<C_FLOAT64>(tmpcpl)));
 
   C_FLOAT64 tmp_product = 1.0;
 
   for (i = 0; (C_INT) i < mN; i++)
     {
       tmp_product *= mR[i] / (1 - 0.99 * exp(-mI[i]));
     }
 
   mBifurcationIndicator_Hopf_BDT = tmp_product;
 
   mOscillationIndicator = 0.0; //not used at the moment.
 }

Friends And Related Function Documentation

std::ostream& operator<<	(	std::ostream &	os,
		const CEigen &	A
	)

friend

Definition at line 693 of file CEigen.cpp.

 {
   os << std::endl;
   os << "KINETIC STABILITY ANALYSIS";
   os << std::endl;
   os << "The linear stability analysis based on the eigenvalues" << std::endl;
   os << "of the Jacobian matrix is only valid for steady states." << std::endl;
   os << std::endl;
   os << "Summary:" << std::endl;
   os << "This state ";
 
   // Output statistics
 
   if (A.mMaxrealpart > A.mResolution)
     os << "is unstable";
   else if (A.mMaxrealpart < -A.mResolution)
     os << "is asymptotically stable";
   else
     os << "'s stability is undetermined";
 
   if (A.mMaximagpart > A.mResolution)
     {
       os << "," << std::endl;
       os << "transient states in its vicinity have oscillatory components";
     }
 
   os << "." << std::endl;
   os << std::endl;
 
   os << "Eigenvalue statistics:" << std::endl;
   // Output Max Real Part
   os << " Largest real part: ";
   os << std::setprecision(6) << A.mMaxrealpart << std::endl;
   // Output Max imaginary Part
   os << " Largest absolute imaginary part:  ";
   os << std::setprecision(6) << A.mMaximagpart << std::endl;
 
   if (A.mImagOfMaxComplex > A.mResolution)
     os << " The complex eigenvalues with the largest real part are:  "
        << A.mMaxRealOfComplex << " +|- " << A.mImagOfMaxComplex << "i" << std::endl;
 
   // Output Eigen-nreal
   os.unsetf(std::ios_base::scientific);
   os.unsetf(std::ios_base::showpoint);
   os << " " << A.mNreal;
   os << " are purely real" << std::endl;
   // Output Eigen-nimage
   os << " " << A.mNimag;
   os << " are purely imaginary" << std::endl;
   // Output Eigen-ncplxconj
   os << " " << A.mNcplxconj;
   os << " are complex" << std::endl;
   // Output Eigen-nzero
   os << " " << A.mNzero;
   os << " are equal to zero" << std::endl;
   // Output Eigen-nposreal
   os << " " << A.mNposreal;
   os << " have positive real part" << std::endl;
   // Output Eigen-nnegreal
   os << " " << A.mNnegreal;
   os << " have negative real part" << std::endl;
 
   // Set point manipulators
   os.setf(std::ios_base::showpoint);
   // Output Eigne-stiffness
   os << " stiffness = " << A.mStiffness << std::endl;
   os << " time hierarchy = " << A.mHierarchy << std::endl;
 
   return os;
 }

Member Data Documentation

CMatrix< C_FLOAT64 > CEigen::mA

private

#5: (input/output) The double precision array, dimension (LDA,N) On entry, the N-by-N matrix A On exit, A has been overwritten by its real Schur form T Use a pointer variable here instead of array since we don't know the dimension yet.

Definition at line 180 of file CEigen.h.

Referenced by calcEigenValues().

C_FLOAT64 CEigen::mBifurcationIndicator_Fold

private

A bifurcation test function, it is 0 for a Fold-type bifurcation.

Definition at line 126 of file CEigen.h.

Referenced by initObjects(), and stabilityAnalysis().

C_FLOAT64 CEigen::mBifurcationIndicator_Fold_BDT

private

A bifurcation test function, defined according to the "Bifurcation Discovery Tool", Chickarmane et al 2005

Definition at line 136 of file CEigen.h.

Referenced by initObjects(), and stabilityAnalysis().

C_FLOAT64 CEigen::mBifurcationIndicator_Hopf

private

A bifurcation test function, it is 0 for a Hopf-bifurcation. Definition according to Kuznetsov. This is only a necessary, not a sufficient condition for a H.-bifurcation

Definition at line 121 of file CEigen.h.

Referenced by initObjects(), and stabilityAnalysis().

C_FLOAT64 CEigen::mBifurcationIndicator_Hopf_BDT

private

A bifurcation test function, defined according to the "Bifurcation Discovery Tool", Chickarmane et al 2005

Definition at line 131 of file CEigen.h.

Referenced by initObjects(), and stabilityAnalysis().

C_FLOAT64 CEigen::mFreqOfMaxComplex

private

frequency of complex with largest real part

Definition at line 105 of file CEigen.h.

Referenced by initObjects(), and stabilityAnalysis().

C_FLOAT64 CEigen::mHierarchy

private

the hierary of the eigenvalues

Definition at line 89 of file CEigen.h.

Referenced by getHierarchy(), initObjects(), operator<<(), and stabilityAnalysis().

CVector< C_FLOAT64 > CEigen::mI

private

#9: array with dimension (mN)

Definition at line 203 of file CEigen.h.

Referenced by calcEigenValues(), getI(), initialize(), initObjects(), and stabilityAnalysis().

C_FLOAT64 CEigen::mImagOfMaxComplex

private

imaginary part of complex with largest real part

Definition at line 100 of file CEigen.h.

Referenced by initObjects(), operator<<(), and stabilityAnalysis().

C_INT CEigen::mInfo

private

#15: (output) an integer =0: successful exit <0: if mInfo=-i, the ith argument had an illegal value >0: if mInfo=i, and i is <=N: the QR algorithm failed to compute all the eigenvalues; =N+1: the eigenvalues could not be reordered because some eigenvalues were too close to separate (ill-conditioned) =N+2: after reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues in the Schur form no longer satisfy mSelect=.True. This could caused by underflow due to scaling

Definition at line 249 of file CEigen.h.

Referenced by calcEigenValues().

char CEigen::mJobvs

private

#1: (input) characer*1 = 'N': Schur vectors are not computed = 'V': Schur vectors are computed

Definition at line 150 of file CEigen.h.

Referenced by calcEigenValues().

C_INT CEigen::mLDA

private

#6: (input) The leading dimension of the array A. LDA >= max(1,N)

Definition at line 185 of file CEigen.h.

Referenced by calcEigenValues(), and initialize().

C_INT CEigen::mLdvs

private

#11: an integer, the leading dimension of the array VS. mLdvs >= 1; if mJobvs='V', mLdvs >= N.

Definition at line 216 of file CEigen.h.

Referenced by calcEigenValues().

C_INT CEigen::mLWork

private

#13: (input) Dimension of array Work, its value >= max(1,3*mN). For good performance, it must generally be larger

Definition at line 228 of file CEigen.h.

Referenced by calcEigenValues().

C_FLOAT64 CEigen::mMaximagpart

private

the imaginary part of the maximum eigenvalue

Definition at line 49 of file CEigen.h.

Referenced by getMaximagpart(), initObjects(), operator<<(), and stabilityAnalysis().

C_FLOAT64 CEigen::mMaxRealOfComplex

private

largest real part of a complex eigenvalue. If there is no complex eigenvalue it will be the smallest ev.

Definition at line 95 of file CEigen.h.

Referenced by initObjects(), operator<<(), and stabilityAnalysis().

C_FLOAT64 CEigen::mMaxrealpart

private

the real part of the maximum eigenvalue

Definition at line 44 of file CEigen.h.

Referenced by getMaxrealpart(), initObjects(), operator<<(), and stabilityAnalysis().

C_INT CEigen::mN

private

#4: (input) The order of the matrix A

Definition at line 171 of file CEigen.h.

Referenced by calcEigenValues(), initialize(), and stabilityAnalysis().

size_t CEigen::mNcplxconj

private

Definition at line 74 of file CEigen.h.

Referenced by getNcplxconj(), initialize(), initObjects(), operator<<(), and stabilityAnalysis().

size_t CEigen::mNimag

private

the number of imaginary eigenvalue numbers

Definition at line 69 of file CEigen.h.

Referenced by getNimag(), initialize(), initObjects(), operator<<(), and stabilityAnalysis().

size_t CEigen::mNnegreal

private

the number of eigenvalues with negative real part

Definition at line 59 of file CEigen.h.

Referenced by getNnegreal(), initialize(), initObjects(), operator<<(), and stabilityAnalysis().

size_t CEigen::mNposreal

private

the number of eigenvalues with positive real part

Definition at line 54 of file CEigen.h.

Referenced by getNposreal(), initialize(), initObjects(), operator<<(), and stabilityAnalysis().

size_t CEigen::mNreal

private

the number of real eigenvalues

Definition at line 64 of file CEigen.h.

Referenced by getNreal(), initialize(), initObjects(), operator<<(), and stabilityAnalysis().

size_t CEigen::mNzero

private

the number of eigenvalues with value of zero

Definition at line 79 of file CEigen.h.

Referenced by getNzero(), initialize(), initObjects(), operator<<(), and stabilityAnalysis().

C_FLOAT64 CEigen::mOscillationIndicator

private

An index that is supposed to indicate the presence of oscillations

Definition at line 110 of file CEigen.h.

Referenced by initObjects(), and stabilityAnalysis().

C_FLOAT64 CEigen::mOscillationIndicator_EV

private

A heuristical index that is supposed to indicate the presence of oscillations, based on Eigenvalues

Definition at line 115 of file CEigen.h.

Referenced by initObjects(), and stabilityAnalysis().

C_LOGICAL* CEigen::mpBWork

private

#14: (workspace) Logical array, dimension (N) Not referenced if mSort = 'N'

Definition at line 234 of file CEigen.h.

Referenced by calcEigenValues().

C_INT32* CEigen::mpSelect

private

#3: (input) Logical function of two double precision arguments It must be declared external = 'S': Select is used to select eigenvalues to sort to the top left of the Schur form = 'N': Eigenvalues are ordered Select is not refereced.

Definition at line 166 of file CEigen.h.

C_FLOAT64* CEigen::mpVS

private

#10: (output) array with dimension (mLdvs, mN) If mJobvs='V', mVS contains the orthogoanl matrix Z of Schur vectors If mJobvs='N', mVS is not referenced

Definition at line 210 of file CEigen.h.

Referenced by calcEigenValues().

CVector< C_FLOAT64 > CEigen::mR

private

#8: array with dimension (mN)

Definition at line 198 of file CEigen.h.

Referenced by calcEigenValues(), getR(), initialize(), initObjects(), and stabilityAnalysis().