Excitation-contraction coupling/mitochondrial energetics (ECME) model (Cortassa et al. 2006)

 Download zip file 
Help downloading and running models
Accession:105383
"An intricate network of reactions is involved in matching energy supply with demand in the heart. This complexity arises because energy production both modulates and is modulated by the electrophysiological and contractile activity of the cardiac myocyte. Here, we present an integrated mathematical model of the cardiac cell that links excitation-contraction coupling with mitochondrial energy generation. The dynamics of the model are described by a system of 50 ordinary differential equations. The formulation explicitly incorporates cytoplasmic ATP-consuming processes associated with force generation and ion transport, as well as the creatine kinase reaction. Changes in the electrical and contractile activity of the myocyte are coupled to mitochondrial energetics through the ATP, Ca21, and Na1 concentrations in the myoplasmic and mitochondrial matrix compartments. ..."
Reference:
1 . Cortassa S, Aon MA, Marbán E, Winslow RL, O'Rourke B (2003) An integrated model of cardiac mitochondrial energy metabolism and calcium dynamics. Biophys J 84:2734-55 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell; Electrogenic pump;
Brain Region(s)/Organism:
Cell Type(s): Heart cell;
Channel(s): I L high threshold; I Sodium; I Potassium; Na/Ca exchanger; I_SERCA;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Activity Patterns; Temporal Pattern Generation; Signaling pathways; Calcium dynamics;
Implementer(s):
Search NeuronDB for information about:  I L high threshold; I Sodium; I Potassium; Na/Ca exchanger; I_SERCA;
/*******************************************************************
 *                                                                 *
 * File          : cvspgmr.c                                       *
 * Programmers   : Scott D. Cohen, Alan C. Hindmarsh and           * 
 *                 Radu Serban @ LLNL                              *
 * Version of    : 26 June 2002                                    *
 *-----------------------------------------------------------------*
 * Copyright (c) 2002, The Regents of the University of California * 
 * Produced at the Lawrence Livermore National Laboratory          *
 * All rights reserved                                             *
 * For details, see sundials/cvode/LICENSE                         *
 *-----------------------------------------------------------------*
 * This is the implementation file for the CVODE scaled,           *
 * preconditioned GMRES linear solver, CVSPGMR.                    *
 *                                                                 *
 *******************************************************************/


#include <stdio.h>
#include <stdlib.h>
#include "cvspgmr.h"
#include "cvode.h"
#include "sundialstypes.h"
#include "nvector.h"
#include "sundialsmath.h"
#include "iterativ.h"
#include "spgmr.h"


/* Error Messages */

#define CVSPGMR            "CVSpgmr/CVReInitSpgmr-- "

#define MSG_CVMEM_NULL     CVSPGMR "CVode Memory is NULL.\n\n"

#define MSG_MEM_FAIL       CVSPGMR "A memory request failed.\n\n"

#define MSG_BAD_PRETYPE_1  CVSPGMR "pretype=%d illegal.\n"
#define MSG_BAD_PRETYPE_2  "The legal values are NONE=%d, LEFT=%d, "
#define MSG_BAD_PRETYPE_3  "RIGHT=%d, and BOTH=%d.\n\n"
#define MSG_BAD_PRETYPE    MSG_BAD_PRETYPE_1 MSG_BAD_PRETYPE_2 MSG_BAD_PRETYPE_3

#define MSG_PSOLVE_REQ_1   CVSPGMR "pretype!=NONE, but PSOLVE=NULL is "
#define MSG_PSOLVE_REQ_2   "illegal.\n\n"
#define MSG_PSOLVE_REQ     MSG_PSOLVE_REQ_1 MSG_PSOLVE_REQ_2

#define MSG_BAD_GSTYPE_1   CVSPGMR "gstype=%d illegal.\n"
#define MSG_BAD_GSTYPE_2   "The legal values are MODIFIED_GS=%d and "
#define MSG_BAD_GSTYPE_3   "CLASSICAL_GS=%d.\n\n"
#define MSG_BAD_GSTYPE     MSG_BAD_GSTYPE_1 MSG_BAD_GSTYPE_2 MSG_BAD_GSTYPE_3

/* Other Constants */

