research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Remove old version of CVODE

Browse source 
We now use the current version from the Sundials git submodule
This commit is contained in:
Ray Speth 2015-11-13 14:57:19 -05:00
parent b4a1fb2db1
commit 4e551f69ab
25 changed files with 0 additions and 10843 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

609
ext/cvode/include/band.h
View file

@ -1,609 +0,0 @@
/******************************************************************
* *
* File : band.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 5 May 1998 *
*----------------------------------------------------------------*
* This is the header file for a generic BAND linear solver *
* package. There are two sets of band solver routines listed in *
* this file: one set uses type BandMat defined below and the *
* other set uses the type real ** for band matrix arguments. *
* The two sets of band solver routines make it easy to work *
* with two types of band matrices: *
* *
* (1) The BandMat type is intended for use with large *
* band matrices whose elements/columns may be stored in *
* non-contiguous memory locations or even distributed into *
* different processor memories. This type may be modified to *
* include such distribution information. If this is done, *
* then all the routines that use BandMat must be modified to *
* reflect the new data structure. *
* *
* (2) The set of routines that use real ** (and NOT the BandMat *
* type) is intended for use with small matrices which can *
* easily be allocated within a contiguous block of memory *
* on a single processor. *
* *
* Routines that work with the type BandMat begin with "Band". *
* The BandAllocMat function allocates a band matrix for use in *
* the other matrix routines listed in this file. Matrix storage *
* details are given in the documentation for the type BandMat. *
* The BandAllocPiv function allocates memory for pivot *
* information. The storage allocated by BandAllocMat and *
* BandAllocPiv is deallocated by the routines BandFreeMat and *
* BandFreePiv, respectively. The BandFactor and BandBacksolve *
* routines perform the actual solution of a band linear system. *
* Note that the BandBacksolve routine has a parameter b of type *
* N_Vector. The current implementation makes use of a machine *
* environment specific macro (N_VDATA) which may not exist for *
* other implementations of the type N_Vector. Thus, the *
* implementation of BandBacksolve may need to change if the *
* type N_Vector is changed. *
* *
* Routines that work with real ** begin with "band" (except for *
* the factor and solve routines which are called gbfa and gbsl, *
* respectively). The underlying matrix storage is described in *
* the documentation for bandalloc. *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _band_h
#define _band_h
#include "llnltyps.h"
#include "nvector.h"
/******************************************************************
* *
* Type: BandMat *
*----------------------------------------------------------------*
* The type BandMat is the type of a large (possibly distributed) *
* band matrix. It is defined to be a pointer to a structure *
* with the following fields: *
* *
* size is the number of columns (== number of rows) *
* *
* mu is the upper bandwidth, 0 <= mu <= size-1 *
* *
* ml is the lower bandwidth, 0 <= ml <= size-1 *
* *
* smu is the storage upper bandwidth, mu <= smu <= size-1. *
* The BandFactor routine writes the LU factors *
* into the storage for A. The upper triangular factor U, *
* however, may have an upper bandwidth as big as *
* MIN(size-1,mu+ml) because of partial pivoting. The smu *
* field holds the upper bandwidth allocated for A. *
* *
* data is a two dimensional array used for component storage. *
* The elements of a band matrix of type BandMat are *
* stored columnwise (i.e. columns are stored one on top *
* of the other in memory). Only elements within the *
* specified bandwidths are stored. *
* *
* If we number rows and columns in the band matrix starting *
* from 0, then *
* *
* data[0] is a pointer to (smu+ml+1)*size contiguous locations *
* which hold the elements within the band of A *
* *
* data[j] is a pointer to the uppermost element within the band *
* in the jth column. This pointer may be treated as *
* an array indexed from smu-mu (to access the *
* uppermost element within the band in the jth *
* column) to smu+ml (to access the lowest element *
* within the band in the jth column). (Indices from 0 *
* to smu-mu-1 give access to extra storage elements *
* required by BandFactor.) *
* *
* data[j][i-j+smu] is the (i,j)th element, j-mu <= i <= j+ml. *
* *
* The macros below allow a user to access individual matrix *
* elements without writing out explicit data structure *
* references and without knowing too much about the underlying *
* element storage. The only storage assumption needed is that *
* elements are stored columnwise and that a pointer into the jth *
* column of elements can be obtained via the BAND_COL macro. The *
* BAND_COL_ELEM macro selects an element from a column which has *
* already been isolated via BAND_COL. BAND_COL_ELEM allows the *
* user to avoid the translation from the matrix location (i,j) *
* to the index in the array returned by BAND_COL at which the *
* (i,j)th element is stored. See the documentation for BAND_COL *
* and BAND_COL_ELEM for usage details. Users should use these *
* macros whenever possible. *
* *
******************************************************************/
typedef struct {
integer size;
integer mu, ml, smu;
real **data;
} *BandMat;
/* BandMat accessor macros */
/******************************************************************
* *
* Macro : BAND_ELEM *
* Usage : BAND_ELEM(A,i,j) = a_ij; OR *
* a_ij = BAND_ELEM(A,i,j); *
*----------------------------------------------------------------*
* BAND_ELEM(A,i,j) references the (i,j)th element of the *
* N by N band matrix A, where 0 <= i,j <= N-1. The location *
* (i,j) should further satisfy j-(A->mu) <= i <= j+(A->ml). *
* *
******************************************************************/
#define BAND_ELEM(A,i,j) ((A->data)[j][i-j+(A->smu)])
/******************************************************************
* *
* Macro : BAND_COL *
* Usage : col_j = BAND_COL(A,j); *
*----------------------------------------------------------------*
* BAND_COL(A,j) references the diagonal element of the jth *
* column of the N by N band matrix A, 0 <= j <= N-1. The type of *
* the expression BAND_COL(A,j) is real *. The pointer returned *
* by the call BAND_COL(A,j) can be treated as an array which is *
* indexed from -(A->mu) to (A->ml). *
* *
******************************************************************/
#define BAND_COL(A,j) (((A->data)[j])+(A->smu))
/******************************************************************
* *
* Macro : BAND_COL_ELEM *
* Usage : col_j = BAND_COL(A,j); *
* BAND_COL_ELEM(col_j,i,j) = a_ij; OR *
* a_ij = BAND_COL_ELEM(col_j,i,j); *
*----------------------------------------------------------------*
* This macro references the (i,j)th entry of the band matrix A *
* when used in conjunction with BAND_COL as shown above. The *
* index (i,j) should satisfy j-(A->mu) <= i <= j+(A->ml). *
* *
******************************************************************/
#define BAND_COL_ELEM(col_j,i,j) (col_j[i-j])
/* Functions that use the BandMat representation for a band matrix */
/******************************************************************
* *
* Function : BandAllocMat *
* Usage : A = BandAllocMat(N, mu, ml, smu); *
* if (A == NULL) ... memory request failed *
*----------------------------------------------------------------*
* BandAllocMat allocates memory for an N by N band matrix with *
* upper bandwidth mu, lower bandwidth ml, and storage upper *
* bandwidth smu. Pass smu as follows depending on whether A will *
* be factored by BandFactor: *
* *
* (1) Pass smu = mu if A will not be factored. *
* *
* (2) Pass smu = MIN(N-1,mu+ml) if A will be factored. *
* *
* BandAllocMat returns the storage allocated (type BandMat) or *
* NULL if the request for matrix storage cannot be satisfied. *
* See the documentation for the type BandMat for matrix storage *
* details. *
* *
******************************************************************/
BandMat BandAllocMat(integer N, integer mu, integer ml, integer smu);
/******************************************************************
* *
* Function : BandAllocPiv *
* Usage : p = BandAllocPiv(N); *
* if (p == NULL) ... memory request failed *
*----------------------------------------------------------------*
* BandAllocPiv allocates memory for pivot information to be *
* filled in by the BandFactor routine during the factorization *
* of an N by N band matrix. The underlying type for pivot *
* information is an array of N integers and this routine returns *
* the pointer to the memory it allocates. If the request for *
* pivot storage cannot be satisfied, BandAllocPiv returns NULL. *
* *
******************************************************************/
integer *BandAllocPiv(integer N);
/******************************************************************
* *
* Function : BandFactor *
* Usage : ier = BandFactor(A, p); *
* if (ier != 0) ... A is singular *
*----------------------------------------------------------------*
* BandFactor performs the LU factorization of the N by N band *
* matrix A. This is done using standard Gaussian elimination *
* with partial pivoting. *
* *
* A successful LU factorization leaves the "matrix" A and the *
* pivot array p with the following information: *
* *
* (1) p[k] contains the row number of the pivot element chosen *
* at the beginning of elimination step k, k=0, 1, ..., N-1. *
* *
* (2) If the unique LU factorization of A is given by PA = LU, *
* where P is a permutation matrix, L is a lower triangular *
* matrix with all 1's on the diagonal, and U is an upper *
* triangular matrix, then the upper triangular part of A *
* (including its diagonal) contains U and the strictly lower *
* triangular part of A contains the multipliers, I-L. *
* *
* BandFactor returns 0 if successful. Otherwise it encountered *
* a zero diagonal element during the factorization. In this case *
* it returns the column index (numbered from one) at which *
* it encountered the zero. *
* *
* Important Note: A must be allocated to accommodate the increase*
* in upper bandwidth that occurs during factorization. If, *
* mathematically, A is a band matrix with upper bandwidth mu and *
* lower bandwidth ml, then the upper triangular factor U can *
* have upper bandwidth as big as smu=MIN(n-1,mu+ml). The lower *
* triangular factor L has lower bandwidth ml. Allocate A with *
* call A = BandAllocMat(N,mu,ml,smu), where mu, ml, and smu are *
* as defined above. The user does not have to zero the "extra" *
* storage allocated for the purpose of factorization. This will *
* handled by the BandFactor routine. *
* *
******************************************************************/
integer BandFactor(BandMat A, integer *p);
/******************************************************************
* *
* Function : BandBacksolve *
* Usage : BandBacksolve(A, p, b); *
*----------------------------------------------------------------*
* BandBacksolve solves the N-dimensional system A x = b using *
* the LU factorization in A and the pivot information in p *
* computed in BandFactor. The solution x is returned in b. This *
* routine cannot fail if the corresponding call to BandFactor *
* did not fail. *
* *
******************************************************************/
void BandBacksolve(BandMat A, integer *p, N_Vector b);
/******************************************************************
* *
* Function : BandZero *
* Usage : BandZero(A); *
*----------------------------------------------------------------*
* A(i,j) <- 0.0, j-(A->mu) <= i <= j+(A->ml). *
* *
******************************************************************/
void BandZero(BandMat A);
/******************************************************************
* *
* Function : BandCopy *
* Usage : BandCopy(A, B, copymu, copyml); *
*----------------------------------------------------------------*
* BandCopy copies the submatrix with upper and lower bandwidths *
* copymu, copyml of the N by N band matrix A into the N by N *
* band matrix B. *
* *
******************************************************************/
void BandCopy(BandMat A, BandMat B, integer copymu, integer copyml);
/******************************************************************
* *
* Function: BandScale *
* Usage : BandScale(c, A); *
*----------------------------------------------------------------*
* A(i,j) <- c*A(i,j), j-(A->mu) <= i <= j+(A->ml). *
* *
******************************************************************/
void BandScale(real c, BandMat A);
/******************************************************************
* *
* Function : BandAddI *
* Usage : BandAddI(A); *
*----------------------------------------------------------------*
* A(j,j) <- A(j,j)+1.0, 0 <= j <= (A->size)-1. *
* *
******************************************************************/
void BandAddI(BandMat A);
/******************************************************************
* *
* Function : BandFreeMat *
* Usage : BandFreeMat(A); *
*----------------------------------------------------------------*
* BandFreeMat frees the memory allocated by BandAllocMat for *
* the band matrix A. *
* *
******************************************************************/
void BandFreeMat(BandMat A);
/******************************************************************
* *
* Function : BandFreePiv *
* Usage : BandFreePiv(p); *
*----------------------------------------------------------------*
* BandFreePiv frees the memory allocated by BandAllocPiv for *
* the pivot information array p. *
* *
******************************************************************/
void BandFreePiv(integer *p);
/******************************************************************
* *
* Function : BandPrint *
* Usage : BandPrint(A); *
*----------------------------------------------------------------*
* This routine prints the N by N band matrix A (upper and lower *
* bandwidths A->mu and A->ml, respectively) to standard output *
* as it would normally appear on paper. It is intended as a *
* debugging tool with small values of N. The elements are *
* printed using the %g option. A blank line is printed before *
* and after the matrix. *
* *
******************************************************************/
void BandPrint(BandMat A);
/* Functions that use the real ** representation for a band matrix */
/******************************************************************
* *
* Function : bandalloc *
* Usage : real **a; *
* a = bandalloc(n, smu, ml); *
* if (a == NULL) ... memory request failed *
*----------------------------------------------------------------*
* bandalloc(n, smu, ml) allocates storage for an n by n band *
* matrix A with storage upper bandwidth smu and lower bandwidth *
* ml. It returns a pointer to the newly allocated storage if *
* successful. If the memory request cannot be satisfied, then *
* bandalloc returns NULL. If, mathematically, A has upper and *
* lower bandwidths mu and ml, respectively, then the value *
* passed to bandalloc for smu may need to be greater than mu. *
* The gbfa routine writes the LU factors into the storage (named *
* "a" in the above usage documentation) for A (thus destroying *
* the original elements of A). The upper triangular factor U, *
* however, may have a larger upper bandwidth than the upper *
* bandwidth mu of A. Thus some "extra" storage for A must be *
* allocated if A is to be factored by gbfa. Pass smu as follows: *
* *
* (1) Pass smu = mu if A will not be factored. *
* *
* (2) Pass smu = MIN(n-1,mu+ml) if A will be factored. *
* *
* The underlying type of the band matrix returned is real **. If *
* we allocate a band matrix A in real **a by *
* a = bandalloc(n,smu,ml), then a[0] is a pointer to *
* n * (smu + ml + 1) contiguous storage locations and a[j] is a *
* pointer to the uppermost element in the storage for the jth *
* column. The expression a[j][i-j+smu] references the (i,j)th *
* element of A, where 0 <= i,j <= n-1 and j-mu <= i <= j+ml. *
* (The elements a[j][0], a[j][1], ..., a[j][smu-mu-1] are used *
* by gbfa and gbsl.) *
* *
******************************************************************/
real **bandalloc(integer n, integer smu, integer ml);
/******************************************************************
* *
* Function : bandallocpiv *
* Usage : integer *pivot; *
* pivot = bandallocpiv(n); *
* if (pivot == NULL) ... memory request failed *
*----------------------------------------------------------------*
* bandallocpiv(n) allocates an array of n integers. It returns a *
* pointer to the first element in the array if successful. It *
* returns NULL if the memory request could not be satisfied. *
* *
******************************************************************/
integer *bandallocpiv(integer n);
/******************************************************************
* *
* Function : gbfa *
* Usage : integer ier; *
* ier = gbfa(a,n,mu,ml,smu,p); *
* if (ier > 0) ... zero element encountered during *
* the factorization *
*----------------------------------------------------------------*
* gbfa(a,n,mu,ml,smu,p) factors the n by n band matrix A (upper *
* and lower bandwidths mu and ml, storage upper bandwidth smu) *
* stored in "a". It overwrites the elements of A with the LU *
* factors and it keeps track of the pivot rows chosen in the *
* pivot array p. *
* *
* A successful LU factorization leaves a and pivot array p with *
* the following information: *
* *
* (1) p[k] contains the row number of the pivot element chosen *
* at the beginning of elimination step k, k=0, 1, ..., n-1. *
* *
* (2) If the unique LU factorization of A is given by PA = LU, *
* where P is a permutation matrix, L is a lower triangular *
* matrix with all 1's on the diagonal, and U is an upper *
* triangular matrix, then the upper triangular part of A *
* (including its diagonal) contains U and the strictly lower *
* triangular part of A contains the multipliers, I-L. *
* *
* gbfa returns 0 if successful. Otherwise it encountered a zero *
* diagonal element during the factorization. In this case it *
* returns the column index (numbered from one) at which it *
* encountered the zero. *
* *
* IMPORTANT NOTE: Suppose A is a band matrix with upper *
* bandwidth mu and lower bandwidth ml, then the upper triangular *
* factor U can have upper bandwidth as big as MIN(n-1,mu+ml) *
* because of partial pivoting. The lower triangular factor L has *
* lower bandwidth ml. Thus, if A is to be factored and *
* backsolved using gbfa and gbsl, then it should be allocated *
* as a = bandalloc(n,smu,ml), where smu = MIN(n-1,mu+ml). The *
* call to gbfa is ier = gbfa(a,n,mu,ml,smu,p). The corresponding *
* call to gbsl is gbsl(a,n,smu,ml,p,b). The user does not need *
* to zero the "extra" storage allocated for the purpose of *
* factorization. This is handled by the gbfa routine. If A is *
* not going to be factored and backsolved, then it can be *
* allocated as a = bandalloc(n,smu,ml). In either case, all *
* routines in this section use the parameter name smu for a *
* parameter which must be the "storage upper bandwidth" which *
* was passed to bandalloc. *
* *
******************************************************************/
integer gbfa(real **a, integer n, integer mu, integer ml, integer smu,
integer *p);
/******************************************************************
* *
* Function : gbsl *
* Usage : real *b; *
* ier = gbfa(a,n,mu,ml,smu,p); *
* if (ier == 0) gbsl(a,n,smu,ml,p,b); *
*----------------------------------------------------------------*
* gbsl(a,n,smu,ml,p,b) solves the n by n linear system *
* Ax = b, where A is band matrix stored in "a" with storage *
* upper bandwidth smu and lower bandwidth ml. It assumes that A *
* has been LU factored and the pivot array p has been set by a *
* successful call gbfa(a,n,mu,ml,smu,p). The solution x is *
* written into the b array. *
* *
******************************************************************/
void gbsl(real **a, integer n, integer smu, integer ml, integer *p, real *b);
/******************************************************************
* *
* Function : bandzero *
* Usage : bandzero(a,n,mu,ml,smu); *
*----------------------------------------------------------------*
* a(i,j) <- 0.0, 0 <= i,j <= n-1, j-mu <= i <= j+ml. *
* *
******************************************************************/
void bandzero(real **a, integer n, integer mu, integer ml, integer smu);
/******************************************************************
* *
* Function : bandcopy *
* Usage : bandcopy(a,b,n,a_smu,b_smu,copymu,copyml); *
*----------------------------------------------------------------*
* b(i,j) <- a(i,j), 0 <= i,j <= n-1, j-copymu <= i <= j+copyml. *
* *
******************************************************************/
void bandcopy(real **a, real **b, integer n, integer a_smu, integer b_smu,
integer copymu, integer copyml);
/******************************************************************
* *
* Function : bandscale *
* Usage : bandscale(c,a,n,mu,ml); *
*----------------------------------------------------------------*
* a(i,j) <- c*a(i,j), 0 <= i,j <= n-1, j-mu <= i <= j+ml. *
* *
******************************************************************/
void bandscale(real c, real **a, integer n, integer mu, integer ml,
integer smu);
/******************************************************************
* *
* Function : bandaddI *
* Usage : bandaddI(a,n,smu); *
*----------------------------------------------------------------*
* a(j,j) <- a(j,j)+1.0, 0 <= j <= n-1. *
* *
******************************************************************/
void bandaddI(real **a, integer n, integer smu);
/******************************************************************
* *
* Function : bandfreepiv *
* Usage : bandfreepiv(p); *
*----------------------------------------------------------------*
* bandfreepiv(p) frees the pivot array p allocated by *
* bandallocpiv. *
* *
******************************************************************/
void bandfreepiv(integer *p);
/******************************************************************
* *
* Function : bandfree *
* Usage : bandfree(a); *
*----------------------------------------------------------------*
* bandfree(a) frees the band matrix a allocated by bandalloc. *
* *
******************************************************************/
void bandfree(real **a);
/******************************************************************
* *
* Function : bandprint *
* Usage : bandprint(a,n,mu,ml,smu); *
*----------------------------------------------------------------*
* bandprint(a,n,mu,ml,smu) prints the n by n band matrix stored *
* in a (with upper bandwidth mu and lower bandwidth ml) to *
* standard output as it would normally appear on paper. It is *
* intended as a debugging tool with small values of n. The *
* elements are printed using the %g option. A blank line is *
* printed before and after the matrix. *
* *
******************************************************************/
void bandprint(real **a, integer n, integer mu, integer ml, integer smu);
#endif
#ifdef __cplusplus
}
#endif

221
ext/cvode/include/cvband.h
View file

@ -1,221 +0,0 @@
/******************************************************************
* *
* File : cvband.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 24 February 2000 *
*----------------------------------------------------------------*
* This is the header file for the CVODE band linear solver, *
* CVBAND. *
* *
* Note: The type integer must be large enough to store the value *
* N + mupper + mlower, where N is the linear system size and *
* mupper and mlower are the upper and lower bandwidths, *
* respectively, passed to CVBand. *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _cvband_h
#define _cvband_h
#include <stdio.h>
#include "cvode.h"
#include "llnltyps.h"
#include "band.h"
#include "nvector.h"
/******************************************************************
* *
* CVBAND solver statistics indices *
*----------------------------------------------------------------*
* The following enumeration gives a symbolic name to each *
* CVBAND statistic. The symbolic names are used as indices into *
* the iopt and ropt arrays passed to CVodeMalloc. *
* The CVBAND statistics are: *
* *
* iopt[BAND_NJE] : number of Jacobian evaluations, i.e. of *
* calls made to the band Jacobian routine *
* (default or user-supplied). *
* *
* iopt[BAND_LRW] : size (in real words) of real workspace *
* matrices and vectors used by this solver. *
* *
* iopt[BAND_LIW] : size (in integer words) of integer *
* workspace vectors used by this solver. *
* *
******************************************************************/
enum { BAND_NJE=CVODE_IOPT_SIZE, BAND_LRW, BAND_LIW };
/******************************************************************
* *
* CVBAND solver constants *
*----------------------------------------------------------------*
* CVB_MSBJ : maximum number of steps between band Jacobian *
* evaluations *
* *
* CVB_DGMAX : maximum change in gamma between band Jacobian *
* evaluations *
* *
******************************************************************/
#define CVB_MSBJ 50
#define CVB_DGMAX RCONST(0.2)
/******************************************************************
* *
* Type : CVBandJacFn *
*----------------------------------------------------------------*
* A band Jacobian approximation function Jac must have the *
* prototype given below. Its parameters are: *
* *
* N is the length of all vector arguments. *
* *
* mupper is the upper half-bandwidth of the approximate banded *
* Jacobian. This parameter is the same as the mupper parameter *
* passed by the user to the CVBand function. *
* *
* mlower is the lower half-bandwidth of the approximate banded *
* Jacobian. This parameter is the same as the mlower parameter *
* passed by the user to the CVBand function. *
* *
* J is the band matrix (of type BandMat) that will be loaded *
* by a CVBandJacFn with an approximation to the Jacobian matrix *
* J = (df_i/dy_j) at the point (t,y). *
* J is preset to zero, so only the nonzero elements need to be *
* loaded. Three efficient ways to load J are: *
* *
* (1) (with macros - no explicit data structure references) *
* for (j=0; j < N; j++) { *
* col_j = BAND_COL(J,j); *
* for (i=j-mupper; i <= j+mlower; i++) { *
* generate J_ij = the (i,j)th Jacobian element *
* BAND_COL_ELEM(col_j,i,j) = J_ij; *
* } *
* } *
* *
* (2) (with BAND_COL macro, but without BAND_COL_ELEM macro) *
* for (j=0; j < N; j++) { *
* col_j = BAND_COL(J,j); *
* for (k=-mupper; k <= mlower; k++) { *
* generate J_ij = the (i,j)th Jacobian element, i=j+k *
* col_j[k] = J_ij; *
* } *
* } *
* *
* (3) (without macros - explicit data structure references) *
* offset = J->smu; *
* for (j=0; j < N; j++) { *
* col_j = ((J->data)[j])+offset; *
* for (k=-mupper; k <= mlower; k++) { *
* generate J_ij = the (i,j)th Jacobian element, i=j+k *
* col_j[k] = J_ij; *
* } *
* } *
* Caution: J->smu is generally NOT the same as mupper. *
* *
* The BAND_ELEM(A,i,j) macro is appropriate for use in small *
* problems in which efficiency of access is NOT a major concern. *
* *
* f is the right hand side function for the ODE problem. *
* *
* f_data is a pointer to user data to be passed to f, the same *
* as the F_data parameter passed to CVodeMalloc. *
* *
* t is the current value of the independent variable. *
* *
* y is the current value of the dependent variable vector, *
* namely the predicted value of y(t). *
* *
* fy is the vector f(t,y). *
* *
* ewt is the error weight vector. *
* *
* h is a tentative step size in t. *
* *
* uround is the machine unit roundoff. *
* *
* jac_data is a pointer to user data - the same as the jac_data *
* parameter passed to CVBand. *
* *
* nfePtr is a pointer to the memory location containing the *
* CVODE problem data nfe = number of calls to f. The Jacobian *
* routine should update this counter by adding on the number *
* of f calls made in order to approximate the Jacobian, if any. *
* For example, if the routine calls f a total of N times, then *
* the update is *nfePtr += N. *
* *
* vtemp1, vtemp2, and vtemp3 are pointers to memory allocated *
* for vectors of length N which can be used by a CVBandJacFn *
* as temporary storage or work space. *
* *
******************************************************************/
typedef void (*CVBandJacFn)(integer N, integer mupper, integer mlower,
BandMat J, RhsFn f, void *f_data, real t,
N_Vector y, N_Vector fy, N_Vector ewt, real h,
real uround, void *jac_data, long int *nfePtr,
N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3);
/******************************************************************
* *
* Function : CVBand *
*----------------------------------------------------------------*
* A call to the CVBand function links the main CVODE integrator *
* with the CVBAND linear solver. *
* *
* cvode_mem is the pointer to CVODE memory returned by *
* CVodeMalloc. *
* *
* mupper is the upper bandwidth of the band Jacobian *
* approximation. *
* *
* mlower is the lower bandwidth of the band Jacobian *
* approximation. *
* *
* *
* bjac is the band Jacobian approximation routine to be used. *
* A user-supplied bjac routine must be of type *
* CVBandJacFn. Pass NULL for bjac to use the default *
* difference quotient routine CVBandDQJac supplied *
* with this solver. *
* *
* jac_data is a pointer to user data which is passed to the *
* bjac routine every time it is called. *
* *
******************************************************************/
void CVBand(void *cvode_mem, integer mupper, integer mlower, CVBandJacFn bjac,
void *jac_data);
/******************************************************************
* *
* Function : CVBandDQJac *
*----------------------------------------------------------------*
* This routine generates a banded difference quotient *
* approximation to the Jacobian of f(t,y). *
* *
******************************************************************/
void CVBandDQJac(integer N, integer mupper, integer mlower, BandMat J,
RhsFn f, void *f_data, real t, N_Vector y, N_Vector fy,
N_Vector ewt, real h, real uround, void *jac_data,
long int *nfePtr, N_Vector vtemp1, N_Vector vtemp2,
N_Vector vtemp3);
#endif
#ifdef __cplusplus
}
#endif

151
ext/cvode/include/cvbandpre.h
View file

