newton.c

GESPI 2.0动态系统模拟工具

#include "copyleft.h"

/*
    GEPASI - a simulator of metabolic pathways and other dynamical systems
    Copyright (C) 1989, 1992  Pedro Mendes
*/

/*************************************/
/*                                   */
/*      steady-state solution by     */
/*      the damped Newton method     */
/*                                   */
/*        Zortech C/C++ 3.0 r4       */
/*          MICROSOFT C 6.00         */
/*          Visual C/C++ 1.0         */
/*           QuickC/WIN 1.0          */
/*             ULTRIX cc             */
/*              GNU gcc              */
/*                                   */
/*   (include here compilers that    */
/*   compiled GEPASI successfully)   */
/*                                   */
/*************************************/


#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "globals.h"                       /* global symbols                 */
#include "globvar.h"                       /* extern global variable         */
#include "matrix.h"                        /* matrix manipulation            */
#include "rates.h"                         /* rates and jacobian calculation */
#include "datab.h"						   /* database structures			 */

#define N_OK 0
#define N_LMIN 1
#define N_NOCONV 2
#define N_JSING 3

int newton( void )
{
 register int m, i, j, k;
 unsigned long df;                         /* dumping factor (power of two)  */
 static double xpss[MAX_MET],              /* x at iteration m               */
               ma[MAX_MET][MAX_MET],       /* inverse of jacobian matrix     */
               mb[MAX_MET][MAX_MET],       /* copy of jacobian matrix        */
               h[MAX_MET], h2[MAX_MET],    /* increment and damped increment */
               re,nre;                     /* residual error & new res.err.  */


 for(j=0;j<totmet;j++)                     /* make x[i] the first guess      */
  xpss[j] = xss[j] = x[j];
                                           /* initial values of f(x) and re  */
 calcrates( xpss );                        /* f(x)                           */
 for(j=0, re = (double) 0; j<nmetab;j++)   /* ||f(x)||2 residual error       */
  re += rate[j] * rate[j];
 re = sqrt( re );
 if ( re < options.hrcz ) return N_OK;     /* if already in s.s. job is done!*/
 for(m=0;m<newtlim;m++)                    /* iterate "until satisfied"      */
 {
  if (options.debug) printf("\nnewton() - iteration %2d, ", m);
  calcjacob( xpss );                       /* f'(x) jacobian                 */
  memcpy( mb, jacob, MAX_MET * MAX_MET * sizeof(double) );
  lu_inverse( (double(*)[MAX_MET][MAX_MET])mb,
              (double(*)[MAX_MET][MAX_MET])ma,
              indmet );                    /* inv( f'(x) )                   */
  multmtxv( ma, rate, h,
            indmet, indmet, indmet );      /* h = inv(f'(x)) . f(x)          */
  nre = 2 * re;                            /* just to start with nre > re    */
  for( i=0, df=1; (i<32) && (nre>re) ; i++ )
  {
   if (options.debug) printf(".");
   if (i) df *= 2;                         /* 2 raised to the i power        */
   for( j=0; j<indmet; j++ )
    h2[j] = h[j]/df;
   addvctsm( xpss, h2, -1.0, xss, indmet );/* xss = xpss - h2                */
   /* update the conserved moieties			*/
   for( j=indmet; j<nmetab; j++)
   {
    xss[j] = moiety[j];
    for( k=0; k<indmet; k++)
     xss[j] -= ld[j][k] * xss[k];
   }
   calcrates( xss );						/* f(x)                           */
   for(j=0, nre=0.0; j<indmet; j++)			/* ||f(x)||2 residual error       */
    nre += rate[j] * rate[j];
   nre = sqrt( nre );
  }
  re = nre;									/* keep nre for next step         */
  if (i==32)
  {
   if ( re < options.hrcz )  return N_OK;   /* this is a zero, after all!     */
   else										/* stuck in a local minima		  */
   {
    if (options.debug) printf(" local minima ");
    return N_LMIN;
   }
  }
  if (options.debug) printf(" re=% .4le", re );
  for(j=0; j<nmetab; j++)
   xpss[j] = xss[j];
  if ( re < options.hrcz )
   return N_OK;
 }
 if (options.debug) printf(" no solution ");
 return N_NOCONV;                          /* no convergence                  */
}