#define ZERO RCONST(0.0)
#define ONE  RCONST(1.0)

/******************************************************************
 *                                                                *           
 * Types : CVSpgmrMemRec, CVSpgmrMem                              *
 *----------------------------------------------------------------*
 * The type CVSpgmrMem is pointer to a CVSpgmrMemRec. This        *
 * structure contains CVSpgmr solver-specific data.               *
 *                                                                *
 ******************************************************************/

typedef struct {

  int  g_pretype;     /* type of preconditioning                      */
  int  g_gstype;      /* type of Gram-Schmidt orthogonalization       */
  realtype g_sqrtN;   /* sqrt(N)                                      */
  realtype g_delt;    /* delt = user specified or DELT_DEFAULT        */
  realtype g_deltar;  /* deltar = delt * tq4                          */
  realtype g_delta;   /* delta = deltar * sqrtN                       */
  int  g_maxl;        /* maxl = maximum dimension of the Krylov space */

  long int g_nstlpre;  /* value of nst at the last precond call       */     
  long int g_npe;      /* npe = total number of precond calls         */   
  long int g_nli;      /* nli = total number of linear iterations     */
  long int g_nps;      /* nps = total number of psolve calls          */
  long int g_ncfl;     /* ncfl = total number of convergence failures */

  N_Vector g_ytemp;    /* temp vector passed to jtimes and psolve     */ 
  N_Vector g_x;        /* temp vector used by CVSpgmrSolve            */
  N_Vector g_ycur;     /* CVODE current y vector in Newton Iteration  */
  N_Vector g_fcur;     /* fcur = f(tn, ycur)                          */

  CVSpgmrPrecondFn g_precond; /* precond = user-supplied routine to   */
                              /* compute a preconditioner             */

  CVSpgmrPSolveFn g_psolve;   /* psolve = user-supplied routine to    */
                              /* solve preconditioner linear system   */

  void *g_P_data;            /* P_data passed to psolve and precond   */
  SpgmrMem g_spgmr_mem;      /* spgmr_mem is memory used by the       */
                             /* generic Spgmr solver                  */

  CVSpgmrJtimesFn g_jtimes;  /* jtimes = Jacobian * vector routine    */
                             /*          to be called                 */
  void *g_j_data;            /* j_data is passed to jtimes            */

} CVSpgmrMemRec, *CVSpgmrMem;


/* CVSPGMR linit, lsetup, lsolve, and lfree routines */

static int CVSpgmrInit(CVodeMem cv_mem);

static int CVSpgmrSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
                        N_Vector fpred, booleantype *jcurPtr, N_Vector vtemp1,
                        N_Vector vtemp2, N_Vector vtemp3);

static int CVSpgmrSolve(CVodeMem cv_mem, N_Vector b, N_Vector ycur,
                        N_Vector fcur);

static void CVSpgmrFree(CVodeMem cv_mem);

/* CVSPGMR Atimes and PSolve routines called by generic SPGMR solver */

static int CVSpgmrAtimes(void *cv_mem, N_Vector v, N_Vector z);

static int CVSpgmrPSolve(void *cv_mem, N_Vector r, N_Vector z, int lr);

/* CVSPGMR difference quotient routine for J*v */

static int CVSpgmrDQJtimes(integertype N, N_Vector v, N_Vector Jv, RhsFn f,
                       void *f_data, realtype tn, N_Vector y, N_Vector fy,
                       realtype vnrm, N_Vector ewt, realtype h, realtype uround,
                       void *jac_data, long int *nfePtr, N_Vector ytemp);


/*************** CVSpgmrDQJtimes *************************************

 This routine generates a difference quotient approximation to
 the Jacobian times vector f_y(t,y) * v. The approximation is 
 Jv = vnrm[f(y + v/vnrm) - f(y)], where vnrm = (WRMS norm of v) is
 input, i.e. the WRMS norm of v/vnrm is 1.

**********************************************************************/