@ -1,151 +0,0 @@
/******************************************************************
* *
* File : cvbandpre.h *
* Programmers : Michael Wittman and Alan C. Hindmarsh @ LLNL *
* Version of : 23 March 2000 *
*----------------------------------------------------------------*
* This is the header file for the CVBANDPRE module, which *
* provides a banded difference quotient Jacobian-based *
* preconditioner and solver routines for use with CVSPGMR. *
* *
* Summary: *
* These routines provide a band matrix preconditioner based on *
* difference quotients of the ODE right-hand side function f. *
* The user supplies parameters *
* mu = upper half-bandwidth (number of super-diagonals) *
* ml = lower half-bandwidth (number of sub-diagonals) *
* The routines generate a band matrix of bandwidth ml + mu + 1 *
* and use this to form a preconditioner for use with the Krylov *
* linear solver in CVSPGMR. Although this matrix is intended *
* to approximate the Jacobian df/dy, it may be a very crude *
* approximation. The true Jacobian need not be banded, or its *
* true bandwith may be larger than ml + mu + 1, as long as the *
* banded approximation generated here is sufficiently accurate *
* to speed convergence as a preconditioner. *
* *
* Usage: *
* The following is a summary of the usage of this module. *
* Details of the calls to CVodeMalloc, CVSpgmr, and CVode are *
* available in the CVODE User Guide. *
* To use these routines, the sequence of calls in the user *
* main program should be as follows: *
* *
* CVBandPreData bp_data; *
* ... *
* bp_data = CVBandPreAlloc(N, f, f_data, mu, ml); *
* ... *
* cvode_mem = CVodeMalloc(...); *
* ... *
* CVSpgmr(cvode_mem, pretype, gstype, maxl, delt, *
* CVBandPrecond, CVBandPSolve, bp_data); *
* ... *
* flag = CVode(...); *
* ... *
* CVBandPreFree(bp_data); *
* ... *
* CVodeFree(cvode_mem); *
* *
* Notes: *
* (1) Include this file for the CVBandPreData type definition. *
* (2) In the CVBandPreAlloc call, the arguments N, f, and f_data *
* are the same as in the call to CVodeMalloc. *
* (3) In the CVSpgmr call, the user is free to specify the inputs*
* pretype and gstype, and the optional inputs maxl and delt. *
* But the last three arguments must be as shown, with the *
* last argument being the pointer returned by CVBandPreAlloc.*
* (4) The CVBandPrecond and CVBandPSolve functions are never *
* called by the user explicitly; they are simply passed to *
* the CVSpgmr function. *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _cvbandpre_h
#define _cvbandpre_h
#include "cvode.h"
#include "llnltyps.h"
#include "nvector.h"
#include "band.h"
/************* CVBandPreData type definition ************/
typedef struct {
/* Data set by user in CVBandPreAlloc: */
RhsFn f;
void *f_data;
integer ml, mu;
/* Data set by CVBandPrecond: */
BandMat savedJ;
BandMat savedP;
integer *pivots;
} *CVBandPreData;
/******************************************************************
* *
* Function : CVBandPreAlloc *
*----------------------------------------------------------------*
* CVBandPreAlloc allocates and initializes a CVBandPreData *
* structure to be passed to CVSpgmr (and subsequently used by *
* CVBandPrecond and CVBandPSolve). *
* *
* The parameters of CVBandPreAlloc are as follows: *
* *
* N is the length of all vector arguments. *
* *
* f is the right hand side function. *
* *
* f_data is a pointer to the optional extra data for f. *
* *
* mu is the upper half bandwidth. *
* *
* ml is the lower half bandwidth. *
* *
* CVBandPreAlloc returns the storage pointer (type CVBandPreData)*
* or NULL if the request for storage cannot be satisfied. *
* *
******************************************************************/
CVBandPreData CVBandPreAlloc(integer N, RhsFn f, void *f_data,
integer mu, integer ml);
/******************************************************************
* *
* Function : CVBandPreFree *
*----------------------------------------------------------------*
* CVBandPreFree frees the memory allocated by CVBandPreAlloc in *
* the argument pdata. *
* *
******************************************************************/
void CVBandPreFree(CVBandPreData pdata);
/* Prototypes of CVBandPrecond and CVBandPSolve */
int CVBandPrecond(integer N, real t, N_Vector y, N_Vector fy, boole jok,
boole *jcurPtr, real gamma, N_Vector ewt, real h,
real uround, long int *nfePtr, void *bp_data,
N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3);
int CVBandPSolve(integer N, real t, N_Vector y, N_Vector fy, N_Vector vtemp,
real gamma, N_Vector ewt, real delta, long int *nfePtr,
N_Vector r, int lr, void *bp_data, N_Vector z);
#endif
#ifdef __cplusplus
}
#endif

191
ext/cvode/include/cvdense.h
View file

@ -1,191 +0,0 @@
/******************************************************************
* *
* File : cvdense.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 24 February 2000 *
*----------------------------------------------------------------*
* This is the header file for the CVODE dense linear solver, *
* CVDENSE. *
* *
* Note: The type integer must be large enough to store the value *
* of the linear system size N. *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _cvdense_h
#define _cvdense_h
#include <stdio.h>
#include "cvode.h"
#include "llnltyps.h"
#include "dense.h"
#include "nvector.h"
/******************************************************************
* *
* CVDENSE solver statistics indices *
*----------------------------------------------------------------*
* The following enumeration gives a symbolic name to each *
* CVDENSE statistic. The symbolic names are used as indices into *
* the iopt and ropt arrays passed to CVodeMalloc. *
* The CVDENSE statistics are: *
* *
* iopt[DENSE_NJE] : number of Jacobian evaluations, i.e. of *
* calls made to the dense Jacobian routine *
* (default or user-supplied). *
* *
* iopt[DENSE_LRW] : size (in real words) of real workspace *
* matrices and vectors used by this solver. *
* *
* iopt[DENSE_LIW] : size (in integer words) of integer *
* workspace vectors used by this solver. *
* *
******************************************************************/
enum { DENSE_NJE=CVODE_IOPT_SIZE, DENSE_LRW, DENSE_LIW };
/******************************************************************
* *
* CVDENSE solver constants *
*----------------------------------------------------------------*
* CVD_MSBJ : maximum number of steps between dense Jacobian *
* evaluations *
* *
* CVD_DGMAX : maximum change in gamma between dense Jacobian *
* evaluations *
* *
******************************************************************/
#define CVD_MSBJ 50
#define CVD_DGMAX RCONST(0.2)
/******************************************************************
* *
* Type : CVDenseJacFn *
*----------------------------------------------------------------*
* A dense Jacobian approximation function Jac must have the *
* prototype given below. Its parameters are: *
* *
* N is the length of all vector arguments. *
* *
* J is the dense matrix (of type DenseMat) that will be loaded *
* by a CVDenseJacFn with an approximation to the Jacobian matrix *
* J = (df_i/dy_j) at the point (t,y). *
* J is preset to zero, so only the nonzero elements need to be *
* loaded. Two efficient ways to load J are: *
* *
* (1) (with macros - no explicit data structure references) *
* for (j=0; j < N; j++) { *
* col_j = DENSE_COL(J,j); *
* for (i=0; i < N; i++) { *
* generate J_ij = the (i,j)th Jacobian element *
* col_j[i] = J_ij; *
* } *
* } *
* *
* (2) (without macros - explicit data structure references) *
* for (j=0; j < N; j++) { *
* col_j = (J->data)[j]; *
* for (i=0; i < N; i++) { *
* generate J_ij = the (i,j)th Jacobian element *
* col_j[i] = J_ij; *
* } *
* } *
* *
* The DENSE_ELEM(A,i,j) macro is appropriate for use in small *
* problems in which efficiency of access is NOT a major concern. *
* *
* f is the right hand side function for the ODE problem. *
* *
* f_data is a pointer to user data to be passed to f, the same *
* as the F_data parameter passed to CVodeMalloc. *
* *
* t is the current value of the independent variable. *
* *
* y is the current value of the dependent variable vector, *
* namely the predicted value of y(t). *
* *
* fy is the vector f(t,y). *
* *
* ewt is the error weight vector. *
* *
* h is a tentative step size in t. *
* *
* uround is the machine unit roundoff. *
* *
* jac_data is a pointer to user data - the same as the jac_data *
* parameter passed to CVDense. *
* *
* nfePtr is a pointer to the memory location containing the *
* CVODE problem data nfe = number of calls to f. The Jacobian *
* routine should update this counter by adding on the number *
* of f calls made in order to approximate the Jacobian, if any. *
* For example, if the routine calls f a total of N times, then *
* the update is *nfePtr += N. *
* *
* vtemp1, vtemp2, and vtemp3 are pointers to memory allocated *
* for vectors of length N which can be used by a CVDenseJacFn *
* as temporary storage or work space. *
* *
******************************************************************/
typedef void (*CVDenseJacFn)(integer N, DenseMat J, RhsFn f, void *f_data,
real t, N_Vector y, N_Vector fy, N_Vector ewt,
real h, real uround, void *jac_data,
long int *nfePtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3);
/******************************************************************
* *
* Function : CVDense *
*----------------------------------------------------------------*
* A call to the CVDense function links the main CVODE integrator *
* with the CVDENSE linear solver. *
* *
* cvode_mem is the pointer to CVODE memory returned by *
* CVodeMalloc. *
* *
* djac is the dense Jacobian approximation routine to be used. *
* A user-supplied djac routine must be of type *
* CVDenseJacFn. Pass NULL for djac to use the default *
* difference quotient routine CVDenseDQJac supplied *
* with this solver. *
* *
* jac_data is a pointer to user data which is passed to the *
* djac routine every time it is called. *
* *
******************************************************************/
void CVDense(void *cvode_mem, CVDenseJacFn djac, void *jac_data);
/******************************************************************
* *
* Function : CVDenseDQJac *
*----------------------------------------------------------------*
* This routine generates a dense difference quotient *
* approximation to the Jacobian of f(t,y). *
* *
******************************************************************/
void CVDenseDQJac(integer N, DenseMat J, RhsFn f, void *f_data, real t,
N_Vector y, N_Vector fy, N_Vector ewt, real h, real uround,
void *jac_data, long int *nfePtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3);
#endif
#ifdef __cplusplus
}
#endif

71
ext/cvode/include/cvdiag.h
View file

@ -1,71 +0,0 @@
/******************************************************************
* *
* File : cvdiag.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 4 May 1998 *
*----------------------------------------------------------------*
* This is the header file for the CVODE diagonal linear solver, *
* CVDIAG. *
* *
* Note: The type integer must be large enough to store the value *
* of the linear system size N. *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _cvdiag_h
#define _cvdiag_h
#include <stdio.h>
#include "cvode.h"
#include "llnltyps.h"
#include "nvector.h"
/******************************************************************
* *
* CVDIAG solver statistics indices *
*----------------------------------------------------------------*
* The following enumeration gives a symbolic name to each *
* CVDIAG statistic. The symbolic names are used as indices into *
* the iopt and ropt arrays passed to CVodeMalloc. *
* The CVDIAG statistics are: *
* *
* iopt[DIAG_LRW] : size (in real words) of real workspace *
* vectors used by this solver. *
* *
* iopt[DIAG_LIW] : size (in integer words) of integer *
* workspace vectors used by this solver. *
* *
* The number of diagonal approximate Jacobians formed is equal *
* to the number of CVDiagSetup calls. This number is available *
* in cv_iopt[NSETUPS]. *
* *
******************************************************************/
enum { DIAG_LRW=CVODE_IOPT_SIZE, DIAG_LIW };
/******************************************************************
* *
* Function : CVDiag *
*----------------------------------------------------------------*
* A call to the CVDiag function links the main CVODE integrator *
* with the CVDIAG linear solver. *
* *
* cvode_mem is the pointer to CVODE memory returned by *
* CVodeMalloc. *
* *
******************************************************************/
void CVDiag(void *cvode_mem);
#endif
#ifdef __cplusplus
}
#endif

823
ext/cvode/include/cvode.h
View file

@ -1,823 +0,0 @@
/******************************************************************
* *
* File : cvode.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 29 February 2000 *
*----------------------------------------------------------------*
* This is the interface file for the main CVODE integrator. *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _cvode_h
#define _cvode_h
#include <stdlib.h>
#include <stdio.h>
#include "llnltyps.h"
#include "nvector.h"
/******************************************************************
* *
* CVODE is used to solve numerically the ordinary initial value *
* problem : *
* *
* y' = f(t,y), *
* y(t0) = y0, *
* *
* where t0, y0 in R^N, and f: R x R^N -> R^N are given. *
* *
******************************************************************/
/******************************************************************
* *
* Enumerations for inputs to CVodeMalloc, CVReInit, and CVode. *
*----------------------------------------------------------------*
* Symbolic constants for the lmm, iter, and itol input *
* parameters to CVodeMalloc and CVReInit, as well as the input *
* parameter itask to CVode, are given below. *
* *
* lmm : The user of the CVODE package specifies whether to use *
* the ADAMS or BDF (backward differentiation formula) *
* linear multistep method. The BDF method is recommended *
* for stiff problems, and the ADAMS method is recommended *
* for nonstiff problems. *
* *
* iter : At each internal time step, a nonlinear equation must *
* be solved. The user can specify either FUNCTIONAL *
* iteration, which does not require linear algebra, or a *
* NEWTON iteration, which requires the solution of linear *
* systems. In the NEWTON case, the user also specifies a *
* CVODE linear solver. NEWTON is recommended in case of *
* stiff problems. *
* *
* itol : This parameter specifies the relative and absolute *
* tolerance types to be used. The SS tolerance type means *
* a scalar relative and absolute tolerance, while the SV *
* tolerance type means a scalar relative tolerance and a *
* vector absolute tolerance (a potentially different *
* absolute tolerance for each vector component). *
* *
* itask : The itask input parameter to CVode indicates the job *
* of the solver for the next user step. The NORMAL *
* itask is to have the solver take internal steps until *
* it has reached or just passed the user specified tout *
* parameter. The solver then interpolates in order to *
* return an approximate value of y(tout). The ONE_STEP *
* option tells the solver to just take one internal step *
* and return the solution at the point reached by that *
* step. *
* *
******************************************************************/
enum { ADAMS, BDF }; /* lmm */
enum { FUNCTIONAL, NEWTON }; /* iter */
enum { SS, SV }; /* itol */
enum { NORMAL, ONE_STEP }; /* itask */
/******************************************************************
* *
* Type : RhsFn *
*----------------------------------------------------------------*
* The f function which defines the right hand side of the ODE *
* system y' = f(t,y) must have type RhsFn. *
* f takes as input the problem size N, the independent variable *
* value t, and the dependent variable vector y. It stores the *
* result of f(t,y) in the vector ydot. The y and ydot arguments *
* are of type N_Vector. *
* (Allocation of memory for ydot is handled within CVODE.) *
* The f_data parameter is the same as the f_data *
* parameter passed by the user to the CVodeMalloc routine. This *
* user-supplied pointer is passed to the user's f function *
* every time it is called. *
* A RhsFn f does not have a return value. *
* *
******************************************************************/
typedef void (*RhsFn)(integer N, real t, N_Vector y, N_Vector ydot,
void *f_data);
/******************************************************************
* *
* Function : CVodeMalloc *
*----------------------------------------------------------------*
* CVodeMalloc allocates and initializes memory for a problem to *
* to be solved by CVODE. *
* *
* N is the number of equations in the ODE system. *
* *
* f is the right hand side function in y' = f(t,y). *
* *
* t0 is the initial value of t. *
* *
* y0 is the initial condition vector y(t0). *
* *
* lmm is the type of linear multistep method to be used. *
* The legal values are ADAMS and BDF (see previous *
* description). *
* *
* iter is the type of iteration used to solve the nonlinear *
* system that arises during each internal time step. *
* The legal values are FUNCTIONAL and NEWTON. *
* *
* itol is the type of tolerances to be used. *
* The legal values are: *
* SS (scalar relative and absolute tolerances), *
* SV (scalar relative tolerance and vector *
* absolute tolerance). *
* *
* reltol is a pointer to the relative tolerance scalar. *
* *
* abstol is a pointer to the absolute tolerance scalar or *
* an N_Vector of absolute tolerances. *
* *
* The parameters itol, reltol, and abstol define a vector of *
* error weights, ewt, with components *
* ewt[i] = 1/(reltol*abs(y[i]) + abstol) (if itol = SS), or *
* ewt[i] = 1/(reltol*abs(y[i]) + abstol[i]) (if itol = SV). *
* This vector is used in all error and convergence tests, which *
* use a weighted RMS norm on all error-like vectors v: *
* WRMSnorm(v) = sqrt( (1/N) sum(i=1..N) (v[i]*ewt[i])^2 ). *
* *
* f_data is a pointer to user data that will be passed to the *
* user's f function every time f is called. *
* *
* errfp is the file pointer for an error file where all CVODE *
* warning and error messages will be written. This *
* parameter can be stdout (standard output), stderr *
* (standard error), a file pointer (corresponding to *
* a user error file opened for writing) returned by *
* fopen, or NULL. If the user passes NULL, then all *
* messages will be written to standard output. *
* *
* optIn is a flag indicating whether there are any optional *
* inputs from the user in the arrays iopt and ropt. *
* Pass FALSE to indicate no optional inputs and TRUE *
* to indicate that optional inputs are present. *
* *
* iopt is the user-allocated array (of size OPT_SIZE given *
* later) that will hold optional integer inputs and *
* outputs. The user can pass NULL if he/she does not *
* wish to use optional integer inputs or outputs. *
* If optIn is TRUE, the user should preset to 0 those *
* locations for which default values are to be used. *
* *
* ropt is the user-allocated array (of size OPT_SIZE given *
* later) that will hold optional real inputs and *
* outputs. The user can pass NULL if he/she does not *
* wish to use optional real inputs or outputs. *
* If optIn is TRUE, the user should preset to 0.0 the *
* locations for which default values are to be used. *
* *
* machEnv is a pointer to machine environment-specific *
* information. *
* *
* Note: The tolerance values may be changed in between calls to *
* CVode for the same problem. These values refer to *
* (*reltol) and either (*abstol), for a scalar absolute *
* tolerance, or the components of abstol, for a vector *
* absolute tolerance. *
* *
* If successful, CVodeMalloc returns a pointer to initialized *
* problem memory. This pointer should be passed to CVode. If *
* an initialization error occurs, CVodeMalloc prints an error *
* message to the file specified by errfp and returns NULL. *
* *
******************************************************************/
void *CVodeMalloc(integer N, RhsFn f, real t0, N_Vector y0, int lmm, int iter,
int itol, real *reltol, void *abstol, void *f_data,
FILE *errfp, boole optIn, long int iopt[], real ropt[],
void *machEnv);
/******************************************************************
* *
* Function : CVReInit *
*----------------------------------------------------------------*
* CVReInit re-initializes CVode for the solution of a problem, *
* where a prior call to CVodeMalloc has been made with the same *
* problem size N. CVReInit performs the same input checking *
* and initializations that CVodeMalloc does (except for N). *
* But it does no memory allocation, assuming that the existing *
* internal memory is sufficient for the new problem. *
* *
* The use of CVReInit requires that the maximum method order, *
* maxord, is no larger for the new problem than for the problem *
* specified in the last call to CVodeMalloc. This condition is *
* automatically fulfilled if the multistep method parameter lmm *
* is unchanged (or changed from ADAMS to BDF) and the default *
* value for maxord is specified. *
* *
* The first argument to CVReInit is: *
* *
* cvode_mem = pointer to CVODE memory returned by CVodeMalloc. *
* *
* All the remaining arguments to CVReInit have names and *
* meanings identical to those of CVodeMalloc. Note that the *
* problem size N is not passed as an argument to CVReInit, *
* as that is assumed to unchanged since the CVodeMalloc call. *
* *
* The return value of CVReInit is equal to SUCCESS = 0 if there *
* were no errors; otherwise it is a negative int equal to: *
* CVREI_NO_MEM indicating cvode_mem was NULL, or *
* CVREI_ILL_INPUT indicating an input argument was illegal *
* (including an attempt to increase maxord). *
* In case of an error return, an error message is also printed. *
* *
* Note: the reported workspace sizes iopt[LENRW] and iopt[LENIW] *
* are left unchanged from the values computed by CVodeMalloc, and*
* so may be larger than would be computed for the new problem. *
******************************************************************/
int CVReInit(void *cvode_mem, RhsFn f, real t0, N_Vector y0,
int lmm, int iter, int itol, real *reltol, void *abstol,
void *f_data, FILE *errfp, boole optIn, long int iopt[],
real ropt[], void *machEnv);
/* CVReInit return values: */
/* SUCCESS = 0 (Defined under CVode return values, but listed
here also for completeness) */
enum {CVREI_NO_MEM = -1, CVREI_ILL_INPUT = -2};
/******************************************************************
* *
* Function : CVode *
*----------------------------------------------------------------*
* CVode integrates the ODE over an interval in t. *
* If itask is NORMAL, then the solver integrates from its *
* current internal t value to a point at or beyond tout, then *
* interpolates to t = tout and returns y(tout) in the user- *
* allocated vector yout. If itask is ONE_STEP, then the solver *
* takes one internal time step and returns in yout the value of *
* y at the new internal time. In this case, tout is used only *
* during the first call to CVode to determine the direction of *
* integration and the rough scale of the problem. In either *
* case, the time reached by the solver is placed in (*t). The *
* user is responsible for allocating the memory for this value. *
* *
* cvode_mem is the pointer to CVODE memory returned by *
* CVodeMalloc. *
* *
* tout is the next time at which a computed solution is desired *
* *
* yout is the computed solution vector. In NORMAL mode with no *
* errors, yout=y(tout). *
* *
* t is a pointer to a real location. CVode sets (*t) to the *
* time reached by the solver and returns yout=y(*t). *
* *
* itask is either NORMAL or ONE_STEP mode. These two modes have *
* described above. *
* *
* The return values for CVode are defined later in this file. *
* Here is a brief description of each return value: *
* *
* SUCCESS : CVode succeeded. *
* *
* CVODE_NO_MEM : The cvode_mem argument was NULL. *
* *
* ILL_INPUT : One of the inputs to CVode is illegal. This *
* includes the situation when a component of the *
* error weight vectors becomes < 0 during *
* internal time-stepping. The ILL_INPUT flag *
* will also be returned if the linear solver *
* routine CV--- (called by the user after *
* calling CVodeMalloc) failed to set one of the *
* linear solver-related fields in cvode_mem or *
* if the linear solver's init routine failed. In *
* any case, the user should see the printed *
* error message for more details. *
* *
* TOO_MUCH_WORK : The solver took mxstep internal steps but *
* could not reach tout. The default value for *
* mxstep is MXSTEP_DEFAULT = 500. *
* *
* TOO_MUCH_ACC : The solver could not satisfy the accuracy *
* demanded by the user for some internal step. *
* *
* ERR_FAILURE : Error test failures occurred too many times *
* (= MXNEF = 7) during one internal time step or *
* occurred with |h| = hmin. *
* *
* CONV_FAILURE : Convergence test failures occurred too many *
* times (= MXNCF = 10) during one internal time *
* step or occurred with |h| = hmin. *
* *
* SETUP_FAILURE : The linear solver's setup routine failed in an *
* unrecoverable manner. *
* *
* SOLVE_FAILURE : The linear solver's solve routine failed in an *
* unrecoverable manner. *
* *
******************************************************************/
int CVode(void *cvode_mem, real tout, N_Vector yout, real *t, int itask);
/* CVode return values */
enum { SUCCESS=0, CVODE_NO_MEM=-1, ILL_INPUT=-2, TOO_MUCH_WORK=-3,
TOO_MUCH_ACC=-4, ERR_FAILURE=-5, CONV_FAILURE=-6,
SETUP_FAILURE=-7, SOLVE_FAILURE=-8 };
/******************************************************************
* *
* Function : CVodeDky *
*----------------------------------------------------------------*
* CVodeDky computes the kth derivative of the y function at *
* time t, where tn-hu <= t <= tn, tn denotes the current *
* internal time reached, and hu is the last internal step size *
* successfully used by the solver. The user may request *
* k=0, 1, ..., qu, where qu is the current order. The *
* derivative vector is returned in dky. This vector must be *
* allocated by the caller. It is only legal to call this *
* function after a successful return from CVode. *
* *
* cvode_mem is the pointer to CVODE memory returned by *
* CVodeMalloc. *
* *
* t is the time at which the kth derivative of y is evaluated. *
* The legal range for t is [tn-hu,tn] as described above. *
* *
* k is the order of the derivative of y to be computed. The *
* legal range for k is [0,qu] as described above. *
* *
* dky is the output derivative vector [(D_k)y](t). *
* *
* The return values for CVodeDky are defined later in this file. *
* Here is a brief description of each return value: *
* *
* OKAY : CVodeDky succeeded. *
* *
* BAD_K : k is not in the range 0, 1, ..., qu. *
* *
* BAD_T : t is not in the interval [tn-hu,tn]. *
* *
* BAD_DKY : The dky argument was NULL. *
* *
* DKY_NO_MEM : The cvode_mem argument was NULL. *
* *
******************************************************************/
int CVodeDky(void *cvode_mem, real t, int k, N_Vector dky);
/* CVodeDky return values */
enum { OKAY=0, BAD_K=-1, BAD_T=-2, BAD_DKY=-3, DKY_NO_MEM=-4 };
/******************************************************************
* *
* Function : CVodeFree *
*----------------------------------------------------------------*
* CVodeFree frees the problem memory cvode_mem allocated by *
* CVodeMalloc. Its only argument is the pointer cvode_mem *
* returned by CVodeMalloc. *
* *
******************************************************************/
void CVodeFree(void *cvode_mem);
/******************************************************************
* *
* Optional Inputs and Outputs *
*----------------------------------------------------------------*
* The user should declare two arrays for optional input and *
* output, an iopt array for optional integer input and output *
* and an ropt array for optional real input and output. The *
* size of both these arrays should be OPT_SIZE. *
* So the user's declaration should look like: *
* *
* long int iopt[OPT_SIZE]; *
* real ropt[OPT_SIZE]; *
* *
* The enumerations below the OPT_SIZE definition *
* are indices into the iopt and ropt arrays. Here is a brief *
* description of the contents of these positions: *
* *
* iopt[MAXORD] : maximum lmm order to be used by the solver. *
* Optional input. (Default = 12 for ADAMS, 5 for *
* BDF). *
* *
* iopt[MXSTEP] : maximum number of internal steps to be taken by *
* the solver in its attempt to reach tout. *
* Optional input. (Default = 500). *
* *
* iopt[MXHNIL] : maximum number of warning messages issued *
* by the solver that t+h==t on the next internal *
* step. Optional input. (Default = 10). *
* *
* iopt[NST] : cumulative number of internal steps taken by *
* the solver (total so far). Optional output. *
* *
* iopt[NFE] : number of calls to the user's f function. *
* Optional output. *
* *
* iopt[NSETUPS] : number of calls made to the linear solver's *
* setup routine. Optional output. *
* *
* iopt[NNI] : number of NEWTON iterations performed. *
* Optional output. *
* *
* iopt[NCFN] : number of nonlinear convergence failures *
* that have occurred. Optional output. *
* *
* iopt[NETF] : number of local error test failures that *
* have occurred. Optional output. *
* *
* iopt[QU] : order used during the last internal step. *
* Optional output. *
* *
* iopt[QCUR] : order to be used on the next internal step. *
* Optional output. *
* *
* iopt[LENRW] : size of required CVODE internal real work *
* space, in real words. Optional output. *
* *
* iopt[LENIW] : size of required CVODE internal integer work *
* space, in integer words. Optional output. *
* *
* ropt[H0] : initial step size. Optional input. *
* *
* ropt[HMAX] : maximum absolute value of step size allowed. *
* Optional input. (Default is infinity). *
* *
* ropt[HMIN] : minimum absolute value of step size allowed. *
* Optional input. (Default is 0.0). *
* *
* ropt[HU] : step size for the last internal step. *
* Optional output. *
* *
* ropt[HCUR] : step size to be attempted on the next internal *
* step. Optional output. *
* *
* ropt[TCUR] : current internal time reached by the solver. *
* Optional output. *
* *
* ropt[TOLSF] : a suggested factor by which the user's *
* tolerances should be scaled when too much *
* accuracy has been requested for some internal *
* step. Optional output. *
* *
******************************************************************/
/* iopt, ropt array sizes */
#define OPT_SIZE 40
/* iopt and ropt offsets *
* The constants CVODE_IOPT_SIZE and CVODE_ROPT_SIZE are equal to *
* the number of integer and real optional inputs and outputs *
* actually accessed in cvode.c. The locations beyond these *
* values are used by the linear solvers. */
#define CVODE_IOPT_SIZE 13
#define CVODE_ROPT_SIZE 7
/* iopt indices */
enum { MAXORD, MXSTEP, MXHNIL,
NST, NFE, NSETUPS, NNI, NCFN, NETF, QU, QCUR,
LENRW, LENIW };
/* ropt indices */
enum { H0, HMAX, HMIN,
HU, HCUR, TCUR, TOLSF };
/* Basic CVODE constants */
#define ADAMS_Q_MAX 12 /* max value of q for lmm == ADAMS */
#define BDF_Q_MAX 5 /* max value of q for lmm == BDF */
#define Q_MAX ADAMS_Q_MAX /* max value of q for either lmm */
#define L_MAX (Q_MAX+1) /* max value of L for either lmm */
#define NUM_TESTS 5 /* number of error test quantities */
/******************************************************************
* *
* Types : struct CVodeMemRec, CVodeMem *
*----------------------------------------------------------------*
* The type CVodeMem is type pointer to struct CVodeMemRec. This *
* structure contains fields to keep track of problem state. *
* *
******************************************************************/
typedef struct CVodeMemRec {
real cv_uround; /* machine unit roundoff */
/* Problem Specification Data */
integer cv_N; /* ODE system size */
RhsFn cv_f; /* y' = f(t,y(t)) */
void *cv_f_data; /* user pointer passed to f */
int cv_lmm; /* lmm = ADAMS or BDF */
int cv_iter; /* iter = FUNCTIONAL or NEWTON */
int cv_itol; /* itol = SS or SV */
real *cv_reltol; /* ptr to relative tolerance */
void *cv_abstol; /* ptr to absolute tolerance */
/* Nordsieck History Array */
N_Vector cv_zn[L_MAX]; /* Nordsieck array N x (q+1), */
/* zn[j] is a vector of length N, j=0, ... , q */
/* zn[j] = h^j * jth derivative of the */
/* interpolating polynomial */
/* Vectors of length N */
N_Vector cv_ewt; /* error weight vector */
N_Vector cv_y; /* y is used as temporary storage by the solver */
/* The memory is provided by the user to CVode */
/* where the vector is named yout. */
N_Vector cv_acor; /* In the context of the solution of the */
/* nonlinear equation, acor = y_n(m) - y_n(0). */
/* On return, this vector is scaled to give */
/* the estimated local error in y. */
N_Vector cv_tempv; /* temporary storage vector */
N_Vector cv_ftemp; /* temporary storage vector */
/* Step Data */
int cv_q; /* current order */
int cv_qprime; /* order to be used on the next step */
/* = q-1, q, or q+1 */
int cv_qwait; /* number of internal steps to wait before */
/* considering a change in q */
int cv_L; /* L = q + 1 */
real cv_h; /* current step size */
real cv_hprime; /* step size to be used on the next step */
real cv_eta; /* eta = hprime / h */
real cv_hscale; /* value of h used in zn */
real cv_tn; /* current internal value of t */
real cv_tau[L_MAX+1]; /* array of previous q+1 successful step */
/* sizes indexed from 1 to q+1 */
real cv_tq[NUM_TESTS+1]; /* array of test quantities indexed from */
/* 1 to NUM_TESTS(=5) */
real cv_l[L_MAX]; /* coefficients of l(x) (degree q poly) */
real cv_rl1; /* 1 / l[1] */
real cv_gamma; /* gamma = h * rl1 */
real cv_gammap; /* gamma at the last setup call */
real cv_gamrat; /* gamma / gammap */
real cv_crate; /* estimated corrector convergence rate */
real cv_acnrm; /* | acor | wrms */
int cv_mnewt; /* Newton iteration counter */
/* Limits */
int cv_qmax; /* q <= qmax */
int cv_mxstep; /* maximum number of internal steps for one user call */
int cv_maxcor; /* maximum number of corrector iterations for the */
/* solution of the nonlinear equation */
int cv_mxhnil; /* maximum number of warning messages issued to the */
/* user that t + h == t for the next internal step */
real cv_hmin; /* |h| >= hmin */
real cv_hmax_inv; /* |h| <= 1/hmax_inv */
real cv_etamax; /* eta <= etamax */
/* Counters */
long int cv_nst; /* number of internal steps taken */
long int cv_nfe; /* number of f calls */
long int cv_ncfn; /* number of corrector convergence failures */
long int cv_netf; /* number of error test failures */
long int cv_nni; /* number of Newton iterations performed */
long int cv_nsetups; /* number of setup calls */
int cv_nhnil; /* number of messages issued to the user that */
/* t + h == t for the next iternal step */
long int cv_lrw; /* number of real words in CVODE work vectors */
long int cv_liw; /* no. of integer words in CVODE work vectors */
/* Linear Solver Data */
/* Linear Solver functions to be called */
int (*cv_linit)(struct CVodeMemRec *cv_mem, boole *setupNonNull);
int (*cv_lsetup)(struct CVodeMemRec *cv_mem, int convfail, N_Vector ypred,
N_Vector fpred, boole *jcurPtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3);
int (*cv_lsolve)(struct CVodeMemRec *cv_mem, N_Vector b, N_Vector ycur,
N_Vector fcur);
void (*cv_lfree)(struct CVodeMemRec *cv_mem);
/* Linear Solver specific memory */
void *cv_lmem;
/* Flag to indicate successful cv_linit call */
boole cv_linitOK;
/* Saved Values */
int cv_qu; /* last successful q value used */
long int cv_nstlp; /* step number of last setup call */
real cv_hu; /* last successful h value used */
real cv_saved_tq5; /* saved value of tq[5] */
boole cv_jcur; /* Is the Jacobian info used by */
/* linear solver current? */
real cv_tolsf; /* tolerance scale factor */
boole cv_setupNonNull;/* Does setup do something? */
/* Arrays for Optional Input and Optional Output */
long int *cv_iopt; /* long int optional input, output */
real *cv_ropt; /* real optional input, output */
/* Error File */
FILE *cv_errfp; /* CVODE error messages are sent to errfp */
/* Pointer to Machine Environment-Specific Information */
void *cv_machenv;
} *CVodeMem;
/******************************************************************
* *
* Communication between cvode.c and a CVODE Linear Solver *
*----------------------------------------------------------------*
* (1) cv_linit return values *
* *
* LINIT_OK : The cv_linit routine succeeded. *
* *
* LINIT_ERR : The cv_linit routine failed. Each linear solver *
* init routine should print an appropriate error *
* message to (cv_mem->errfp). *
* *
* (2) convfail (input to cv_lsetup) *
* *
* NO_FAILURES : Either this is the first cv_setup call for this *
* step, or the local error test failed on the *
* previous attempt at this step (but the Newton *
* iteration converged). *
* *
* FAIL_BAD_J : This value is passed to cv_lsetup if *
* *
* (1) The previous Newton corrector iteration *
* did not converge and the linear solver's *
* setup routine indicated that its Jacobian- *
* related data is not current. *
* or *
* (2) During the previous Newton corrector *
* iteration, the linear solver's solve routine *
* failed in a recoverable manner and the *
* linear solver's setup routine indicated that *
* its Jacobian-related data is not current. *
* *
* FAIL_OTHER : During the current internal step try, the *
* previous Newton iteration failed to converge *
* even though the linear solver was using current *
* Jacobian-related data. *
* *
* (3) Parameter documentation, as well as a brief description *
* of purpose, for each CVODE linear solver routine to be *
* called in cvode.c is given below the constant declarations *
* that follow. *
* *
******************************************************************/
/* cv_linit return values */
#define LINIT_OK 0
#define LINIT_ERR -1
/* Constants for convfail (input to cv_lsetup) */
#define NO_FAILURES 0
#define FAIL_BAD_J 1
#define FAIL_OTHER 2
/*******************************************************************
* *
* int (*cv_linit)(CVodeMem cv_mem, boole *setupNonNull); *
*-----------------------------------------------------------------*
* The purpose of cv_linit is to allocate memory for the *
* solver-specific fields in the structure *(cv_mem->cv_lmem) and *
* perform any needed initializations of solver-specific memory, *
* such as counters/statistics. The cv_linit routine should set *
* *setupNonNull to be TRUE if the setup operation for the linear *
* solver is non-empty and FALSE if the setup operation does *
* nothing. An LInitFn should return LINIT_OK (== 0) if it has *
* successfully initialized the CVODE linear solver and LINIT_ERR *
* (== -1) otherwise. These constants are defined above. If an *
* error does occur, an appropriate message should be sent to *
* (cv_mem->errfp). *
* *
*******************************************************************/
/*******************************************************************
* *
* int (*cv_lsetup)(CVodeMem cv_mem, int convfail, N_Vector ypred, *
* N_Vector fpred, boole *jcurPtr, N_Vector vtemp1,*
* N_Vector vtemp2, N_Vector vtemp3); *
*-----------------------------------------------------------------*
* The job of cv_lsetup is to prepare the linear solver for *
* subsequent calls to cv_lsolve. It may re-compute Jacobian- *
* related data is it deems necessary. Its parameters are as *
* follows: *
* *
* cv_mem - problem memory pointer of type CVodeMem. See the big *
* typedef earlier in this file. *
* *
* convfail - a flag to indicate any problem that occurred during *
* the solution of the nonlinear equation on the *
* current time step for which the linear solver is *
* being used. This flag can be used to help decide *
* whether the Jacobian data kept by a CVODE linear *
* solver needs to be updated or not. *
* Its possible values have been documented above. *
* *
* ypred - the predicted y vector for the current CVODE internal *
* step. *
* *
* fpred - f(tn, ypred). *
* *
* jcurPtr - a pointer to a boolean to be filled in by cv_lsetup. *
* The function should set *jcurPtr=TRUE if its Jacobian *
* data is current after the call and should set *
* *jcurPtr=FALSE if its Jacobian data is not current. *
* Note: If cv_lsetup calls for re-evaluation of *
* Jacobian data (based on convfail and CVODE state *
* data), it should return *jcurPtr=TRUE unconditionally;*
* otherwise an infinite loop can result. *
* *
* vtemp1 - temporary N_Vector provided for use by cv_lsetup. *
* *
* vtemp3 - temporary N_Vector provided for use by cv_lsetup. *
* *
* vtemp3 - temporary N_Vector provided for use by cv_lsetup. *
* *
* The cv_lsetup routine should return 0 if successful, *
* a positive value for a recoverable error, and a negative value *
* for an unrecoverable error. *
* *
*******************************************************************/
/*******************************************************************
* *
* int (*cv_lsolve)(CVodeMem cv_mem, N_Vector b, N_Vector ycur, *
* N_Vector fcur); *
*-----------------------------------------------------------------*
* cv_lsolve must solve the linear equation P x = b, where *
* P is some approximation to (I - gamma J), J = (df/dy)(tn,ycur) *
* and the RHS vector b is input. The N-vector ycur contains *
* the solver's current approximation to y(tn) and the vector *
* fcur contains the N-vector f(tn,ycur). The solution is to be *
* returned in the vector b. cv_lsolve returns a positive value *
* for a recoverable error and a negative value for an *
* unrecoverable error. Success is indicated by a 0 return value. *
* *
*******************************************************************/
/*******************************************************************
* *
* void (*cv_lfree)(CVodeMem cv_mem); *
*-----------------------------------------------------------------*
* cv_lfree should free up any memory allocated by the linear *
* solver. This routine is called once a problem has been *
* completed and the linear solver is no longer needed. *
* *
*******************************************************************/
#endif
#ifdef __cplusplus
}
#endif

336
ext/cvode/include/cvspgmr.h
View file

@ -1,336 +0,0 @@
/******************************************************************
* *
* File : cvspgmr.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 4 May 1998 *
*----------------------------------------------------------------*
* This is the header file for the CVODE scaled, preconditioned *
* GMRES linear solver, CVSPGMR. *
* *
* Note: The type integer must be large enough to store the value *
* of the linear system size N. *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _cvspgmr_h
#define _cvspgmr_h
#include <stdio.h>
#include "cvode.h"
#include "spgmr.h"
#include "llnltyps.h"
#include "nvector.h"
/******************************************************************
* *
* CVSPGMR solver statistics indices *
*----------------------------------------------------------------*
* The following enumeration gives a symbolic name to each *
* CVSPGMR statistic. The symbolic names are used as indices into *
* the iopt and ropt arrays passed to CVodeMalloc. *
* The CVSPGMR statistics are: *
* *
* iopt[SPGMR_NPE] : number of preconditioner evaluations, *
* i.e. of calls made to user's precond *
* function with jok == FALSE. *
* *
* iopt[SPGMR_NLI] : number of linear iterations. *
* *
* iopt[SPGMR_NPS] : number of calls made to user's psolve *
* function. *
* *
* iopt[SPGMR_NCFL] : number of linear convergence failures. *
* *
* iopt[SPGMR_LRW] : size (in real words) of real workspace *
* vectors and small matrices used by this *
* solver. *
* *
* iopt[SPGMR_LIW] : size (in integer words) of integer *
* workspace vectors used by this solver. *
* *
******************************************************************/
enum { SPGMR_NPE = CVODE_IOPT_SIZE,
SPGMR_NLI, SPGMR_NPS, SPGMR_NCFL, SPGMR_LRW, SPGMR_LIW };
/******************************************************************
* *
* CVSPGMR solver constants *
*----------------------------------------------------------------*
* CVSPGMR_MAXL : default value for the maximum Krylov *
* dimension is MIN(N, CVSPGMR_MAXL) *
* *
* CVSPGMR_MSBPRE : maximum number of steps between *
* preconditioner evaluations *
* *
* CVSPGMR_DGMAX : maximum change in gamma between *
* preconditioner evaluations *
* *
* CVSPGMR_DELT : default value for factor by which the *
* tolerance on the nonlinear iteration is *
* multiplied to get a tolerance on the linear *
* iteration *
* *
******************************************************************/
#define CVSPGMR_MAXL 5
#define CVSPGMR_MSBPRE 50
#define CVSPGMR_DGMAX RCONST(0.2)
#define CVSPGMR_DELT RCONST(0.05)
/******************************************************************
* *
* Type : CVSpgmrPrecondFn *
*----------------------------------------------------------------*
* The user-supplied preconditioner setup function Precond and *
* the user-supplied preconditioner solve function PSolve *
* together must define left and right preconditoner matrices *
* P1 and P2 (either of which may be trivial), such that the *
* product P1*P2 is an approximation to the Newton matrix *
* M = I - gamma*J. Here J is the system Jacobian J = df/dy, *
* and gamma is a scalar proportional to the integration step *
* size h. The solution of systems P z = r, with P = P1 or P2, *
* is to be carried out by the PSolve function, and Precond is *
* to do any necessary setup operations. *
* *
* The user-supplied preconditioner setup function Precond *
* is to evaluate and preprocess any Jacobian-related data *
* needed by the preconditioner solve function PSolve. *
* This might include forming a crude approximate Jacobian, *
* and performing an LU factorization on the resulting *
* approximation to M. This function will not be called in *
* advance of every call to PSolve, but instead will be called *
* only as often as necessary to achieve convergence within the *
* Newton iteration in CVODE. If the PSolve function needs no *
* preparation, the Precond function can be NULL. *
* *
* For greater efficiency, the Precond function may save *
* Jacobian-related data and reuse it, rather than generating it *
* from scratch. In this case, it should use the input flag jok *
* to decide whether to recompute the data, and set the output *
* flag *jcurPtr accordingly. *
* *
* Each call to the Precond function is preceded by a call to *
* the RhsFn f with the same (t,y) arguments. Thus the Precond *
* function can use any auxiliary data that is computed and *
* saved by the f function and made accessible to Precond. *
* *
* The error weight vector ewt, step size h, and unit roundoff *
* uround are provided to the Precond function for possible use *
* in approximating Jacobian data, e.g. by difference quotients. *
* *
* A function Precond must have the prototype given below. *
* Its parameters are as follows: *
* *
* N is the length of all vector arguments. *
* *
* t is the current value of the independent variable. *
* *
* y is the current value of the dependent variable vector, *
* namely the predicted value of y(t). *
* *
* fy is the vector f(t,y). *
* *
* jok is an input flag indicating whether Jacobian-related *
* data needs to be recomputed, as follows: *
* jok == FALSE means recompute Jacobian-related data *
* from scratch. *
* jok == TRUE means that Jacobian data, if saved from *
* the previous Precond call, can be reused *
* (with the current value of gamma). *
* A Precond call with jok == TRUE can only occur after *
* a call with jok == FALSE. *
* *
* jcurPtr is a pointer to an output integer flag which is *
* to be set by Precond as follows: *
* Set *jcurPtr = TRUE if Jacobian data was recomputed. *
* Set *jcurPtr = FALSE if Jacobian data was not *
* recomputed, but saved data was reused. *
* *
* gamma is the scalar appearing in the Newton matrix. *
* *
* ewt is the error weight vector. *
* *
* h is a tentative step size in t. *
* *
* uround is the machine unit roundoff. *
* *
* nfePtr is a pointer to the memory location containing the *
* CVODE problem data nfe = number of calls to f. *
* The Precond routine should update this counter by *
* adding on the number of f calls made in order to *
* approximate the Jacobian, if any. For example, if *
* the routine calls f a total of W times, then the *
* update is *nfePtr += W. *
* *
* P_data is a pointer to user data - the same as the P_data *
* parameter passed to CVSpgmr. *
* *
* vtemp1, vtemp2, and vtemp3 are pointers to memory allocated *
* for vectors of length N which can be used by *
* CVSpgmrPrecondFn as temporary storage or work space. *
* *
* *
* Returned value: *
* The value to be returned by the Precond function is a flag *
* indicating whether it was successful. This value should be *
* 0 if successful, *
* > 0 for a recoverable error (step will be retried), *
* < 0 for an unrecoverable error (integration is halted). *
* *
******************************************************************/
typedef int (*CVSpgmrPrecondFn)(integer N, real t, N_Vector y, N_Vector fy,
boole jok, boole *jcurPtr, real gamma,
N_Vector ewt, real h, real uround,
long int *nfePtr, void *P_data,
N_Vector vtemp1, N_Vector vtemp2,
N_Vector vtemp3);
/******************************************************************
* *
* Type : CVSpgmrPSolveFn *
*----------------------------------------------------------------*
* The user-supplied preconditioner solve function PSolve *
* is to solve a linear system P z = r in which the matrix P is *
* one of the preconditioner matrices P1 or P2, depending on the *
* type of preconditioning chosen. *
* *
* A function PSolve must have the prototype given below. *
* Its parameters are as follows: *
* *
* N is the length of all vector arguments. *
* *
* t is the current value of the independent variable. *
* *
* y is the current value of the dependent variable vector. *
* *
* fy is the vector f(t,y). *
* *
* vtemp is a pointer to memory allocated for a vector of *
* length N which can be used by PSolve for work space. *
* *
* gamma is the scalar appearing in the Newton matrix. *
* *
* ewt is the error weight vector (input). See delta below. *
* *
* delta is an input tolerance for use by PSolve if it uses *
* an iterative method in its solution. In that case, *
* the residual vector Res = r - P z of the system *
* should be made less than delta in weighted L2 norm, *
* i.e., sqrt [ Sum (Res[i]*ewt[i])^2 ] < delta . *
* *
* nfePtr is a pointer to the memory location containing the *
* CVODE problem data nfe = number of calls to f. The *
* PSolve routine should update this counter by adding *
* on the number of f calls made in order to carry out *
* the solution, if any. For example, if the routine *
* calls f a total of W times, then the update is *
* *nfePtr += W. *
* *
* r is the right-hand side vector of the linear system. *
* *
* lr is an input flag indicating whether PSolve is to use *
* the left preconditioner P1 or right preconditioner *
* P2: lr = 1 means use P1, and lr = 2 means use P2. *
* *
* P_data is a pointer to user data - the same as the P_data *
* parameter passed to CVSpgmr. *
* *
* z is the output vector computed by PSolve. *
* *
* Returned value: *
* The value to be returned by the PSolve function is a flag *
* indicating whether it was successful. This value should be *
* 0 if successful, *
* positive for a recoverable error (step will be retried), *
* negative for an unrecoverable error (integration is halted). *
* *
******************************************************************/
typedef int (*CVSpgmrPSolveFn)(integer N, real t, N_Vector y, N_Vector fy,
N_Vector vtemp, real gamma, N_Vector ewt,
real delta, long int *nfePtr, N_Vector r,
int lr, void *P_data, N_Vector z);
/******************************************************************
* *
* Function : CVSpgmr *
*----------------------------------------------------------------*
* A call to the CVSpgmr function links the main CVODE integrator *
* with the CVSPGMR linear solver. *
* *
* cvode_mem is the pointer to CVODE memory returned by *
* CVodeMalloc. *
* *
* pretype is the type of user preconditioning to be done. *
* This must be one of the four enumeration constants *
* NONE, LEFT, RIGHT, or BOTH defined in iterativ.h. *
* These correspond to no preconditioning, *
* left preconditioning only, right preconditioning *
* only, and both left and right preconditioning, *
* respectively. *
* *
* gstype is the type of Gram-Schmidt orthogonalization to be *
* used. This must be one of the two enumeration *
* constants MODIFIED_GS or CLASSICAL_GS defined in *
* iterativ.h. These correspond to using modified *
* Gram-Schmidt and classical Gram-Schmidt, *
* respectively. *
* *
* maxl is the maximum Krylov dimension. This is an *
* optional input to the CVSPGMR solver. Pass 0 to *
* use the default value MIN(N, CVSPGMR_MAXL=5). *
* *
* delt is the factor by which the tolerance on the *
* nonlinear iteration is multiplied to get a *
* tolerance on the linear iteration. This is an *
* optional input to the CVSPGMR solver. Pass 0 to *
* use the default value CVSPGMR_DELT = 0.05. *
* *
* precond is the user's preconditioner routine. It is used to *
* evaluate and preprocess any Jacobian-related data *
* needed by the psolve routine. See the *
* documentation for the type CVSpgmrPrecondFn for *
* full details. Pass NULL if no such setup of *
* Jacobian data is required. A precond routine is *
* NOT required for any of the four possible values *
* of pretype. *
* *
* psolve is the user's preconditioner solve routine. It is *
* used to solve Pz=r, where P is a preconditioner *
* matrix. See the documentation for the type *
* CVSpgmrPSolveFn for full details. The only case *
* in which psolve is allowed to be NULL is when *
* pretype is NONE. A valid psolve function must be *
* supplied when any preconditioning is to be done. *
* *
* P_data is a pointer to user preconditioner data. This *
* pointer is passed to precond and psolve every time *
* these routines are called. *
* *
******************************************************************/
void CVSpgmr(void *cvode_mem, int pretype, int gstype, int maxl, real delt,
CVSpgmrPrecondFn precond, CVSpgmrPSolveFn psolve, void *P_data);
#endif
#ifdef __cplusplus
}
#endif

494
ext/cvode/include/dense.h
View file

@ -1,494 +0,0 @@
/******************************************************************
* *
* File : dense.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 6 May 1998 *
*----------------------------------------------------------------*
* This is the header file for a generic DENSE linear solver *
* package. There are two sets of dense solver routines listed in *
* this file: one set uses type DenseMat defined below and the *
* other set uses the type real ** for dense matrix arguments. *
* The two sets of dense solver routines make it easy to work *
* with two types of dense matrices: *
* *
* (1) The DenseMat type is intended for use with large dense *
* matrices whose elements/columns may be stored in *
* non-contiguous memory locations or even distributed into *
* different processor memories. This type may be modified to *
* include such distribution information. If this is done, *
* then all the routines that use DenseMat must be modified *
* to reflect the new data structure. *
* *
* (2) The set of routines that use real ** (and NOT the DenseMat *
* type) is intended for use with small matrices which can *
* easily be allocated within a contiguous block of memory *
* on a single processor. *
* *
* Routines that work with the type DenseMat begin with "Dense". *
* The DenseAllocMat function allocates a dense matrix for use in *
* the other DenseMat routines listed in this file. Matrix *
* storage details are given in the documentation for the type *
* DenseMat. The DenseAllocPiv function allocates memory for *
* pivot information. The storage allocated by DenseAllocMat and *
* DenseAllocPiv is deallocated by the routines DenseFreeMat and *
* DenseFreePiv, respectively. The DenseFactor and DenseBacksolve *
* routines perform the actual solution of a dense linear system. *
* Note that the DenseBacksolve routine has a parameter b of type *
* N_Vector. The current implementation makes use of a machine *
* environment-specific macro (N_VDATA) which may not exist for *
* other implementations of the type N_Vector. Thus, the *
* implementation of DenseBacksolve may need to change if the *
* type N_Vector is changed. *
* *
* Routines that work with real ** begin with "den" (except for *
* the factor and solve routines which are called gefa and gesl, *
* respectively). The underlying matrix storage is described in *
* the documentation for denalloc. *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _dense_h
#define _dense_h
#include "llnltyps.h"
#include "nvector.h"
/******************************************************************
* *
* Type: DenseMat *
*----------------------------------------------------------------*
* The type DenseMat is defined to be a pointer to a structure *
* with a size and a data field. The size field indicates the *
* number of columns (== number of rows) of a dense matrix, while *
* the data field is a two dimensional array used for component *
* storage. The elements of a dense matrix are stored columnwise *
* (i.e columns are stored one on top of the other in memory). If *
* A is of type DenseMat, then the (i,j)th element of A (with *
* 0 <= i,j <= size-1) is given by the expression (A->data)[j][i] *
* or by the expression (A->data)[0][j*n+i]. The macros below *
* allow a user to access efficiently individual matrix *
* elements without writing out explicit data structure *
* references and without knowing too much about the underlying *
* element storage. The only storage assumption needed is that *
* elements are stored columnwise and that a pointer to the jth *
* column of elements can be obtained via the DENSE_COL macro. *
* Users should use these macros whenever possible. *
* *
******************************************************************/
typedef struct {
integer size;
real **data;
} *DenseMat;
/* DenseMat accessor macros */
/******************************************************************
* *
* Macro : DENSE_ELEM *
* Usage : DENSE_ELEM(A,i,j) = a_ij; OR *
* a_ij = DENSE_ELEM(A,i,j); *
*----------------------------------------------------------------*
* DENSE_ELEM(A,i,j) references the (i,j)th element of the N by N *
* DenseMat A, 0 <= i,j <= N-1. *
* *
******************************************************************/
#define DENSE_ELEM(A,i,j) ((A->data)[j][i])
/******************************************************************
* *
* Macro : DENSE_COL *
* Usage : col_j = DENSE_COL(A,j); *
*----------------------------------------------------------------*
* DENSE_COL(A,j) references the jth column of the N by N *
* DenseMat A, 0 <= j <= N-1. The type of the expression *
* DENSE_COL(A,j) is real *. After the assignment in the usage *
* above, col_j may be treated as an array indexed from 0 to N-1. *
* The (i,j)th element of A is referenced by col_j[i]. *
* *
******************************************************************/
#define DENSE_COL(A,j) ((A->data)[j])
/* Functions that use the DenseMat representation for a dense matrix */
/******************************************************************
* *
* Function : DenseAllocMat *
* Usage : A = DenseAllocMat(N); *
* if (A == NULL) ... memory request failed *
*----------------------------------------------------------------*
* DenseAllocMat allocates memory for an N by N dense matrix and *
* returns the storage allocated (type DenseMat). DenseAllocMat *
* returns NULL if the request for matrix storage cannot be *
* satisfied. See the above documentation for the type DenseMat *
* for matrix storage details. *
* *
******************************************************************/
DenseMat DenseAllocMat(integer N);
/******************************************************************
* *
* Function : DenseAllocPiv *
* Usage : p = DenseAllocPiv(N); *
* if (p == NULL) ... memory request failed *
*----------------------------------------------------------------*
* DenseAllocPiv allocates memory for pivot information to be *
* filled in by the DenseFactor routine during the factorization *
* of an N by N dense matrix. The underlying type for pivot *
* information is an array of N integers and this routine returns *
* the pointer to the memory it allocates. If the request for *
* pivot storage cannot be satisfied, DenseAllocPiv returns NULL. *
* *
******************************************************************/
integer *DenseAllocPiv(integer N);
/******************************************************************
* *
* Function : DenseFactor *
* Usage : ier = DenseFactor(A, p); *
* if (ier != 0) ... A is singular *
*----------------------------------------------------------------*
* DenseFactor performs the LU factorization of the N by N dense *
* matrix A. This is done using standard Gaussian elimination *
* with partial pivoting. *
* *
* A successful LU factorization leaves the matrix A and the *
* pivot array p with the following information: *
* *
* (1) p[k] contains the row number of the pivot element chosen *
* at the beginning of elimination step k, k=0, 1, ..., N-1. *
* *
* (2) If the unique LU factorization of A is given by PA = LU, *
* where P is a permutation matrix, L is a lower triangular *
* matrix with all 1's on the diagonal, and U is an upper *
* triangular matrix, then the upper triangular part of A *
* (including its diagonal) contains U and the strictly lower *
* triangular part of A contains the multipliers, I-L. *
* *
* DenseFactor returns 0 if successful. Otherwise it encountered *
* a zero diagonal element during the factorization. In this case *
* it returns the column index (numbered from one) at which *
* it encountered the zero. *
* *
******************************************************************/
integer DenseFactor(DenseMat A, integer *p);
/******************************************************************
* *
* Function : DenseBacksolve *
* Usage : DenseBacksolve(A, p, b); *
*----------------------------------------------------------------*
* DenseBacksolve solves the N-dimensional system A x = b using *
* the LU factorization in A and the pivot information in p *
* computed in DenseFactor. The solution x is returned in b. This *
* routine cannot fail if the corresponding call to DenseFactor *
* did not fail. *
* *
******************************************************************/
void DenseBacksolve(DenseMat A, integer *p, N_Vector b);
/******************************************************************
* *
* Function : DenseZero *
* Usage : DenseZero(A); *
*----------------------------------------------------------------*
* DenseZero sets all the elements of the N by N matrix A to 0.0. *
* *
******************************************************************/
void DenseZero(DenseMat A);
/******************************************************************
* *
* Function : DenseCopy *
* Usage : DenseCopy(A, B); *
*----------------------------------------------------------------*
* DenseCopy copies the contents of the N by N matrix A into the *
* N by N matrix B. *
* *
******************************************************************/
void DenseCopy(DenseMat A, DenseMat B);
/******************************************************************
* *
* Function: DenseScale *
* Usage : DenseScale(c, A); *
*----------------------------------------------------------------*
* DenseScale scales the elements of the N by N matrix A by the *
* constant c and stores the result back in A. *
* *
******************************************************************/
void DenseScale(real c, DenseMat A);
/******************************************************************
* *
* Function : DenseAddI *
* Usage : DenseAddI(A); *
*----------------------------------------------------------------*
* DenseAddI adds the identity matrix to A and stores the result *
* back in A. *
* *
******************************************************************/
void DenseAddI(DenseMat A);
/******************************************************************
* *
* Function : DenseFreeMat *
* Usage : DenseFreeMat(A); *
*----------------------------------------------------------------*
* DenseFreeMat frees the memory allocated by DenseAllocMat for *
* the N by N matrix A. *
* *
******************************************************************/
void DenseFreeMat(DenseMat A);
/******************************************************************
* *
* Function : DenseFreePiv *
* Usage : DenseFreePiv(p); *
*----------------------------------------------------------------*
* DenseFreePiv frees the memory allocated by DenseAllocPiv for *
* the pivot information array p. *
* *
******************************************************************/
void DenseFreePiv(integer *p);
/******************************************************************
* *
* Function : DensePrint *
* Usage : DensePrint(A); *
*----------------------------------------------------------------*
* This routine prints the N by N dense matrix A to standard *
* output as it would normally appear on paper. It is intended *
* as a debugging tool with small values of N. The elements are *
* printed using the %g option. A blank line is printed before *
* and after the matrix. *
* *
******************************************************************/
void DensePrint(DenseMat A);
/* Functions that use the real ** representation for a dense matrix */
/******************************************************************
* *
* Function : denalloc *
* Usage : real **a; *
* a = denalloc(n); *
* if (a == NULL) ... memory request failed *
*----------------------------------------------------------------*
* denalloc(n) allocates storage for an n by n dense matrix. It *
* returns a pointer to the newly allocated storage if *
* successful. If the memory request cannot be satisfied, then *
* denalloc returns NULL. The underlying type of the dense matrix *
* returned is real **. If we allocate a dense matrix real **a by *
* a = denalloc(n), then a[j][i] references the (i,j)th element *
* of the matrix a, 0 <= i,j <= n-1, and a[j] is a pointer to the *
* first element in the jth column of a. The location a[0] *
* contains a pointer to n^2 contiguous locations which contain *
* the elements of a. *
* *
******************************************************************/
real **denalloc(integer n);
/******************************************************************
* *
* Function : denallocpiv *
* Usage : integer *pivot; *
* pivot = denallocpiv(n); *
* if (pivot == NULL) ... memory request failed *
*----------------------------------------------------------------*
* denallocpiv(n) allocates an array of n integers. It returns a *
* pointer to the first element in the array if successful. It *
* returns NULL if the memory request could not be satisfied. *
* *
******************************************************************/
integer *denallocpiv(integer n);
/******************************************************************
* *
* Function : gefa *
* Usage : integer ier; *
* ier = gefa(a,n,p); *
* if (ier > 0) ... zero element encountered during *
* the factorization *
*----------------------------------------------------------------*
* gefa(a,n,p) factors the n by n dense matrix a. It overwrites *
* the elements of a with its LU factors and keeps track of the *
* pivot rows chosen in the pivot array p. *
* *
* A successful LU factorization leaves the matrix a and the *
* pivot array p with the following information: *
* *
* (1) p[k] contains the row number of the pivot element chosen *
* at the beginning of elimination step k, k=0, 1, ..., n-1. *
* *
* (2) If the unique LU factorization of a is given by Pa = LU, *
* where P is a permutation matrix, L is a lower triangular *
* matrix with all 1's on the diagonal, and U is an upper *
* triangular matrix, then the upper triangular part of a *
* (including its diagonal) contains U and the strictly lower *
* triangular part of a contains the multipliers, I-L. *
* *
* gefa returns 0 if successful. Otherwise it encountered a zero *
* diagonal element during the factorization. In this case it *
* returns the column index (numbered from one) at which it *
* encountered the zero. *
* *
******************************************************************/
integer gefa(real **a, integer n, integer *p);
/******************************************************************
* *
* Function : gesl *
* Usage : real *b; *
* ier = gefa(a,n,p); *
* if (ier == 0) gesl(a,n,p,b); *
*----------------------------------------------------------------*
* gesl(a,n,p,b) solves the n by n linear system ax = b. It *
* assumes that a has been LU factored and the pivot array p has *
* been set by a successful call to gefa(a,n,p). The solution x *
* is written into the b array. *
* *
******************************************************************/
void gesl(real **a, integer n, integer *p, real *b);
/******************************************************************
* *
* Function : denzero *
* Usage : denzero(a,n); *
*----------------------------------------------------------------*
* denzero(a,n) sets all the elements of the n by n dense matrix *
* a to be 0.0. *
* *
******************************************************************/
void denzero(real **a, integer n);
/******************************************************************
* *
* Function : dencopy *
* Usage : dencopy(a,b,n); *
*----------------------------------------------------------------*
* dencopy(a,b,n) copies the n by n dense matrix a into the *
* n by n dense matrix b. *
* *
******************************************************************/
void dencopy(real **a, real **b, integer n);
/******************************************************************
* *
* Function : denscale *
* Usage : denscale(c,a,n); *
*----------------------------------------------------------------*
* denscale(c,a,n) scales every element in the n by n dense *
* matrix a by c. *
* *
******************************************************************/
void denscale(real c, real **a, integer n);
/******************************************************************
* *
* Function : denaddI *
* Usage : denaddI(a,n); *
*----------------------------------------------------------------*
* denaddI(a,n) increments the n by n dense matrix a by the *
* identity matrix. *
* *
******************************************************************/
void denaddI(real **a, integer n);
/******************************************************************
* *
* Function : denfreepiv *
* Usage : denfreepiv(p); *
*----------------------------------------------------------------*
* denfreepiv(p) frees the pivot array p allocated by *
* denallocpiv. *
* *
******************************************************************/
void denfreepiv(integer *p);
/******************************************************************
* *
* Function : denfree *
* Usage : denfree(a); *
*----------------------------------------------------------------*
* denfree(a) frees the dense matrix a allocated by denalloc. *
* *
******************************************************************/
void denfree(real **a);
/******************************************************************
* *
* Function : denprint *
* Usage : denprint(a,n); *
*----------------------------------------------------------------*
* denprint(a,n) prints the n by n dense matrix a to standard *
* output as it would normally appear on paper. It is intended as *
* a debugging tool with small values of n. The elements are *
* printed using the %g option. A blank line is printed before *
* and after the matrix. *
* *
******************************************************************/
void denprint(real **a, integer n);
#endif
#ifdef __cplusplus
}
#endif

243
ext/cvode/include/iterativ.h
View file

@ -1,243 +0,0 @@
/******************************************************************
* *
* File : iterativ.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 4 May 1998 *
*----------------------------------------------------------------*
* This header file contains declarations intended for use by *
* generic iterative solvers of Ax = b. The enumeration gives *
* symbolic names for the type of preconditioning to be used. *
* The function type declarations give the prototypes for the *
* functions to be called within an iterative linear solver, that *
* are responsible for *
* multiplying A by a given vector v (ATimesFn), and *
* solving the preconditioner equation Pz = r (PSolveFn). *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _iterativ_h
#define _iterativ_h
#include "llnltyps.h"
#include "nvector.h"
/******************************************************************
* *
* enum : types of preconditioning *
*----------------------------------------------------------------*
* NONE : The iterative linear solver should not use *
* preconditioning. *
* *
* LEFT : The iterative linear solver uses preconditioning on *
* the left only. *
* *
* RIGHT : The iterative linear solver uses preconditioning on *
* the right only. *
* *
* BOTH : The iterative linear solver uses preconditioning on *
* both the left and the right. *
* *
******************************************************************/
enum { NONE, LEFT, RIGHT, BOTH };
/******************************************************************
* *
* enum : types of Gram-Schmidt routines *
*----------------------------------------------------------------*
* MODIFIED_GS : The iterative solver uses the modified *
* Gram-Schmidt routine ModifiedGS listed in this *
* file. *
* *
* CLASSICAL_GS : The iterative solver uses the classical *
* Gram-Schmidt routine ClassicalGS listed in this *
* file. *
* *
******************************************************************/
enum { MODIFIED_GS, CLASSICAL_GS };
/******************************************************************
* *
* Type: ATimesFn *
*----------------------------------------------------------------*
* An ATimesFn multiplies Av and stores the result in z. The *
* caller is responsible for allocating memory for the z vector. *
* The parameter A_data is a pointer to any information about A *
* which the function needs in order to do its job. The vector v *
* is unchanged. An ATimesFn returns 0 if successful and a *
* non-zero value if unsuccessful. *
* *
******************************************************************/
typedef int (*ATimesFn)(void *A_data, N_Vector v, N_Vector z);
/******************************************************************
* *
* Type: PSolveFn *
*----------------------------------------------------------------*
* A PSolveFn solves the preconditioner equation Pz = r for the *
* vector z. The caller is responsible for allocating memory for *
* the z vector. The parameter P_data is a pointer to any *
* information about P which the function needs in order to do *
* its job. The parameter lr is input, and indicates whether P *
* is to be taken as the left preconditioner or the right *
* preconditioner: lr = 1 for left and lr = 2 for right. *
* If preconditioning is on one side only, lr can be ignored. *
* The vector r is unchanged. *
* A PSolveFn returns 0 if successful and a non-zero value if *
* unsuccessful. On a failure, a negative return value indicates *
* an unrecoverable condition, while a positive value indicates *
* a recoverable one, in which the calling routine may reattempt *
* the solution after updating preconditioner data. *
* *
******************************************************************/
typedef int (*PSolveFn)(void *P_data, N_Vector r, N_Vector z, int lr);
/******************************************************************
* *
* Function: ModifiedGS *
*----------------------------------------------------------------*
* ModifiedGS performs a modified Gram-Schmidt orthogonalization *
* of the N_Vector v[k] against the p unit N_Vectors at *
* v[k-1], v[k-2], ..., v[k-p]. *
* *
* v is an array of (k+1) N_Vectors v[i], i=0, 1, ..., k. *
* v[k-1], v[k-2], ..., v[k-p] are assumed to have L2-norm *
* equal to 1. *
* *
* h is the output k by k Hessenberg matrix of inner products. *
* This matrix must be allocated row-wise so that the (i,j)th *
* entry is h[i][j]. The inner products (v[i],v[k]), *
* i=i0, i0+1, ..., k-1, are stored at h[i][k-1]. Here *
* i0=MAX(0,k-p). *
* *
* k is the index of the vector in the v array that needs to be *
* orthogonalized against previous vectors in the v array. *
* *
* p is the number of previous vectors in the v array against *
* which v[k] is to be orthogonalized. *
* *
* new_vk_norm is a pointer to memory allocated by the caller to *
* hold the Euclidean norm of the orthogonalized vector v[k]. *
* *
* If (k-p) < 0, then ModifiedGS uses p=k. The orthogonalized *
* v[k] is NOT normalized and is stored over the old v[k]. Once *
* the orthogonalization has been performed, the Euclidean norm *
* of v[k] is stored in (*new_vk_norm). *
* *
* ModifiedGS returns 0 to indicate success. It cannot fail. *
* *
******************************************************************/
int ModifiedGS(N_Vector *v, real **h, int k, int p, real *new_vk_norm);
/******************************************************************
* *
* Function: ClassicalGS *
*----------------------------------------------------------------*
* ClassicalGS performs a classical Gram-Schmidt *
* orthogonalization of the N_Vector v[k] against the p unit *
* N_Vectors at v[k-1], v[k-2], ..., v[k-p]. The parameters v, h, *
* k, p, and new_vk_norm are as described in the documentation *
* for ModifiedGS. *
* *
* temp is an N_Vector which can be used as workspace by the *
* ClassicalGS routine. *
* *
* s is a length k array of reals which can be used as workspace *
* by the ClassicalGS routine. *
* *
* ClassicalGS returns 0 to indicate success. It cannot fail. *
* *
******************************************************************/
int ClassicalGS(N_Vector *v, real **h, int k, int p, real *new_vk_norm,
N_Vector temp, real *s);
/******************************************************************
* *
* Function: QRfact *
*----------------------------------------------------------------*
* QRfact performs a QR factorization of the Hessenberg matrix H. *
* *
* n is the problem size; the matrix H is (n+1) by n. *
* *
* h is the (n+1) by n Hessenberg matrix H to be factored. It is *
* stored row-wise. *
* *
* q is an array of length 2*n containing the Givens rotations *
* computed by this function. A Givens rotation has the form: *
* | c -s | *
* | s c |. *
* The components of the Givens rotations are stored in q as *
* (c, s, c, s, ..., c, s). *
* *
* job is a control flag. If job==0, then a new QR factorization *
* is performed. If job!=0, then it is assumed that the first *
* n-1 columns of h have already been factored and only the last *
* column needs to be updated. *
* *
* QRfact returns 0 if successful. If a zero is encountered on *
* the diagonal of the triangular factor R, then QRfact returns *
* the equation number of the zero entry, where the equations are *
* numbered from 1, not 0. If QRsol is subsequently called in *
* this situation, it will return an error because it could not *
* divide by the zero diagonal entry. *
* *
******************************************************************/
int QRfact(int n, real **h, real *q, int job);
/******************************************************************
* *
* Function: QRsol *
*----------------------------------------------------------------*
* QRsol solves the linear least squares problem *
* *
* min (b - H*x, b - H*x), x in R^n, *
* *
* where H is a Hessenberg matrix, and b is in R^(n+1). *
* It uses the QR factors of H computed by QRfact. *
* *
* n is the problem size; the matrix H is (n+1) by n. *
* *
* h is a matrix (computed by QRfact) containing the upper *
* triangular factor R of the original Hessenberg matrix H. *
* *
* q is an array of length 2*n (computed by QRfact) containing *
* the Givens rotations used to factor H. *
* *
* b is the (n+1)-vector appearing in the least squares problem *
* above. *
* *
* On return, b contains the solution x of the least squares *
* problem, if QRsol was successful. *
* *
* QRsol returns a 0 if successful. Otherwise, a zero was *
* encountered on the diagonal of the triangular factor R. *
* In this case, QRsol returns the equation number (numbered *
* from 1, not 0) of the zero entry. *
* *
******************************************************************/
int QRsol(int n, real **h, real *q, real *b);
#endif
#ifdef __cplusplus
}
#endif

121
ext/cvode/include/llnlmath.h
View file

@ -1,121 +0,0 @@
/******************************************************************
* *
* File : llnlmath.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 4 May 1998 *
*----------------------------------------------------------------*
* This is the header file for a C math library. The routines *
* listed here work with the type real as defined in llnltyps.h. *
* To do single precision floating point arithmetic, set the type *
* real to be float. To do double precision arithmetic, set the *
* type real to be double. The default implementations for *
* RPowerR and RSqrt call standard math library functions which *
* do double precision arithmetic. If this is unacceptable when *
* real is float, then the user should re-implement these two *
* routines by calling single precision routines available on *
* his/her machine. *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _llnlmath_h
#define _llnlmath_h
#include "llnltyps.h"
/******************************************************************
* *
* Macros : MIN, MAX, ABS, SQR *
*----------------------------------------------------------------*
* MIN(A, B) returns the minimum of A and B. *
* *
* MAX(A, B) returns the maximum of A and B. *
* *
* ABS(A) returns the absolute value of A. *
* *
* SQR(A) returns the square of A. *
* *
******************************************************************/
#ifndef MIN
#define MIN(A, B) ((A) < (B) ? (A) : (B))
#endif
#ifndef MAX
#define MAX(A, B) ((A) > (B) ? (A) : (B))
#endif
#ifndef ABS
#define ABS(A) ((A < 0) ? -(A) : (A))
#endif
#ifndef SQR
#define SQR(A) ((A) * (A))
#endif
/******************************************************************
* *
* Function : UnitRoundoff *
* Usage : real uround; *
* uround = UnitRoundoff(); *
*----------------------------------------------------------------*
* UnitRoundoff returns the unit roundoff u for real floating *
* point arithmetic, where u is defined to be the smallest *
* positive real such that 1.0 + u != 1.0. *
* *
******************************************************************/
real UnitRoundoff(void);
/******************************************************************
* *
* Function : RPowerI *
* Usage : int exponent; *
* real base, ans; *
* ans = RPowerI(base,exponent); *
*----------------------------------------------------------------*
* RPowerI returns the value base^exponent, where base is a real *
* and exponent is an int. *
* *
******************************************************************/
real RPowerI(real base, int exponent);
/******************************************************************
* *
* Function : RPowerR *
* Usage : real base, exponent, ans; *
* ans = RPowerR(base,exponent); *
*----------------------------------------------------------------*
* RPowerR returns the value base^exponent, where both base and *
* exponent are reals. If base < 0.0, then RPowerR returns 0.0. *
* *
******************************************************************/
real RPowerR(real base, real exponent);
/******************************************************************
* *
* Function : RSqrt *
* Usage : real sqrt_x; *
* sqrt_x = RSqrt(x); *
*----------------------------------------------------------------*
* RSqrt(x) returns the square root of x. If x < 0.0, then RSqrt *
* returns 0.0. *
* *
******************************************************************/
real RSqrt(real x);
#endif
#ifdef __cplusplus
}
#endif

130
ext/cvode/include/llnltyps.h
View file

@ -1,130 +0,0 @@
/******************************************************************
* *
* File : llnltyps.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 4 May 1998 *
*----------------------------------------------------------------*
* This header file exports three types: real, integer, and boole *
* (short for boolean), as well as the constants TRUE and FALSE. *
* *
* Users should #include "llnltyps.h" in any file that should *
* be easily modifiable to work with different real or integer *
* types and use the exported names real and integer within such *
* a file. The types for real and integer below have been set to *
* double and int, respectively. A user should modify these *
* type declarations as he/she sees fit. For example, if a user *
* wants the work with type float because double precision *
* floating point arithmetic is too expensive on the user's *
* machine, then the definition below should be changed to: *
* *
* typedef float real; *
* *
* Similarly, if a user needs to work with extremely large *
* integers (see the system header file <limits.h> for the limits *
* on type int and long int on your machine), then the user *
* should change the definition below to: *
* *
* typedef long int integer; *
* *
* The constants LLNL_FLOAT, LLNL_DOUBLE, LLNL_INT, LLNL_LONG *
* indicate the underlying types for real and integer. *
* They should be set as follows: *
* *
* (1) #define LLNL_FLOAT 1 *
* #define LLNL_DOUBLE 0 (real is float) *
* *
* (2) #define LLNL_FLOAT 0 *
* #define LLNL_DOUBLE 1 (real is double) *
* *
* (3) #define LLNL_INT 1 *
* #define LLNL_LONG 0 (integer is int) *
* *
* (4) #define LLNL_INT 0 *
* #define LLNL_LONG 1 (integer is long int) *
* *
* Thus the legal types for real are float and double, while *
* the legal types for integer are int and long int. The macro *
* RCONST gives a user a convenient way to define real *
* constants. To use the real constant 1.0, for example, the *
* user should write *
* *
* #define ONE RCONST(1.0) *
* *
* If real is double, then RCONST(1.0) expands to 1.0. If real is *
* float, then RCONST(1.0) expands to 1.0F. There is never a *
* need to explicitly cast 1.0 to (real). *
* *
******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _llnltyps_h
#define _llnltyps_h
/******************************************************************
* *
* Types : real, integer *
*----------------------------------------------------------------*
* The types real and integer are currently set to double and *
* int, respectively. See the documentation at the top for *
* usage details and a description of associated constants and *
* macros. *
* *
******************************************************************/
typedef double real;
//typedef long int integer;
typedef int integer;
#define LLNL_FLOAT 0
#define LLNL_DOUBLE 1
#define LLNL_INT 1
#define LLNL_LONG 0
#if LLNL_FLOAT
#define RCONST(x) x##F
#elif LLNL_DOUBLE
#define RCONST(x) x
#endif
/******************************************************************
* *
* Type : boole *
* Constants : FALSE, TRUE *
*----------------------------------------------------------------*
* ANSI C does not have a built-in boolean type. Below is the *
* definition for a new type boole. The advantage of using the *
* name boole (instead of int) is an increase in code readability.*
* It allows the programmer to make a distinction between int and *
* boolean data. Variables of type boole are intended to have only*
* the two values FALSE and TRUE which are defined below to be *
* equal to 0 and 1, respectively. *
* *
******************************************************************/
#ifndef boole
#define boole int
#endif
#ifndef FALSE
#define FALSE 0
#endif
#ifndef TRUE
#define TRUE 1
#endif
#endif
#ifdef __cplusplus
}
#endif

535
ext/cvode/include/nvector.h
View file

@ -1,535 +0,0 @@
/****************************************************************
* *
* File : nvector.h *
* Programmers : Scott D. Cohen, Alan C. Hindmarsh, and *
* : Allan G. Taylor, LLNL *
* Version of : 17 December 1999 *
*--------------------------------------------------------------*
* *
* This is the header file for a generic serial NVECTOR package.*
* It exports the type N_Vector. *
* *
* Part I of this file contains declarations which are specific *
* to the particular machine environment in which this version *
* of the vector package is to be used. This includes the *
* typedef for the type N_Vector, as well as accessor macros *
* that allow the user to use efficiently the type N_Vector *
* without making explicit references to its underlying *
* representation. The underlying type of N_Vector will always *
* be some pointer type. *
* *
* Part II of this file contains the prototypes for the vector *
* kernels which operate on the type N_Vector. These prototypes *
* are fixed for all implementations of the vector package. The *
* definitions of the types real and integer are in the header *
* file llnltyps.h and these may be changed according to the *
* user's needs. The llnltyps.h file also contains the *
* definition for the type boole (short for boolean) that is the*
* return type for the routine N_VInvTest. *
* *
* Important Note: N_Vector arguments to arithmetic kernels *
* need not be distinct. Thus, for example, the call *
* N_VLinearSum(a,x,b,y,y); y <- ax+by *
* is legal. *
* *
* This version of nvector.h is for the ordinary sequential *
* machine environment. In the documentation given below, N is *
* the length of all N_Vector parameters and x[i] denotes the *
* ith component of the N_Vector x, where 0 <= i <= N-1. *
* *
****************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef nvector_h
#define nvector_h
#include "llnltyps.h"
/* Part I: Machine Environment-Dependent Declarations */
/* Environment: Sequential */
typedef struct {
int dummy; /* dummy element */
} *machEnvType; /* dummy machEnvType definition */
/***************************************************************
* *
* Type: N_Vector *
*-------------------------------------------------------------*
* The type N_Vector is an abstract vector type. The fields of *
* its concrete representation should not be accessed *
* directly, but rather through the macros given below. *
* *
* A user may assume that the N components of an N_Vector *
* are stored contiguously. A pointer to the first component *
* can be obtained via the macro N_VDATA. *
* *
***************************************************************/
typedef struct {
integer length;
real *data;
} *N_Vector;
/***************************************************************
* *
* Macros: N_VMAKE, N_VDISPOSE, N_VDATA, N_VLENGTH, N_VIth *
*-------------------------------------------------------------*
* In the descriptions below, the following user *
* declarations are assumed: *
* *
* N_Vector v; real *v_data, r; integer v_len, i; *
* *
* (1) N_VMAKE, N_VDISPOSE *
* *
* These companion routines are used to create and *
* destroy an N_Vector with a component array v_data *
* allocated by the user. *
* *
* The call N_VMAKE(v, v_data, v_len) makes v an *
* N_Vector with component array v_data and length v_len. *
* N_VMAKE stores the pointer v_data so that changes *
* made by the user to the elements of v_data are *
* simultaneously reflected in v. There is no copying of *
* elements. *
* *
* The call N_VDISPOSE(v) frees all memory associated *
* with v except for its component array. This memory was *
* allocated by the user and, therefore, should be *
* deallocated by the user. *
* *
* (2) N_VDATA, N_VLENGTH *
* *
* These routines give individual access to the parts of *
* an N_Vector. *
* *
* The assignment v_data=N_VDATA(v) sets v_data to be *
* a pointer to the first component of v. The assignment *
* N_VDATA(v)=v_data sets the component array of v to *
* be v_data by storing the pointer v_data. *
* *
* The assignment v_len=N_VLENGTH(v) sets v_len to be *
* the length of v. The call N_VLENGTH(v)=len_v sets *
* the length of v to be len_v. *
* *
* (3) N_VIth *
* *
* In the following description, the components of an *
* N_Vector are numbered 0..N-1, where N is the length of *
* v. *
* *
* The assignment r=N_VIth(v,i) sets r to be the value of *
* the ith component of v. The assignment N_VIth(v,i)=r *
* sets the value of the ith component of v to be r. *
* *
* Notes.. *
* *
* Users who use the macros (1) must #include<stdlib.h> *
* since these macros expand to calls to malloc and free. *
* *
* When looping over the components of an N_Vector v, it is *
* more efficient to first obtain the component array via *
* v_data=N_VDATA(v) and then access v_data[i] within the *
* loop than it is to use N_VIth(v,i) within the loop. *
* *
* N_VMAKE and N_VDISPOSE are similar to N_VNew and N_VFree. *
* The difference is one of responsibility for component *
* memory allocation and deallocation. N_VNew allocates memory *
* for the N_Vector components and N_VFree frees the component *
* memory allocated by N_VNew. For N_VMAKE and N_VDISPOSE, the *
* component memory is allocated and freed by the user of *
* this package. *
* *
***************************************************************/
#define N_VMAKE(v, v_data, v_len) v = (N_Vector) malloc(sizeof(*v)); \
v->data = v_data; \
v->length = v_len
#define N_VDISPOSE(v) free(v)
#define N_VDATA(v) (v->data)
#define N_VLENGTH(v) (v->length)
#define N_VIth(v,i) ((v->data)[i])
/* Part II: N_Vector Kernel Prototypes (Machine Environment-Independent) */
/***************************************************************
* *
* Memory Allocation and Deallocation: N_VNew, N_VFree *
* *
***************************************************************/
/***************************************************************
* *
* Function : N_VNew *
* Usage : x = N_VNew(N, machEnv); *
*-------------------------------------------------------------*
* *
* Returns a new N_Vector of length N. The parameter machEnv *
* is a pointer to machine environment-specific information. *
* It is ignored in the sequential machine environment and the *
* user in this environment should simply pass NULL for this *
* argument. If there is not enough memory for a new N_Vector, *
* then N_VNew returns NULL. *
* *
***************************************************************/
N_Vector N_VNew(integer n, void *machEnv);
/***************************************************************
* *
* Function : N_VFree *
* Usage : N_VFree(x); *
*-------------------------------------------------------------*
* *
* Frees the N_Vector x. It is illegal to use x after the call *
* N_VFree(x). *
* *
***************************************************************/
void N_VFree(N_Vector x);
/***************************************************************
* *
* N_Vector Arithmetic: N_VLinearSum, N_VConst, N_VProd, *
* N_VDiv, N_VScale, N_VAbs, N_VInv, *
* N_VAddConst *
* *
***************************************************************/
/***************************************************************
* *
* Function : N_VLinearSum *
* Operation : z = a x + b y *
* *
***************************************************************/
void N_VLinearSum(real a, N_Vector x, real b, N_Vector y, N_Vector z);
/***************************************************************
* *
* Function : N_VConst *
* Operation : z[i] = c for i=0, 1, ..., N-1 *
* *
***************************************************************/
void N_VConst(real c, N_Vector z);
/***************************************************************
* *
* Function : N_VProd *
* Operation : z[i] = x[i] * y[i] for i=0, 1, ..., N-1 *
* *
***************************************************************/
void N_VProd(N_Vector x, N_Vector y, N_Vector z);
/***************************************************************
* *
* Function : N_VDiv *
* Operation : z[i] = x[i] / y[i] for i=0, 1, ..., N-1 *
* *
***************************************************************/
void N_VDiv(N_Vector x, N_Vector y, N_Vector z);
/***************************************************************
* *
* Function : N_VScale *
* Operation : z = c x *
* *
***************************************************************/
void N_VScale(real c, N_Vector x, N_Vector z);
/***************************************************************
* *
* Function : N_VAbs *
* Operation : z[i] = |x[i]|, for i=0, 1, ..., N-1 *
* *
***************************************************************/
void N_VAbs(N_Vector x, N_Vector z);
/***************************************************************
* *
* Function : N_VInv *
* Operation : z[i] = 1.0 / x[i] for i = 0, 1, ..., N-1 *
*-------------------------------------------------------------*
* *
* This routine does not check for division by 0. It should be *
* called only with an N_Vector x which is guaranteed to have *
* all non-zero components. *
* *
***************************************************************/
void N_VInv(N_Vector x, N_Vector z);
/***************************************************************
* *
* Function : N_VAddConst *
* Operation : z[i] = x[i] + b for i = 0, 1, ..., N-1 *
* *
***************************************************************/
void N_VAddConst(N_Vector x, real b, N_Vector z);
/***************************************************************
* *
* N_Vector Measures: N_VDotProd, N_VMaxNorm, VWrmsNorm, *
* N_VMin, N_VWL2Norm, N_VL1Norm *
* *
* *
***************************************************************/
/***************************************************************
* *
* Function : N_VDotProd *
* Usage : dotprod = N_VDotProd(x, y); *
*-------------------------------------------------------------*
* *
* Returns the value of the ordinary dot product of x and y: *
* *
* -> sum (i=0 to N-1) {x[i] * y[i]} *
* *
* Returns 0.0 if N <= 0. *
* *
***************************************************************/
real N_VDotProd(N_Vector x, N_Vector y);
/***************************************************************
* *
* Function : N_VMaxNorm *
* Usage : maxnorm = N_VMaxNorm(x); *
*-------------------------------------------------------------*
* *
* Returns the maximum norm of x: *
* *
* -> max (i=0 to N-1) |x[i]| *
* *
* Returns 0.0 if N <= 0. *
* *
***************************************************************/
real N_VMaxNorm(N_Vector x);
/***************************************************************
* *
* Function : N_VWrmsNorm *
* Usage : wrmsnorm = N_VWrmsNorm(x, w); *
*-------------------------------------------------------------*
* *
* Returns the weighted root mean square norm of x with *
* weight vector w: *
* *
* -> sqrt [(sum (i=0 to N-1) {(x[i] * w[i])^2}) / N] *
* *
* Returns 0.0 if N <= 0. *
* *
***************************************************************/
real N_VWrmsNorm(N_Vector x, N_Vector w);
/***************************************************************
* *
* Function : N_VMin *
* Usage : min = N_VMin(x); *
*-------------------------------------------------------------*
* *
* Returns min x[i] if N > 0 and returns 0.0 if N <= 0. *
* i *
* *
***************************************************************/
real N_VMin(N_Vector x);
/***************************************************************
* *
* Function : N_VWL2Norm *
* Usage : wl2norm = N_VWL2Norm(x, w); *
*-------------------------------------------------------------*
* *
* Returns the weighted Euclidean L2 norm of x with *
* weight vector w: *
* *
* -> sqrt [(sum (i=0 to N-1) {(x[i] * w[i])^2}) ] *
* *
* Returns 0.0 if N <= 0. *
* *
***************************************************************/
real N_VWL2Norm(N_Vector x, N_Vector w);
/***************************************************************
* *
* Function : N_VL1Norm *
* Usage : l1norm = N_VL1Norm(x); *
*-------------------------------------------------------------*
* *
* Returns sum of ABS(x[i]) if N > 0 and returns 0.0 if N <= 0.*
* i *
*
* i.e., calculates and returns the L1 norm of x *
* *
***************************************************************/
real N_VL1Norm(N_Vector x);
/***************************************************************
* *
* Miscellaneous : N_VOneMask, N_VCompare, N_VInvTest, *
* N_VConstrProdPos, N_VConstrMask, and N_VMinQuotient *
* *
***************************************************************/
/***************************************************************
* *
* Function : N_VOneMask *
* Operation : x[i] = 1.0 if |x[i]| != 0. i = 0, 1, ..., N-1 *
* 0.0 otherwise *
* *
***************************************************************/
void N_VOneMask(N_Vector x);
/***************************************************************
* *
* Function : N_VCompare *
* Operation : z[i] = 1.0 if |x[i]| >= c i = 0, 1, ..., N-1 *
* 0.0 otherwise *
* *
***************************************************************/
void N_VCompare(real c, N_Vector x, N_Vector z);
/***************************************************************
* *
* Function : N_VInvTest *
* Operation : z[i] = 1.0 / x[i] with a test for x[i]==0.0 *
* before inverting x[i]. *
*-------------------------------------------------------------*
* *
* This routine returns TRUE if all components of x are *
* non-zero (successful inversion) and returns FALSE *
* otherwise. *
* *
***************************************************************/
boole N_VInvTest(N_Vector x, N_Vector z);
/***************************************************************
* *
* Function : N_VConstrProdPos *
* Usage : booltest = N_VConstrProdPos(c,x); *
*-------------------------------------------------------------*
* *
* Returns a boolean FALSE if some c[i]!=0.0 and x[i]*c[i]<=0.0*
* and TRUE otherwise *
* *
* This routine is used for constraint checking. *
* *
***************************************************************/
boole N_VConstrProdPos(N_Vector c, N_Vector x);
/***************************************************************
* *
* Function : N_VConstrMask *
* Operation : m[i] = 1.0 , if constraint test fails, for i *
* m[i] = 0.0 , if constraint test passes, for i *
* where the constraint tests parallel those *
* of routine N_VConstrProdPos *
*-------------------------------------------------------------*
* This routine returns a boole FALSE if any element failed *
* the constraint test, TRUE if all passed. It also creates a *
* mask vector, m, which has all elements whose corresponding *
* constraint test failed, marked with 1.0, passed with 0.0 *
* This routine is specialized in that it is used only for *
* constraint checking. *
***************************************************************/
boole N_VConstrMask(N_Vector c, N_Vector x, N_Vector m);
/***************************************************************
* *
* Function : N_VMinQuotient *
* Operation : minq = min ( num[i]/denom[i]) over all i such *
* that denom[i] != 0. *
*-------------------------------------------------------------*
* *
* This routine returns the minimum of the quotients obtained *
* by term-wise dividing num[i] by denom[i]. A zero element *
* in denom will be skipped. If no such quotients are found, *
* then the large value 1.e99 is returned. *
* *
***************************************************************/
real N_VMinQuotient(N_Vector num, N_Vector denom);
/***************************************************************
* *
* Debugging Tools : N_VPrint *
* *
***************************************************************/
/***************************************************************
* *
* Function : N_VPrint *
* Usage : N_VPrint(x); *
*-------------------------------------------------------------*
* *
* Prints the N_Vector x to stdout. Each component of x is *
* printed on a separate line using the %g specification. This *
* routine is provided as an aid in debugging code which uses *
* this vector package. *
* *
***************************************************************/
void N_VPrint(N_Vector x);
#endif
#ifdef __cplusplus
}
#endif

295
ext/cvode/include/spgmr.h
View file

@ -1,295 +0,0 @@
/*****************************************************************************
* File : spgmr.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 17 December 1999 *
*---------------------------------------------------------------------------*
* This is the header file for the implementation of SPGMR Krylov *
* iterative linear solver. The SPGMR algorithm is based on the Scaled *
* Preconditioned GMRES (Generalized Minimal Residual) method. *
* *
* The SPGMR algorithm solves a N by N linear system A x = b. *
* Preconditioning is allowed on the left, right, or both. *
* Scaling is allowed on both sides, and restarts are also allowed. *
* We denote the preconditioner and scaling matrices as follows: *
* P1 = left preconditioner *
* P2 = right preconditioner *
* S1 = diagonal matrix of scale factors for P1-inverse b *
* S2 = diagonal matrix of scale factors for P2 x *
* The matrices A, P1, and P2 are not required explicitly; only routines *
* that provide A, P1-inverse, and P2-inverse as operators are required. *
* *
* In this notation, SPGMR applies the underlying GMRES method to the *
* equivalent transformed system *
* Abar xbar = bbar , where *
* Abar = S1 (P1-inverse) A (P2-inverse) (S2-inverse) , *
* bbar = S1 (P1-inverse) b , and xbar = S2 P2 x . *
* *
* The scaling matrices must be chosen so that vectors S1 P1-inverse b *
* and S2 P2 x have dimensionless components. If preconditioning is done *
* on the left only (P2 = I), by a matrix P, then S2 must be a scaling *
* for x, while S1 is a scaling for P-inverse b, and so may also be taken *
* as a scaling for x. Similarly, if preconditioning is done on the *
* right only (P1 = I, P2 = P), then S1 must be a scaling for b, while S2 *
* is a scaling for P x, and may also be taken as a scaling for b. *
* *
* The stopping test for the SPGMR iterations is on the L2 norm of the *
* scaled preconditioned residual: *
* || bbar - Abar xbar ||_2 < delta *
* with an input test constant delta. *
* *
* The usage of this SPGMR solver involves supplying two routines and *
* making three calls. The user-supplied routines are *
* atimes (A_data, x, y) to compute the product y = A x, given x, *
* and *
* psolve (P_data, x, y, lr) to solve P1 x = y or P2 x = y for x, given y.*
* The three user calls are: *
* mem = SpgmrMalloc(N, lmax, machEnv); *
* to initialize memory, *
* flag = SpgmrSolve(mem,A_data,x,b,...,P_data,s1,s2,atimes,psolve,...); *
* to solve the system, and *
* SpgmrFree(mem); *
* to free the memory created by SpgmrMalloc. *
* Complete details for specifying atimes and psolve and for the usage calls *
* are given in the paragraphs below and in iterativ.h. *
* *
*****************************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _spgmr_h
#define _spgmr_h
#include "llnltyps.h"
#include "iterativ.h"
#include "nvector.h"
/******************************************************************
* *
* Types: SpgmrMemRec, SpgmrMem *
*----------------------------------------------------------------*
* SpgmrMem is a pointer to an SpgmrMemRec which contains *
* the memory needed by SpgmrSolve. The SpgmrMalloc routine *
* returns a pointer of type SpgmrMem which should then be passed *
* in subsequent calls to SpgmrSolve. The SpgmrFree routine frees *
* the memory allocated by SpgmrMalloc. *
* *
* N is the linear system size. *
* *
* l_max is the maximum Krylov dimension that SpgmrSolve will be *
* permitted to use. *
* *
* V is the array of Krylov basis vectors v_1, ..., v_(l_max+1), *
* stored in V[0], ..., V[l_max], where l_max is the second *
* parameter to SpgmrMalloc. Each v_i is a length N vector of *
* type N_Vector. (N is the first parameter to SpgmrMalloc and *
* represents the size of the linear system.) *
* *
* Hes is the (l_max+1) x l_max Hessenberg matrix. It is stored *
* row-wise so that the (i,j)th element is given by Hes[i][j]. *
* *
* givens is a length 2*l_max array which represents the *
* Givens rotation matrices that arise in the algorithm. The *
* Givens rotation matrices F_0, F_1, ..., F_j, where F_i is *
* *
* 1 *
* 1 *
* c_i -s_i <--- row i *
* s_i c_i *
* 1 *
* 1 *
* *
* are represented in the givens vector as *
* givens[0]=c_0, givens[1]=s_0, givens[2]=c_1, givens[3]=s_1, *
* ..., givens[2j]=c_j, givens[2j+1]=s_j. *
* *
* xcor is a length N vector (type N_Vector) which holds the *
* scaled, preconditioned correction to the initial guess. *
* *
* yg is a length (l_max+1) array of reals used to hold "short" *
* vectors (e.g. y and g). *
* *
* vtemp is a length N vector (type N_Vector) used as temporary *
* vector storage during calculations. *
* *
******************************************************************/
typedef struct {
integer N;
int l_max;
N_Vector *V;
real **Hes;
real *givens;
N_Vector xcor;
real *yg;
N_Vector vtemp;
} SpgmrMemRec, *SpgmrMem;
/******************************************************************
* *
* Function : SpgmrMalloc *
*----------------------------------------------------------------*
* SpgmrMalloc allocates the memory used by SpgmrSolve. It *
* returns a pointer of type SpgmrMem which the user of the *
* SPGMR package should pass to SpgmrSolve. The parameter N *
* is the size of the system to be solved by SpgmrSolve and l_max *
* is the maximum Krylov dimension that SpgmrSolve will be *
* permitted to use. The parameter machEnv is a pointer to *
* machine environment-specific information. Pass NULL in the *
* ordinary sequential case (see nvector.h). This routine returns *
* NULL if there is a memory request failure. *
* *
******************************************************************/
SpgmrMem SpgmrMalloc(integer N, int l_max, void *machEnv);
/******************************************************************
* *
* Function : SpgmrSolve *
*----------------------------------------------------------------*
* SpgmrSolve solves the linear system Ax = b using the SPGMR *
* method. The return values are given by the symbolic constants *
* below. The first SpgmrSolve parameter is a pointer to memory *
* allocated by a prior call to SpgmrMalloc. The system size N *
* passed in the call to SpgmrMalloc should be the same as the *
* length of all N_Vector arguments passed to SpgmrSolve. *
* *
* mem is the pointer returned by SpgmrMalloc to the structure *
* containing the memory needed by SpgmrSolve. *
* *
* A_data is a pointer to information about the coefficient *
* matrix A. This pointer is passed to the user-supplied function *
* atimes. *
* *
* x is the initial guess x_0 upon entry and the solution *
* N_Vector upon exit with return value SPGMR_SUCCESS or *
* SPGMR_RES_REDUCED. For all other return values, the output x *
* is undefined. *
* *
* b is the right hand side N_Vector. It is undisturbed by this *
* function. *
* *
* pretype is the type of preconditioning to be used. Its *
* legal possible values are enumerated in iterativ.h. These *
* values are NONE=0, LEFT=1, RIGHT=2, and BOTH=3. *
* *
* gstype is the type of Gram-Schmidt orthogonalization to be *
* used. Its legal values are enumerated in iterativ.h. These *
* values are MODIFIED_GS=0 and CLASSICAL_GS=1. *
* *
* delta is the tolerance on the L2 norm of the scaled, *
* preconditioned residual. On return with value SPGMR_SUCCESS, *
* this residual satisfies || s1 P1_inv (b - Ax) ||_2 <= delta. *
* *
* max_restarts is the maximum number of times the algorithm is *
* allowed to restart. *
* *
* P_data is a pointer to preconditioner information. This *
* pointer is passed to the user-supplied function psolve. *
* *
* s1 is an N_Vector of positive scale factors for P1-inv b, where*
* P1 is the left preconditioner. (Not tested for positivity.) *
* Pass NULL if no scaling on P1-inv b is required. *
* *
* s2 is an N_Vector of positive scale factors for P2 x, where *
* P2 is the right preconditioner. (Not tested for positivity.) *
* Pass NULL if no scaling on P2 x is required. *
* *
* atimes is the user-supplied function which performs the *
* operation of multiplying A by a given vector. Its description *
* is given in iterativ.h. *
* *
* psolve is the user-supplied function which solves a *
* preconditioner system Pz = r, where P is P1 or P2. Its full *
* description is given in iterativ.h. The psolve function will *
* not be called if pretype is NONE; in that case, the user *
* should pass NULL for psolve. *
* *
* res_norm is a pointer to the L2 norm of the scaled, *
* preconditioned residual. On return with value SPGMR_SUCCESS or *
* SPGMR_RES_REDUCED, (*res_norm) contains the value *
* || s1 P1_inv (b - Ax) ||_2 for the computed solution x. *
* For all other return values, (*res_norm) is undefined. The *
* caller is responsible for allocating the memory (*res_norm) *
* to be filled in by SpgmrSolve. *
* *
* nli is a pointer to the number of linear iterations done in *
* the execution of SpgmrSolve. The caller is responsible for *
* allocating the memory (*nli) to be filled in by SpgmrSolve. *
* *
* nps is a pointer to the number of calls made to psolve during *
* the execution of SpgmrSolve. The caller is responsible for *
* allocating the memory (*nps) to be filled in by SpgmrSolve. *
* *
* Note.. Repeated calls can be made to SpgmrSolve with varying *
* input arguments. If, however, the problem size N or the *
* maximum Krylov dimension l_max changes, then a call to *
* SpgmrMalloc must be made to obtain new memory for SpgmrSolve *
* to use. *
* *
******************************************************************/
int SpgmrSolve(SpgmrMem mem, void *A_data, N_Vector x, N_Vector b,
int pretype, int gstype, real delta, int max_restarts,
void *P_data, N_Vector s1, N_Vector s2, ATimesFn atimes,
PSolveFn psolve, real *res_norm, int *nli, int *nps);
/* Return values for SpgmrSolve */
#define SPGMR_SUCCESS 0 /* Converged */
#define SPGMR_RES_REDUCED 1 /* Did not converge, but reduced
norm of residual */
#define SPGMR_CONV_FAIL 2 /* Failed to converge */
#define SPGMR_QRFACT_FAIL 3 /* QRfact found singular matrix */
#define SPGMR_PSOLVE_FAIL_REC 4 /* psolve failed recoverably */
#define SPGMR_MEM_NULL -1 /* mem argument is NULL */
#define SPGMR_ATIMES_FAIL -2 /* atimes returned failure flag */
#define SPGMR_PSOLVE_FAIL_UNREC -3 /* psolve failed unrecoverably */
#define SPGMR_GS_FAIL -4 /* Gram-Schmidt routine
returned failure flag */
#define SPGMR_QRSOL_FAIL -5 /* QRsol found singular R */
/******************************************************************
* *
* Function : SpgmrFree *
*----------------------------------------------------------------*
* SpgmrMalloc frees the memory allocated by SpgmrMalloc. It is *
* illegal to use the pointer mem after a call to SpgmrFree. *
* *
******************************************************************/
void SpgmrFree(SpgmrMem mem);
/******************************************************************
* Macro: SPGMR_VTEMP *
* *
*----------------------------------------------------------------*
* This macro provides access to the work vector vtemp in the *
* memory block of the SPGMR module. The argument mem is the *
* memory pointer returned by SpgmrMalloc, of type SpgmrMem, *
* and the macro value is of type N_Vector. *
* On a return from SpgmrSolve with *nli = 0, this vector *
* contains the scaled preconditioned initial residual, *
* s1 * P1_inverse * (b - A x_0). *
******************************************************************/
#define SPGMR_VTEMP(mem) (mem->vtemp)
#endif
#ifdef __cplusplus
}
#endif

356
ext/cvode/source/band.c
View file

@ -1,356 +0,0 @@
/******************************************************************
* *
* File : band.c *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 25 February 2000 *
*----------------------------------------------------------------*
* This is the implementation file for a generic BAND linear *
* solver package. *
* *
******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "band.h"
#include "llnltyps.h"
#include "nvector.h"
#include "llnlmath.h"
#define ZERO RCONST(0.0)
#define ONE RCONST(1.0)
#define ROW(i,j,smu) (i-j+smu)
/* Implementation */
BandMat BandAllocMat(integer N, integer mu, integer ml, integer smu)
{
BandMat A;
if (N <= 0) return(NULL);
A = (BandMat) malloc(sizeof *A);
if (A == NULL) return (NULL);
A->data = bandalloc(N, smu, ml);
if (A->data == NULL) {
free(A);
return(NULL);
}
A->size = N;
A->mu = mu;
A->ml = ml;
A->smu = smu;
return(A);
}
integer *BandAllocPiv(integer N)
{
if (N <= 0) return(NULL);
return((integer *) malloc(N * sizeof(integer)));
}
integer BandFactor(BandMat A, integer *p)
{
return(gbfa(A->data, A->size, A->mu, A->ml, A->smu, p));
}
void BandBacksolve(BandMat A, integer *p, N_Vector b)
{
gbsl(A->data, A->size, A->smu, A->ml, p, N_VDATA(b));
}
void BandZero(BandMat A)
{
bandzero(A->data, A->size, A->mu, A->ml, A->smu);
}
void BandCopy(BandMat A, BandMat B, integer copymu, integer copyml)
{
bandcopy(A->data, B->data, A->size, A->smu, B->smu, copymu, copyml);
}
void BandScale(real c, BandMat A)
{
bandscale(c, A->data, A->size, A->mu, A->ml, A->smu);
}
void BandAddI(BandMat A)
{
bandaddI(A->data, A->size, A->smu);
}
void BandFreeMat(BandMat A)
{
bandfree(A->data);
free(A);
}
void BandFreePiv(integer *p)
{
free(p);
}
void BandPrint(BandMat A)
{
bandprint(A->data, A->size, A->mu, A->ml, A->smu);
}
real **bandalloc(integer n, integer smu, integer ml)
{
real **a;
integer j, colSize;
if (n <= 0) return(NULL);
a = (real **) malloc(n * sizeof(real *));
if (a == NULL) return(NULL);
colSize = smu + ml + 1;
a[0] = (real *) malloc(n * colSize * sizeof(real));
if (a[0] == NULL) {
free(a);
return(NULL);
}
for (j=1; j < n; j++) a[j] = a[0] + j * colSize;
return(a);
}
integer *bandallocpiv(integer n)
{
if (n <= 0) return(NULL);
return((integer *) malloc(n * sizeof(integer)));
}
integer gbfa(real **a, integer n, integer mu, integer ml, integer smu,
integer *p)
{
integer c, r, num_rows;
integer i, j, k, l, storage_l, storage_k, last_col_k, last_row_k;
real *a_c, *col_k, *diag_k, *sub_diag_k, *col_j, *kptr, *jptr;
real max, temp, mult, a_kj;
boole swap;
/* zero out the first smu - mu rows of the rectangular array a */
num_rows = smu - mu;
if (num_rows > 0) {
for (c=0; c < n; c++) {
a_c = a[c];
for (r=0; r < num_rows; r++) {
a_c[r] = ZERO;
}
}
}
/* k = elimination step number */
for (k=0; k < n-1; k++, p++) {
col_k = a[k];
diag_k = col_k + smu;
sub_diag_k = diag_k + 1;
last_row_k = MIN(n-1,k+ml);
/* find l = pivot row number */
l=k;
max = ABS(*diag_k);
for (i=k+1, kptr=sub_diag_k; i <= last_row_k; i++, kptr++) {
if (ABS(*kptr) > max) {
l=i;
max = ABS(*kptr);
}
}
storage_l = ROW(l, k, smu);
*p = l;
/* check for zero pivot element */
if (col_k[storage_l] == ZERO) return(k+1);
/* swap a(l,k) and a(k,k) if necessary */
if ( (swap = (l != k) )) {
temp = col_k[storage_l];
col_k[storage_l] = *diag_k;
*diag_k = temp;
}
/* Scale the elements below the diagonal in */
/* column k by -1.0 / a(k,k). After the above swap, */
/* a(k,k) holds the pivot element. This scaling */
/* stores the pivot row multipliers -a(i,k)/a(k,k) */
/* in a(i,k), i=k+1, ..., MIN(n-1,k+ml). */
mult = -ONE / (*diag_k);
for (i=k+1, kptr = sub_diag_k; i <= last_row_k; i++, kptr++)
(*kptr) *= mult;
/* row_i = row_i - [a(i,k)/a(k,k)] row_k, i=k+1, ..., MIN(n-1,k+ml) */
/* row k is the pivot row after swapping with row l. */
/* The computation is done one column at a time, */
/* column j=k+1, ..., MIN(k+smu,n-1). */
last_col_k = MIN(k+smu,n-1);
for (j=k+1; j <= last_col_k; j++) {
col_j = a[j];
storage_l = ROW(l,j,smu);
storage_k = ROW(k,j,smu);
a_kj = col_j[storage_l];
/* Swap the elements a(k,j) and a(k,l) if l!=k. */
if (swap) {
col_j[storage_l] = col_j[storage_k];
col_j[storage_k] = a_kj;
}
/* a(i,j) = a(i,j) - [a(i,k)/a(k,k)]*a(k,j) */
/* a_kj = a(k,j), *kptr = - a(i,k)/a(k,k), *jptr = a(i,j) */
if (a_kj != ZERO) {
for (i=k+1, kptr=sub_diag_k, jptr=col_j+ROW(k+1,j,smu);
i <= last_row_k;
i++, kptr++, jptr++)
(*jptr) += a_kj * (*kptr);
}
}
}
/* set the last pivot row to be n-1 and check for a zero pivot */
*p = n-1;
if (a[n-1][smu] == ZERO) return(n);
/* return 0 to indicate success */
return(0);
}
void gbsl(real **a, integer n, integer smu, integer ml, integer *p,
real *b)
{
integer k, l, i, first_row_k, last_row_k;
real mult, *diag_k;
/* Solve Ly = Pb, store solution y in b */
for (k=0; k < n-1; k++) {
l = p[k];
mult = b[l];
if (l != k) {
b[l] = b[k];
b[k] = mult;
}
diag_k = a[k]+smu;
last_row_k = MIN(n-1,k+ml);
for (i=k+1; i <= last_row_k; i++)
b[i] += mult * diag_k[i-k];
}
/* Solve Ux = y, store solution x in b */
for (k=n-1; k >= 0; k--) {
diag_k = a[k]+smu;
first_row_k = MAX(0,k-smu);
b[k] /= (*diag_k);
mult = -b[k];
for (i=first_row_k; i <= k-1; i++)
b[i] += mult*diag_k[i-k];
}
}
void bandzero(real **a, integer n, integer mu, integer ml, integer smu)
{
integer i, j, colSize;
real *col_j;
colSize = mu + ml + 1;
for (j=0; j < n; j++) {
col_j = a[j]+smu-mu;
for (i=0; i < colSize; i++)
col_j[i] = ZERO;
}
}
void bandcopy(real **a, real **b, integer n, integer a_smu, integer b_smu,
integer copymu, integer copyml)
{
integer i, j, copySize;
real *a_col_j, *b_col_j;
copySize = copymu + copyml + 1;
for (j=0; j < n; j++) {
a_col_j = a[j]+a_smu-copymu;
b_col_j = b[j]+b_smu-copymu;
for (i=0; i < copySize; i++)
b_col_j[i] = a_col_j[i];
}
}
void bandscale(real c, real **a, integer n, integer mu, integer ml,
integer smu)
{
integer i, j, colSize;
real *col_j;
colSize = mu + ml + 1;
for(j=0; j < n; j++) {
col_j = a[j]+smu-mu;
for (i=0; i < colSize; i++)
col_j[i] *= c;
}
}
void bandaddI(real **a, integer n, integer smu)
{
integer j;
for(j=0; j < n; j++)
a[j][smu] += ONE;
}
void bandfreepiv(integer *p)
{
free(p);
}
void bandfree(real **a)
{
free(a[0]);
free(a);
}
void bandprint(real **a, integer n, integer mu, integer ml, integer smu)
{
integer i, j, start, finish;
printf("\n");
for (i=0; i < n; i++) {
start = MAX(0,i-ml);
finish = MIN(n-1,i+mu);
for (j=0; j < start; j++) printf("%10s","");
for (j=start; j <= finish; j++) {
printf("%10g", a[j][i-j+smu]);
}
printf("\n");
}
printf("\n");
}

415
ext/cvode/source/cvband.c
View file

@ -1,415 +0,0 @@
/******************************************************************
* *
* File : cvband.c *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 24 February 2000 *
*----------------------------------------------------------------*
* This is the implementation file for the CVODE band linear *
* solver, CVBAND. *
* *
******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "cvband.h"
#include "cvode.h"
#include "band.h"
#include "llnltyps.h"
#include "nvector.h"
#include "llnlmath.h"
/* Error Messages */
#define CVBAND_INIT "CVBandInit-- "
#define MSG_MEM_FAIL CVBAND_INIT "A memory request failed.\n\n"
#define MSG_BAD_SIZES_1 CVBAND_INIT "Illegal bandwidth parameter(s) "
#define MSG_BAD_SIZES_2 "ml = %ld, mu = %ld.\n"
#define MSG_BAD_SIZES_3 "Must have 0 <= ml, mu <= N-1=%ld.\n\n"
#define MSG_BAD_SIZES MSG_BAD_SIZES_1 MSG_BAD_SIZES_2 MSG_BAD_SIZES_3
/* Other Constants */
#define MIN_INC_MULT RCONST(1000.0)
#define ZERO RCONST(0.0)
#define ONE RCONST(1.0)
#define TWO RCONST(2.0)
/******************************************************************
* *
* Types : CVBandMemRec, CVBandMem *
*----------------------------------------------------------------*
* The type CVBandMem is pointer to a CVBandMemRec. This *
* structure contains CVBand solver-specific data. *
* *
******************************************************************/
typedef struct {
CVBandJacFn b_jac; /* jac = Jacobian routine to be called */
integer b_ml; /* b_ml = lower bandwidth of savedJ */
integer b_mu; /* b_mu = upper bandwidth of savedJ */
integer b_storage_mu; /* upper bandwith of M = MIN(N-1,b_mu+b_ml) */
BandMat b_M; /* M = I - gamma J, gamma = h / l1 */
integer *b_pivots; /* pivots = pivot array for PM = LU */
BandMat b_savedJ; /* savedJ = old Jacobian */
long int b_nstlj; /* nstlj = nst at last Jacobian eval. */
long int b_nje; /* nje = no. of calls to jac */
void *b_J_data; /* J_data is passed to jac */
} CVBandMemRec, *CVBandMem;
/* CVBAND linit, lsetup, lsolve, and lfree routines */
static int CVBandInit(CVodeMem cv_mem, boole *setupNonNull);
static int CVBandSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
N_Vector fpred, boole *jcurPtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3);
static int CVBandSolve(CVodeMem cv_mem, N_Vector b, N_Vector ycur,
N_Vector fcur);
static void CVBandFree(CVodeMem cv_mem);
/*************** CVBandDQJac *****************************************
This routine generates a banded difference quotient approximation to
the Jacobian of f(t,y). It assumes that a band matrix of type
BandMat is stored column-wise, and that elements within each column
are contiguous. This makes it possible to get the address of a column
of J via the macro BAND_COL and to write a simple for loop to set
each of the elements of a column in succession.
**********************************************************************/
void CVBandDQJac(integer N, integer mupper, integer mlower, BandMat J,
RhsFn f, void *f_data, real tn, N_Vector y,
N_Vector fy, N_Vector ewt, real h, real uround,
void *jac_data, long int *nfePtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3)
{
real fnorm, minInc, inc, inc_inv, srur;
N_Vector ftemp, ytemp;
integer group, i, j, width, ngroups, i1, i2;
real *col_j, *ewt_data, *fy_data, *ftemp_data, *y_data, *ytemp_data;
/* Rename work vectors for use as temporary values of y and f */
ftemp = vtemp1;
ytemp = vtemp2;
/* Obtain pointers to the data for ewt, fy, ftemp, y, ytemp */
ewt_data = N_VDATA(ewt);
fy_data = N_VDATA(fy);
ftemp_data = N_VDATA(ftemp);
y_data = N_VDATA(y);
ytemp_data = N_VDATA(ytemp);
/* Load ytemp with y = predicted y vector */
N_VScale(ONE, y, ytemp);
/* Set minimum increment based on uround and norm of f */
srur = RSqrt(uround);
fnorm = N_VWrmsNorm(fy, ewt);
minInc = (fnorm != ZERO) ?
(MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;
/* Set bandwidth and number of column groups for band differencing */
width = mlower + mupper + 1;
ngroups = MIN(width, N);
for (group=1; group <= ngroups; group++) {
/* Increment all y_j in group */
for(j=group-1; j < N; j+=width) {
inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
ytemp_data[j] += inc;
}
/* Evaluate f with incremented y */
f(N, tn, ytemp, ftemp, f_data);
/* Restore ytemp, then form and load difference quotients */
for (j=group-1; j < N; j+=width) {
ytemp_data[j] = y_data[j];
col_j = BAND_COL(J,j);
inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
inc_inv = ONE/inc;
i1 = MAX(0, j-mupper);
i2 = MIN(j+mlower, N-1);
for (i=i1; i <= i2; i++)
BAND_COL_ELEM(col_j,i,j) =
inc_inv * (ftemp_data[i] - fy_data[i]);
}
}
/* Increment counter nfe = *nfePtr */
*nfePtr += ngroups;
}
/* Readability Replacements */
#define N (cv_mem->cv_N)
#define lmm (cv_mem->cv_lmm)
#define f (cv_mem->cv_f)
#define f_data (cv_mem->cv_f_data)
#define uround (cv_mem->cv_uround)
#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 gamrat (cv_mem->cv_gamrat)
#define ewt (cv_mem->cv_ewt)
#define nfe (cv_mem->cv_nfe)
#define errfp (cv_mem->cv_errfp)
#define iopt (cv_mem->cv_iopt)
#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 jac (cvband_mem->b_jac)
#define M (cvband_mem->b_M)
#define mu (cvband_mem->b_mu)
#define ml (cvband_mem->b_ml)
#define storage_mu (cvband_mem->b_storage_mu)
#define pivots (cvband_mem->b_pivots)
#define savedJ (cvband_mem->b_savedJ)
#define nstlj (cvband_mem->b_nstlj)
#define nje (cvband_mem->b_nje)
#define J_data (cvband_mem->b_J_data)
/*************** CVBand **********************************************
This routine initializes the memory record and sets various function
fields specific to the band linear solver module. CVBand sets the
cv_linit, cv_lsetup, cv_lsolve, and cv_lfree fields in (*cvode_mem)
to be CVBandInit, CVBandSetup, CVBandSolve, and CVBandFree,
respectively. It allocates memory for a structure of type
CVBandMemRec and sets the cv_lmem field in (*cvode_mem) to the
address of this structure. Finally, it sets b_J_data field in the
CVBandMemRec structure to be the input parameter jac_data, b_mu to
be mupper, b_ml to be mlower, and the b_jac field to be:
(1) the input parameter bjac if bjac != NULL or
(2) CVBandDQJac if bjac == NULL.
**********************************************************************/
void CVBand(void *cvode_mem, integer mupper, integer mlower, CVBandJacFn bjac,
void *jac_data)
{
CVodeMem cv_mem;
CVBandMem cvband_mem;
/* Return immediately if cvode_mem is NULL */
cv_mem = (CVodeMem) cvode_mem;
if (cv_mem == NULL) return; /* CVode reports this error */
/* Set four main function fields in cv_mem */
linit = CVBandInit;
lsetup = CVBandSetup;
lsolve = CVBandSolve;
lfree = CVBandFree;
/* Get memory for CVBandMemRec */
lmem = cvband_mem = (CVBandMem) malloc(sizeof(CVBandMemRec));
if (cvband_mem == NULL) return; /* CVBandInit reports this error */
/* Set Jacobian routine field to user's bjac or CVBandDQJac */
if (bjac == NULL) {
jac = CVBandDQJac;
} else {
jac = bjac;
}
J_data = jac_data;
/* Load half-bandwiths in cvband_mem */
ml = mlower;
mu = mupper;
}
/*************** CVBandInit ******************************************
This routine initializes remaining memory specific to the band
linear solver. If any memory request fails, all memory previously
allocated is freed, and an error message printed, before returning.
**********************************************************************/
static int CVBandInit(CVodeMem cv_mem, boole *setupNonNull)
{
CVBandMem cvband_mem;
cvband_mem = (CVBandMem) lmem;
/* Print error message and return if cvband_mem is NULL */
if (cvband_mem == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
return(LINIT_ERR);
}
/* Set flag setupNonNull = TRUE */
*setupNonNull = TRUE;
/* Test ml and mu for legality */
if ((ml < 0) || (mu < 0) || (ml >= N) || (mu >= N)) {
fprintf(errfp, MSG_BAD_SIZES, ml, mu, N-1);
return(LINIT_ERR);
}
/* Set extended upper half-bandwith for M (required for pivoting) */
storage_mu = MIN(N-1, mu + ml);
/* Allocate memory for M, savedJ, and pivot arrays */
M = BandAllocMat(N, mu, ml, storage_mu);
if (M == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
return(LINIT_ERR);
}
savedJ = BandAllocMat(N, mu, ml, mu);
if (savedJ == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
BandFreeMat(M);
return(LINIT_ERR);
}
pivots = BandAllocPiv(N);
if (pivots == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
BandFreeMat(M);
BandFreeMat(savedJ);
return(LINIT_ERR);
}
/* Initialize nje and nstlj, and set workspace lengths */
nje = 0;
if (iopt != NULL) {
iopt[BAND_NJE] = nje;
iopt[BAND_LRW] = N*(storage_mu + mu + 2*ml + 2);
iopt[BAND_LIW] = N;
}
nstlj = 0;
return(LINIT_OK);
}
/*************** CVBandSetup *****************************************
This routine does the setup operations for the band linear solver.
It makes a decision whether or not to call the Jacobian evaluation
routine based on various state variables, and if not it uses the
saved copy. In any case, it constructs the Newton matrix
M = I - gamma*J, updates counters, and calls the band LU
factorization routine.
**********************************************************************/
static int CVBandSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
N_Vector fpred, boole *jcurPtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3)
{
boole jbad, jok;
real dgamma;
integer ier;
CVBandMem cvband_mem;
cvband_mem = (CVBandMem) lmem;
/* Use nst, gamma/gammap, and convfail to set J eval. flag jok */
dgamma = ABS((gamma/gammap) - ONE);
jbad = (nst == 0) || (nst > nstlj + CVB_MSBJ) ||
((convfail == FAIL_BAD_J) && (dgamma < CVB_DGMAX)) ||
(convfail == FAIL_OTHER);
jok = !jbad;
if (jok) {
/* If jok = TRUE, use saved copy of J */
*jcurPtr = FALSE;
BandCopy(savedJ, M, mu, ml);
} else {
/* If jok = FALSE, call jac routine for new J value */
nje++;
if (iopt != NULL) iopt[BAND_NJE] = nje;
nstlj = nst;
*jcurPtr = TRUE;
BandZero(M);
jac(N, mu, ml, M, f, f_data, tn, ypred, fpred, ewt,
h, uround, J_data, &nfe, vtemp1, vtemp2, vtemp3);
BandCopy(M, savedJ, mu, ml);
}
/* Scale and add I to get M = I - gamma*J */
BandScale(-gamma, M);
BandAddI(M);
/* Do LU factorization of M */
ier = BandFactor(M, pivots);
/* Return 0 if the LU was complete; otherwise return 1 */
if (ier > 0) return(1);
return(0);
}
/*************** CVBandSolve *****************************************
This routine handles the solve operation for the band linear solver
by calling the band backsolve routine. The return value is 0.
**********************************************************************/
static int CVBandSolve(CVodeMem cv_mem, N_Vector b, N_Vector ycur,
N_Vector fcur)
{
CVBandMem cvband_mem;
cvband_mem = (CVBandMem) lmem;
BandBacksolve(M, pivots, b);
/* If BDF, scale the correction to account for change in gamma */
if ((lmm == BDF) && (gamrat != ONE)) {
N_VScale(TWO/(ONE + gamrat), b, b);
}
return(0);
}
/*************** CVBandFree ******************************************
This routine frees memory specific to the band linear solver.
**********************************************************************/
static void CVBandFree(CVodeMem cv_mem)
{
CVBandMem cvband_mem;
cvband_mem = (CVBandMem) lmem;
BandFreeMat(M);
BandFreeMat(savedJ);
BandFreePiv(pivots);
free(lmem);
}

318
ext/cvode/source/cvbandpre.c
View file

@ -1,318 +0,0 @@
/******************************************************************
* *
* File : cvbandpre.c *
* Programmers : Michael Wittman and Alan C. Hindmarsh @ LLNL *
* Version of : 23 March 2000 *
*----------------------------------------------------------------*
* This file contains implementations of the banded difference *
* quotient Jacobian-based preconditioner and solver routines for *
* use with CVSpgmr. *
* *
******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "cvbandpre.h"
#include "cvode.h"
#include "llnltyps.h"
#include "nvector.h"
#include "llnlmath.h"
#include "band.h"
#define MIN_INC_MULT RCONST(1000.0)
#define ZERO RCONST(0.0)
#define ONE RCONST(1.0)
/* Prototype for difference quotient Jacobian calculation routine */
static void CVBandPDQJac(integer N, integer mupper, integer mlower, BandMat J,
RhsFn f, void *f_data, real tn, N_Vector y,
N_Vector fy, N_Vector ewt, real h, real uround,
N_Vector ftemp, N_Vector ytemp);
/********************** Malloc and Free Functions **********************/
CVBandPreData CVBandPreAlloc(integer N, RhsFn f, void *f_data,
integer mu, integer ml)
{
CVBandPreData pdata;
integer mup, mlp, storagemu;
pdata = (CVBandPreData) malloc(sizeof *pdata); /* Allocate data memory */
if (pdata == NULL) return(NULL);
/* Load pointers and bandwidths into pdata block. */
pdata->f = f;
pdata->f_data = f_data;
pdata->mu = mup = MIN( N-1, MAX(0,mu) );
pdata->ml = mlp = MIN( N-1, MAX(0,ml) );
/* Allocate memory for saved banded Jacobian approximation. */
pdata->savedJ = BandAllocMat(N, mup, mlp, mup);
if (pdata->savedJ == NULL) {
free(pdata);
return(NULL);
}
/* Allocate memory for banded preconditioner. */
storagemu = MIN( N-1, mup + mlp);
pdata->savedP = BandAllocMat(N, mup, mlp, storagemu);
if (pdata->savedP == NULL) {
BandFreeMat(pdata->savedJ);
free(pdata);
return(NULL);
}
/* Allocate memory for pivot array. */
pdata->pivots = BandAllocPiv(N);
if (pdata->savedJ == NULL) {
BandFreeMat(pdata->savedP);
BandFreeMat(pdata->savedJ);
free(pdata);
return(NULL);
}
return(pdata);
}
void CVBandPreFree(CVBandPreData pdata)
{
BandFreeMat(pdata->savedJ);
BandFreeMat(pdata->savedP);
BandFreePiv(pdata->pivots);
free(pdata);
}
/***************** Preconditioner setup and solve functions *******/
/* Readability Replacements */
#define f (pdata->f)
#define f_data (pdata->f_data)
#define mu (pdata->mu)
#define ml (pdata->ml)
#define pivots (pdata->pivots)
#define savedJ (pdata->savedJ)
#define savedP (pdata->savedP)
/* Preconditioner setup routine CVBandPrecond. */
/******************************************************************
* Together CVBandPrecond and CVBandPSolve use a banded *
* difference quotient Jacobian to create a preconditioner. *
* CVBandPrecond calculates a new J, if necessary, then *
* calculates P = I - gamma*J, and does an LU factorization of P. *
* *
* The parameters of CVBandPrecond are as follows: *
* *
* N is the length of all vector arguments. *
* *
* t is the current value of the independent variable. *
* *
* y is the current value of the dependent variable vector, *
* namely the predicted value of y(t). *
* *
* fy is the vector f(t,y). *
* *
* jok is an input flag indicating whether Jacobian-related *
* data needs to be recomputed, as follows: *
* jok == FALSE means recompute Jacobian-related data *
* from scratch. *
* jok == TRUE means that Jacobian data from the *
* previous Precond call will be reused *
* (with the current value of gamma). *
* A CVBandPrecond call with jok == TRUE should only *
* occur after a call with jok == FALSE. *
* *
* jcurPtr is a pointer to an output integer flag which is *
* set by CVBandPrecond as follows: *
* *jcurPtr = TRUE if Jacobian data was recomputed. *
* *jcurPtr = FALSE if Jacobian data was not recomputed,*
* but saved data was reused. *
* *
* gamma is the scalar appearing in the Newton matrix. *
* *
* ewt is the error weight vector. *
* *
* h is a tentative step size in t. *
* *
* uround is the machine unit roundoff. *
* *
* nfePtr is a pointer to the memory location containing the *
* CVODE problem data nfe = number of calls to f. *
* The routine calls f a total of ml+mu+1 times, so *
* it increments *nfePtr by ml+mu+1. *
* *
* bp_data is a pointer to preconditoner data - the same as the *
* bp_data parameter passed to CVSpgmr. *
* *
* vtemp1, vtemp2, and vtemp3 are pointers to memory allocated *
* for vectors of length N for work space. This *
* routine uses only vtemp1 and vtemp2. *
* *
* *
* The value to be returned by the CVBandPrecond function is *
* 0 if successful, or *
* 1 if the band factorization failed. *
******************************************************************/
int CVBandPrecond(integer N, real t, N_Vector y, N_Vector fy,
boole jok, boole *jcurPtr, real gamma,
N_Vector ewt, real h, real uround,
long int *nfePtr, void *bp_data,
N_Vector vtemp1, N_Vector vtemp2,
N_Vector vtemp3)
{
integer ier;
CVBandPreData pdata;
/* Assume matrix and pivots have already been allocated. */
pdata = (CVBandPreData) bp_data;
if (jok) {
/* If jok = TRUE, use saved copy of J. */
*jcurPtr = FALSE;
BandCopy(savedJ, savedP, mu, ml);
} else {
/* If jok = FALSE, call CVBandPDQJac for new J value. */
*jcurPtr = TRUE;
BandZero(savedJ);
CVBandPDQJac(N, mu, ml, savedJ, f, f_data, t, y, fy, ewt,
h, uround, vtemp1, vtemp2);
BandCopy(savedJ, savedP, mu, ml);
*nfePtr += MIN( N, ml + mu + 1 );
}
/* Scale and add I to get savedP = I - gamma*J. */
BandScale(-gamma, savedP);
BandAddI(savedP);
/* Do LU factorization of matrix. */
ier = BandFactor(savedP, pivots);
/* Return 0 if the LU was complete; otherwise return 1. */
if (ier > 0) return(1);
return(0);
}
/* Preconditioner solve routine CVBandPSolve */
/******************************************************************
* CVBandPSolve solves a linear system P z = r, where P is the *
* matrix computed by CVBandPrecond. *
* *
* The parameters of CVBandPSolve used here are as follows: *
* *
* r is the right-hand side vector of the linear system. *
* *
* bp_data is a pointer to preconditioner data - the same as the *
* bp_data parameter passed to CVSpgmr. *
* *
* z is the output vector computed by CVBandPSolve. *
* *
* The value returned by the CVBandPSolve function is always 0, *
* indicating success. *
* *
******************************************************************/
int CVBandPSolve(integer N, real t, N_Vector y, N_Vector fy,
N_Vector vtemp, real gamma, N_Vector ewt,
real delta, long int *nfePtr, N_Vector r,
int lr, void *bp_data, N_Vector z)
{
CVBandPreData pdata;
/* Assume matrix and pivots have already been allocated. */
pdata = (CVBandPreData) bp_data;
/* Copy r to z. */
N_VScale(ONE, r, z);
/* Do band backsolve on the vector z. */
BandBacksolve(savedP, pivots, z);
return(0);
}
#undef f
#undef f_data
#undef mu
#undef ml
#undef pivots
#undef savedJ
#undef savedP
/*************** CVBandPDQJac ****************************************
This routine generates a banded difference quotient approximation to
the Jacobian of f(t,y). It assumes that a band matrix of type
BandMat is stored column-wise, and that elements within each column
are contiguous. This makes it possible to get the address of a column
of J via the macro BAND_COL and to write a simple for loop to set
each of the elements of a column in succession.
**********************************************************************/
static void CVBandPDQJac(integer N, integer mupper, integer mlower, BandMat J,
RhsFn f, void *f_data, real tn, N_Vector y,
N_Vector fy, N_Vector ewt, real h, real uround,
N_Vector ftemp, N_Vector ytemp)
{
real fnorm, minInc, inc, inc_inv, srur;
integer group, i, j, width, ngroups, i1, i2;
real *col_j, *ewt_data, *fy_data, *ftemp_data, *y_data, *ytemp_data;
/* Obtain pointers to the data for ewt, fy, ftemp, y, ytemp. */
ewt_data = N_VDATA(ewt);
fy_data = N_VDATA(fy);
ftemp_data = N_VDATA(ftemp);
y_data = N_VDATA(y);
ytemp_data = N_VDATA(ytemp);
/* Load ytemp with y = predicted y vector. */
N_VScale(ONE, y, ytemp);
/* Set minimum increment based on uround and norm of f. */
srur = RSqrt(uround);
fnorm = N_VWrmsNorm(fy, ewt);
minInc = (fnorm != ZERO) ?
(MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;
/* Set bandwidth and number of column groups for band differencing. */
width = mlower + mupper + 1;
ngroups = MIN(width, N);
for (group = 1; group <= ngroups; group++) {
/* Increment all y_j in group. */
for(j = group-1; j < N; j += width) {
inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
ytemp_data[j] += inc;
}
/* Evaluate f with incremented y. */
f(N, tn, ytemp, ftemp, f_data);
/* Restore ytemp, then form and load difference quotients. */
for (j = group-1; j < N; j += width) {
ytemp_data[j] = y_data[j];
col_j = BAND_COL(J,j);
inc = MAX(srur*ABS(y_data[j]), minInc/ewt_data[j]);
inc_inv = ONE/inc;
i1 = MAX(0, j-mupper);
i2 = MIN(j+mlower, N-1);
for (i=i1; i <= i2; i++)
BAND_COL_ELEM(col_j,i,j) =
inc_inv * (ftemp_data[i] - fy_data[i]);
}
}
}

372
ext/cvode/source/cvdense.c
View file

@ -1,372 +0,0 @@
/******************************************************************
* *
* File : cvdense.c *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 25 February 2000 *
*----------------------------------------------------------------*
* This is the implementation file for the CVODE dense linear *
* solver, CVDENSE. *
* *
******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "cvdense.h"
#include "cvode.h"
#include "dense.h"
#include "llnltyps.h"
#include "nvector.h"
#include "llnlmath.h"
/* Error Messages */
#define CVDENSE_INIT "CVDenseInit-- "
#define MSG_MEM_FAIL CVDENSE_INIT "A memory request failed.\n\n"
/* Other Constants */
#define MIN_INC_MULT RCONST(1000.0)
#define ZERO RCONST(0.0)
#define ONE RCONST(1.0)
#define TWO RCONST(2.0)
/******************************************************************
* *
* Types : CVDenseMemRec, CVDenseMem *
*----------------------------------------------------------------*
* The type CVDenseMem is pointer to a CVDenseMemRec. This *
* structure contains CVDense solver-specific data. *
* *
******************************************************************/
typedef struct {
CVDenseJacFn d_jac; /* jac = Jacobian routine to be called */
DenseMat d_M; /* M = I - gamma J, gamma = h / l1 */
integer *d_pivots; /* pivots = pivot array for PM = LU */
DenseMat d_savedJ; /* savedJ = old Jacobian */
long int d_nstlj; /* nstlj = nst at last Jacobian eval. */
long int d_nje; /* nje = no. of calls to jac */
void *d_J_data; /* J_data is passed to jac */
} CVDenseMemRec, *CVDenseMem;
/* CVDENSE linit, lsetup, lsolve, and lfree routines */
static int CVDenseInit(CVodeMem cv_mem, boole *setupNonNull);
static int CVDenseSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
N_Vector fpred, boole *jcurPtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3);
static int CVDenseSolve(CVodeMem cv_mem, N_Vector b, N_Vector ycur,
N_Vector fcur);
static void CVDenseFree(CVodeMem cv_mem);
/*************** CVDenseDQJac ****************************************
This routine generates a dense difference quotient approximation to
the Jacobian of f(t,y). It assumes that a dense matrix of type
DenseMat is stored column-wise, and that elements within each column
are contiguous. The address of the jth column of J is obtained via
the macro DENSE_COL and an N_Vector with the jth column as the
component array is created using N_VMAKE and N_VDATA. Finally, the
actual computation of the jth column of the Jacobian is done with a
call to N_VLinearSum.
**********************************************************************/
void CVDenseDQJac(integer N, DenseMat J, RhsFn f, void *f_data, real tn,
N_Vector y, N_Vector fy, N_Vector ewt, real h, real uround,
void *jac_data, long int *nfePtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3)
{
real fnorm, minInc, inc, inc_inv, yjsaved, srur;
real *y_data, *ewt_data;
N_Vector ftemp, jthCol;
integer j;
ftemp = vtemp1; /* Rename work vector for use as f vector value */
/* Obtain pointers to the data for ewt, y */
ewt_data = N_VDATA(ewt);
y_data = N_VDATA(y);
/* Set minimum increment based on uround and norm of f */
srur = RSqrt(uround);
fnorm = N_VWrmsNorm(fy, ewt);
minInc = (fnorm != ZERO) ?
(MIN_INC_MULT * ABS(h) * uround * N * fnorm) : ONE;
N_VMAKE(jthCol, y_data, N); /* j loop overwrites this data address */
/* This is the only for loop for 0..N-1 in CVODE */
for (j = 0; j < N; j++) {
/* Generate the jth col of J(tn,y) */
N_VDATA(jthCol) = DENSE_COL(J,j);
yjsaved = y_data[j];
inc = MAX(srur*ABS(yjsaved), minInc/ewt_data[j]);
y_data[j] += inc;
f(N, tn, y, ftemp, f_data);
inc_inv = ONE/inc;
N_VLinearSum(inc_inv, ftemp, -inc_inv, fy, jthCol);
y_data[j] = yjsaved;
}
N_VDISPOSE(jthCol);
/* Increment counter nfe = *nfePtr */
*nfePtr += N;
}
/* Readability Replacements */
#define N (cv_mem->cv_N)
#define lmm (cv_mem->cv_lmm)
#define f (cv_mem->cv_f)
#define f_data (cv_mem->cv_f_data)
#define uround (cv_mem->cv_uround)
#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 gamrat (cv_mem->cv_gamrat)
#define ewt (cv_mem->cv_ewt)
#define nfe (cv_mem->cv_nfe)
#define errfp (cv_mem->cv_errfp)
#define iopt (cv_mem->cv_iopt)
#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 jac (cvdense_mem->d_jac)
#define M (cvdense_mem->d_M)
#define pivots (cvdense_mem->d_pivots)
#define savedJ (cvdense_mem->d_savedJ)
#define nstlj (cvdense_mem->d_nstlj)
#define nje (cvdense_mem->d_nje)
#define J_data (cvdense_mem->d_J_data)
/*************** CVDense *********************************************
This routine initializes the memory record and sets various function
fields specific to the dense linear solver module. CVDense sets the
cv_linit, cv_lsetup, cv_lsolve, and cv_lfree fields in (*cvode_mem)
to be CVDenseInit, CVDenseSetup, CVDenseSolve, and CVDenseFree,
respectively. It allocates memory for a structure of type
CVDenseMemRec and sets the cv_lmem field in (*cvode_mem) to the
address of this structure. Finally, it sets d_J_data field in the
CVDenseMemRec structure to be the input parameter jac_data and the
d_jac field to be:
(1) the input parameter djac if djac != NULL or
(2) CVDenseDQJac if djac == NULL.
**********************************************************************/
void CVDense(void *cvode_mem, CVDenseJacFn djac, void *jac_data)
{
CVodeMem cv_mem;
CVDenseMem cvdense_mem;
/* Return immediately if cvode_mem is NULL */
cv_mem = (CVodeMem) cvode_mem;
if (cv_mem == NULL) return; /* CVode reports this error */
/* Set four main function fields in cv_mem */
linit = CVDenseInit;
lsetup = CVDenseSetup;
lsolve = CVDenseSolve;
lfree = CVDenseFree;
/* Get memory for CVDenseMemRec */
lmem = cvdense_mem = (CVDenseMem) malloc(sizeof(CVDenseMemRec));
if (cvdense_mem == NULL) return; /* CVDenseInit reports this error */
/* Set Jacobian routine field to user's djac or CVDenseDQJac */
if (djac == NULL) {
jac = CVDenseDQJac;
} else {
jac = djac;
}
J_data = jac_data;
}
/*************** CVDenseInit *****************************************
This routine initializes remaining memory specific to the dense
linear solver. If any memory request fails, all memory previously
allocated is freed, and an error message printed, before returning.
**********************************************************************/
static int CVDenseInit(CVodeMem cv_mem, boole *setupNonNull)
{
CVDenseMem cvdense_mem;
cvdense_mem = (CVDenseMem) lmem;
/* Print error message and return if cvdense_mem is NULL */
if (cvdense_mem == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
return(LINIT_ERR);
}
/* Set flag setupNonNull = TRUE */
*setupNonNull = TRUE;
/* Allocate memory for M, savedJ, and pivot array */
M = DenseAllocMat(N);
if (M == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
return(LINIT_ERR);
}
savedJ = DenseAllocMat(N);
if (savedJ == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
DenseFreeMat(M);
return(LINIT_ERR);
}
pivots = DenseAllocPiv(N);
if (pivots == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
DenseFreeMat(M);
DenseFreeMat(savedJ);
return(LINIT_ERR);
}
/* Initialize nje and nstlj, and set workspace lengths */
nje = 0;
if (iopt != NULL) {
iopt[DENSE_NJE] = nje;
iopt[DENSE_LRW] = 2*N*N;
iopt[DENSE_LIW] = N;
}
nstlj = 0;
return(LINIT_OK);
}
/*************** CVDenseSetup ****************************************
This routine does the setup operations for the dense linear solver.
It makes a decision whether or not to call the Jacobian evaluation
routine based on various state variables, and if not it uses the
saved copy. In any case, it constructs the Newton matrix
M = I - gamma*J, updates counters, and calls the dense LU
factorization routine.
**********************************************************************/
static int CVDenseSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
N_Vector fpred, boole *jcurPtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3)
{
boole jbad, jok;
real dgamma;
integer ier;
CVDenseMem cvdense_mem;
cvdense_mem = (CVDenseMem) lmem;
/* Use nst, gamma/gammap, and convfail to set J eval. flag jok */
dgamma = ABS((gamma/gammap) - ONE);
jbad = (nst == 0) || (nst > nstlj + CVD_MSBJ) ||
((convfail == FAIL_BAD_J) && (dgamma < CVD_DGMAX)) ||
(convfail == FAIL_OTHER);
jok = !jbad;
if (jok) {
/* If jok = TRUE, use saved copy of J */
*jcurPtr = FALSE;
DenseCopy(savedJ, M);
} else {
/* If jok = FALSE, call jac routine for new J value */
nje++;
if (iopt != NULL) iopt[DENSE_NJE] = nje;
nstlj = nst;
*jcurPtr = TRUE;
DenseZero(M);
jac(N, M, f, f_data, tn, ypred, fpred, ewt, h,
uround, J_data, &nfe, vtemp1, vtemp2, vtemp3);
DenseCopy(M, savedJ);
}
/* Scale and add I to get M = I - gamma*J */
DenseScale(-gamma, M);
DenseAddI(M);
/* Do LU factorization of M */
ier = DenseFactor(M, pivots);
/* Return 0 if the LU was complete; otherwise return 1 */
if (ier > 0) return(1);
return(0);
}
/*************** CVDenseSolve ****************************************
This routine handles the solve operation for the dense linear solver
by calling the dense backsolve routine. The returned value is 0.
**********************************************************************/
static int CVDenseSolve(CVodeMem cv_mem, N_Vector b, N_Vector ycur,
N_Vector fcur)
{
CVDenseMem cvdense_mem;
cvdense_mem = (CVDenseMem) lmem;
DenseBacksolve(M, pivots, b);
/* If BDF, scale the correction to account for change in gamma */
if ((lmm == BDF) && (gamrat != ONE)) {
N_VScale(TWO/(ONE + gamrat), b, b);
}
return(0);
}
/*************** CVDenseFree *****************************************
This routine frees memory specific to the dense linear solver.
**********************************************************************/
static void CVDenseFree(CVodeMem cv_mem)
{
CVDenseMem cvdense_mem;
cvdense_mem = (CVDenseMem) lmem;
DenseFreeMat(M);
DenseFreeMat(savedJ);
DenseFreePiv(pivots);
free(lmem);
}

292
ext/cvode/source/cvdiag.c
View file

@ -1,292 +0,0 @@
/******************************************************************
* *
* File : cvdiag.c *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 4 May 1998 *
*----------------------------------------------------------------*
* This is the implementation file for the CVODE diagonal linear *
* solver, CVDIAG. *
* *
******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "cvdiag.h"
#include "cvode.h"
#include "llnltyps.h"
#include "nvector.h"
/* Error Messages */
#define CVDIAG_INIT "CVDiagInit-- "
#define MSG_MEM_FAIL CVDIAG_INIT "A memory request failed.\n\n"
/* Other Constants */
#define FRACT RCONST(0.1)
#define ONE RCONST(1.0)
/******************************************************************
* *
* Types : CVDiagMemRec, CVDiagMem *
*----------------------------------------------------------------*
* The type CVDiagMem is pointer to a CVDiagMemRec. This *
* structure contains CVDiag solver-specific data. *
* *
******************************************************************/
typedef struct {
real di_gammasv; /* gammasv = gamma at the last call to setup */
/* or solve */
N_Vector di_M; /* M = (I - gamma J)^{-1} , gamma = h / l1 */
N_Vector di_bit; /* temporary storage vector */
N_Vector di_bitcomp; /* temporary storage vector */
} CVDiagMemRec, *CVDiagMem;
/* CVDIAG linit, lsetup, lsolve, and lfree routines */
static int CVDiagInit(CVodeMem cv_mem, boole *setupNonNull);
static int CVDiagSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
N_Vector fpred, boole *jcurPtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3);
static int CVDiagSolve(CVodeMem cv_mem, N_Vector b, N_Vector ycur,
N_Vector fcur);
static void CVDiagFree(CVodeMem cv_mem);
/* Readability Replacements */
#define N (cv_mem->cv_N)
#define f (cv_mem->cv_f)
#define f_data (cv_mem->cv_f_data)
#define uround (cv_mem->cv_uround)
#define tn (cv_mem->cv_tn)
#define h (cv_mem->cv_h)
#define rl1 (cv_mem->cv_rl1)
#define gamma (cv_mem->cv_gamma)
#define ewt (cv_mem->cv_ewt)
#define nfe (cv_mem->cv_nfe)
#define errfp (cv_mem->cv_errfp)
#define iopt (cv_mem->cv_iopt)
#define zn (cv_mem->cv_zn)
#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 gammasv (cvdiag_mem->di_gammasv)
#define M (cvdiag_mem->di_M)
#define bit (cvdiag_mem->di_bit)
#define bitcomp (cvdiag_mem->di_bitcomp)
/*************** CVDiag **********************************************
This routine initializes the memory record and sets various function
fields specific to the diagonal linear solver module. CVDiag sets the
cv_linit, cv_lsetup, cv_lsolve, and cv_lfree fields in (*cvode_mem)
to be CVDiagInit, CVDiagSetup, CVDiagSolve, and CVDiagFree,
respectively. It allocates memory for a structure of type
CVDiagMemRec and sets the cv_lmem field in (*cvode_mem) to the
address of this structure.
**********************************************************************/
void CVDiag(void *cvode_mem)
{
CVodeMem cv_mem;
CVDiagMem cvdiag_mem;
/* Return immediately if cvode_mem is NULL */
cv_mem = (CVodeMem) cvode_mem;
if (cv_mem == NULL) return; /* CVode reports this error */
/* Set four main function fields in cv_mem */
linit = CVDiagInit;
lsetup = CVDiagSetup;
lsolve = CVDiagSolve;
lfree = CVDiagFree;
/* Get memory for CVDiagMemRec */
lmem = cvdiag_mem = (CVDiagMem) malloc(sizeof(CVDiagMemRec));
if (cvdiag_mem == NULL) return; /* CVDiagInit reports this error */
}
/*************** CVDiagInit ******************************************
This routine initializes remaining memory specific to the diagonal
linear solver. If any memory request fails, all memory previously
allocated is freed, and an error message printed, before returning.
**********************************************************************/
static int CVDiagInit(CVodeMem cv_mem, boole *setupNonNull)
{
CVDiagMem cvdiag_mem;
cvdiag_mem = (CVDiagMem) lmem;
/* Print error message and return if cvdiag_mem is NULL */
if (cvdiag_mem == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
return(LINIT_ERR);
}
/* Set flag setupNonNull = TRUE */
*setupNonNull = TRUE;
/* Allocate memory for M, bit, and bitcomp */
M = N_VNew(N, machenv);
if (M == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
return(LINIT_ERR);
}
bit = N_VNew(N, machenv);
if (bit == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
N_VFree(M);
return(LINIT_ERR);
}
bitcomp = N_VNew(N, machenv);
if (bitcomp == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
N_VFree(M);
N_VFree(bit);
return(LINIT_ERR);
}
/* Set workspace lengths */
if (iopt != NULL) {
iopt[DIAG_LRW] = N*3;
iopt[DIAG_LIW] = 0;
}
return(LINIT_OK);
}
/*************** CVDiagSetup *****************************************
This routine does the setup operations for the diagonal linear
solver. It constructs a diagonal approximation to the Newton matrix
M = I - gamma*J, updates counters, and inverts M.
**********************************************************************/
static int CVDiagSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
N_Vector fpred, boole *jcurPtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3)
{
real r;
N_Vector ftemp, y;
boole invOK;
CVDiagMem cvdiag_mem;
cvdiag_mem = (CVDiagMem) lmem;
/* Rename work vectors for use as temporary values of y and f */
ftemp = vtemp1;
y = vtemp2;
/* Form y with perturbation = FRACT*(func. iter. correction) */
r = FRACT * rl1;
N_VLinearSum(h, fpred, -ONE, zn[1], ftemp);
N_VLinearSum(r, ftemp, ONE, ypred, y);
/* Evaluate f at perturbed y */
f(N, tn, y, M, f_data);
nfe++;
/* Construct M = I - gamma*J with J = diag(deltaf_i/deltay_i) */
N_VLinearSum(ONE, M, -ONE, fpred, M);
N_VLinearSum(FRACT, ftemp, -h, M, M);
N_VProd(ftemp, ewt, y);
/* Protect against deltay_i being at roundoff level */
N_VCompare(uround, y, bit);
N_VAddConst(bit, -ONE, bitcomp);
N_VProd(ftemp, bit, y);
N_VLinearSum(FRACT, y, -ONE, bitcomp, y);
N_VDiv(M, y, M);
N_VProd(M, bit, M);
N_VLinearSum(ONE, M, -ONE, bitcomp, M);
/* Invert M with test for zero components */
invOK = N_VInvTest(M, M);
if (!invOK) return(1);
/* Set jcur = TRUE, save gamma in gammasv, and return */
*jcurPtr = TRUE;
gammasv = gamma;
return(0);
}
/*************** CVDiagSolve *****************************************
This routine performs the solve operation for the diagonal linear
solver. If necessary it first updates gamma in M = I - gamma*J.
**********************************************************************/
static int CVDiagSolve(CVodeMem cv_mem, N_Vector b, N_Vector ycur,
N_Vector fcur)
{
boole invOK;
real r;
CVDiagMem cvdiag_mem;
cvdiag_mem = (CVDiagMem) lmem;
/* If gamma has changed, update factor in M, and save gamma value */
if (gammasv != gamma) {
r = gamma / gammasv;
N_VInv(M, M);
N_VAddConst(M, -ONE, M);
N_VScale(r, M, M);
N_VAddConst(M, ONE, M);
invOK = N_VInvTest(M, M);
if (!invOK) return (1);
gammasv = gamma;
}
/* Apply M-inverse to b */
N_VProd(b, M, b);
return(0);
}
/*************** CVDiagFree ******************************************
This routine frees memory specific to the diagonal linear solver.
**********************************************************************/
static void CVDiagFree(CVodeMem cv_mem)
{
CVDiagMem cvdiag_mem;
cvdiag_mem = (CVDiagMem) lmem;
N_VFree(M);
N_VFree(bit);
N_VFree(bitcomp);
free(lmem);
}

2634
ext/cvode/source/cvode.c
View file

File diff suppressed because it is too large Load diff

499
ext/cvode/source/cvspgmr.c
View file

@ -1,499 +0,0 @@
/******************************************************************
* *
* File : cvspgmr.c *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 25 February 2000 *
*----------------------------------------------------------------*
* 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 "llnltyps.h"
#include "nvector.h"
#include "llnlmath.h"
#include "iterativ.h"
#include "spgmr.h"
/* Error Messages */
#define CVSPGMR_INIT "CVSpgmrInit-- "
#define MSG_MEM_FAIL CVSPGMR_INIT "A memory request failed.\n\n"
#define MSG_BAD_PRETYPE_1 CVSPGMR_INIT "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_INIT "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_INIT "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 */
real g_sqrtN; /* sqrt(N) */
real g_delt; /* delt = user specified or DELT_DEFAULT */
real g_deltar; /* deltar = delt * tq4 */
real 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 used by CVAtimesDQ */
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 */
} CVSpgmrMemRec, *CVSpgmrMem;
/* CVSPGMR linit, lsetup, lsolve, and lfree routines */
static int CVSpgmrInit(CVodeMem cv_mem, boole *setupNonNull);
static int CVSpgmrSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
N_Vector fpred, boole *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 CVSpgmrAtimesDQ(void *cv_mem, N_Vector v, N_Vector z);
static int CVSpgmrPSolve(void *cv_mem, N_Vector r, N_Vector z, int lr);
/* 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 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 sets the
cv_linit, cv_lsetup, cv_lsolve, and 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. CVSpgmr sets the following fields in the
CVSpgmrMemRec structure:
g_pretype = pretype
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
**********************************************************************/
void CVSpgmr(void *cvode_mem, int pretype, int gstype, int maxl, real delt,
CVSpgmrPrecondFn precond, CVSpgmrPSolveFn psolve, void *P_data)
{
CVodeMem cv_mem;
CVSpgmrMem cvspgmr_mem;
/* Return immediately if cvode_mem is NULL */
cv_mem = (CVodeMem) cvode_mem;
if (cv_mem == NULL) return; /* CVode reports this error */
/* 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) return; /* CVSpgmrInit reports this error */
/* Set Spgmr parameters that have been passed in call sequence */
cvspgmr_mem->g_pretype = pretype;
cvspgmr_mem->g_gstype = gstype;
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;
}
/* 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)
/*************** CVSpgmrInit *****************************************
This routine initializes remaining memory specific to the Spgmr
linear solver. If any memory request fails, all memory previously
allocated is freed, and an error message printed, before returning.
**********************************************************************/
static int CVSpgmrInit(CVodeMem cv_mem, boole *setupNonNull)
{
CVSpgmrMem cvspgmr_mem;
cvspgmr_mem = (CVSpgmrMem) lmem;
/* Print error message and return if cvspgmr_mem is NULL */
if (cvspgmr_mem == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
return(LINIT_ERR);
}
/* 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(LINIT_ERR);
}
if ((pretype != NONE) && (psolve == NULL)) {
fprintf(errfp, MSG_PSOLVE_REQ);
return(LINIT_ERR);
}
/* Check for legal gstype */
if ((gstype != MODIFIED_GS) && (gstype != CLASSICAL_GS)) {
fprintf(errfp, MSG_BAD_GSTYPE, gstype, MODIFIED_GS, CLASSICAL_GS);
return(LINIT_ERR);
}
/* Allocate memory for ytemp and x */
ytemp = N_VNew(N, machenv);
if (ytemp == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
return(LINIT_ERR);
}
x = N_VNew(N, machenv);
if (x == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
N_VFree(ytemp);
return(LINIT_ERR);
}
/* Call SpgmrMalloc to allocate workspace for Spgmr */
spgmr_mem = SpgmrMalloc(N, maxl, machenv);
if (spgmr_mem == NULL) {
fprintf(errfp, MSG_MEM_FAIL);
N_VFree(ytemp);
N_VFree(x);
return(LINIT_ERR);
}
/* 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;
}
/* Set setupNonNull to TRUE iff there is preconditioning */
/* (pretype != NONE) and there is a preconditioning setup phase */
/* (precond != NULL) */
*setupNonNull = (pretype != NONE) && (precond != NULL);
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, boole *jcurPtr, N_Vector vtemp1,
N_Vector vtemp2, N_Vector vtemp3)
{
boole jbad, jok;
real 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)
{
real 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, CVSpgmrAtimesDQ, 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(lmem);
}
/*************** CVSpgmrAtimesDQ *************************************
This routine generates the matrix-vector product z = Mv, where
M = I - gamma*J, by using a difference quotient approximation to
the product Jv. The approximation is Jv = rho[f(y + v/rho) - f(y)],
where rho = (WRMS norm of v), i.e. the WRMS norm of v/rho is 1.
**********************************************************************/
static int CVSpgmrAtimesDQ(void *cvode_mem, N_Vector v, N_Vector z)
{
real rho;
CVodeMem cv_mem;
CVSpgmrMem cvspgmr_mem;
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);
}
/* Set ytemp = ycur + (1/rho) v */
N_VLinearSum(ONE/rho, v, ONE, ycur, ytemp);
/* Set z = f(tn, ytemp) */
f(N, tn, ytemp, z, f_data);
nfe++;
/* Replace z by v - (gamma*rho)(z - fcur) */
N_VLinearSum(ONE, z, -ONE, fcur, z);
N_VLinearSum(-gamma*rho, z, ONE, v, 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);
}

311
ext/cvode/source/dense.c
View file

@ -1,311 +0,0 @@
/******************************************************************
* *
* File : dense.c *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 25 February 2000 *
*----------------------------------------------------------------*
* This is the implementation file for a generic DENSE linear *
* solver package. *
* *
******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "dense.h"
#include "llnltyps.h"
#include "nvector.h"
#include "llnlmath.h"
#define ZERO RCONST(0.0)
#define ONE RCONST(1.0)
/* Implementation */
DenseMat DenseAllocMat(integer N)
{
DenseMat A;
if (N <= 0) return(NULL);
A = (DenseMat) malloc(sizeof *A);
if (A==NULL) return (NULL);
A->data = denalloc(N);
if (A->data == NULL) {
free(A);
return(NULL);
}
A->size = N;
return(A);
}
integer *DenseAllocPiv(integer N)
{
if (N <= 0) return(NULL);
return((integer *) malloc(N * sizeof(integer)));
}
integer DenseFactor(DenseMat A, integer *p)
{
return(gefa(A->data, A->size, p));
}
void DenseBacksolve(DenseMat A, integer *p, N_Vector b)
{
gesl(A->data, A->size, p, N_VDATA(b));
}
void DenseZero(DenseMat A)
{
denzero(A->data, A->size);
}
void DenseCopy(DenseMat A, DenseMat B)
{
dencopy(A->data, B->data, A->size);
}
void DenseScale(real c, DenseMat A)
{
denscale(c, A->data, A->size);
}
void DenseAddI(DenseMat A)
{
denaddI(A->data, A->size);
}
void DenseFreeMat(DenseMat A)
{
denfree(A->data);
free(A);
}
void DenseFreePiv(integer *p)
{
free(p);
}
void DensePrint(DenseMat A)
{
denprint(A->data, A->size);
}
real **denalloc(integer n)
{
integer j;
real **a;
if (n <= 0) return(NULL);
a = (real **) malloc(n * sizeof(real *));
if (a == NULL) return(NULL);
a[0] = (real *) malloc(n * n * sizeof(real));
if (a[0] == NULL) {
free(a);
return(NULL);
}
for (j=1; j < n; j++) a[j] = a[0] + j * n;
return(a);
}
integer *denallocpiv(integer n)
{
if (n <= 0) return(NULL);
return((integer *) malloc(n * sizeof(integer)));
}
integer gefa(real **a, integer n, integer *p)
{
integer i, j, k, l;
real *col_j, *col_k, *diag_k;
real temp, mult, a_kj;
boole swap;
/* k = elimination step number */
for (k=0; k < n-1; k++, p++) {
col_k = a[k];
diag_k = col_k + k;
/* find l = pivot row number */
l=k;
for (i=k+1; i < n; i++)
if (ABS(col_k[i]) > ABS(col_k[l])) l=i;
*p = l;
/* check for zero pivot element */
if (col_k[l] == ZERO) return(k+1);
/* swap a(l,k) and a(k,k) if necessary */
if ( (swap = (l != k) )) {
temp = col_k[l];
col_k[l] = *diag_k;
*diag_k = temp;
}
/* Scale the elements below the diagonal in */
/* column k by -1.0 / a(k,k). After the above swap, */
/* a(k,k) holds the pivot element. This scaling */
/* stores the pivot row multipliers -a(i,k)/a(k,k) */
/* in a(i,k), i=k+1, ..., n-1. */
mult = -ONE / (*diag_k);
for(i=k+1; i < n; i++)
col_k[i] *= mult;
/* row_i = row_i - [a(i,k)/a(k,k)] row_k, i=k+1, ..., n-1 */
/* row k is the pivot row after swapping with row l. */
/* The computation is done one column at a time, */
/* column j=k+1, ..., n-1. */
for (j=k+1; j < n; j++) {
col_j = a[j];
a_kj = col_j[l];
/* Swap the elements a(k,j) and a(k,l) if l!=k. */
if (swap) {
col_j[l] = col_j[k];
col_j[k] = a_kj;
}
/* a(i,j) = a(i,j) - [a(i,k)/a(k,k)]*a(k,j) */
/* a_kj = a(k,j), col_k[i] = - a(i,k)/a(k,k) */
if (a_kj != ZERO) {
for (i=k+1; i < n; i++)
col_j[i] += a_kj * col_k[i];
}
}
}
/* set the last pivot row to be n-1 and check for a zero pivot */
*p = n-1;
if (a[n-1][n-1] == ZERO) return(n);
/* return 0 to indicate success */
return(0);
}
void gesl(real **a, integer n, integer *p, real *b)
{
integer k, l, i;
real mult, *col_k;
/* Solve Ly = Pb, store solution y in b */
for (k=0; k < n-1; k++) {
l = p[k];
mult = b[l];
if (l != k) {
b[l] = b[k];
b[k] = mult;
}
col_k = a[k];
for (i=k+1; i < n; i++)
b[i] += mult*col_k[i];
}
/* Solve Ux = y, store solution x in b */
for (k=n-1; k >= 0; k--) {
col_k = a[k];
b[k] /= col_k[k];
mult = -b[k];
for (i=0; i < k; i++)
b[i] += mult*col_k[i];
}
}
void denzero(real **a, integer n)
{
integer i, j;
real *col_j;
for (j=0; j < n; j++) {
col_j = a[j];
for (i=0; i < n; i++)
col_j[i] = ZERO;
}
}
void dencopy(real **a, real **b, integer n)
{
integer i, j;
real *a_col_j, *b_col_j;
for (j=0; j < n; j++) {
a_col_j = a[j];
b_col_j = b[j];
for (i=0; i < n; i++)
b_col_j[i] = a_col_j[i];
}
}
void denscale(real c, real **a, integer n)
{
integer i, j;
real *col_j;
for (j=0; j < n; j++) {
col_j = a[j];
for (i=0; i < n; i++)
col_j[i] *= c;
}
}
void denaddI(real **a, integer n)
{
integer i;
for (i=0; i < n; i++) a[i][i] += ONE;
}
void denfreepiv(integer *p)
{
free(p);
}
void denfree(real **a)
{
free(a[0]);
free(a);
}
void denprint(real **a, integer n)
{
integer i, j;
printf("\n");
for (i=0; i < n; i++) {
for (j=0; j < n; j++) {
printf("%10g", a[j][i]);
}
printf("\n");
}
printf("\n");
}

258
ext/cvode/source/iterativ.c
View file

@ -1,258 +0,0 @@
/******************************************************************
* *
* File : iterativ.c *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 16 January 1998 *
*----------------------------------------------------------------*
* This is the implementation file for the iterativ.h header *
* file. It contains the implementation of functions that may be *
* useful for many different iterative solvers of A x = b. *
* *
******************************************************************/
#include "iterativ.h"
#include "llnltyps.h"
#include "nvector.h"
#include "llnlmath.h"
#define FACTOR RCONST(1000.0)
#define ZERO RCONST(0.0)
#define ONE RCONST(1.0)
/************************* ModifiedGS ***********************************
This implementation of ModifiedGS is a slight modification of a previous
modified Gram-Schmidt routine (called mgs) written by Milo Dorr.
*************************************************************************/
int ModifiedGS(N_Vector *v, real **h, int k, int p, real *new_vk_norm)
{
int i, k_minus_1, i0;
real new_norm_2, new_product, vk_norm, temp;
vk_norm = RSqrt(N_VDotProd(v[k],v[k]));
k_minus_1 = k - 1;
i0 = MAX(k-p, 0);
/* Perform modified Gram-Schmidt */
for (i=i0; i < k; i++) {
h[i][k_minus_1] = N_VDotProd(v[i], v[k]);
N_VLinearSum(ONE, v[k], -h[i][k_minus_1], v[i], v[k]);
}
/* Compute the norm of the new vector at v[k]. */
*new_vk_norm = RSqrt(N_VDotProd(v[k], v[k]));
/* If the norm of the new vector at v[k] is less than
FACTOR (== 1000) times unit roundoff times the norm of the
input vector v[k], then the vector will be reorthogonalized
in order to ensure that nonorthogonality is not being masked
by a very small vector length. */
temp = FACTOR * vk_norm;
if ((temp + (*new_vk_norm)) != temp) return(0);
new_norm_2 = ZERO;
for (i=i0; i < k; i++) {
new_product = N_VDotProd(v[i], v[k]);
temp = FACTOR * h[i][k_minus_1];
if ((temp + new_product) == temp) continue;
h[i][k_minus_1] += new_product;
N_VLinearSum(ONE, v[k],-new_product, v[i], v[k]);
new_norm_2 += SQR(new_product);
}
if (new_norm_2 != ZERO) {
new_product = SQR(*new_vk_norm) - new_norm_2;
*new_vk_norm = (new_product > ZERO) ? RSqrt(new_product) : ZERO;
}
return(0);
}
/************************ ClassicalGS ********************************
This implementation of ClassicalGS was contributed to by Homer Walker
and Peter Brown.
**********************************************************************/
int ClassicalGS(N_Vector *v, real **h, int k, int p, real *new_vk_norm,
N_Vector temp, real *s)
{
int i, k_minus_1, i0;
real vk_norm;
k_minus_1 = k - 1;
/* Perform Classical Gram-Schmidt */
vk_norm = RSqrt(N_VDotProd(v[k], v[k]));
i0 = MAX(k-p, 0);
for (i=i0; i < k; i++) {
h[i][k_minus_1] = N_VDotProd(v[i], v[k]);
}
for (i=i0; i < k; i++) {
N_VLinearSum(ONE, v[k], -h[i][k_minus_1], v[i], v[k]);
}
/* Compute the norm of the new vector at v[k]. */
*new_vk_norm = RSqrt(N_VDotProd(v[k], v[k]));
/* Reorthogonalize if necessary */
if ((FACTOR * (*new_vk_norm)) < vk_norm) {
for (i=i0; i < k; i++) {
s[i] = N_VDotProd(v[i], v[k]);
}
if (i0 < k) {
N_VScale(s[i0], v[i0], temp);
h[i0][k_minus_1] += s[i0];
}
for (i=i0+1; i < k; i++) {
N_VLinearSum(s[i], v[i], ONE, temp, temp);
h[i][k_minus_1] += s[i];
}
N_VLinearSum(ONE, v[k], -ONE, temp, v[k]);
*new_vk_norm = RSqrt(N_VDotProd(v[k],v[k]));
}
return(0);
}
/*************** QRfact **********************************************
This implementation of QRfact is a slight modification of a previous
routine (called qrfact) written by Milo Dorr.
**********************************************************************/
int QRfact(int n, real **h, real *q, int job)
{
real c, s, temp1, temp2, temp3;
int i, j, k, q_ptr, n_minus_1, code=0;
switch (job) {
case 0:
/* Compute a new factorization of H. */
code = 0;
for (k=0; k < n; k++) {
/* Multiply column k by the previous k-1 Givens rotations. */
for (j=0; j < k-1; j++) {
i = 2*j;
temp1 = h[j][k];
temp2 = h[j+1][k];
c = q[i];
s = q[i+1];
h[j][k] = c*temp1 - s*temp2;
h[j+1][k] = s*temp1 + c*temp2;
}
/* Compute the Givens rotation components c and s */
q_ptr = 2*k;
temp1 = h[k][k];
temp2 = h[k+1][k];
if( temp2 == ZERO) {
c = ONE;
s = ZERO;
} else if (ABS(temp2) >= ABS(temp1)) {
temp3 = temp1/temp2;
s = -ONE/RSqrt(ONE+SQR(temp3));
c = -s*temp3;
} else {
temp3 = temp2/temp1;
c = ONE/RSqrt(ONE+SQR(temp3));
s = -c*temp3;
}
q[q_ptr] = c;
q[q_ptr+1] = s;
if( (h[k][k] = c*temp1 - s*temp2) == ZERO) code = k+1;
}
break;
default:
/* Update the factored H to which a new column has been added. */
n_minus_1 = n - 1;
code = 0;
/* Multiply the new column by the previous n-1 Givens rotations. */
for (k=0; k < n_minus_1; k++) {
i = 2*k;
temp1 = h[k][n_minus_1];
temp2 = h[k+1][n_minus_1];
c = q[i];
s = q[i+1];
h[k][n_minus_1] = c*temp1 - s*temp2;
h[k+1][n_minus_1] = s*temp1 + c*temp2;
}
/* Compute new Givens rotation and multiply it times the last two
entries in the new column of H. Note that the second entry of
this product will be 0, so it is not necessary to compute it. */
temp1 = h[n_minus_1][n_minus_1];
temp2 = h[n][n_minus_1];
if (temp2 == ZERO) {
c = ONE;
s = ZERO;
} else if (ABS(temp2) >= ABS(temp1)) {
temp3 = temp1/temp2;
s = -ONE/RSqrt(ONE+SQR(temp3));
c = -s*temp3;
} else {
temp3 = temp2/temp1;
c = ONE/RSqrt(ONE+SQR(temp3));
s = -c*temp3;
}
q_ptr = 2*n_minus_1;
q[q_ptr] = c;
q[q_ptr+1] = s;
if ((h[n_minus_1][n_minus_1] = c*temp1 - s*temp2) == ZERO)
code = n;
}
return (code);
}
/*************** QRsol ************************************************
This implementation of QRsol is a slight modification of a previous
routine (called qrsol) written by Milo Dorr.
**********************************************************************/
int QRsol(int n, real **h, real *q, real *b)
{
real c, s, temp1, temp2;
int i, k, q_ptr, code=0;
/* Compute Q*b. */
for (k=0; k < n; k++) {
q_ptr = 2*k;
c = q[q_ptr];
s = q[q_ptr+1];
temp1 = b[k];
temp2 = b[k+1];
b[k] = c*temp1 - s*temp2;
b[k+1] = s*temp1 + c*temp2;
}
/* Solve R*x = Q*b. */
for (k=n-1; k >= 0; k--) {
if (h[k][k] == ZERO) {
code = k + 1;
break;
}
b[k] /= h[k][k];
for (i=0; i < k; i++) b[i] -= b[k]*h[i][k];
}
return (code);
}

