COPASI API  4.16.103
CTSSAMethod.cpp
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) 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /**
16  * CTSSAMethod class.
17  * This class describes the interface to all time scale separation analysis methods.
18  * The various method like ILDM or CSP have to be derived from
19  * this class.
20  *
21  */
22 
23 #include <sstream>
24 
25 #include "copasi.h"
26 
27 #include "CTSSAMethod.h"
28 #include "CTSSAProblem.h"
29 #include "CILDMMethod.h"
30 #include "CILDMModifiedMethod.h"
31 #include "CCSPMethod.h"
32 
33 #include "CTSSAProblem.h"
34 #include "model/CState.h"
35 #include "model/CCompartment.h"
38 #include "model/CModel.h"
39 
40 #include "lapack/lapackwrap.h" // CLAPACK
41 #include "lapack/blaswrap.h" // BLAS
42 
45 {
46  CTSSAMethod * pMethod = NULL;
47 
48  switch (subType)
49  {
50  case unset:
51  case tssILDM:
52  pMethod = new CILDMMethod();
53  break;
54 
55  case tssILDMModified:
56  pMethod = new CILDMModifiedMethod();
57  break;
58 
59  case tssCSP:
60  pMethod = new CCSPMethod();
61  break;
62 
63  default:
64  fatalError();
65  break;
66  }
67 
68  return pMethod;
69 }
70 
71 /**
72  * Default constructor.
73  */
75  const CCopasiContainer * pParent) :
76  CCopasiMethod(CCopasiTask::tssAnalysis, subType, pParent),
77  mpCurrentState(NULL),
78  mpProblem(NULL),
79  mpState(NULL),
80  mData(),
81  mY(NULL),
82  mYdot(),
83  mY_initial(),
84  mTime(0.0),
85  mJacobian(),
86  mJacobian_initial(),
87  mQ(),
88  mQ_desc(),
89  mR(),
90  mR_desc(),
91  mTd(),
92  mTdInverse(),
93  mQz(),
94  mTd_save(),
95  mTdInverse_save(),
96  mCfast(),
97  mY_cons(),
98  mVslow(),
99  mVslow_metab(),
100  mVslow_space(),
101  mVfast_space(),
102  mSlow(0),
103  mLsodaStatus(1),
104  mReducedModel(false),
105  mRtol(1e-5),
106  mAtol(),
107  mErrorMsg(),
108  mLSODA(),
109  mState(0),
110  mDWork(),
111  mIWork(),
112  mJType(0),
113  mpModel(NULL),
114  mDtol(1e-6),
115  mEPS(0.01),
116  mVec_SlowModes(),
117  mCurrentTime(),
118  mVec_TimeScale(),
119  mCurrentStep(0)
120 {}
121 
122 /**
123  * Copy constructor.
124  * @param "const CTSSAMethod &" src
125  */
127  const CCopasiContainer * pParent):
128  CCopasiMethod(src, pParent),
129  mpCurrentState(src.mpCurrentState),
130  mpProblem(src.mpProblem),
131  mpState(NULL),
132  mData(),
133  mY(NULL),
134  mYdot(),
135  mY_initial(),
136  mTime(src.mTime),
137  mJacobian(),
138  mJacobian_initial(),
139  mQ(),
140  mQ_desc(),
141  mR(),
142  mR_desc(),
143  mTd(),
144  mTdInverse(),
145  mQz(),
146  mTd_save(),
147  mTdInverse_save(),
148  mCfast(),
149  mY_cons(),
150  mVslow(),
151  mVslow_metab(),
152  mVslow_space(),
153  mVfast_space(),
154  mSlow(0),
155  mLsodaStatus(1),
156  mReducedModel(src.mReducedModel),
157  mRtol(src.mRtol),
158  mAtol(src.mAtol),
159  mErrorMsg(),
160  mLSODA(),
161  mState(0),
162  mDWork(),
163  mIWork(),
164  mJType(0),
165  mpModel(NULL),
166  mDtol(src.mDtol),
167  mEPS(src.mEPS),
168  mVec_SlowModes(),
169  mCurrentTime(),
170  mVec_TimeScale(),
171  mCurrentStep(0)
172 {}
173 
174 /**
175  * Destructor.
176  */
178 {}
179 
181 {
182  mpCurrentState = currentState;
183 }
184 
185 /**
186  * Set a pointer to the problem.
187  * This method is used by CTSSA
188  * @param "CTSSAProblem *" problem
189  */
191 {mpProblem = problem;}
192 
193 /**
194  * This instructs the method to calculate a a time step of deltaT
195  * starting with the current state, i.e., the result of the previous
196  * step.
197  * The new state (after deltaT) is expected in the current state.
198  * The return value is the actual time step taken.
199  * @param "const double &" deltaT
200  */
201 void CTSSAMethod::step(const double & C_UNUSED(deltaT))
202 {return ;}
203 
204 /**
205  * This instructs the method to calculate a a time step of deltaT
206  * starting with the initialState given.
207  * The new state (after deltaT) is expected in the current state.
208  * The return value is the actual time step taken.
209  * @param "double &" deltaT
210  * @param "const CState *" initialState
211  * @return "const double &" actualDeltaT
212  */
213 void CTSSAMethod::start(const CState * C_UNUSED(initialState))
214 {return;}
215 
217 {
218  mpModel = model;
219 }
220 
222 {
223  return;
224 }
225 
226 //virtual
228 {
229 
230  if (!CCopasiMethod::isValidProblem(pProblem)) return false;
231 
232  const CTSSAProblem * pTP = dynamic_cast<const CTSSAProblem *>(pProblem);
233 
234  if (!pTP)
235  {
237  return false;
238  }
239 
240  CModel * pModel = pTP->getModel();
241 
242  if (pModel == NULL)
243  return false;
244 
245  if (pModel->getMetabolites().size() == 0)
246  {
248  return false;
249  }
250 
251  if (pModel->getCompartments().size() != 1)
252  {
253  CCopasiMethod::SubType subType;
254 
255  subType = mData.pMethod->getSubType();
256 
257  switch (subType)
258  {
259  case tssILDM:
260  case tssILDMModified:
261 
263  return false;
264 
265  case tssCSP:
266  return true;
267 
268  default:
269  fatalError();
270  break;
271  }
272  }
273 
274 // Check if the model has a species with an assigments or an ODE
275  if (pModel->getNumODEMetabs() != 0 || pModel->getNumAssignmentMetabs() != 0)
276  {
277  CCopasiMessage(CCopasiMessage::ERROR, "TSSA can not be applyed for systems with species determined by assigments or ODE.");
278  return false;
279  }
280 
281 // Check if the model has a compartment with an ODE
284 
285  for (; it != end; ++it)
286  if ((*it)->getStatus() == CModelEntity::ODE || (*it)->getStatus() == CModelEntity::ASSIGNMENT)
287 
288  {
289  CCopasiMessage(CCopasiMessage::ERROR, " TSSA can not be applyed for systems with non-constant volumes");
290  return false;
291  }
292 
293 // Check if the model has a model parameter with an ODE
294 
295  size_t i;
296 
297  for (i = 0; i < pModel->getModelValues().size(); i++)
298  {
299  if (pModel->getModelValues()[i]->getStatus() == CModelEntity::ODE)
300  {
301  CCopasiMessage(CCopasiMessage::ERROR, "TSSA can not be applyed for systems with parameters defined by ODE.");
302  return false;
303  }
304  }
305 
306 // Check if the model contains events
307  if (pModel->getEvents().size() != 0)
308  {
309  CCopasiMessage(CCopasiMessage::ERROR, "TSSA can not be applyed for systems with events");
310  return false;
311  }
312 
313  return true;
314 }
315 
317 {
318  if ((int) mCurrentTime.size() > step)
319  return mCurrentTime[step];
320  else
321  return std::numeric_limits<C_FLOAT64>::quiet_NaN();
322 };
323 
325 {return;}
326 
328 {
329  CCopasiParameter *pParm;
330 
331  assertParameter("Integrate Reduced Model", CCopasiParameter::BOOL, (bool) true);
332  assertParameter("Relative Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-6);
333  assertParameter("Absolute Tolerance", CCopasiParameter::UDOUBLE, (C_FLOAT64) 1.0e-12);
334  assertParameter("Max Internal Steps", CCopasiParameter::UINT, (unsigned C_INT32) 10000);
335 
336  // Check whether we have a method with the old parameter names
337  if ((pParm = getParameter("LSODA.RelativeTolerance")) != NULL)
338  {
339  setValue("Relative Tolerance", *pParm->getValue().pUDOUBLE);
340  removeParameter("LSODA.RelativeTolerance");
341 
342  if ((pParm = getParameter("LSODA.AbsoluteTolerance")) != NULL)
343  {
344  setValue("Absolute Tolerance", *pParm->getValue().pUDOUBLE);
345  removeParameter("LSODA.AbsoluteTolerance");
346  }
347 
348  if ((pParm = getParameter("LSODA.AdamsMaxOrder")) != NULL)
349  {
350  removeParameter("LSODA.AdamsMaxOrder");
351  }
352 
353  if ((pParm = getParameter("LSODA.BDFMaxOrder")) != NULL)
354  {
355  removeParameter("LSODA.BDFMaxOrder");
356  }
357 
358  if ((pParm = getParameter("LSODA.MaxStepsInternal")) != NULL)
359  {
360  setValue("Max Internal Steps", *pParm->getValue().pUINT);
361  removeParameter("LSODA.MaxStepsInternal");
362  }
363  }
364 
365  // Check whether we have a method with "Use Default Absolute Tolerance"
366  if ((pParm = getParameter("Use Default Absolute Tolerance")) != NULL)
367  {
368  C_FLOAT64 NewValue;
369 
370  if (*pParm->getValue().pBOOL)
371  {
372  // The default
373  NewValue = 1.e-12;
374  }
375  else
376  {
377  C_FLOAT64 OldValue = *getValue("Absolute Tolerance").pUDOUBLE;
378 
379  CCopasiDataModel* pDataModel = getObjectDataModel();
380  assert(pDataModel != NULL);
381  CModel * pModel = pDataModel->getModel();
382 
383  if (pModel == NULL)
384  // The default
385  NewValue = 1.e-12;
386  else
387  {
388  const CCopasiVectorNS< CCompartment > & Compartment = pModel->getCompartments();
389  size_t i, imax;
391 
392  for (i = 0, imax = Compartment.size(); i < imax; i++)
393  if (Compartment[i]->getValue() < Volume)
394  Volume = Compartment[i]->getValue();
395 
397  // The default
398  NewValue = 1.e-12;
399  else
400  // Invert the scaling as best as we can
401  NewValue = OldValue / (Volume * pModel->getQuantity2NumberFactor());
402  }
403  }
404 
405  setValue("Absolute Tolerance", NewValue);
406  removeParameter("Use Default Absolute Tolerance");
407  }
408 
409  // These parameters are no longer supported.
410  removeParameter("Adams Max Order");
411  removeParameter("BDF Max Order");
412 
413  //These parametera are no longer defined by user
414  removeParameter("Relative Tolerance");
415  removeParameter("Absolute Tolerance");
416  removeParameter("Integrate Reduced Model");
417  removeParameter("Max Internal Steps");
418 }
419 
421 {
423  return true;
424 }
425 
426 void CTSSAMethod::integrationStep(const double & deltaT)
427 {
428  if (mData.dim == 0) //just do nothing if there are no variables
429  {
430  mTime = mTime + deltaT;
433 
434  return;
435  }
436 
437  C_FLOAT64 EndTime = mTime + deltaT;
438 
439  C_INT ITOL = 2; // mRtol scalar, mAtol vector
440  C_INT one = 1;
441  C_INT DSize = (C_INT) mDWork.size();
442  C_INT ISize = (C_INT) mIWork.size();
443 
444  mLSODA(&EvalF, // 1. evaluate F
445  &mData.dim, // 2. number of variables
446  mY, // 3. the array of current concentrations
447  &mTime, // 4. the current time
448  &EndTime, // 5. the final time
449  &ITOL, // 6. error control
450  &mRtol, // 7. relative tolerance array
451  mAtol.array(), // 8. absolute tolerance array
452  &mState, // 9. output by overshoot & interpolation
453  &mLsodaStatus, // 10. the state control variable
454  &one, // 11. further options (one)
455  mDWork.array(), // 12. the double work array
456  &DSize, // 13. the double work array size
457  mIWork.array(), // 14. the int work array
458  &ISize, // 15. the int work array size
459  NULL, // 16. evaluate J (not given)
460  &mJType); // 17. the type of jacobian calculate (2)
461 
462  // Why did we ignore this error?
463  // if (mLsodaStatus == -1) mLsodaStatus = 2;
464 
465  if ((mLsodaStatus <= 0))
466  {
468  }
469 
470  if (mLsodaStatus == 3)
471  {
472  // It is sufficient to switch to 2. Eventual state changes due to events
473  // are indicated via the method stateChanged()
474  mLsodaStatus = 2;
475  }
476 
479 
480  return;
481 }
482 /**
483 MAT_ANAL_MOD: mathematical analysis of matrices mTdInverse for post-analysis
484  */
486 {
487 
488  C_INT i, j, dim;
489 
490  dim = mData.dim;
491 
492  CVector<C_FLOAT64> denom;
493  denom.resize(dim);
494 
495  CMatrix<C_FLOAT64> Matrix_anal;
496  Matrix_anal.resize(dim, dim);
497 
498  // norm of mTd
499 
500  if (slow < dim)
501  {
502  for (j = 0; j < dim; j++)
503  denom[j] = 0;
504 
505  for (i = 0; i < dim; i++)
506  {
507  for (j = 0; j < dim; j++)
508  denom[i] = denom[i] + fabs(mTdInverse(i, j));
509  }
510 
511  for (i = 0; i < dim; i++)
512  {
513  for (j = 0; j < dim; j++)
514  mVslow(i, j) = fabs(mTdInverse(i, j)) / denom[i] * 100;
515 
516  // mVslow(i, j) = fabs(mTd(i, j)) / denom[i] * 100;
517  }
518  }
519  else
520  for (i = 0; i < dim; i++)
521  for (j = 0; j < dim; j++)
522  mVslow(i, j) = 0;
523 
524  return;
525 }
526 
527 /**
528 MAT_ANAL_METAB: mathematical analysis of matrices mTd for post-analysis
529  */
531 {
532  C_INT i, j, dim;
533 
534  dim = mData.dim;
535 
536  CVector<C_FLOAT64> denom;
537  denom.resize(dim);
538 
539  if (slow < dim)
540  {
541  for (j = 0; j < dim; j++)
542  denom[j] = 0;
543 
544  for (i = 0; i < dim; i++)
545  for (j = 0; j < dim; j++)
546  denom[i] = denom[i] + fabs(mTd(i, j));
547 
548  for (i = 0; i < dim; i++)
549  for (j = 0; j < dim; j++)
550  mVslow_metab(i, j) = fabs(mTd(i, j)) / denom[i] * 100;
551  }
552  else
553  for (i = 0; i < dim; i++)
554  for (j = 0; j < dim; j++)
555  mVslow_metab(i, j) = 0;
556 
557  return;
558 }
559 
560 /**
561 MAT_ANAL_MOD_space: mathematical analysis of matrices mTdInverse for post-analysis
562  */
563 
565 {
566  C_FLOAT64 denom, length;
567  C_INT i, j, dim;
568 
569  dim = mData.dim;
570  C_INT fast;
571  fast = dim - slow;
572 
573  CMatrix<C_FLOAT64> Matrix_anal;
574  Matrix_anal.resize(dim, dim);
575 
576  for (j = 0; j < dim; j++)
577  {
578  length = 0;
579 
580  for (i = 0; i < dim; i++)
581  length = length + mTdInverse(i, j) * mTdInverse(i, j);
582 
583  length = sqrt(length);
584  length = 1;
585 
586  for (i = 0; i < dim; i++)
587  Matrix_anal(i, j) = mTdInverse(i, j) / length;
588  }
589 
590  if ((slow < dim) & (slow > 0))
591  {
592  denom = 0.0;
593 
594  for (i = 0; i < dim; i++)
595  {
596  for (j = 0; j < slow; j++)
597  denom = denom + fabs(Matrix_anal(j, i));
598  }
599 
600  for (i = 0; i < dim; i++)
601  mVslow_space[i] = 0.0;
602 
603  for (j = 0; j < dim; j++)
604  {
605  for (i = 0; i < slow; i++)
606  mVslow_space[j] = mVslow_space[j] + fabs(Matrix_anal(i, j));
607 
608  mVslow_space[j] = (mVslow_space[j] / denom) * 100;
609  }
610  }
611  else
612  for (i = 0; i < dim; i++)
613  mVslow_space[i] = 0;
614 
615  return;
616 }
617 
618 /**
619 MAT_ANAL_fast_space: mathematical analysis of matrices mTdInverse for post-analysis
620  */
621 
623 {
624  C_FLOAT64 denom, length;
625  C_INT i, j, dim;
626 
627  dim = mData.dim;
628  C_INT fast;
629  fast = dim - slow;
630 
631  CMatrix<C_FLOAT64> Matrix_anal;
632  Matrix_anal.resize(dim, dim);
633 
634  for (j = 0; j < dim; j++)
635  {
636  length = 0;
637 
638  for (i = 0; i < dim; i++)
639  length = length + mTdInverse(i, j) * mTdInverse(i, j);
640 
641  length = sqrt(length);
642  length = 1;
643 
644  for (i = 0; i < dim; i++)
645  Matrix_anal(i, j) = mTdInverse(i, j) / length;
646  }
647 
648  if (slow < dim)
649  {
650  denom = 0.0;
651 
652  for (i = 0; i < dim; i++)
653  {
654  for (j = slow; j < dim; j++)
655  denom = denom + fabs(Matrix_anal(j, i));
656  }
657 
658  for (i = 0; i < dim; i++)
659  mVfast_space[i] = 0.0;
660 
661  for (j = 0; j < dim; j++)
662  {
663  for (i = slow; i < dim; i++)
664  mVfast_space[j] = mVfast_space[j] + fabs(Matrix_anal(i, j));
665 
666  mVfast_space[j] = (mVfast_space[j] / denom) * 100;
667  }
668  }
669  else
670  for (i = 0; i < dim; i++)
671  mVfast_space[i] = 0;
672 
673  return;
674 }
675 
676 /**
677 MAT_ANAL_fast_space: mathematical analysis of matrices mTdInverse for post-analysis
678  */
679 
681 {
682  C_FLOAT64 scalar_product, absolute_value_1;
683  C_INT i, j, k, dim;
684 
685  dim = mData.dim;
686  C_INT fast;
687  fast = dim - slow;
688 
689  C_FLOAT64 number2conc = mpModel->getNumber2QuantityFactor() / mpModel->getCompartments()[0]->getInitialValue();
690  //C_FLOAT64 number2conc = 1.;
691 
692  //this is an ugly hack that only makes sense if all metabs are in the same compartment
693  //at the moment is is the only case the algorithm deals with
694 
695  CVector<C_FLOAT64> Xconc; //current state converted to concentrations
696  Xconc.resize(dim);
697 
698  for (i = 0; i < dim; ++i)
699  Xconc[i] = mY_initial[i] * number2conc;
700 
701  CVector<C_FLOAT64> d_full;
702  d_full.resize(dim);
703  CVector<C_FLOAT64> d_transformed;
704  d_transformed.resize(dim);
705 
706  for (j = 0; j < dim; j++)
707  {
708  for (i = 0; i < slow; i++)
709  d_full[i] = 0.0;
710 
711  for (i = slow; i < dim; i++)
712  d_full[i] = mQ(j, i) * Xconc[j];
713 
714  for (i = 0; i < dim; i ++)
715  {
716  d_transformed[i] = 0.0;
717 
718  for (k = 0; k < dim; k++)
719  d_transformed[i] += mQ(i, k) * d_full[k];
720  }
721 
722  scalar_product = d_transformed[j];
723  absolute_value_1 = 0.0;
724 
725  for (i = 0; i < dim; i++)
726  absolute_value_1 += d_transformed[i] * d_transformed[i];
727 
728  absolute_value_1 = sqrt(absolute_value_1);
729 
730  if (absolute_value_1 * Xconc[j] > 0.0)
731  scalar_product = scalar_product / absolute_value_1;
732  else
733  scalar_product = -2.0;
734 
735  mVfast_space[j] = scalar_product;
736  }
737 
738  return;
739 }
740 
741 double CTSSAMethod::orthog(C_INT & number1, C_INT & number2)
742 {
743  C_FLOAT64 product = 0;
744  C_INT k, dim;
745 
746  dim = mData.dim;
747 
748  for (k = 0; k < dim; k++)
749  product = product + mTdInverse(k, number1) * mTdInverse(k, number2);
750 
751  return product;
752 }
753 
754 /**
755  SCHUR: Schur Decomposition of Jacobian (reordered).
756  Output: mQ - transformation matrix
757  mR - block upper triangular matrix (with ordered eigenvalues)
758  */
759 
761 {
762 
763  /* int dgees_(char *jobvs, char *sort, L_fp select, integer *n,
764  * doublereal *a, integer *lda, integer *sdim, doublereal *wr,
765  * doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work,
766  * integer *lwork, logical *bwork, integer *info)
767  *
768  * -- LAPACK driver routine (version 3.0) --
769  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
770  * Courant Institute, Argonne National Lab, and Rice University
771  * June 30, 1999
772  *
773  *
774  * Purpose
775  * =======
776  *
777  * DGEES computes for an N-by-N real nonsymmetric matrix A, the
778  * eigenvalues, the real Schur form T, and, optionally, the matrix of
779  * Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T).
780  *
781  * Optionally, it also orders the eigenvalues on the diagonal of the
782  * real Schur form so that selected eigenvalues are at the top left.
783  * The leading columns of Z then form an orthonormal basis for the
784  * invariant subspace corresponding to the selected eigenvalues.
785  *
786  * A matrix is in real Schur form if it is upper quasi-triangular with
787  * 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the
788  * form
789  * [ a b ]
790  * [ c a ]
791  *
792  * where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
793  *
794  * Arguments
795  * =========
796  *
797  * JOBVS (input) CHARACTER*1
798  * = 'N': Schur vectors are not computed;
799  * = 'V': Schur vectors are computed.
800  *
801  * SORT (input) CHARACTER*1
802  * Specifies whether or not to order the eigenvalues on the
803  * diagonal of the Schur form.
804  * = 'N': Eigenvalues are not ordered;
805  * = 'S': Eigenvalues are ordered (see SELECT).
806  *
807  * SELECT (input) LOGICAL FUNCTION of two DOUBLE PRECISION arguments
808  * SELECT must be declared EXTERNAL in the calling subroutine.
809  * If SORT = 'S', SELECT is used to select eigenvalues to sort
810  * to the top left of the Schur form.
811  * If SORT = 'N', SELECT is not referenced.
812  * An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
813  * SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
814  * conjugate pair of eigenvalues is selected, then both complex
815  * eigenvalues are selected.
816  * Note that a selected complex eigenvalue may no longer
817  * satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
818  * ordering may change the value of complex eigenvalues
819  * (especially if the eigenvalue is ill-conditioned); in this
820  * case INFO is set to N+2 (see INFO below).
821  *
822  * N (input) INTEGER
823  * The order of the matrix A. N >= 0.
824  *
825  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
826  * On entry, the N-by-N matrix A.
827  * On exit, A has been overwritten by its real Schur form T.
828  *
829  * LDA (input) INTEGER
830  * The leading dimension of the array A. LDA >= max(1,N).
831  *
832  * SDIM (output) INTEGER
833  * If SORT = 'N', SDIM = 0.
834  * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
835  * for which SELECT is true. (Complex conjugate
836  * pairs for which SELECT is true for either
837  * eigenvalue count as 2.)
838  *
839  * WR (output) DOUBLE PRECISION array, dimension (N)
840  * WI (output) DOUBLE PRECISION array, dimension (N)
841  * WR and WI contain the real and imaginary parts,
842  * respectively, of the computed eigenvalues in the same order
843  * that they appear on the diagonal of the output Schur form T.
844  * Complex conjugate pairs of eigenvalues will appear
845  * consecutively with the eigenvalue having the positive
846  * imaginary part first.
847  *
848  * VS (output) DOUBLE PRECISION array, dimension (LDVS,N)
849  * If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
850  * vectors.
851  * If JOBVS = 'N', VS is not referenced.
852  *
853  * LDVS (input) INTEGER
854  * The leading dimension of the array VS. LDVS >= 1; if
855  * JOBVS = 'V', LDVS >= N.
856  *
857  * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
858  * On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
859  *
860  * LWORK (input) INTEGER
861  * The dimension of the array WORK. LWORK >= max(1,3*N).
862  * For good performance, LWORK must generally be larger.
863  *
864  * If LWORK = -1, then a workspace query is assumed; the routine
865  * only calculates the optimal size of the WORK array, returns
866  * this value as the first entry of the WORK array, and no error
867  * message related to LWORK is issued by XERBLA.
868  *
869  * BWORK (workspace) LOGICAL array, dimension (N)
870  * Not referenced if SORT = 'N'.
871  *
872  * INFO (output) INTEGER
873  * = 0: successful exit
874  * < 0: if INFO = -i, the i-th argument had an illegal value.
875  * > 0: if INFO = i, and i is
876  * <= N: the QR algorithm failed to compute all the
877  * eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
878  * contain those eigenvalues which have converged; if
879  * JOBVS = 'V', VS contains the matrix which reduces A
880  * to its partially converged Schur form.
881  * = N+1: the eigenvalues could not be reordered because some
882  * eigenvalues were too close to separate (the problem
883  * is very ill-conditioned);
884  * = N+2: after reordering, roundoff changed values of some
885  * complex eigenvalues so that leading eigenvalues in
886  * the Schur form no longer satisfy SELECT=.TRUE. This
887  * could also be caused by underflow due to scaling.
888  *
889  */
890 
891  char V = 'V';
892  char N = 'N';
893  // TO REMOVE : L_fp select;
894  C_INT dim = mData.dim;
895  C_INT SDIM = 0;
896 
898  R.resize(dim * dim);
899  C_INT i, j;
900 
901  for (i = 0; i < dim; i++)
902  for (j = 0; j < dim; j++)
903  R[j + dim * i] = mJacobian_initial(j, i);
904 
905  CVector<C_FLOAT64> eval_r;
906  CVector<C_FLOAT64> eval_i;
907  eval_r.resize(dim);
908  eval_i.resize(dim);
909 
911  Q.resize(dim * dim);
912 
913  C_INT lwork = 10 * dim;
915  work.resize(10 * dim);
916 
917  CVector< C_LOGICAL > Bwork;
918  Bwork.resize(dim);
919 
920  dgees_(&V, &N, NULL, &dim, R.array(), &dim, &SDIM, eval_r.array(), eval_i.array(), Q.array(), &dim, work.array(), & lwork, Bwork.array(), &info);
921 
922  if (info)
923  {
924  return;
925  }
926 
927  /* Sorting of Schurmatix */
928 
929  CVector<C_FLOAT64> eval_reor;
930  eval_reor.resize(dim);
931 
932  CVector<C_INT> index;
933  index.resize(dim);
934 
935  for (i = 0; i < dim; i++)
936  eval_reor[i] = eval_r [i];
937 
938  for (i = 0; i < dim; i++) index[i] = 0;
939 
940  map_index(eval_reor.array(), index.array(), dim);
941 
942  CVector<C_INT> nid;
943  CVector<C_INT> pid;
944  nid.resize(dim);
945  pid.resize(dim);
946 
947  if (dim > 2)
948  {
949  update_nid(index.array(), nid.array(), dim);
950  update_pid(index.array(), pid.array(), dim);
951  }
952  else
953  {
954  for (i = 0; i < dim; i++)
955  {
956  nid[i] = 0;
957  pid[i] = 0;
958  }
959  }
960 
961  bool changed = true;
962  C_INT count;
963 
964  while (changed == true)
965  {
966  changed = false;
967  count = 0;
968 
969  while (count < dim - 1)
970  {
971  C_INT first;
972  C_INT second;
973 
974  bool diagorder = false;
975 
976  CCopasiMethod::SubType subType;
977 
978  subType = mData.pMethod->getSubType();
979 
980  switch (subType)
981  {
982  case tssILDM:
983  diagorder = (index[count + 1] < index[count]);
984  break;
985 
986  case tssILDMModified:
987  diagorder = (index[count + 1] < index[count]);
988  break;
989 
990  case tssCSP:
991  diagorder = (index[count + 1] > index[count]);
992  break;
993 
994  default:
995  fatalError();
996  break;
997  }
998 
999  //ILDM : if (index[count + 1] < index[count])
1000  //CSP: if (index[count + 1] > index[count])
1001  if (diagorder)
1002  {
1003  changed = true;
1004  first = count + 2;
1005  second = count + 1;
1006 
1007  CVector<C_FLOAT64> j_diag;
1008  j_diag.resize(dim);
1009 
1010  for (i = 0; i < dim; i++)
1011  j_diag[i] = R.array()[i * dim + i];
1012 
1013  map_index(j_diag.array(), index.array(), dim);
1014  update_nid(index.array(), nid.array(), dim);
1015  update_pid(index.array(), pid.array(), dim);
1016 
1017  /* int dtrexc_(char *compq, integer *n, doublereal *t, integer *
1018  * ldt, doublereal *q, integer *ldq, integer *ifst, integer *ilst,
1019  * doublereal *work, integer *info)
1020  * -- LAPACK routine (version 3.0) --
1021  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1022  * Courant Institute, Argonne National Lab, and Rice University
1023  * March 31, 1993
1024  *
1025  *
1026  * Purpose
1027  * =======
1028  *
1029  * DTREXC reorders the real Schur factorization of a real matrix
1030  * A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
1031  * moved to row ILST.
1032  *
1033  * The real Schur form T is reordered by an orthogonal similarity
1034  * transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
1035  * is updated by postmultiplying it with Z.
1036  *
1037  * T must be in Schur canonical form (as returned by DHSEQR), that is,
1038  * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
1039  * 2-by-2 diagonal block has its diagonal elements equal and its
1040  * off-diagonal elements of opposite sign.
1041  *
1042  * Arguments
1043  * =========
1044  *
1045  * COMPQ (input) CHARACTER*1
1046  * = 'V': update the matrix Q of Schur vectors;
1047  * = 'N': do not update Q.
1048  *
1049  * N (input) INTEGER
1050  * The order of the matrix T. N >= 0.
1051  *
1052  * T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
1053  * On entry, the upper quasi-triangular matrix T, in Schur
1054  * Schur canonical form.
1055  * On exit, the reordered upper quasi-triangular matrix, again
1056  * in Schur canonical form.
1057  *
1058  * LDT (input) INTEGER
1059  * The leading dimension of the array T. LDT >= max(1,N).
1060  *
1061  * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
1062  * On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
1063  * On exit, if COMPQ = 'V', Q has been postmultiplied by the
1064  * orthogonal transformation matrix Z which reorders T.
1065  * If COMPQ = 'N', Q is not referenced.
1066  *
1067  * LDQ (input) INTEGER
1068  * The leading dimension of the array Q. LDQ >= max(1,N).
1069  *
1070  * IFST (input/output) INTEGER
1071  * ILST (input/output) INTEGER
1072  * Specify the reordering of the diagonal blocks of T.
1073  * The block with row index IFST is moved to row ILST, by a
1074  * sequence of transpositions between adjacent blocks.
1075  * On exit, if IFST pointed on entry to the second row of a
1076  * 2-by-2 block, it is changed to point to the first row; ILST
1077  * always points to the first row of the block in its final
1078  * position (which may differ from its input value by +1 or -1).
1079  * 1 <= IFST <= N; 1 <= ILST <= N.
1080  *
1081  * WORK (workspace) DOUBLE PRECISION array, dimension (N)
1082  *
1083  * INFO (output) INTEGER
1084  * = 0: successful exit
1085  * < 0: if INFO = -i, the i-th argument had an illegal value
1086  * = 1: two adjacent blocks were too close to swap (the problem
1087  * is very ill-conditioned); T may have been partially
1088  * reordered, and ILST points to the first row of the
1089  * current position of the block being moved.
1090  */
1091 
1092  CVector< C_FLOAT64 > work1;
1093  work1.resize(dim);
1094 
1095  dtrexc_(&V, &dim, R.array(), &dim, Q.array(), &dim, &first, &second, work1.array(), &info);
1096 
1097  if (info)
1098  {
1099  /* TODO */
1100  return;
1101  }
1102 
1103  C_INT tmp;
1104 
1105  if (nid[count] == 0)
1106  {
1107  if (pid[count] == 0)
1108  {
1109  /* both are not degenerate */
1110 
1111  tmp = index[count];
1112  index[count] = index[count + 1];
1113  index[count + 1] = tmp;
1114 
1115  update_nid(index.array(), nid.array(), dim);
1116  update_pid(index.array(), pid.array(), dim);
1117 
1118  count = count + 1;
1119  }
1120  else
1121  {
1122  /* the first is degenerate, the second not */
1123 
1124  tmp = index[count - 1];
1125  index[count - 1] = index[count + 1];
1126  index[count + 1] = tmp;
1127 
1128  update_nid(index.array(), nid.array(), dim);
1129  update_pid(index.array(), pid.array(), dim);
1130 
1131  count = count + 1;
1132  }
1133  }
1134  else
1135  {
1136  if (pid[count] == 0)
1137  {
1138  /* the next is deg, prev is not */
1139 
1140  tmp = index[count];
1141  index[count] = index[count + 2];
1142  index[count + 2] = tmp;
1143 
1144  update_nid(index.array(), nid.array(), dim);
1145  update_pid(index.array(), pid.array(), dim);
1146 
1147  count = count + 1;
1148  }
1149  else
1150  {
1151  /* both are deg. */
1152 
1153  tmp = index[count];
1154  index[count - 1] = index[count + 1];
1155  index[count] = index[count + 1];
1156  index[count + 1] = tmp;
1157  index[count + 2] = tmp;
1158 
1159  update_nid(index.array(), nid.array(), dim);
1160  update_pid(index.array(), pid.array(), dim);
1161 
1162  count = count + 1;
1163  }
1164  }
1165  }
1166  else
1167  count = count + 1;
1168  }
1169  }
1170 
1171  for (i = 0; i < dim; i++)
1172  for (j = 0; j < dim; j++)
1173  {
1174  mQ(j, i) = Q.array()[j + dim * i];
1175  mR(j, i) = R.array()[j + dim * i];
1176  }
1177 
1178  return;
1179 }
1180 
1181 /**
1182  SCHUR_desc: Schur Decomposition of Jacobian (reordered).
1183  Output: mQ_desc - transformation matrix
1184  mR_desc - block upper triangular matrix (with ordered eigenvalues)
1185  */
1186 
1188 {
1189 
1190  /* int dgees_(char *jobvs, char *sort, L_fp select, integer *n,
1191  * doublereal *a, integer *lda, integer *sdim, doublereal *wr,
1192  * doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work,
1193  * integer *lwork, logical *bwork, integer *info)
1194  *
1195  * -- LAPACK driver routine (version 3.0) --
1196  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1197  * Courant Institute, Argonne National Lab, and Rice University
1198  * June 30, 1999
1199  *
1200  *
1201  * Purpose
1202  * =======
1203  *
1204  * DGEES computes for an N-by-N real nonsymmetric matrix A, the
1205  * eigenvalues, the real Schur form T, and, optionally, the matrix of
1206  * Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T).
1207  *
1208  * Optionally, it also orders the eigenvalues on the diagonal of the
1209  * real Schur form so that selected eigenvalues are at the top left.
1210  * The leading columns of Z then form an orthonormal basis for the
1211  * invariant subspace corresponding to the selected eigenvalues.
1212  *
1213  * A matrix is in real Schur form if it is upper quasi-triangular with
1214  * 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the
1215  * form
1216  * [ a b ]
1217  * [ c a ]
1218  *
1219  * where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
1220  *
1221  * Arguments
1222  * =========
1223  *
1224  * JOBVS (input) CHARACTER*1
1225  * = 'N': Schur vectors are not computed;
1226  * = 'V': Schur vectors are computed.
1227  *
1228  * SORT (input) CHARACTER*1
1229  * Specifies whether or not to order the eigenvalues on the
1230  * diagonal of the Schur form.
1231  * = 'N': Eigenvalues are not ordered;
1232  * = 'S': Eigenvalues are ordered (see SELECT).
1233  *
1234  * SELECT (input) LOGICAL FUNCTION of two DOUBLE PRECISION arguments
1235  * SELECT must be declared EXTERNAL in the calling subroutine.
1236  * If SORT = 'S', SELECT is used to select eigenvalues to sort
1237  * to the top left of the Schur form.
1238  * If SORT = 'N', SELECT is not referenced.
1239  * An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
1240  * SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
1241  * conjugate pair of eigenvalues is selected, then both complex
1242  * eigenvalues are selected.
1243  * Note that a selected complex eigenvalue may no longer
1244  * satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
1245  * ordering may change the value of complex eigenvalues
1246  * (especially if the eigenvalue is ill-conditioned); in this
1247  * case INFO is set to N+2 (see INFO below).
1248  *
1249  * N (input) INTEGER
1250  * The order of the matrix A. N >= 0.
1251  *
1252  * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
1253  * On entry, the N-by-N matrix A.
1254  * On exit, A has been overwritten by its real Schur form T.
1255  *
1256  * LDA (input) INTEGER
1257  * The leading dimension of the array A. LDA >= max(1,N).
1258  *
1259  * SDIM (output) INTEGER
1260  * If SORT = 'N', SDIM = 0.
1261  * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
1262  * for which SELECT is true. (Complex conjugate
1263  * pairs for which SELECT is true for either
1264  * eigenvalue count as 2.)
1265  *
1266  * WR (output) DOUBLE PRECISION array, dimension (N)
1267  * WI (output) DOUBLE PRECISION array, dimension (N)
1268  * WR and WI contain the real and imaginary parts,
1269  * respectively, of the computed eigenvalues in the same order
1270  * that they appear on the diagonal of the output Schur form T.
1271  * Complex conjugate pairs of eigenvalues will appear
1272  * consecutively with the eigenvalue having the positive
1273  * imaginary part first.
1274  *
1275  * VS (output) DOUBLE PRECISION array, dimension (LDVS,N)
1276  * If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
1277  * vectors.
1278  * If JOBVS = 'N', VS is not referenced.
1279  *
1280  * LDVS (input) INTEGER
1281  * The leading dimension of the array VS. LDVS >= 1; if
1282  * JOBVS = 'V', LDVS >= N.
1283  *
1284  * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
1285  * On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
1286  *
1287  * LWORK (input) INTEGER
1288  * The dimension of the array WORK. LWORK >= max(1,3*N).
1289  * For good performance, LWORK must generally be larger.
1290  *
1291  * If LWORK = -1, then a workspace query is assumed; the routine
1292  * only calculates the optimal size of the WORK array, returns
1293  * this value as the first entry of the WORK array, and no error
1294  * message related to LWORK is issued by XERBLA.
1295  *
1296  * BWORK (workspace) LOGICAL array, dimension (N)
1297  * Not referenced if SORT = 'N'.
1298  *
1299  * INFO (output) INTEGER
1300  * = 0: successful exit
1301  * < 0: if INFO = -i, the i-th argument had an illegal value.
1302  * > 0: if INFO = i, and i is
1303  * <= N: the QR algorithm failed to compute all the
1304  * eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
1305  * contain those eigenvalues which have converged; if
1306  * JOBVS = 'V', VS contains the matrix which reduces A
1307  * to its partially converged Schur form.
1308  * = N+1: the eigenvalues could not be reordered because some
1309  * eigenvalues were too close to separate (the problem
1310  * is very ill-conditioned);
1311  * = N+2: after reordering, roundoff changed values of some
1312  * complex eigenvalues so that leading eigenvalues in
1313  * the Schur form no longer satisfy SELECT=.TRUE. This
1314  * could also be caused by underflow due to scaling.
1315  *
1316  */
1317 
1318  char V = 'V';
1319  char N = 'N';
1320  // TO REMOVE : L_fp select;
1321  C_INT dim = mData.dim;
1322  C_INT SDIM = 0;
1323 
1325  R.resize(dim * dim);
1326  C_INT i, j;
1327 
1328  for (i = 0; i < dim; i++)
1329  for (j = 0; j < dim; j++)
1330  R[j + dim * i] = mJacobian_initial(j, i);
1331 
1332  CVector<C_FLOAT64> eval_r;
1333  CVector<C_FLOAT64> eval_i;
1334  eval_r.resize(dim);
1335  eval_i.resize(dim);
1336 
1338  Q.resize(dim * dim);
1339 
1340  C_INT lwork = 10 * dim;
1341  CVector< C_FLOAT64 > work;
1342  work.resize(10 * dim);
1343 
1344  CVector< C_LOGICAL > Bwork;
1345  Bwork.resize(dim);
1346 
1347  dgees_(&V, &N, NULL, &dim, R.array(), &dim, &SDIM, eval_r.array(), eval_i.array(), Q.array(), &dim, work.array(), & lwork, Bwork.array(), &info);
1348 
1349  if (info)
1350  {
1351  return;
1352  }
1353 
1354  /* Sorting of Schurmatix */
1355 
1356  CVector<C_FLOAT64> eval_reor;
1357  eval_reor.resize(dim);
1358 
1359  CVector<C_INT> index;
1360  index.resize(dim);
1361 
1362  for (i = 0; i < dim; i++)
1363  eval_reor[i] = eval_r [i];
1364 
1365  for (i = 0; i < dim; i++) index[i] = 0;
1366 
1367  map_index_desc(eval_reor.array(), index.array(), dim);
1368 
1369  CVector<C_INT> nid;
1370  CVector<C_INT> pid;
1371  nid.resize(dim);
1372  pid.resize(dim);
1373 
1374  if (dim > 2)
1375  {
1376  update_nid(index.array(), nid.array(), dim);
1377  update_pid(index.array(), pid.array(), dim);
1378  }
1379  else
1380  {
1381  for (i = 0; i < dim; i++)
1382  {
1383  nid[i] = 0;
1384  pid[i] = 0;
1385  }
1386  }
1387 
1388  bool changed = true;
1389  C_INT count;
1390 
1391  while (changed == true)
1392  {
1393  changed = false;
1394  count = 0;
1395 
1396  while (count < dim - 1)
1397  {
1398  C_INT first;
1399  C_INT second;
1400 
1401  if (index[count + 1] < index[count])
1402  {
1403  changed = true;
1404  first = count + 2;
1405  second = count + 1;
1406 
1407  CVector<C_FLOAT64> j_diag;
1408  j_diag.resize(dim);
1409 
1410  for (i = 0; i < dim; i++)
1411  j_diag[i] = R.array()[i * dim + i];
1412 
1413  map_index_desc(j_diag.array(), index.array(), dim);
1414  update_nid(index.array(), nid.array(), dim);
1415  update_pid(index.array(), pid.array(), dim);
1416 
1417  /* int dtrexc_(char *compq, integer *n, doublereal *t, integer *
1418  * ldt, doublereal *q, integer *ldq, integer *ifst, integer *ilst,
1419  * doublereal *work, integer *info)
1420  * -- LAPACK routine (version 3.0) --
1421  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1422  * Courant Institute, Argonne National Lab, and Rice University
1423  * March 31, 1993
1424  *
1425  *
1426  * Purpose
1427  * =======
1428  *
1429  * DTREXC reorders the real Schur factorization of a real matrix
1430  * A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
1431  * moved to row ILST.
1432  *
1433  * The real Schur form T is reordered by an orthogonal similarity
1434  * transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
1435  * is updated by postmultiplying it with Z.
1436  *
1437  * T must be in Schur canonical form (as returned by DHSEQR), that is,
1438  * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
1439  * 2-by-2 diagonal block has its diagonal elements equal and its
1440  * off-diagonal elements of opposite sign.
1441  *
1442  * Arguments
1443  * =========
1444  *
1445  * COMPQ (input) CHARACTER*1
1446  * = 'V': update the matrix Q of Schur vectors;
1447  * = 'N': do not update Q.
1448  *
1449  * N (input) INTEGER
1450  * The order of the matrix T. N >= 0.
1451  *
1452  * T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
1453  * On entry, the upper quasi-triangular matrix T, in Schur
1454  * Schur canonical form.
1455  * On exit, the reordered upper quasi-triangular matrix, again
1456  * in Schur canonical form.
1457  *
1458  * LDT (input) INTEGER
1459  * The leading dimension of the array T. LDT >= max(1,N).
1460  *
1461  * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
1462  * On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
1463  * On exit, if COMPQ = 'V', Q has been postmultiplied by the
1464  * orthogonal transformation matrix Z which reorders T.
1465  * If COMPQ = 'N', Q is not referenced.
1466  *
1467  * LDQ (input) INTEGER
1468  * The leading dimension of the array Q. LDQ >= max(1,N).
1469  *
1470  * IFST (input/output) INTEGER
1471  * ILST (input/output) INTEGER
1472  * Specify the reordering of the diagonal blocks of T.
1473  * The block with row index IFST is moved to row ILST, by a
1474  * sequence of transpositions between adjacent blocks.
1475  * On exit, if IFST pointed on entry to the second row of a
1476  * 2-by-2 block, it is changed to point to the first row; ILST
1477  * always points to the first row of the block in its final
1478  * position (which may differ from its input value by +1 or -1).
1479  * 1 <= IFST <= N; 1 <= ILST <= N.
1480  *
1481  * WORK (workspace) DOUBLE PRECISION array, dimension (N)
1482  *
1483  * INFO (output) INTEGER
1484  * = 0: successful exit
1485  * < 0: if INFO = -i, the i-th argument had an illegal value
1486  * = 1: two adjacent blocks were too close to swap (the problem
1487  * is very ill-conditioned); T may have been partially
1488  * reordered, and ILST points to the first row of the
1489  * current position of the block being moved.
1490  */
1491 
1492  CVector< C_FLOAT64 > work1;
1493  work1.resize(dim);
1494 
1495  dtrexc_(&V, &dim, R.array(), &dim, Q.array(), &dim, &first, &second, work1.array(), &info);
1496 
1497  if (info)
1498  {
1499  /* TODO */
1500  return;
1501  }
1502 
1503  C_INT tmp;
1504 
1505  if (nid[count] == 0)
1506  {
1507  if (pid[count] == 0)
1508  {
1509  /* both are not degenerate */
1510 
1511  tmp = index[count];
1512  index[count] = index[count + 1];
1513  index[count + 1] = tmp;
1514 
1515  update_nid(index.array(), nid.array(), dim);
1516  update_pid(index.array(), pid.array(), dim);
1517 
1518  count = count + 1;
1519  }
1520  else
1521  {
1522  /* the first is degenerate, the second not */
1523 
1524  tmp = index[count - 1];
1525  index[count - 1] = index[count + 1];
1526  index[count + 1] = tmp;
1527 
1528  update_nid(index.array(), nid.array(), dim);
1529  update_pid(index.array(), pid.array(), dim);
1530 
1531  count = count + 1;
1532  }
1533  }
1534  else
1535  {
1536  if (pid[count] == 0)
1537  {
1538  /* the next is deg, prev is not */
1539 
1540  tmp = index[count];
1541  index[count] = index[count + 2];
1542  index[count + 2] = tmp;
1543 
1544  update_nid(index.array(), nid.array(), dim);
1545  update_pid(index.array(), pid.array(), dim);
1546 
1547  count = count + 1;
1548  }
1549  else
1550  {
1551  /* both are deg. */
1552 
1553  tmp = index[count];
1554  index[count - 1] = index[count + 1];
1555  index[count] = index[count + 1];
1556  index[count + 1] = tmp;
1557  index[count + 2] = tmp;
1558 
1559  update_nid(index.array(), nid.array(), dim);
1560  update_pid(index.array(), pid.array(), dim);
1561 
1562  count = count + 1;
1563  }
1564  }
1565  }
1566  else
1567  count = count + 1;
1568  }
1569  }
1570 
1571  for (i = 0; i < dim; i++)
1572  for (j = 0; j < dim; j++)
1573  {
1574  mQ_desc(j, i) = Q.array()[j + dim * i];
1575  mR_desc(j, i) = R.array()[j + dim * i];
1576  }
1577 
1578  return;
1579 }
1580 
1581 /**
1582 SYLVESTER:
1583 Solution of Sylvester equation for given slow, mQ,mR
1584 Output: mTd, mTdinverse, mQz (is used later for Newton iterations)
1585  */
1586 
1588 {
1589  char N1 = 'N';
1590  char N2 = 'N';
1591  C_INT isgn = -1;
1592 
1593  C_INT dim = mData.dim;
1594 
1595  C_INT fast = dim - slow;
1596  // C_INT info;
1597 
1598  C_INT i, j, k;
1599 
1600  C_FLOAT64 scale = -1;
1601 
1602  CVector<C_FLOAT64> st_slow;
1603  st_slow.resize(slow * slow);
1604 
1605  CVector<C_FLOAT64> st_fast;
1606  st_fast.resize(fast * fast);
1607 
1608  CVector<C_FLOAT64> st_coup;
1609  st_coup.resize(slow * fast);
1610 
1611  CMatrix<C_FLOAT64> S_coup;
1612  S_coup.resize(slow, fast);
1613 
1614  for (i = 0; i < slow; i++)
1615  for (j = 0; j < slow; j++)
1616  st_slow[j + slow * i] = mR(j, i);
1617 
1618  for (i = 0; i < fast; i++)
1619  for (j = 0; j < fast; j++)
1620  st_fast[j + fast * i] = mR(j + slow, i + slow);
1621 
1622  for (i = 0; i < slow; i++)
1623  for (j = 0; j < fast; j++)
1624  S_coup(i, j) = mR(i, j + slow);
1625 
1626  for (j = 0; j < fast; j++)
1627  for (i = 0; i < slow; i++)
1628  st_coup[i + slow * j] = -S_coup(i, j);
1629 
1630  /* int dtrsyl_(char *trana, char *tranb, integer *isgn, integer
1631  * *m, integer *n, doublereal *a, integer *lda, doublereal *b, integer *
1632  * ldb, doublereal *c__, integer *ldc, doublereal *scale, integer *info)
1633  *
1634  * -- LAPACK routine (version 3.0) --
1635  * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1636  * Courant Institute, Argonne National Lab, and Rice University
1637  * March 31, 1993
1638  *
1639  *
1640  * Purpose
1641  * =======
1642  *
1643  * DTRSYL solves the real Sylvester matrix equation:
1644  *
1645  * op(A)*X + X*op(B) = scale*C or
1646  * op(A)*X - X*op(B) = scale*C,
1647  *
1648  * where op(A) = A or A**T, and A and B are both upper quasi-
1649  * triangular. A is M-by-M and B is N-by-N; the right hand side C and
1650  * the solution X are M-by-N; and scale is an output scale factor, set
1651  * <= 1 to avoid overflow in X.
1652  *
1653  * A and B must be in Schur canonical form (as returned by DHSEQR), that
1654  * is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
1655  * each 2-by-2 diagonal block has its diagonal elements equal and its
1656  * off-diagonal elements of opposite sign.
1657  *
1658  * Arguments
1659  * =========
1660  *
1661  * TRANA (input) CHARACTER*1
1662  * Specifies the option op(A):
1663  * = 'N': op(A) = A (No transpose)
1664  * = 'T': op(A) = A**T (Transpose)
1665  * = 'C': op(A) = A**H (Conjugate transpose = Transpose)
1666  *
1667  * TRANB (input) CHARACTER*1
1668  * Specifies the option op(B):
1669  * = 'N': op(B) = B (No transpose)
1670  * = 'T': op(B) = B**T (Transpose)
1671  * = 'C': op(B) = B**H (Conjugate transpose = Transpose)
1672  *
1673  * ISGN (input) INTEGER
1674  * Specifies the sign in the equation:
1675  * = +1: solve op(A)*X + X*op(B) = scale*C
1676  * = -1: solve op(A)*X - X*op(B) = scale*C
1677  *
1678  * M (input) INTEGER
1679  * The order of the matrix A, and the number of rows in the
1680  * matrices X and C. M >= 0.
1681  *
1682  * N (input) INTEGER
1683  * The order of the matrix B, and the number of columns in the
1684  * matrices X and C. N >= 0.
1685  *
1686  * A (input) DOUBLE PRECISION array, dimension (LDA,M)
1687  * The upper quasi-triangular matrix A, in Schur canonical form.
1688  *
1689  * LDA (input) INTEGER
1690  * The leading dimension of the array A. LDA >= max(1,M).
1691  *
1692  * B (input) DOUBLE PRECISION array, dimension (LDB,N)
1693  * The upper quasi-triangular matrix B, in Schur canonical form.
1694  *
1695  * LDB (input) INTEGER
1696  * The leading dimension of the array B. LDB >= max(1,N).
1697  *
1698  * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
1699  * On entry, the M-by-N right hand side matrix C.
1700  * On exit, C is overwritten by the solution matrix X.
1701  *
1702  * LDC (input) INTEGER
1703  * The leading dimension of the array C. LDC >= max(1,M)
1704  *
1705  * SCALE (output) DOUBLE PRECISION
1706  * The scale factor, scale, set <= 1 to avoid overflow in X.
1707  *
1708  * INFO (output) INTEGER
1709  * = 0: successful exit
1710  * < 0: if INFO = -i, the i-th argument had an illegal value
1711  * = 1: A and B have common or very close eigenvalues; perturbed
1712  * values were used to solve the equation (but the matrices
1713  * A and B are unchanged).
1714  */
1715 
1716  dtrsyl_(&N1, &N2, &isgn, &slow, &fast, st_slow.array(), &slow, st_fast.array(), &fast, st_coup.array(), &slow, &scale, &info);
1717 
1718  /* if (info) TODO*/
1719  if (info)
1720  {
1721  return;
1722  }
1723 
1724  CMatrix<C_FLOAT64> Cmat;
1725  Cmat.resize(dim, dim);
1726 
1728  C.resize(dim, dim);
1729 
1730  for (i = 0; i < dim; i++)
1731  for (j = 0; j < dim; j++)
1732  Cmat(i, j) = 0.;
1733 
1734  for (j = 0; j < fast; j++)
1735  for (i = 0; i < slow; i++)
1736  Cmat(i, j + slow) = st_coup[i + slow * j];
1737 
1738  for (i = 0; i < dim; i++)
1739  for (j = 0; j < dim; j++)
1740  if (i == j)
1741  C(i, j) = 1. + Cmat(i, j);
1742  else
1743  C(i, j) = Cmat(i, j);
1744 
1745  for (i = 0; i < dim; i++)
1746  {
1747  for (j = 0; j < dim; j++)
1748  {
1749  mTd(i, j) = 0.;
1750 
1751  for (k = 0; k < dim; k++)
1752  mTd(i, j) = mTd(i, j) + mQ(i, k) * C(k, j);
1753  }
1754  }
1755 
1756  for (i = 0; i < dim; i++)
1757  for (j = 0; j < dim; j++)
1758  if (i == j)
1759  C(i, j) = 1. - Cmat(i, j);
1760  else
1761  C(i, j) = - Cmat(i, j);
1762 
1763  CMatrix<C_FLOAT64> transp_Qt;
1764  transp_Qt.resize(dim, dim);
1765 
1766  for (i = 0; i < dim; i++)
1767  for (j = 0; j < dim; j++)
1768  transp_Qt(i, j) = mQ(j, i);
1769 
1770  for (i = 0; i < dim; i++)
1771  {
1772  for (j = 0; j < dim; j++)
1773  {
1774  mTdInverse(i, j) = 0.0;
1775 
1776  for (k = 0; k < dim; k++)
1777  mTdInverse(i, j) = mTdInverse(i, j) + C(i, k) * transp_Qt(k, j);
1778  }
1779  }
1780 
1782  E.resize(dim, dim);
1783 
1785  S.resize(dim, dim);
1786 
1787  for (i = 0; i < dim; i++)
1788  {
1789  for (j = 0; j < dim; j++)
1790  {
1791  E(i, j) = 0.;
1792 
1793  for (k = 0; k < dim; k++)
1794  E(i, j) = E(i, j) + mJacobian_initial(i, k) * mTd(k, j);
1795  }
1796  }
1797 
1798  for (i = 0; i < dim; i++)
1799  {
1800  for (j = 0; j < dim; j++)
1801  {
1802  S(i, j) = 0.;
1803 
1804  for (k = 0; k < dim; k++)
1805  S(i, j) = S(i, j) + mTdInverse(i, k) * E(k, j);
1806  }
1807  }
1808 
1809  C_INT flag_sylvester;
1810 
1811  flag_sylvester = 1;
1812 
1813  for (i = 0; i < dim; i++)
1814  for (j = 0; j < dim; j++)
1815  mQz(i, j) = 0.;
1816 
1817  for (i = 0; i < fast; i++)
1818  for (j = 0; j < fast; j++)
1819  mQz(i, j) = S(i + slow, j + slow);
1820 
1821  return;
1822 }
1823 
1825 {
1826  size_t i, imax;
1827  size_t indep;
1828 
1829  // indep = mpModel->getNumIndependentMetabs();
1831 
1832  CVector<C_FLOAT64> tmp;
1833  tmp.resize(indep);
1834 
1835  /* make copy of the current state concentrations */
1836  for (i = 0, imax = indep; i < imax; i++)
1837  tmp[i] = mpModel->getMetabolitesX()[i]->getValue();
1838 
1839  C_FLOAT64 conc2number = mpModel->getQuantity2NumberFactor() * mpModel->getCompartments()[0]->getInitialValue();
1840  //C_FLOAT64 conc2number = 1.;
1841 
1842  /* write new concentrations in the current state */
1843  for (i = 0, imax = indep; i < imax; i++)
1844  mpModel->getMetabolitesX()[i]->setValue(X1[i]*conc2number);
1845 
1846  //mpState->setUpdateDependentRequired(true);
1848  // TO REMOVE: mpModel->applyAssignments();
1850 
1851  C_FLOAT64 number2conc = mpModel->getNumber2QuantityFactor() / mpModel->getCompartments()[0]->getInitialValue();
1852  //C_FLOAT64 number2conc = 1.;
1853 
1854  for (i = 0; i < imax; ++i)
1855  Y1[i] *= number2conc;
1856 
1857  /* write back concentrations of the current state*/
1858  for (i = 0, imax = indep; i < imax; i++)
1859  mpModel->getMetabolitesX()[i]->setValue(tmp[i]);
1860 
1861  //mpState->setUpdateDependentRequired(true);
1863 
1864  return;
1865 }
1866 
1868 {
1869  C_INT i, imax;
1870  C_INT nmetab;
1871 
1872  nmetab = mData.dim;
1873 
1874  CVector<C_FLOAT64> tmp;
1875  tmp.resize(nmetab);
1876 
1877  /* make copy of the current state concentrations */
1878  for (i = 0, imax = nmetab; i < imax; i++)
1879  tmp[i] = mpModel->getMetabolites()[i]->getValue();
1880 
1881  C_FLOAT64 conc2number = mpModel->getQuantity2NumberFactor() * mpModel->getCompartments()[0]->getInitialValue();
1882  //C_FLOAT64 conc2number = 1.;
1883 
1884  /* write new concentrations in the current state */
1885  for (i = 0, imax = nmetab; i < imax; i++)
1886  mpModel->getMetabolites()[i]->setValue(X1[i]*conc2number);
1887 
1890 
1891  C_FLOAT64 number2conc = mpModel->getNumber2QuantityFactor() / mpModel->getCompartments()[0]->getInitialValue();
1892  //C_FLOAT64 number2conc = 1.;
1893 
1894  for (i = 0; i < imax; ++i)
1895  Y1[i] *= number2conc;
1896 
1897  /* write back concentrations of the current state*/
1898  for (i = 0, imax = nmetab; i < imax; i++)
1899  mpModel->getMetabolites()[i]->setValue(tmp[i]);
1900 
1902 
1903  return;
1904 }
1905 
1906 /**
1907 MAP_INDEX used for sorting of SchurDecompostion
1908  */
1909 
1910 void CTSSAMethod::map_index(C_FLOAT64 *eval_r, C_INT *index, const C_INT & dim)
1911 {
1912  C_INT i, j, count;
1913  C_INT max;
1914  C_FLOAT64 max_value, factor;
1915 
1916  CVector< C_FLOAT64 > abs_eval_r(dim);
1917 
1918  max_value = eval_r[0];
1919 
1920  for (i = 1; i < dim; i++)
1921  if (eval_r[i] > max_value)
1922  max_value = eval_r[i];
1923 
1924  if (max_value > 0)
1925  factor = 1.1;
1926  else
1927  {
1928  if (max_value == 0)
1929  {
1930  max_value = 10;
1931  factor = 1;
1932  }
1933  else
1934  factor = 0.;
1935  }
1936 
1937  for (i = 0; i < dim; i++)
1938  {
1939  index[i] = 0;
1940  //abs_eval_r[i] = fabs(eval_r[i]);
1941  abs_eval_r[i] = (eval_r[i]);
1942  }
1943 
1944  count = dim;
1945 
1946  for (i = 0; i < dim; i++)
1947  {
1948  max = i;
1949 
1950  for (j = 0; j < dim; j++)
1951  {
1952  //if (abs_eval_r[j] > abs_eval_r[max])
1953  if (abs_eval_r[j] < abs_eval_r[max])
1954  max = j;
1955  }
1956 
1957  index[max] = count;
1958  abs_eval_r[max] = factor * max_value;
1959  count --;
1960  }
1961 
1962  for (i = 0; i < dim - 1; i++)
1963  if (eval_r[i] == eval_r[i + 1])
1964  index[i + 1] = index[i];
1965 
1966  return;
1967 }
1968 
1969 /**
1970 This is the test only. We try to reorder the Schur matrix from slowest mode (on the bottom) to the fastest (on the top)
1971 The function map_index_desc() is used on the end of schur() in order to produce the orthogonal slow space
1972  */
1973 void CTSSAMethod::map_index_desc(C_FLOAT64 *eval_r, C_INT *index, const C_INT & dim)
1974 {
1975  C_INT i, j, count;
1976  C_INT min;
1977  C_FLOAT64 min_value, factor;
1978 
1979  CVector< C_FLOAT64 > abs_eval_r(dim);
1980 
1981  factor = 1.1;
1982 
1983  min_value = eval_r[0];
1984 
1985  for (i = 1; i < dim; i++)
1986  if (eval_r[i] < min_value)
1987  min_value = eval_r[i];
1988 
1989  for (i = 0; i < dim; i++)
1990  {
1991  index[i] = 0;
1992  //abs_eval_r[i] = fabs(eval_r[i]);
1993  abs_eval_r[i] = (eval_r[i]);
1994  }
1995 
1996  count = dim;
1997 
1998  for (i = 0; i < dim; i++)
1999  {
2000  min = i;
2001 
2002  for (j = 0; j < dim; j++)
2003  {
2004  //if (abs_eval_r[j] > abs_eval_r[max])
2005  if (abs_eval_r[j] >= abs_eval_r[min])
2006  min = j;
2007  }
2008 
2009  index[min] = count;
2010  abs_eval_r[min] = factor * min_value;
2011  count --;
2012  }
2013 
2014  for (i = 0; i < dim - 1; i++)
2015  if (eval_r[i] == eval_r[i + 1])
2016  index[i + 1] = index[i];
2017 
2018  return;
2019 }
2020 
2021 /**
2022 UPDATE_NID: used for sorting of SchurDecompostion
2023  */
2024 void CTSSAMethod::update_nid(C_INT *index, C_INT *nid, const C_INT & dim)
2025 {
2026  C_INT k;
2027 
2028  for (k = 0; k < dim; k++)
2029  nid[k] = 0;
2030 
2031  for (k = 1; k < dim - 1; k++)
2032  if (index[k] == index[k + 1])
2033  nid[k - 1] = k;
2034 
2035  return;
2036 }
2037 
2038 /**
2039 UPDATE_PID: used for sorting of SchurDecompostion
2040  */
2041 void CTSSAMethod::update_pid(C_INT *index, C_INT *pid, const C_INT & dim)
2042 {
2043  C_INT k;
2044 
2045  for (k = 0; k < dim; k++)
2046  pid[k] = 0;
2047 
2048  for (k = 1; k < dim; k++)
2049  if (index[k] == index[k - 1])
2050  pid[k] = k;
2051 
2052  return;
2053 }
2054 
2056 {
2057  /* Retrieve the model to calculate */
2058  mpModel = mpProblem->getModel();
2059 
2060  /* Reset lsoda */
2061 
2062  mLsodaStatus = 1;
2063  mState = 1;
2064  mJType = 2;
2065  mErrorMsg.str("");
2067 
2068  /* Release previous state and make the initialState the current */
2069  pdelete(mpState);
2070  mpState = new CState(*initialState);
2072 
2073  mTime = mpState->getTime();
2074 
2075  if (mReducedModel)
2076  {
2078  }
2079  else
2080  {
2082  }
2083 
2084  mYdot.resize(mData.dim);
2085  // mY_initial.resize(mData.dim);
2087 
2088  /* Configure lsoda */
2089  mRtol = 1.e-6; // * getValue("Relative Tolerance").pUDOUBLE;
2090  initializeAtol();
2091 
2092  mDWork.resize(22 + mData.dim * std::max<C_INT>(16, mData.dim + 9));
2093  mDWork[4] = mDWork[5] = mDWork[6] = mDWork[7] = mDWork[8] = mDWork[9] = 0.0;
2094  mIWork.resize(20 + mData.dim);
2095  mIWork[4] = mIWork[6] = mIWork[9] = 0;
2096 
2097  mIWork[5] = 10000 ; // * getValue("Max Internal Steps").pUINT;
2098  mIWork[7] = 12;
2099  mIWork[8] = 5;
2100 
2101  return;
2102 }
2103 
2105 {
2106  //C_FLOAT64 * pTolerance = getValue("Absolute Tolerance").pUDOUBLE;
2107 
2108  mAtol.resize(mData.dim);
2109  C_FLOAT64 * pAtol = mAtol.array();
2110  C_FLOAT64 * pEnd = pAtol + mAtol.size();
2111 
2112  CModelEntity *const* ppEntity = mpModel->getStateTemplate().beginIndependent();
2113  const CMetab * pMetab;
2114 
2115  for (; pAtol != pEnd; ++pAtol, ++ppEntity)
2116  {
2117  *pAtol = 1.e-12 ; // *pTolerance;
2118 
2119  // Rescale for metabolites as we are using particle numbers
2120  if ((pMetab = dynamic_cast< const CMetab * >(*ppEntity)) != NULL)
2121  {
2122  *pAtol *=
2124  }
2125  }
2126 }
2127 
2128 void CTSSAMethod::EvalF(const C_INT * n, const C_FLOAT64 * t, const C_FLOAT64 * y, C_FLOAT64 * ydot)
2129 {static_cast<Data *>((void *) n)->pMethod->evalF(t, y, ydot);}
2130 
2131 void CTSSAMethod::evalF(const C_FLOAT64 * t, const C_FLOAT64 * /* y */, C_FLOAT64 * ydot)
2132 {
2133  mpState->setTime(*t);
2134 
2136  // TO REMOVE : mpModel->applyAssignments();
2138 
2139  if (mReducedModel)
2141  else
2143 
2144  return;
2145 }
2146 
2147 /**************** display matrix methods ***********************************/
2148 
2149 // to activate the printing set flag_...=0:
2150 
2151 // flag_jacob=0 to print Jacobian
2152 // flag_schur=0 to print matrices of Schur decomposition
2153 // flag_tab =0 to print the Tabs with slow space Analysis
2154 // flag_deufl=0 to prove the Deuflhard algorithm
2155 // flag_Td =0 to print the transformation matrices mTd and mTdInverse
2156 // flag_sylvester=0 to print the transformed Jacobian: mTdInverse*Jacobian_initial*mTd (should be diagonal)
2157 // flag_norm =0 for printing "norm story"
2158 // flag_orthog =0 to print the matrices proved the orthogonality of transformation
2159 
2160 const int & CTSSAMethod::getCurrentStep() const
2161 {
2162  return mCurrentStep;
2163 }
2164 
2165 /**
2166  * return mVec_TimeScale for visualization in ILDM-tab
2167  * in the CQTSSAResultSubWidget
2168  **/
2170 {
2171  return mVec_TimeScale[step - 1];
2172 }
2173 
2174 /**
2175  * Empty every vector to be able to fill them with new values for a new calculation.
2176  * Also nullify the step counter.
2177  **/
2179 {}
2180 
2181 /**
2182  *upgrade all vectors with values from actually calculation for current step
2183  **/
2184 void CTSSAMethod::setVectors(int /* slowMode */)
2185 {}
2186 
2187 /**
2188  * Create the CArraAnnotations for every ILDM-tab in the CQTSSAResultSubWidget.
2189  * Input for each CArraAnnotations is a separate CMatrix.
2190  **/
2192 {}
void schur(C_INT &info)
CCopasiDataModel * getObjectDataModel()
void update_pid(C_INT *index, C_INT *pid, const C_INT &dim)
#define C_INT
Definition: copasi.h:115
void mat_anal_fast_space_thomas(C_INT &slow)
void sylvester(C_INT slow, C_INT &info)
#define MCTSSAMethod
#define pdelete(p)
Definition: copasi.h:215
void update_nid(C_INT *index, C_INT *nid, const C_INT &dim)
CCopasiVectorN< CEvent > & getEvents()
Definition: CModel.cpp:1110
void setCurrentState(CState *currentState)
std::vector< C_FLOAT64 > mCurrentTime
Definition: CTSSAMethod.h:460
CTSSAProblem * mpProblem
Definition: CTSSAMethod.h:55
int dtrsyl_(char *trana, char *tranb, integer *isgn, integer *m, integer *n, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *c__, integer *ldc, doublereal *scale, integer *info)
const CCopasiVector< CMetab > & getMetabolites() const
Definition: CModel.cpp:1051
const CCopasiVectorN< CModelValue > & getModelValues() const
Definition: CModel.cpp:1060
CVector< C_FLOAT64 > mYdot
Definition: CTSSAMethod.h:190
CMatrix< C_FLOAT64 > mVslow
Definition: CTSSAMethod.h:261
const int & getCurrentStep() const
virtual size_t size() const
CMatrix< C_FLOAT64 > mTdInverse
Definition: CTSSAMethod.h:230
virtual void initializeParameter()
CVector< C_FLOAT64 > mY_initial
Definition: CTSSAMethod.h:195
CState * mpState
Definition: CTSSAMethod.h:175
void setOstream(std::ostream &os)
iterator begin()
void updateSimulatedValues(const bool &updateMoieties)
Definition: CModel.cpp:1851
void createAnnotationsM()
CModel * mpModel
Definition: CTSSAMethod.h:336
#define fatalError()
CState * mpCurrentState
Definition: CTSSAMethod.h:50
CVector< C_FLOAT64 > mVslow_space
Definition: CTSSAMethod.h:271
CVector< C_FLOAT64 > mDWork
Definition: CTSSAMethod.h:321
virtual void predifineAnnotation()
void initializeIntegrationsParameter()
static CTSSAMethod * createMethod(CCopasiMethod::SubType subType=CCopasiMethod::unset)
Definition: CTSSAMethod.cpp:44
Definition: CState.h:305
void integrationStep(const double &deltaT)
CVector< C_FLOAT64 > mAtol
Definition: CTSSAMethod.h:301
void resize(size_t size, const bool &copy=false)
Definition: CVector.h:301
void schur_desc(C_INT &info)
void setVectors(int slowMode)
std::ostringstream mErrorMsg
Definition: CTSSAMethod.h:306
size_t getNumODEMetabs() const
Definition: CModel.cpp:1124
const CCopasiMethod::SubType & getSubType() const
void map_index(C_FLOAT64 *eval_r, C_INT *index, const C_INT &dim)
#define C_UNUSED(p)
Definition: copasi.h:220
#define C_INT32
Definition: copasi.h:90
CVector< C_INT > mIWork
Definition: CTSSAMethod.h:326
double orthog(C_INT &number1, C_INT &number2)
Definition: CMetab.h:178
CMatrix< C_FLOAT64 > mQ
Definition: CTSSAMethod.h:215
void mat_anal_mod_space(C_INT &slow)
void setTime(const C_FLOAT64 &time)
Definition: CState.cpp:326
CVector< C_FLOAT64 > mVfast_space
Definition: CTSSAMethod.h:276
void map_index_desc(C_FLOAT64 *eval_r, C_INT *index, const C_INT &dim)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
static void EvalF(const C_INT *n, const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
CMatrix< C_FLOAT64 > mQ_desc
Definition: CTSSAMethod.h:216
CMatrix< C_FLOAT64 > mJacobian_initial
Definition: CTSSAMethod.h:210
bool removeParameter(const std::string &name)
#define MCTrajectoryMethod
iterator end()
std::vector< CType * >::const_iterator const_iterator
Definition: CCopasiVector.h:57
const C_FLOAT64 & getQuantity2NumberFactor() const
Definition: CModel.cpp:2354
const C_FLOAT64 & getNumber2QuantityFactor() const
Definition: CModel.cpp:2357
size_t getNumDependentReactionMetabs() const
Definition: CModel.cpp:1133
void calculateDerivativesX(C_FLOAT64 *X1, C_FLOAT64 *Y1)
bool elevateChildren()
CMatrix< C_FLOAT64 > mR_desc
Definition: CTSSAMethod.h:221
void setState(const CState &state)
Definition: CModel.cpp:1785
void calculateDerivatives(C_FLOAT64 *X1, C_FLOAT64 *Y1)
void mat_anal_fast_space(C_INT &slow)
int mCurrentStep
Definition: CTSSAMethod.h:466
int dtrexc_(char *compq, integer *n, doublereal *t, integer *ldt, doublereal *q, integer *ldq, integer *ifst, integer *ilst, doublereal *work, integer *info)
const Value & getValue() const
bool mReducedModel
Definition: CTSSAMethod.h:291
#define MCCopasiMethod
bool setValue(const std::string &name, const CType &value)
virtual bool isValidProblem(const CCopasiProblem *pProblem)
unsigned C_INT32 * pUINT
size_t getNumIndependent() const
Definition: CState.cpp:342
void setModel(CModel *model)
CMatrix< C_FLOAT64 > mQz
Definition: CTSSAMethod.h:235
CModelEntity ** beginIndependent()
Definition: CState.cpp:208
CCopasiParameter * getParameter(const std::string &name)
CVector< C_FLOAT64 > getVec_TimeScale(int step)
CLSODA mLSODA
Definition: CTSSAMethod.h:311
virtual void resize(size_t rows, size_t cols, const bool &copy=false)
Definition: CMatrix.h:151
size_t getNumIndependentReactionMetabs() const
Definition: CModel.cpp:1130
void mat_anal_metab(C_INT &slow)
size_t size() const
Definition: CVector.h:100
void calculateDerivativesX(C_FLOAT64 *derivativesX)
Definition: CModel.cpp:1941
CCopasiVectorNS< CCompartment > & getCompartments()
Definition: CModel.cpp:1145
const C_FLOAT64 & getValue() const
const C_FLOAT64 & getTime() const
Definition: CState.cpp:325
CMatrix< C_FLOAT64 > mR
Definition: CTSSAMethod.h:220
C_INT mState
Definition: CTSSAMethod.h:316
C_FLOAT64 mRtol
Definition: CTSSAMethod.h:296
std::vector< CVector< C_FLOAT64 > > mVec_TimeScale
Definition: CTSSAMethod.h:461
CMatrix< C_FLOAT64 > mJacobian
Definition: CTSSAMethod.h:205
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 returnCurrentTime(int step)
void setProblem(CTSSAProblem *problem)
void emptyVectors()
size_t getNumAssignmentMetabs() const
Definition: CModel.cpp:1127
CType * array()
Definition: CVector.h:139
const CStateTemplate & getStateTemplate() const
Definition: CModel.cpp:1172
virtual void start(const CState *initialState)
CMatrix< C_FLOAT64 > mTd
Definition: CTSSAMethod.h:225
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)
void calculateDerivatives(C_FLOAT64 *derivatives)
Definition: CModel.cpp:1903
CCopasiParameter * assertParameter(const std::string &name, const CCopasiParameter::Type type, const CType &defaultValue)
Definition: CModel.h:50
C_INT mJType
Definition: CTSSAMethod.h:331
const CCopasiVector< CMetab > & getMetabolitesX() const
Definition: CModel.cpp:1057
virtual void step(const double &deltaT)
C_FLOAT64 mTime
Definition: CTSSAMethod.h:200
const CCompartment * getCompartment() const
Definition: CMetab.cpp:222
void integrationMethodStart(const CState *initialState)
CModel * getModel() const
#define min(a, b)
Definition: f2c.h:175
CTSSAMethod * pMethod
Definition: CTSSAMethod.h:169
C_INT mLsodaStatus
Definition: CTSSAMethod.h:286
void mat_anal_mod(C_INT &slow)
CMatrix< C_FLOAT64 > mVslow_metab
Definition: CTSSAMethod.h:266
C_FLOAT64 * beginIndependent()
Definition: CState.cpp:328
void initializeAtol()
C_FLOAT64 * mY
Definition: CTSSAMethod.h:185
void evalF(const C_FLOAT64 *t, const C_FLOAT64 *y, C_FLOAT64 *ydot)
#define max(a, b)
Definition: f2c.h:176