Procedures

This routine computes the spatial average value of $c_1$.

This routine computes and loads the vectors of initial values.

This routine computes the matrix-vector product:

This routine checks for constraint violations in the initial approximate solution y.

This subroutine checks for constraint violations in the proposed new approximate solution ynew. If a constraint violation occurs, then a new step length, tau, is calculated, and this value is given to the linesearch routine to calculate a new approximate solution ynew.

This routine solves a nonlinear system of algebraic equations of the form:

This routine solves a nonlinear system of algebraic equations of the form:

The methods in subroutine ddstp use polynomials to approximate the solution. This routine approximates the solution and its derivative at tout by evaluating one of these polynomials, and its derivative, there. Information defining the polynomial is passed from the caller.

This subroutine sets the error weight vector, wt, according to:

This routine solves a system of differential/algebraic equations of the form:

This function computes the weighted root-mean-square norm of the vector v with reciprocal weights rwt:

This routine calculates the scaled preconditioned norm of the nonlinear function used in the nonlinear iteration for obtaining consistent initial conditions. Specifically, it calculates the weighted root-mean-square norm of the vector:

This routine calculates the scaled preconditioned norm of the nonlinear function used in the nonlinear iteration for obtaining consistent initial conditions. Specifically, it calculates the weighted root-mean-square norm of the vector:

DGBFA factors a real band matrix by elimination.

DGBSL solves a real banded system factored by DGBCO or DGBFA.

DGEFA factors a real general matrix.

DGESL solves a real general linear system A * X = B.

This routine solves the least squares problem:

This routine performs a QR decomposition of an upper Hessenberg matrix a using Givens rotations. There are two options available:

This subroutine checks the error weight vector wt for components that are less or equal to 0, and if none are found, it inverts wt in place. This replaces division operations with multiplications in all norm evaluations.

This routine uses a linesearch algorithm to calculate a new $(y,\dot{y})$ pair, denoted $(y_{new},\dot{y}_{new})$, such that:

This routine uses a linesearch algorithm to calculate a new $(y,\dot{y})$ pair, denoted $(y_{new},\dot{y}_{new})$, such that:

This routine computes the iteration matrix:

This routine solves a nonlinear system of algebraic equations of the form:

This routine solves a nonlinear system of algebraic equations of the form:

This routine solves a nonlinear system of algebraic equations of the form:

This routine solves a nonlinear system of algebraic equations of the form:

This routine solves a nonlinear system of algebraic equations of the form:

This routines solves a nonlinear system of algebraic equations of the form:

This routine orthogonalizes the vector vnew against the previous kmp vectors in the v matrix. It uses a modified Gram-Schmidt orthogonalization procedure with conditional reorthogonalization.

This routine manages the solution of the linear system arising in the Newton iteration. Real matrix information and real temporary storage is stored in the array rwm, while integer matrix information is stored in the array iwm. For a dense matrix, the LINPACK routine dgesl is called. For a banded matrix, the LINPACK routine dgbsl is called.

This routine uses a restart algorithm and interfaces to dspigm for the solution of the linear system arising from a Newton iteration.

This routine solves the linear system $A z = r$ using a scaled preconditioned version of the generalized minimum residual method. An initial guess of $z=0$ is assumed.

Jacobian routine.

This routine interfaces to subroutines jac_rbdpre or jac_rbgpre, depending on the flag jbg = ipar(2), to generate and preprocess the block-diagonal Jacobian corresponding to the reaction term.

This routine generates a banded preconditioner matrix $P$ that approximates the Jacobian matrix $J$. The banded matrix $P$ has half-bandwidths $\mathrm{ml}$ and $\mathrm{mu}$, and is computed by making $\mathrm{ml} + \mathrm{mu} + 1$ calls to the user's res routine and forming difference quotients, exactly as in the banded direct method option of daskr. Afterwards, this matrix is LU factorized by dgbfa from LINPACK and the factors are stored in the work arrays rwp and iwp.