static int CVSpgmrDQJtimes(integertype N, N_Vector v, N_Vector Jv, RhsFn f, 
                       void *f_data, realtype tn, N_Vector y, N_Vector fy,
                       realtype vnrm, N_Vector ewt, realtype h, realtype uround,
                       void *jac_data, long int *nfePtr, N_Vector work)
{
  N_Vector ytemp;

  /* Rename work for readability */
  ytemp = work;

  /* Set ytemp = y + (1/vnrm) v */
  N_VLinearSum(ONE/vnrm, v, ONE, y, ytemp);

  /* Set Jv = f(tn, ytemp) */
  f(N, tn, ytemp, Jv, f_data); 
  (*nfePtr)++;
  
  /* Replace Jv by vnrm*(Jv - fy) */
  N_VLinearSum(ONE, Jv, -ONE, fy, Jv);
  N_VScale(vnrm, Jv, Jv);

  return(0);
}

/**********************************************************************/

/* Readability Replacements */

#define N       (cv_mem->cv_N)      
#define uround  (cv_mem->cv_uround)
#define tq      (cv_mem->cv_tq)
#define nst     (cv_mem->cv_nst)
#define tn      (cv_mem->cv_tn)
#define h       (cv_mem->cv_h)
#define gamma   (cv_mem->cv_gamma)
#define gammap  (cv_mem->cv_gammap)   
#define nfe     (cv_mem->cv_nfe)
#define f       (cv_mem->cv_f)
#define f_data  (cv_mem->cv_f_data)
#define ewt     (cv_mem->cv_ewt)
#define errfp   (cv_mem->cv_errfp)
#define mnewt   (cv_mem->cv_mnewt)
#define iopt    (cv_mem->cv_iopt)
#define ropt    (cv_mem->cv_ropt)
#define linit   (cv_mem->cv_linit)
#define lsetup  (cv_mem->cv_lsetup)
#define lsolve  (cv_mem->cv_lsolve)
#define lfree   (cv_mem->cv_lfree)
#define lmem    (cv_mem->cv_lmem)
#define machenv (cv_mem->cv_machenv)
#define setupNonNull (cv_mem->cv_setupNonNull)

#define sqrtN   (cvspgmr_mem->g_sqrtN)   
#define ytemp   (cvspgmr_mem->g_ytemp)
#define x       (cvspgmr_mem->g_x)
#define ycur    (cvspgmr_mem->g_ycur)
#define fcur    (cvspgmr_mem->g_fcur)
#define delta   (cvspgmr_mem->g_delta)
#define deltar  (cvspgmr_mem->g_deltar)
#define npe     (cvspgmr_mem->g_npe)
#define nli     (cvspgmr_mem->g_nli)
#define nps     (cvspgmr_mem->g_nps)
#define ncfl    (cvspgmr_mem->g_ncfl)
#define nstlpre (cvspgmr_mem->g_nstlpre)
#define spgmr_mem (cvspgmr_mem->g_spgmr_mem)


/*************** CVSpgmr *********************************************

 This routine initializes the memory record and sets various function
 fields specific to the Spgmr linear solver module. CVSpgmr first
 calls the existing lfree routine if this is not NULL.  It then sets
 the cv_linit, cv_lsetup, cv_lsolve, cv_lfree fields in (*cvode_mem)
 to be CVSpgmrInit, CVSpgmrSetup, CVSpgmrSolve, and CVSpgmrFree,
 respectively.  It allocates memory for a structure of type
 CVSpgmrMemRec and sets the cv_lmem field in (*cvode_mem) to the
 address of this structure.  It sets setupNonNull in (*cvode_mem),
 and sets the following fields in the CVSpgmrMemRec structure:
   g_pretype = pretype                                       
   g_gstype  = gstype                                       
   g_maxl    = MIN(N,CVSPGMR_MAXL)  if maxl <= 0             
             = maxl                 if maxl > 0              
   g_delt    = CVSPGMR_DELT if delt == 0.0                     
             = delt         if delt != 0.0                     
   g_P_data  = P_data                                        
   g_precond = precond                                       
   g_psolve  = psolve                                        
   g_jtimes  = input parameter jtimes  if jtimes != NULL
             = CVSpgmrDQJtimes         otherwise
   g_j_data  = input parameter jac_data
 Finally, CVSpgmr allocates memory for ytemp and x, and calls
 SpgmrMalloc to allocate memory for the Spgmr solver.  The CVSpgmr
 return value is SUCCESS = 0, LMEM_FAIL = -1, or LIN_ILL_INPUT = -2.

**********************************************************************/

int CVSpgmr(void *cvode_mem, int pretype, int gstype, int maxl, realtype delt,
             CVSpgmrPrecondFn precond, CVSpgmrPSolveFn psolve, void *P_data,
             CVSpgmrJtimesFn jtimes, void *jac_data)

{
  CVodeMem cv_mem;
  CVSpgmrMem cvspgmr_mem;
  int mxl;

  /* Return immediately if cvode_mem is NULL */
  cv_mem = (CVodeMem) cvode_mem;
  if (cv_mem == NULL) {                /* CVode reports this error */
    fprintf(errfp, MSG_CVMEM_NULL);
    return(LMEM_FAIL);
  }

  if (lfree != NULL) lfree(cv_mem);

  /* Set four main function fields in cv_mem */
  linit  = CVSpgmrInit;
  lsetup = CVSpgmrSetup;
  lsolve = CVSpgmrSolve;
  lfree  = CVSpgmrFree;

  /* Get memory for CVSpgmrMemRec */
  lmem = cvspgmr_mem = (CVSpgmrMem) malloc(sizeof(CVSpgmrMemRec));
  if (cvspgmr_mem == NULL) {
    fprintf(errfp, MSG_MEM_FAIL);
    return(LMEM_FAIL);
  }

  /* Set Spgmr parameters that have been passed in call sequence */
  cvspgmr_mem->g_pretype    = pretype;
  cvspgmr_mem->g_gstype     = gstype;
  mxl = cvspgmr_mem->g_maxl = (maxl <= 0) ? MIN(CVSPGMR_MAXL, N) : maxl;
  cvspgmr_mem->g_delt       = (delt == ZERO) ? CVSPGMR_DELT : delt;
  cvspgmr_mem->g_P_data     = P_data;
  cvspgmr_mem->g_precond    = precond;
  cvspgmr_mem->g_psolve     = psolve;

  /* Set Jacobian times vector routine to user's jtimes or CVSpgmrDQJtimes */
  if(jtimes == NULL) {
    cvspgmr_mem->g_jtimes = CVSpgmrDQJtimes;
  } else {
    cvspgmr_mem->g_jtimes = jtimes;
  }

  /* Set Jacobian data */
  cvspgmr_mem->g_j_data = jac_data;

  /* Check for legal pretype, precond, and psolve */ 
  if ((pretype != NONE) && (pretype != LEFT) &&
      (pretype != RIGHT) && (pretype != BOTH)) {
    fprintf(errfp, MSG_BAD_PRETYPE, pretype, NONE, LEFT, RIGHT, BOTH);
    return(LIN_ILL_INPUT);
  }
  if ((pretype != NONE) && (psolve == NULL)) {
    fprintf(errfp, MSG_PSOLVE_REQ);
    return(LIN_ILL_INPUT);
  }

  /* Check for legal gstype */
  if ((gstype != MODIFIED_GS) && (gstype != CLASSICAL_GS)) {
    fprintf(errfp, MSG_BAD_GSTYPE, gstype, MODIFIED_GS, CLASSICAL_GS);
    return(LIN_ILL_INPUT);
  }

  /* Set setupNonNull = TRUE iff there is preconditioning (pretype != NONE)
     and there is a preconditioning setup phase (precond != NULL)          */
  setupNonNull = (pretype != NONE) && (precond != NULL);

  /* Allocate memory for ytemp and x */
  ytemp = N_VNew(N, machenv);
  if (ytemp == NULL) {
    fprintf(errfp, MSG_MEM_FAIL);
    return(LMEM_FAIL);
  }
  x = N_VNew(N, machenv);
  if (x == NULL) {
    fprintf(errfp, MSG_MEM_FAIL);
    N_VFree(ytemp);
    return(LMEM_FAIL);
  }

  /* Call SpgmrMalloc to allocate workspace for Spgmr */
  spgmr_mem = SpgmrMalloc(N, mxl, machenv);
  if (spgmr_mem == NULL) {
    fprintf(errfp, MSG_MEM_FAIL);
    N_VFree(ytemp);
    N_VFree(x);
    return(LMEM_FAIL);
  }

  return(SUCCESS);
}


/*************** CVReInitSpgmr****************************************

 This routine resets the link between the main CVODE module and the
 Spgmr linear solver module CVSPGMR.  No memory freeing or allocation
 operations are done, as the existing linear solver memory is assumed
 sufficient.  All other initializations are the same as in CVSpgmr.
 The return value is SUCCESS=0, LMEM_FAIL=-1, or LIN_ILL_INPUT=-2.

**********************************************************************/

int CVReInitSpgmr(void *cvode_mem, int pretype, int gstype, int maxl,
            realtype delt, CVSpgmrPrecondFn precond, CVSpgmrPSolveFn psolve,
            void *P_data, CVSpgmrJtimesFn jtimes, void *jac_data)

{
  CVodeMem cv_mem;
  CVSpgmrMem cvspgmr_mem;
  int mxl;

  /* Return immediately if cvode_mem is NULL */
  cv_mem = (CVodeMem) cvode_mem;
  if (cv_mem == NULL) {                /* CVode reports this error */
    fprintf(errfp, MSG_CVMEM_NULL);
    return(LMEM_FAIL);
  }

  /* Set four main function fields in cv_mem */
  linit  = CVSpgmrInit;
  lsetup = CVSpgmrSetup;
  lsolve = CVSpgmrSolve;
  lfree  = CVSpgmrFree;

  cvspgmr_mem = lmem;   /* Use existing linear solver memory pointer */

  /* Set Spgmr parameters that have been passed in call sequence */
  cvspgmr_mem->g_pretype    = pretype;
  cvspgmr_mem->g_gstype     = gstype;
  mxl = cvspgmr_mem->g_maxl = (maxl <= 0) ? MIN(CVSPGMR_MAXL, N) : maxl;
  cvspgmr_mem->g_delt       = (delt == ZERO) ? CVSPGMR_DELT : delt;
  cvspgmr_mem->g_P_data     = P_data;
  cvspgmr_mem->g_precond    = precond;
  cvspgmr_mem->g_psolve     = psolve;

  /* Set Jacobian times vector routine to user's jtimes or CVSpgmrDQJtimes */
  if(jtimes == NULL) {
    cvspgmr_mem->g_jtimes = CVSpgmrDQJtimes;
  } else {
    cvspgmr_mem->g_jtimes = jtimes;
  }

  /* Set Jacobian data */
  cvspgmr_mem->g_j_data = jac_data;

  /* Check for legal pretype, precond, and psolve */ 
  if ((pretype != NONE) && (pretype != LEFT) &&
      (pretype != RIGHT) && (pretype != BOTH)) {
    fprintf(errfp, MSG_BAD_PRETYPE, pretype, NONE, LEFT, RIGHT, BOTH);
    return(LIN_ILL_INPUT);
  }
  if ((pretype != NONE) && (psolve == NULL)) {
    fprintf(errfp, MSG_PSOLVE_REQ);
    return(LIN_ILL_INPUT);
  }

  /* Check for legal gstype */
  if ((gstype != MODIFIED_GS) && (gstype != CLASSICAL_GS)) {
    fprintf(errfp, MSG_BAD_GSTYPE, gstype, MODIFIED_GS, CLASSICAL_GS);
    return(LIN_ILL_INPUT);
  }

  /* Set setupNonNull = TRUE iff there is preconditioning (pretype != NONE)
     and there is a preconditioning setup phase (precond != NULL)          */
  setupNonNull = (pretype != NONE) && (precond != NULL);

  return(SUCCESS);
}


/* Additional readability Replacements */

#define pretype (cvspgmr_mem->g_pretype)
#define gstype  (cvspgmr_mem->g_gstype)
#define delt    (cvspgmr_mem->g_delt)
#define maxl    (cvspgmr_mem->g_maxl)
#define psolve  (cvspgmr_mem->g_psolve)
#define precond (cvspgmr_mem->g_precond)
#define P_data  (cvspgmr_mem->g_P_data)
#define jtimes  (cvspgmr_mem->g_jtimes)
#define j_data  (cvspgmr_mem->g_j_data)


/*************** CVSpgmrInit *****************************************

 This routine does remaining initializations specific to the Spgmr 
 linear solver.

**********************************************************************/

static int CVSpgmrInit(CVodeMem cv_mem)
{
  CVSpgmrMem cvspgmr_mem;
  cvspgmr_mem = (CVSpgmrMem) lmem;

  /* Initialize sqrtN and counters, and set workspace lengths */

  sqrtN = RSqrt(N);
  npe = nli = nps = ncfl = nstlpre = 0;
    
  if (iopt != NULL) {
    iopt[SPGMR_NPE] = npe;
    iopt[SPGMR_NLI] = nli;
    iopt[SPGMR_NPS] = nps;
    iopt[SPGMR_NCFL] = ncfl;
    iopt[SPGMR_LRW] = N*(maxl + 5) + maxl*(maxl + 4) + 1;
    iopt[SPGMR_LIW] = 0;
  }

  return(LINIT_OK);
}

/*************** CVSpgmrSetup ****************************************

 This routine does the setup operations for the Spgmr linear solver.
 It makes a decision as to whether or not to signal for re-evaluation
 of Jacobian data in the precond routine, based on various state
 variables, then it calls precond.  If we signal for re-evaluation,
 then we reset jcur = *jcurPtr to TRUE, regardless of the precond output.
 In any case, if jcur == TRUE, we increment npe and save nst in nstlpre.

**********************************************************************/

static int CVSpgmrSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
                        N_Vector fpred, booleantype *jcurPtr, N_Vector vtemp1,
                        N_Vector vtemp2, N_Vector vtemp3)
{
  booleantype jbad, jok;
  realtype dgamma;
  int  ier;
  CVSpgmrMem cvspgmr_mem;

  cvspgmr_mem = (CVSpgmrMem) lmem;

  /* Use nst, gamma/gammap, and convfail to set J eval. flag jok */
  dgamma = ABS((gamma/gammap) - ONE);
  jbad = (nst == 0) || (nst > nstlpre + CVSPGMR_MSBPRE) ||
    ((convfail == FAIL_BAD_J) && (dgamma < CVSPGMR_DGMAX)) ||
    (convfail == FAIL_OTHER);
  *jcurPtr = jbad;
  jok = !jbad;

  /* Call precond routine and possibly reset jcur */
  ier = precond(N, tn, ypred, fpred, jok, jcurPtr, gamma, ewt, h,
                uround, &nfe, P_data, vtemp1, vtemp2, vtemp3);
  if (jbad) *jcurPtr = TRUE;

  /* If jcur = TRUE, increment npe and save nst value */
  if (*jcurPtr) {
    npe++;
    nstlpre = nst;
  }

  /* Set npe, and return the same value ier that precond returned */
  if (iopt != NULL) iopt[SPGMR_NPE] = npe;
  return(ier);
}

/*************** CVSpgmrSolve ****************************************

 This routine handles the call to the generic solver SpgmrSolve
 for the solution of the linear system Ax = b with the SPGMR method,
 without restarts.  The solution x is returned in the vector b.

 If the WRMS norm of b is small, we return x = b (if this is the first
 Newton iteration) or x = 0 (if a later Newton iteration).

 Otherwise, we set the tolerance parameter and initial guess (x = 0),
 call SpgmrSolve, and copy the solution x into b.  The x-scaling and
 b-scaling arrays are both equal to ewt, and no restarts are allowed.

 The counters nli, nps, and ncfl are incremented, and the return value
 is set according to the success of SpgmrSolve.  The success flag is
 returned if SpgmrSolve converged, or if this is the first Newton
 iteration and the residual norm was reduced below its initial value.

**********************************************************************/

static int CVSpgmrSolve(CVodeMem cv_mem, N_Vector b, N_Vector ynow,
                        N_Vector fnow)
{
  realtype bnorm, res_norm;
  CVSpgmrMem cvspgmr_mem;
  int nli_inc, nps_inc, ier;
  
  cvspgmr_mem = (CVSpgmrMem) lmem;

  /* Test norm(b); if small, return x = 0 or x = b */
  deltar = delt*tq[4]; 

  bnorm = N_VWrmsNorm(b, ewt);
  if (bnorm <= deltar) {
    if (mnewt > 0) N_VConst(ZERO, b); 
    return(0);
  }

  /* Set vectors ycur and fcur for use by the Atimes and Psolve routines */
  ycur = ynow;
  fcur = fnow;

  /* Set inputs delta and initial guess x = 0 to SpgmrSolve */  
  delta = deltar * sqrtN;
  N_VConst(ZERO, x);
  
  /* Call SpgmrSolve and copy x to b */
  ier = SpgmrSolve(spgmr_mem, cv_mem, x, b, pretype, gstype, delta, 0,
                   cv_mem, ewt, ewt, CVSpgmrAtimes, CVSpgmrPSolve,
                   &res_norm, &nli_inc, &nps_inc);

  N_VScale(ONE, x, b);
  
  /* Increment counters nli, nps, and ncfl */
  nli += nli_inc;
  nps += nps_inc;
  if (iopt != NULL) {
    iopt[SPGMR_NLI] = nli;
    iopt[SPGMR_NPS] = nps;
  }  
  if (ier != 0) { 
    ncfl++;
    if (iopt != NULL) iopt[SPGMR_NCFL] = ncfl;
  }

  /* Set return value to -1, 0, or 1 */
  if (ier < 0) return(-1);  
  if ((ier == SPGMR_SUCCESS) || 
      ((ier == SPGMR_RES_REDUCED) && (mnewt == 0)))
    return(0);
  return(1);  
}

/*************** CVSpgmrFree *****************************************

 This routine frees memory specific to the Spgmr linear solver.

**********************************************************************/

static void CVSpgmrFree(CVodeMem cv_mem)
{
  CVSpgmrMem cvspgmr_mem;

  cvspgmr_mem = (CVSpgmrMem) lmem;
  
  N_VFree(ytemp);
  N_VFree(x);
  SpgmrFree(spgmr_mem);
  free(cvspgmr_mem);
}

/***************** CVSpgmrAtimes *************************************

 This routine generates the matrix-vector product z = Mv, where
 M = I - gamma*J. The product J*v is obtained by calling the jtimes 
 routine. It is then scaled by -gamma and added to v to obtain M*v.
 The return value is the same as the value returned by jtimes --
 0 if successful, nonzero otherwise.

**********************************************************************/

static int CVSpgmrAtimes(void *cvode_mem, N_Vector v, N_Vector z)
{
  realtype rho;
  CVodeMem   cv_mem;
  CVSpgmrMem cvspgmr_mem;
  int jtflag;

  cv_mem = (CVodeMem) cvode_mem;
  cvspgmr_mem = (CVSpgmrMem) lmem;

  /* If rho = norm(v) is 0, return z = 0 */
  rho = N_VWrmsNorm(v, ewt);
  if (rho == ZERO) {
    N_VConst(ZERO, z);
    return(0);
  }
  
  jtflag = jtimes(N, v, z, f, f_data, tn, ycur, fcur, rho, ewt, h,
                  uround, j_data, &nfe, ytemp);
  if (jtflag != 0) return(jtflag);

  N_VLinearSum(ONE, v, -gamma, z, z);

  return(0);
}

/*************** CVSpgmrPSolve ***************************************

 This routine interfaces between the generic SpgmrSolve routine and
 the user's psolve routine.  It passes to psolve all required state 
 information from cvode_mem.  Its return value is the same as that
 returned by psolve. Note that the generic SPGMR solver guarantees
 that CVSpgmrPSolve will not be called in the case in which
 preconditioning is not done. This is the only case in which the
 user's psolve routine is allowed to be NULL.

**********************************************************************/

static int CVSpgmrPSolve(void *cvode_mem, N_Vector r, N_Vector z, int lr)
{
  CVodeMem   cv_mem;
  CVSpgmrMem cvspgmr_mem;
  int ier;

  cv_mem = (CVodeMem) cvode_mem;
  cvspgmr_mem = (CVSpgmrMem)lmem;

  ier = psolve(N, tn, ycur, fcur, ytemp, gamma, ewt, delta, &nfe, r,
               lr, P_data, z);
  /* This call is counted in nps within the CVSpgmrSolve routine */

  return(ier);     
}

Loading data, please wait...