67
ext/cvode/source/llnlmath.c
View file

@ -1,67 +0,0 @@
/******************************************************************
* *
* File : llnlmath.c *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 1 September 1994 *
*----------------------------------------------------------------*
* This is the implementation file for a C math library. *
* *
******************************************************************/
#include <stdio.h>
#include <math.h>
#include "llnlmath.h"
#include "llnltyps.h"
#define ZERO RCONST(0.0)
#define ONE RCONST(1.0)
#define TWO RCONST(2.0)
real UnitRoundoff(void)
{
real u;
volatile real one_plus_u;
u = ONE;
one_plus_u = ONE + u;
while (one_plus_u != ONE) {
u /= TWO;
one_plus_u = ONE + u;
}
u *= TWO;
return(u);
}
real RPowerI(real base, int exponent)
{
int i, expt;
real prod;
prod = ONE;
expt = ABS(exponent);
for(i=1; i <= expt; i++) prod *= base;
if (exponent < 0) prod = ONE/prod;
return(prod);
}
real RPowerR(real base, real exponent)
{
if (base <= ZERO) return(ZERO);
return((real)pow((double)base,(double)exponent));
}
real RSqrt(real x)
{
if (x <= ZERO) return(ZERO);
return((real) sqrt((double) x));
}

672
ext/cvode/source/nvector.c
View file

@ -1,672 +0,0 @@
/****************************************************************
* *
* File : nvector.c *
* Programmers : Scott D. Cohen, Alan C. Hindmarsh, and *
* : Allan G. Taylor, LLNL *
* Version of : 17 December 1999 *
*--------------------------------------------------------------*
* *
* This is the implementation file for a generic serial NVECTOR *
* package. It contains the implementation of the N_Vector *
* kernels listed in nvector.h. *
* *
****************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "nvector.h"
#include "llnltyps.h"
#include "llnlmath.h"
#define ZERO RCONST(0.0)
#define HALF RCONST(0.5)
#define ONE RCONST(1.0)
#define ONEPT5 RCONST(1.5)
/* Private Helper Prototypes */
static void VCopy(N_Vector x, N_Vector z); /* z=x */
static void VSum(N_Vector x, N_Vector y, N_Vector z); /* z=x+y */
static void VDiff(N_Vector x, N_Vector y, N_Vector z); /* z=x-y */
static void VNeg(N_Vector x, N_Vector z); /* z=-x */
/* z=c(x+y) */
static void VScaleSum(real c, N_Vector x, N_Vector y, N_Vector z);
/* z=c(x-y) */
static void VScaleDiff(real c, N_Vector x, N_Vector y, N_Vector z);
static void VLin1(real a, N_Vector x, N_Vector y, N_Vector z); /* z=ax+y */
static void VLin2(real a, N_Vector x, N_Vector y, N_Vector z); /* z=ax-y */
static void Vaxpy(real a, N_Vector x, N_Vector y); /* y <- ax+y */
static void VScaleBy(real a, N_Vector x); /* x <- ax */
/********************* Exported Functions ************************/
N_Vector N_VNew(integer N, void *machEnv)
{
N_Vector v;
if (N <= 0) return(NULL);
v = (N_Vector) malloc(sizeof *v);
if (v == NULL) return(NULL);
v->data = (real *) malloc(N * sizeof(real));
if (v->data == NULL) {
free(v);
return(NULL);
}
v->length = N;
return(v);
}
void N_VFree(N_Vector x)
{
free(x->data);
free(x);
}
void N_VLinearSum(real a, N_Vector x, real b, N_Vector y, N_Vector z)
{
integer i, N;
real c, *xd, *yd, *zd;
N_Vector v1, v2;
boole test;
if ((b == ONE) && (z == y)) { /* BLAS usage: axpy y <- ax+y */
Vaxpy(a,x,y);
return;
}
if ((a == ONE) && (z == x)) { /* BLAS usage: axpy x <- by+x */
Vaxpy(b,y,x);
return;
}
/* Case: a == b == 1.0 */
if ((a == ONE) && (b == ONE)) {
VSum(x, y, z);
return;
}
/* Cases: (1) a == 1.0, b = -1.0, (2) a == -1.0, b == 1.0 */
if ((test = ((a == ONE) && (b == -ONE))) || ((a == -ONE) && (b == ONE))) {
v1 = test ? y : x;
v2 = test ? x : y;
VDiff(v2, v1, z);
return;
}
/* Cases: (1) a == 1.0, b == other or 0.0, (2) a == other or 0.0, b == 1.0 */
/* if a or b is 0.0, then user should have called N_VScale */
if ((test = (a == ONE)) || (b == ONE)) {
c = test ? b : a;
v1 = test ? y : x;
v2 = test ? x : y;
VLin1(c, v1, v2, z);
return;
}
/* Cases: (1) a == -1.0, b != 1.0, (2) a != 1.0, b == -1.0 */
if ((test = (a == -ONE)) || (b == -ONE)) {
c = test ? b : a;
v1 = test ? y : x;
v2 = test ? x : y;
VLin2(c, v1, v2, z);
return;
}
/* Case: a == b */
/* catches case both a and b are 0.0 - user should have called N_VConst */
if (a == b) {
VScaleSum(a, x, y, z);
return;
}
/* Case: a == -b */
if (a == -b) {
VScaleDiff(a, x, y, z);
return;
}
/* Do all cases not handled above:
(1) a == other, b == 0.0 - user should have called N_VScale
(2) a == 0.0, b == other - user should have called N_VScale
(3) a,b == other, a !=b, a != -b */
N = x->length;
xd = x->data;
yd = y->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = a * (*xd++) + b * (*yd++);
}
void N_VConst(real c, N_Vector z)
{
integer i, N;
real *zd;
N = z->length;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = c;
}
void N_VProd(N_Vector x, N_Vector y, N_Vector z)
{
integer i, N;
real *xd, *yd, *zd;
N = x->length;
xd = x->data;
yd = y->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = (*xd++) * (*yd++);
}
void N_VDiv(N_Vector x, N_Vector y, N_Vector z)
{
integer i, N;
real *xd, *yd, *zd;
N = x->length;
xd = x->data;
yd = y->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = (*xd++) / (*yd++);
}
void N_VScale(real c, N_Vector x, N_Vector z)
{
integer i, N;
real *xd, *zd;
if (z == x) { /* BLAS usage: scale x <- cx */
VScaleBy(c, x);
return;
}
if (c == ONE) {
VCopy(x, z);
} else if (c == -ONE) {
VNeg(x, z);
} else {
N = x->length;
xd = x->data;
zd = z->data;
for (i=0; i < N; i++) *zd++ = c * (*xd++);
}
}
void N_VAbs(N_Vector x, N_Vector z)
{
integer i, N;
real *xd, *zd;
N = x->length;
xd = x->data;
zd = z->data;
for (i=0; i < N; i++, xd++, zd++)
*zd = ABS(*xd);
}
void N_VInv(N_Vector x, N_Vector z)
{
integer i, N;
real *xd, *zd;
N = x->length;
xd = x->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = ONE / (*xd++);
}
void N_VAddConst(N_Vector x, real b, N_Vector z)
{
integer i, N;
real *xd, *zd;
N = x->length;
xd = x->data;
zd = z->data;
for (i=0; i < N; i++) *zd++ = (*xd++) + b;
}
real N_VDotProd(N_Vector x, N_Vector y)
{
integer i, N;
real sum = ZERO, *xd, *yd;
N = x->length;
xd = x->data;
yd = y->data;
for (i=0; i < N; i++)
sum += (*xd++) * (*yd++);
return(sum);
}
real N_VMaxNorm(N_Vector x)
{
integer i, N;
real max = ZERO, *xd;
N = x->length;
xd = x->data;
for (i=0; i < N; i++, xd++) {
if (ABS(*xd) > max) max = ABS(*xd);
}
return(max);
}
real N_VWrmsNorm(N_Vector x, N_Vector w)
{
integer i, N;
real sum = ZERO, prodi, *xd, *wd;
N = x->length;
xd = x->data;
wd = w->data;
for (i=0; i < N; i++) {
prodi = (*xd++) * (*wd++);
sum += prodi * prodi;
}
return(RSqrt(sum / N));
}
real N_VMin(N_Vector x)
{
integer i, N;
real min, *xd;
N = x->length;
xd = x->data;
min = xd[0];
for (i=1; i < N; i++, xd++) {
if ((*xd) < min) min = *xd;
}
return(min);
}
real N_VWL2Norm(N_Vector x, N_Vector w)
{
integer i, N;
real sum = ZERO, prodi, *xd, *wd;
N = x->length;
xd = x->data;
wd = w->data;
for (i=0; i < N; i++) {
prodi = (*xd++) * (*wd++);
sum += prodi * prodi;
}
return(RSqrt(sum));
}
real N_VL1Norm(N_Vector x)
{
integer i, N;
real sum = ZERO, *xd;
N = x->length;
xd = x->data;
for (i=0; i<N; i++)
sum += ABS(xd[i]);
return(sum);
}
void N_VOneMask(N_Vector x)
{
integer i, N;
real *xd;
N = x->length;
xd = x->data;
for (i=0; i<N; i++,xd++) {
if (*xd != ZERO) *xd = ONE;
}
}
void N_VCompare(real c, N_Vector x, N_Vector z)
{
integer i, N;
real *xd, *zd;
N = x->length;
xd = x->data;
zd = z->data;
for (i=0; i < N; i++, xd++, zd++) {
*zd = (ABS(*xd) >= c) ? ONE : ZERO;
}
}
boole N_VInvTest(N_Vector x, N_Vector z)
{
integer i, N;
real *xd, *zd;
N = x->length;
xd = x->data;
zd = z->data;
for (i=0; i < N; i++) {
if (*xd == ZERO) return(FALSE);
*zd++ = ONE / (*xd++);
}
return(TRUE);
}
boole N_VConstrProdPos(N_Vector c, N_Vector x)
{
integer i, N;
real *xd, *cd;
boole test;
N = x->length;
xd = x->data;
cd = c->data;
test = TRUE;
for (i=0; i < N; i++, xd++,cd++) {
if (*cd != ZERO) {
if ((*xd)*(*cd) <= ZERO) {
test = FALSE;
break;
}
}
}
return(test);
}
boole N_VConstrMask(N_Vector c, N_Vector x, N_Vector m)
{
integer i, N;
boole test;
real *cd, *xd, *md;
N = x->length;
cd = c->data;
xd = x->data;
md = m->data;
test = TRUE;
for (i=0; i<N; i++, cd++, xd++, md++) {
if ( *cd == ZERO) *md = ZERO;
else {
if ( *cd > ONEPT5 || (*cd) < -ONEPT5) {
if ( (*xd)*(*cd) <= ZERO) {
test = FALSE;
*md = ONE;
}
else {
*md = ZERO;
}
} else if ( (*cd) > HALF || (*cd) < -HALF) {
if ( (*xd)*(*cd) < ZERO ) {
test = FALSE;
*md = ONE;
} else {
*md = ZERO;
}
}
}
}
return(test);
}
real N_VMinQuotient(N_Vector num, N_Vector denom)
{
boole notEvenOnce;
integer i, N;
real *nd, *dd, min;
N = num->length;
nd = num->data;
dd = denom->data;
notEvenOnce = TRUE;
for (i=0; i<N; i++, nd++, dd++) {
if (*dd == ZERO) continue;
else {
if (notEvenOnce) {
min = *nd / *dd ;
notEvenOnce = FALSE;
}
else
min = MIN(min, (*nd)/(*dd));
}
}
if (notEvenOnce) min = 1.e99;
return(min);
}
void N_VPrint(N_Vector x)
{
integer i, N;
real *xd;
N = x->length;
xd = x->data;
for (i=0; i < N; i++) printf("%11.8g\n", *xd++);
printf("\n");
}
/***************** Private Helper Functions **********************/
static void VCopy(N_Vector x, N_Vector z)
{
integer i, N;
real *xd, *zd;
N = x->length;
xd = x->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = *xd++;
}
static void VSum(N_Vector x, N_Vector y, N_Vector z)
{
integer i, N;
real *xd, *yd, *zd;
N = x->length;
xd = x->data;
yd = y->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = (*xd++) + (*yd++);
}
static void VDiff(N_Vector x, N_Vector y, N_Vector z)
{
integer i, N;
real *xd, *yd, *zd;
N = x->length;
xd = x->data;
yd = y->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = (*xd++) - (*yd++);
}
static void VNeg(N_Vector x, N_Vector z)
{
integer i, N;
real *xd, *zd;
N = x->length;
xd = x->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = -(*xd++);
}
static void VScaleSum(real c, N_Vector x, N_Vector y, N_Vector z)
{
integer i, N;
real *xd, *yd, *zd;
N = x->length;
xd = x->data;
yd = y->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = c * ((*xd++) + (*yd++));
}
static void VScaleDiff(real c, N_Vector x, N_Vector y, N_Vector z)
{
integer i, N;
real *xd, *yd, *zd;
N = x->length;
xd = x->data;
yd = y->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = c * ((*xd++) - (*yd++));
}
static void VLin1(real a, N_Vector x, N_Vector y, N_Vector z)
{
integer i, N;
real *xd, *yd, *zd;
N = x->length;
xd = x->data;
yd = y->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = a * (*xd++) + (*yd++);
}
static void VLin2(real a, N_Vector x, N_Vector y, N_Vector z)
{
integer i, N;
real *xd, *yd, *zd;
N = x->length;
xd = x->data;
yd = y->data;
zd = z->data;
for (i=0; i < N; i++)
*zd++ = a * (*xd++) - (*yd++);
}
static void Vaxpy(real a, N_Vector x, N_Vector y)
{
integer i, N;
real *xd, *yd;
N = x->length;
xd = x->data;
yd = y->data;
if (a == ONE) {
for (i=0; i < N; i++)
*yd++ += (*xd++);
return;
}
if (a == -ONE) {
for (i=0; i < N; i++)
*yd++ -= (*xd++);
return;
}
for (i=0; i < N; i++)
*yd++ += a * (*xd++);
}
static void VScaleBy(real a, N_Vector x)
{
integer i, N;
real *xd;
N = x->length;
xd = x->data;
for (i=0; i < N; i++)
*xd++ *= a;
}

429
ext/cvode/source/spgmr.c
View file

@ -1,429 +0,0 @@
/******************************************************************
* File : spgmr.c *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 17 December 1999 *
*----------------------------------------------------------------*
* This is the implementation file for the scaled preconditioned *
* GMRES (SPGMR) iterative linear solver. *
* *
******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "iterativ.h"
#include "spgmr.h"
#include "llnltyps.h"
#include "nvector.h"
#include "llnlmath.h"
#define ZERO RCONST(0.0)
#define ONE RCONST(1.0)
/*************** Private Helper Function Prototype *******************/
static void FreeVectorArray(N_Vector *A, int indMax);
/* Implementation of SPGMR algorithm */
/*************** SpgmrMalloc *****************************************/
SpgmrMem SpgmrMalloc(integer N, int l_max, void *machEnv)
{
SpgmrMem mem;
N_Vector *V, xcor, vtemp;
real **Hes, *givens, *yg;
int k, i;
/* Check the input parameters. */
if ((N <= 0) || (l_max <= 0)) return(NULL);
/* Get memory for the Krylov basis vectors V[0], ..., V[l_max]. */
V = (N_Vector *) malloc((l_max+1)*sizeof(N_Vector));
if (V == NULL) return(NULL);
for (k = 0; k <= l_max; k++) {
V[k] = N_VNew(N, machEnv);
if (V[k] == NULL) {
FreeVectorArray(V, k-1);
return(NULL);
}
}
/* Get memory for the Hessenberg matrix Hes. */
Hes = (real **) malloc((l_max+1)*sizeof(real *));
if (Hes == NULL) {
FreeVectorArray(V, l_max);
return(NULL);
}
for (k = 0; k <= l_max; k++) {
Hes[k] = (real *) malloc(l_max*sizeof(real));
if (Hes[k] == NULL) {
for (i = 0; i < k; i++) free(Hes[i]);
FreeVectorArray(V, l_max);
return(NULL);
}
}
/* Get memory for Givens rotation components. */
givens = (real *) malloc(2*l_max*sizeof(real));
if (givens == NULL) {
for (i = 0; i <= l_max; i++) free(Hes[i]);
FreeVectorArray(V, l_max);
return(NULL);
}
/* Get memory to hold the correction to z_tilde. */
xcor = N_VNew(N, machEnv);
if (xcor == NULL) {
free(givens);
for (i = 0; i <= l_max; i++) free(Hes[i]);
FreeVectorArray(V, l_max);
return(NULL);
}
/* Get memory to hold SPGMR y and g vectors. */
yg = (real *) malloc((l_max+1)*sizeof(real));
if (yg == NULL) {
N_VFree(xcor);
free(givens);
for (i = 0; i <= l_max; i++) free(Hes[i]);
FreeVectorArray(V, l_max);
return(NULL);
}
/* Get an array to hold a temporary vector. */
vtemp = N_VNew(N, machEnv);
if (vtemp == NULL) {
free(yg);
N_VFree(xcor);
free(givens);
for (i = 0; i <= l_max; i++) free(Hes[i]);
FreeVectorArray(V, l_max);
return(NULL);
}
/* Get memory for an SpgmrMemRec containing SPGMR matrices and vectors. */
mem = (SpgmrMem) malloc(sizeof(SpgmrMemRec));
if (mem == NULL) {
N_VFree(vtemp);
free(yg);
N_VFree(xcor);
free(givens);
for (i = 0; i <= l_max; i++) free(Hes[i]);
FreeVectorArray(V, l_max);
return(NULL);
}
/* Set the fields of mem. */
mem->N = N;
mem->l_max = l_max;
mem->V = V;
mem->Hes = Hes;
mem->givens = givens;
mem->xcor = xcor;
mem->yg = yg;
mem->vtemp = vtemp;
/* Return the pointer to SPGMR memory. */
return(mem);
}
/*************** SpgmrSolve ******************************************/
int SpgmrSolve(SpgmrMem mem, void *A_data, N_Vector x, N_Vector b,
int pretype, int gstype, real delta, int max_restarts,
void *P_data, N_Vector s1, N_Vector s2, ATimesFn atimes,
PSolveFn psolve, real *res_norm, int *nli, int *nps)
{
N_Vector *V, xcor, vtemp;
real **Hes, *givens, *yg;
real beta, rotation_product, r_norm, s_product, rho;
boole preOnLeft, preOnRight, scale2, scale1, converged;
int i, j, k, l, l_plus_1, l_max, krydim, ier, ntries;
if (mem == NULL) return(SPGMR_MEM_NULL);
/* Make local copies of mem variables. */
l_max = mem->l_max;
V = mem->V;
Hes = mem->Hes;
givens = mem->givens;
xcor = mem->xcor;
yg = mem->yg;
vtemp = mem->vtemp;
*nli = *nps = 0; /* Initialize counters */
converged = FALSE; /* Initialize converged flag */
if (max_restarts < 0) max_restarts = 0;
if ((pretype != LEFT) && (pretype != RIGHT) && (pretype != BOTH))
pretype = NONE;
preOnLeft = ((pretype == LEFT) || (pretype == BOTH));
preOnRight = ((pretype == RIGHT) || (pretype == BOTH));
scale1 = (s1 != NULL);
scale2 = (s2 != NULL);
/* Set vtemp and V[0] to initial (unscaled) residual r_0 = b - A*x_0. */
if (N_VDotProd(x, x) == ZERO) {
N_VScale(ONE, b, vtemp);
} else {
if (atimes(A_data, x, vtemp) != 0)
return(SPGMR_ATIMES_FAIL);
N_VLinearSum(ONE, b, -ONE, vtemp, vtemp);
}
N_VScale(ONE, vtemp, V[0]);
/* Apply left preconditioner and left scaling to V[0] = r_0. */
if (preOnLeft) {
ier = psolve(P_data, V[0], vtemp, LEFT);
(*nps)++;
if (ier != 0)
return((ier < 0) ? SPGMR_PSOLVE_FAIL_UNREC : SPGMR_PSOLVE_FAIL_REC);
} else {
N_VScale(ONE, V[0], vtemp);
}
if (scale1) {
N_VProd(s1, vtemp, V[0]);
} else {
N_VScale(ONE, vtemp, V[0]);
}
/* Set r_norm = beta to L2 norm of V[0] = s1 P1_inv r_0, and
return if small. */
*res_norm = r_norm = beta = RSqrt(N_VDotProd(V[0], V[0]));
if (r_norm <= delta)
return(SPGMR_SUCCESS);
/* Set xcor = 0. */
N_VConst(ZERO, xcor);
/* Begin outer iterations: up to (max_restarts + 1) attempts. */
for (ntries = 0; ntries <= max_restarts; ntries++) {
/* Initialize the Hessenberg matrix Hes and Givens rotation
product. Normalize the initial vector V[0]. */
for (i = 0; i <= l_max; i++)
for (j = 0; j < l_max; j++)
Hes[i][j] = ZERO;
rotation_product = ONE;
N_VScale(ONE/r_norm, V[0], V[0]);
/* Inner loop: generate Krylov sequence and Arnoldi basis. */
for (l = 0; l < l_max; l++) {
(*nli)++;
krydim = l_plus_1 = l + 1;
/* Generate A-tilde V[l], where A-tilde = s1 P1_inv A P2_inv s2_inv. */
/* Apply right scaling: vtemp = s2_inv V[l]. */
if (scale2) N_VDiv(V[l], s2, vtemp);
else N_VScale(ONE, V[l], vtemp);
/* Apply right preconditioner: vtemp = P2_inv s2_inv V[l]. */
if (preOnRight) {
N_VScale(ONE, vtemp, V[l_plus_1]);
ier = psolve(P_data, V[l_plus_1], vtemp, RIGHT);
(*nps)++;
if (ier != 0)
return((ier < 0) ? SPGMR_PSOLVE_FAIL_UNREC : SPGMR_PSOLVE_FAIL_REC);
}
/* Apply A: V[l+1] = A P2_inv s2_inv V[l]. */
if (atimes(A_data, vtemp, V[l_plus_1] ) != 0)
return(SPGMR_ATIMES_FAIL);
/* Apply left preconditioning: vtemp = P1_inv A P2_inv s2_inv V[l]. */
if (preOnLeft) {
ier = psolve(P_data, V[l_plus_1], vtemp, LEFT);
(*nps)++;
if (ier != 0)
return((ier < 0) ? SPGMR_PSOLVE_FAIL_UNREC : SPGMR_PSOLVE_FAIL_REC);
} else {
N_VScale(ONE, V[l_plus_1], vtemp);
}
/* Apply left scaling: V[l+1] = s1 P1_inv A P2_inv s2_inv V[l]. */
if (scale1) {
N_VProd(s1, vtemp, V[l_plus_1]);
} else {
N_VScale(ONE, vtemp, V[l_plus_1]);
}
/* Orthogonalize V[l+1] against previous V[i]: V[l+1] = w_tilde. */
if (gstype == CLASSICAL_GS) {
if (ClassicalGS(V, Hes, l_plus_1, l_max, &(Hes[l_plus_1][l]),
vtemp, yg) != 0)
return(SPGMR_GS_FAIL);
} else {
if (ModifiedGS(V, Hes, l_plus_1, l_max, &(Hes[l_plus_1][l])) != 0)
return(SPGMR_GS_FAIL);
}
/* Update the QR factorization of Hes. */
if(QRfact(krydim, Hes, givens, l) != 0 )
return(SPGMR_QRFACT_FAIL);
/* Update residual norm estimate; break if convergence test passes. */
rotation_product *= givens[2*l+1];
*res_norm = rho = ABS(rotation_product*r_norm);
if (rho <= delta) { converged = TRUE; break; }
/* Normalize V[l+1] with norm value from the Gram-Schmidt routine. */
N_VScale(ONE/Hes[l_plus_1][l], V[l_plus_1], V[l_plus_1]);
}
/* Inner loop is done. Compute the new correction vector xcor. */
/* Construct g, then solve for y. */
yg[0] = r_norm;
for (i = 1; i <= krydim; i++) yg[i]=ZERO;
if (QRsol(krydim, Hes, givens, yg) != 0)
return(SPGMR_QRSOL_FAIL);
/* Add correction vector V_l y to xcor. */
for (k = 0; k < krydim; k++)
N_VLinearSum(yg[k], V[k], ONE, xcor, xcor);
/* If converged, construct the final solution vector x and return. */
if (converged) {
/* Apply right scaling and right precond.: vtemp = P2_inv s2_inv xcor. */
if (scale2) N_VDiv(xcor, s2, xcor);
if (preOnRight) {
ier = psolve(P_data, xcor, vtemp, RIGHT);
(*nps)++;
if (ier != 0)
return((ier < 0) ? SPGMR_PSOLVE_FAIL_UNREC : SPGMR_PSOLVE_FAIL_REC);
} else {
N_VScale(ONE, xcor, vtemp);
}
/* Add vtemp to initial x to get final solution x, and return */
N_VLinearSum(ONE, x, ONE, vtemp, x);
return(SPGMR_SUCCESS);
}
/* Not yet converged; if allowed, prepare for restart. */
if (ntries == max_restarts) break;
/* Construct last column of Q in yg. */
s_product = ONE;
for (i = krydim; i > 0; i--) {
yg[i] = s_product*givens[2*i-2];
s_product *= givens[2*i-1];
}
yg[0] = s_product;
/* Scale r_norm and yg. */
r_norm *= s_product;
for (i = 0; i <= krydim; i++)
yg[i] *= r_norm;
r_norm = ABS(r_norm);
/* Multiply yg by V_(krydim+1) to get last residual vector; restart. */
N_VScale(yg[0], V[0], V[0]);
for (k = 1; k <= krydim; k++)
N_VLinearSum(yg[k], V[k], ONE, V[0], V[0]);
}
/* Failed to converge, even after allowed restarts.
If the residual norm was reduced below its initial value, compute
and return x anyway. Otherwise return failure flag. */
if (rho < beta) {
/* Apply right scaling and right precond.: vtemp = P2_inv s2_inv xcor. */
if (scale2) N_VDiv(xcor, s2, xcor);
if (preOnRight) {
ier = psolve(P_data, xcor, vtemp, RIGHT);
(*nps)++;
if (ier != 0)
return((ier < 0) ? SPGMR_PSOLVE_FAIL_UNREC : SPGMR_PSOLVE_FAIL_REC);
} else {
N_VScale(ONE, xcor, vtemp);
}
/* Add vtemp to initial x to get final solution x, and return. */
N_VLinearSum(ONE, x, ONE, vtemp, x);
return(SPGMR_RES_REDUCED);
}
return(SPGMR_CONV_FAIL);
}
/*************** SpgmrFree *******************************************/
void SpgmrFree(SpgmrMem mem)
{
int i, l_max;
real **Hes;
if (mem == NULL) return;
l_max = mem->l_max;
Hes = mem->Hes;
FreeVectorArray(mem->V, l_max);
for (i = 0; i <= l_max; i++) free(Hes[i]);
free(Hes);
free(mem->givens);
N_VFree(mem->xcor);
free(mem->yg);
N_VFree(mem->vtemp);
free(mem);
}
/*************** Private Helper Function: FreeVectorArray ************/
static void FreeVectorArray(N_Vector *A, int indMax)
{
int j;
for (j = 0; j <= indMax; j++) N_VFree(A[j]);
free(A);
}