This subroutine uses finite-differences to calculate the Jacobian matrix in sparse format, and then performs an incomplete LU decomposition using either ilut or ilutp from SPARSKIT.

This routine generates and preprocesses a block-diagonal preconditioner matrix, based on the part of the Jacobian corresponding to the reaction terms $R$ of the problem. It generates a matrix of the form $c_J I_d - \partial R/ \partial u$. It calls dgefa from LINPACK to do the LU decomposition of each diagonal block. The computation of the diagonal blocks uses the mesh information in the module variables. One block per spatial point is computed. The Jacobian elements are generated by difference quotients. This routine calls a user routine rblock to compute the block (jx,jy) of $R$.

This routine generates and preprocesses a block-diagonal preconditioner matrix, based on the part of the Jacobian corresponding to the reaction terms $R$ of the problem, using block-grouping. It generates a matrix of the form $c_J I_d - \partial R/ \partial u$. It calls dgefa from LINPACK to do the LU decomposition of each diagonal block. The computation of the diagonal blocks uses the mesh information in the module variables. One block per group is computed. The Jacobian elements are generated by difference quotients. This routine calls a user routine rblock to compute the block (jx,jy) of $R$.

This routine prints the values of the individual species densities at the current time t, to logical unit lout.

This routine applies the inverse of a product preconditioner matrix to the vector in the array b. Depending on the flag jpre, this involves a call to gauss_seidel, for the inverse of the spatial factor, and/or a call to psol_rbdpre or psol_rbgpre for the inverse of the reaction-based factor $c_J I_d - \partial R / \partial y$. If jbg == 0, the latter factor does not use block-grouping and jac_rbdpre is called. Otherwise, if jbg == 1, block-grouping is used and jac_rbgpre is called. The array b is overwritten with the solution.

This routine solves the linear system $P x = b$ for the banded preconditioner $P$, given a vector $b$, using the LU decomposition produced by jac_banpre. The solution is carried out by dgbsl from LINPACK.

This subroutine solves the linear system $P x = b$ for the sparse preconditioner $P$, given a vector $b$, using the incomplete LU decomposition produced by jac_ilupre. The solution is carried out by lusol from SPARSKIT.

This routine solves a linear system $A_R x = b$, using the LU factors of the diagonal blocks computed in jac_rbdpre and mesh parameters passed as module variables. The solution is carried out by dgesl from LINPACK.

This routine solves a linear system $A_R x = b$, using the LU factors of the diagonal blocks computed in jac_rbgpre and mesh parameters passed as module variables. The solution is carried out by dgesl from LINPACK.

User-defined residuals routine. It computes the residuals for the 2D discretized heat equation, with zero boundary values.

Residuals routine.

Residuals routine.

User-defined residuals routine. It computes the residuals for the 2D discretized heat equation, with zero boundary values.

This routine computes the residuals vector, using subroutine f for the right-hand sides.

User-defined roots routine.

Roots routine.

Roots routine.

User-defined roots routine.

Roots routine.

This routine sets the ID array in iwork, indicating which components are differential and which are algebraic.

This routine sets the basic problem parameters which are passed to the various routines via the module web_par.

Setup routine for the incomplete LU preconditioner. This routine checks the user input and calculates the minimum length needed for the preconditioner workspace arrays.

This routine sets the mesh parameters needed to use the routines jac_rbdpre and psol_rbdpre, assuming a 2D rectangular problem. Additionally, it loads the lengths lrwp and liwp in array iwork.

This routine sets the mesh and block-grouping parameters needed to use the routines jac_rbgpre and psol_rbgpre, assuming a 2D rectangular problem with uniform rectangular block-grouping. Additionally, it loads the lengths lrwp and liwp in array iwork.

This routine computes and loads the vector of initial values. The initial $u$ values are given by the polynomial:

This routine computes and loads the vector of initial values. The initial $u$ values are given by the polynomial: