The FKINSOL interface module is a package of C functions which support the use of the KINSOL solver, for the solution nonlinear systems F(u) = 0, in a mixed FORTRAN/C setting. While KINSOL is written in C, it is assumed here that the user's calling program and user-supplied problem-defining routines are written in FORTRAN. This package provides the necessary interface to KINSOL for both the serial and the parallel NVECTOR implementations.
In this package, the names of the interface functions, and the names of the FORTRAN user routines called by them, appear as dummy names which are mapped to actual values by a series of definitions in the header files fkinsol.hand fkinbbd.h. By default, those mapping definitions depend in turn on the C macro F77_FUNC defined in the header file sundials_config.h by configure. However, the set of flags SUNDIALS_CASE_UPPER, SUNDIALS_CASE_LOWER, SUNDIALS_UNDERSCORE_NONE, SUNDIALS_UNDERSCORE_ONE, and SUNDIALS_UNDERSCORE_TWO can be explicitly defined in the header file sundials_config.h when configuring SUNDIALS via the options -with-f77underscore and -with-f77case to override the default behavior if necessary (see Chapter 2). Either way, the names into which the dummy names are mapped are in upper or lower case and have up to two underscores appended.
The user must also ensure that variables in the user FORTRAN code are declared in a manner consistent with their counterparts in KINSOL. All real variables must be declared as REAL, DOUBLE PRECISION, or perhaps as REAL*n, where n denotes the number of bytes, depending on whether KINSOL was built in single, double or extended precision (see Chapter 2). Moreover, some of the FORTRAN integer variables must be declared as INTEGER*4 or INTEGER*8 according to the C type long int. These integer variables include: the array of integer optional inputs and outputs (IOUT), problem dimensions (NEQ, NLOCAL, NGLOBAL), and Jacobian half-bandwidths (MU and ML). This is particularly important when using KINSOL and the FKINSOL package on 64-bit architectures.
The user-callable functions, with the corresponding KINSOL functions, are as follows:
The user-supplied functions, each listed with the corresponding interface function which calls it (and its type within KINSOL), are as follows:
FKINSOL routine (FORTRAN) | KINSOL function (C) | KINSOL function type |
FKFUN | FKINfunc | KINSysFn |
FKDJAC | FKINDenseJac | KINDenseJacFn |
FKBJAC | FKINBandJac | KINBandJacFn |
FKPSET | FKINPSet | KINSpilsPrecSetupFn |
FKPSOL | FKINPSol | KINSpilsPrecSolveFn |
FKJTIMES | FKINJtimes | KINSpilsJacTimesVecFn |
The usage of FKINSOL requires calls to a few different interface functions, depending on the method options selected, and one or more user-supplied routines which define the problem to be solved. These function calls and user routines are summarized separately below. Some details are omitted, and the user is referred to the description of the corresponding KINSOL functions for information on the arguments of any given user-callable interface routine, or of a given user-supplied function called by an interface function.
Steps marked with [S] in the instructions below apply to the serial NVECTOR implementation (NVECTOR-SERIAL) only, while those marked with [P] apply to NVECTOR-PARALLEL.
The user must in all cases supply the following FORTRAN routine
SUBROUTINE FKFUN (U, FVAL, IER) DIMENSION U(*), FVAL(*)It must set the FVAL array to F(u), the system function, as a function of U= u. IER is an error return flag which should be set to 0 if successful, a positive value if a recoverable error occurred (in which case KINSOL will attempt to corret), or a negative value if it failed unrecoverably (in which case the the solution process is halted).
[S] To initialize the serial NVECTOR module, the user must make the following call:
CALL FNVINITS (KEY, NEQ, IER)where KEY is the solver id (KEY = 3 for KINSOL), NEQ is the size of vectors, and IER is a return completion flag which is 0 on success and -1 if a failure occurred.
[P] To initialize the parallel vector module, the user must make the following call:
CALL FNVINITP (COMM, KEY, NLOCAL, NGLOBAL, IER)in which the arguments are: COMM = MPI communicator, KEY = 3, NLOCAL = the local size of vectors on this processor, and NGLOBAL = the system size (and the global size of all vectors, equal to the sum of all values of NLOCAL). The return flag IER is set to 0 on a successful return and to -1 otherwise.
If the header file sundials_config.h defines SUNDIALS_MPI_COMM_F2C to be 1 (meaning the MPI implementation used to build SUNDIALS includes the MPI_Comm_f2c function), then COMM can be any valid MPI communicator. Otherwise, MPI_COMM_WORLD will be used, so just pass an integer value as a placeholder.
To set various problem and solution parameters and allocate internal memory, make the following call:
CALL FKINMALLOC (IOUT, ROUT, IER) |
The optional outputs associated with the main KINSOL integrator are listed in Table 6.2.
The solution method in KINSOL involves the solution of linear systems related to the Jacobian of the nonlinear system. KINSOL presently includes five choices for the treatment of these systems, and the user of FKINSOL must call a routine with a specific name to make the desired choice.
[S] Dense treatment of the linear system
CALL FKINDENSE (NEQ, IER)where NEQ is the size of the nonlinear system. The argument IER is an error return flag which is 0 for success , -1 if a memory allocation failure occurred, or -2 for illegal input. As an option when using the DENSE linear solver, the user may supply a routine that computes a dense approximation of the system Jacobian J = F/u. If supplied, it must have the following form:
SUBROUTINE FKDJAC (NEQ, U, FVAL, DJAC, & WK1, WK2, IER) DIMENSION U(*), FVAL(*), DJAC(NEQ,*), & WK1(*), WK2(*)Typically this routine will use only NEQ, U, and DJAC. It must compute the Jacobian and store it columnwise in DJAC. The input arguments U and FVAL contain the current values of u, and F(u), respectively. The vectors WK1 and WK2 of length NEQ are provided as work space for use in FKDJAC. IER is an error return flag which should be set to 0 if successful, a positive value if a recoverable error occurred (in which case KINSOL will attempt to correct), or a negative value if FKDJAC failed unrecoverably (in which case the solution process is halted).
If the FKDJAC routine is provided, then, following the call to FKINDENSE, the user must make the call:
CALL FKINDENSESETJAC (FLAG, IER)with FLAG 0 to specify use of the user-supplied Jacobian approximation. The argument IER is an error return flag which is 0 for success or non-zero if an error occurred.
Optional outputs specific to the DENSE case are listed in Table 6.2.
[S] Band treatment of the linear system
CALL FKINBAND (NEQ, MU, ML, IER)The arguments are: MU, the upper half-bandwidth; ML, the lower half-bandwidth; and IER an error return flag which is 0 for success , -1 if a memory allocation failure occurred, or -2 in case an input has an illegal value.
As an option when using the BAND linear solver, the user may supply a routine that computes a band approximation of the system Jacobian J = F/u. If supplied, it must have the following form:
SUBROUTINE FKBJAC (NEQ, MU, ML, MDIM, U, FVAL, BJAC, & WK1, WK2, IER) DIMENSION U(*), FVAL(*), BJAC(MDIM,*), & WK1(*), WK2(*)Typically this routine will use only NEQ, MU, ML, U, and BJAC. It must load the MDIM by N array BJAC with the Jacobian matrix at the current u in band form. Store in BJAC(k, j) the Jacobian element Ji, j with k = i - j + MU +1 ( k = 1 ... ML + MU + 1) and j = 1 ... N. The input arguments U, and FVAL contain the current values of u, and F(u), respectively. The vectors WK1 and WK2 of length NEQ are provided as work space for use in FKBJAC. IER is an error return flag which should be set to 0 if successful, a positive value if a recoverable error occurred (in which case KINSOL will attempt to correct), or a negative value if FKBJAC failed unrecoverably (in which case the solution process is halted).
If the FKBJAC routine is provided, then, following the call to FKINBAND, the user must make the call:
CALL FKINBANDSETJAC (FLAG, IER)with FLAG 0 to specify use of the user-supplied Jacobian approximation. The argument IER is an error return flag which is 0 for success or non-zero if an error occurred.
Optional outputs specific to the BAND case are listed in Table 6.2.
[S][P] SPGMR treatment of the linear systems
For the Scaled Preconditioned GMRES solution of the linear systems, the user must make the call
CALL FKINSPGMR (MAXL, MAXLRST, IER)The arguments are as follows. MAXL is the maximum Krylov subspace dimension. MAXLRST is the maximum number of restarts. IER is an error return flag which is 0 to indicate success, -1 if a memory allocation failure occurred, or -2 to indicate an illegal input.
Optional outputs specific to the SPGMR case are listed in Table 6.2.
For descriptions of the relevant optional user-supplied routines, see User-supplied routines for SPGMR/SPBCG/SPTFQMR below.
[S][P] SPBCG treatment of the linear systems
For the Scaled Preconditioned Bi-CGStab solution of the linear systems, the user must make the call
CALL FKINSPBCG (MAXL, IER)Its arguments are the same as those with the same names for FKINSPGMR.
Optional outputs specific to the SPBCG case are listed in Table 6.2.
For descriptions of the relevant optional user-supplied routines, see User-supplied routines for SPGMR/SPBCG/SPTFQMR below.
[S][P] SPTFQMR treatment of the linear systems
For the Scaled Preconditioned Transpose-Reee Quasi-Minimal Residual solution of the linear systems, the user must make the call
CALL FKINSPTFQMR (MAXL, IER)Its arguments are the same as those with the same names for FKINSPGMR.
Optional outputs specific to the SPTFQMR case are listed in Table 6.2.
For descriptions of the relevant optional user-supplied routines, see below.
[S][P] Functions used by SPGMR/SPBCG/SPTFQMR
An optional user-supplied routine, FKINJTIMES (see below), can be provided for Jacobian-vector products. If it is, then, following the call to FKINSPGMR, FKINSPBCG, or FKINSPTFQMR, the user must make the call:
CALL FKINSPILSSETJAC (FLAG, IER)with FLAG 0 to specify use of the user-supplied Jacobian-times-vector approximation. The argument IER is an error return flag which is 0 for success or non-zero if an error occurred.
If preconditioning is to be done, then the user must call
CALL FKINSPILSSETPREC (FLAG, IER)with FLAG 0. The return flag IER is 0 if successful, or negative if a memory error occurred. In addition, the user program must include preconditioner routines FKPSOL and FKPSET (see below).
[S][P] User-supplied routines for SPGMR/SPBCG/SPTFQMR
With treatment of the linear systems by any of the Krylov iterative solvers, there are three optional user-supplied routines -- FKINJTIMES, FKPSOL, and FKPSET. The specifications for these routines are given below.
As an option when using the SPGMR, SPBCG, or SPTFQMR linear solvers, the user may supply a routine that computes the product of the system Jacobian J = F/u and a given vector v. If supplied, it must have the following form:
SUBROUTINE FKINJTIMES (V, FJV, NEWU, U, IER) DIMENSION V(*), FJV(*), U(*)Typically this routine will use only NEQ, U, V, and FJV. It must compute the product vector Jv, where the vector v is stored in V, and store the product in FJV. The input argument U contains the current value of u. On return, set IER = 0 if FKINJTIMES was successful, and nonzero otherwise.
If preconditioning is to be included, the following routine must be supplied, for solution of the preconditioner linear system:
SUBROUTINE FKPSOL (U, USCALE, FVAL, FSCALE, VTEM, FTEM, IER) DIMENSION U(*), USCALE(*), FVAL(*), FSCALE(*), VTEM(*), FTEM(*)Typically this routine will use only U, FVAL, VTEM and FTEM. It must solve the preconditioned linear system Pz = r, where r = VTEM is input, and store the solution z in VTEM as well. Here P is the right preconditioner. If scaling is being used, the routine supplied must also account for scaling on either coordinate or function value, as given in the arrays USCALE and FSCALE, respectively.
If the user's preconditioner requires that any Jacobian-related data be evaluated or preprocessed, then the following routine can be used for the evaluation and preprocessing of the preconditioner:
SUBROUTINE FKPSET (U, USCALE, FVAL, FSCALE, VTEMP1, VTEMP2, IER) DIMENSION U(*), USCALE(*), FVAL(*), FSCALE(*), VTEMP1(*), VTEMP2(*)It must perform any evaluation of Jacobian-related data and preprocessing needed for the solution of the preconditioned linear systems by FKPSOL. The variables U through FSCALE are for use in the preconditioning setup process. Typically, the system function FKFUN is called before any calls to FKPSET, so that FVAL will have been updated. U is the current solution iterate. The arrays VTEMP1 and VTEMP2 are available for work space. If scaling is being used, USCALE and FSCALE are available for those operations requiring scaling. NEQ is the problem size.
On return, set IER = 0 if FKPSET was successful or set IER = 1 if an error occurred.
If the user calls FKINSPILSSETPREC, the routine FKPSET must be provided, even if it is not needed and must return IER=0.
Solving the nonlinear system is accomplished by making the following call:
CALL FKINSOL (U, GLOBALSTRAT, USCALE, FSCALE, IER)The arguments are as follows. U is an array containing the initial guess on input, and the solution on return. GLOBALSTRAT is an integer (type INTEGER) defining the global strategy choice (1 specifies Inexact Newton, while 2 indicates line search). USCALE is an array of scaling factors for the U vector. FSCALE is an array of scaling factors for the FVAL vector. IER is an integer completion flag and will have one of the following values: 0 to indicate success, 1 to indicate that the initial guess satisfies F(u) = 0 within tolerances, 2 to indicate apparent stalling (small step), or a negative value to indicate an error or failure. These values correspond to the KINSol returns (see §5.5.3 and §10.2). The values of the optional outputs are available in IOPT and ROPT (see Table 6.2).
To free the internal memory created by the call to FKINMALLOC, make the call
CALL FKINFREE
In order to keep the number of user-callable FKINSOL interface routines to a minimum, optional inputs to the KINSOL solver are passed through only two routines: FKINSETIIN for integer optional inputs and FKINSETRIN for real optional inputs. These functions should be called as follows:
CALL FKINSETIIN (KEY, IVAL, IER) CALL FKINSETRIN (KEY, RVAL, IER)where KEY is a quoted string indicating which optional input is set (see Table 6.1), IVAL is the integer input value to the used, RVAL is the real input value to be used, and IER is an integer return flag which is set to 0 on success and a negative value if a filure occurred.
The optional outputs from the KINSOL solver are accessed not through individual functions, but rather through a pair of arrays, IOUT (integer type) of dimension at least 15, and ROUT (real type) of dimension at least 2. These arrays are owned (and allocated) by the user and are passed as arguments to FKINMALLOC. Table 6.2 lists the entries in these two arrays and specifies the optional variable as well as the KINSOL function which is actually called to extract the optional output.
For more details on the optional inputs and outputs, see §5.5.4 and §5.5.5.
The FKINBBD interface sub-module is a package of C functions which, as part of the FKINSOL interface module, support the use of the KINSOL solver with the parallel NVECTOR-PARALLEL module and the KINBBDPRE preconditioner module (see §5.7), for the solution of nonlinear problems in a mixed FORTRAN/C setting.
The user-callable functions in this package, with the corresponding KINSOL and KINBBDPRE functions, are as follows:
In addition to the FORTRAN right-hand side function FKFUN, the user-supplied functions used by this package, are listed below, each with the corresponding interface function which calls it (and its type within KINBBDPRE or KINSOL):
FKINBBD routine (FORTRAN) | KINSOL function (C) | KINSOL function type |
FKLOCFN | FKINgloc | KINLocalFn |
FKCOMMF | FKINgcomm | KINCommFn |
FKJTIMES | FKINJtimes | KINSpilsJacTimesVecFn |
The following is a summary of the usage of this module. Steps that are unchanged from the main program described in §6.3 are grayed-out.
To initialize the KINBBDPRE preconditioner, make the following call:
CALL FKINBBDINIT (NLOCAL, MUDQ, MLDQ, MU, ML, IER)The arguments are as follows. NLOCAL is the local size of vectors for this process. MUDQ and MLDQ are the upper and lower half-bandwidths to be used in the computation of the local Jacobian blocks by difference quotients; these may be smaller than the true half-bandwidths of the Jacobian of the local block of G, when smaller values may provide greater efficiency. MU and ML are the upper and lower half-bandwidths of the band matrix that is retained as an approximation of the local Jacobian block; these may be smaller than MUDQ and MLDQ. IER is a return completion flag. A value of 0 indicates success, while a value of -1 indicates that a memory failure occurred or that an input had an illegal value.
To specify the SPGMR linear system solver and use the KINBBDPRE preconditioner, make the following call:
CALL FKINBBDSPGMR (MAXL, MAXLRST, IER)Its arguments are the same as those of FKINSPGMR (see step 4 in §6.3).
To specify the SPBCG linear system solver and use the KINBBDPRE preconditioner, make the following call:
CALL FKINBBDSPBCG (MAXL, IER)Its arguments are the same as those of FKINSPGMR (see step 4 in §6.3).
To specify the SPTFQMR linear system solver and use the KINBBDPRE preconditioner, make the following call:
CALL FKINBBDSPTFQMR (MAXL, IER)Its arguments are the same as those of FKINSPGMR (see step 4 in §6.3).
Optionally, to specify that SPGMR, SPBCG, or SPTFQMR should use the supplied FKJTIMES, make the call
CALL FKINSPILSSETJAC (FLAG, IER)with FLAG 0. (see step 4 in §6.3).
Optional outputs specific to the SPGMR, SPBCG, or SPTFQMR solver are listed in Table 6.2. To obtain the optional outputs associated with the KINBBDPRE module, make the following call:
CALL FKINBBDOPT (LENRBBD, LENIBBD, NGEBBD)The arguments returned are as follows. LENRBBD is the length of real preconditioner work space, in realtype words. LENIBBD is the length of integer preconditioner work space, in integer words. These sizes are local to the current process. NGEBBD is the cumulative number of G(u) evaluations (calls to FKLOCFN) so far.
To free the internal memory created by the call to FKINBBDINIT, before calling FKINFREE and FNVFREEP, the user must call
CALL FKINBBDFREE
The following two routines must be supplied for use with the KINBBDPRE module:
SUBROUTINE FKLOCFN (NLOC, ULOC, GLOC, IER) DIMENSION ULOC(*), GLOC(*)This routine is to evaluate the function G(u) approximating F (possibly identical to F), in terms of the array ULOC (of length NLOC), which is the sub-vector of u local to this processor. The resulting (local) sub-vector is to be stored in the array GLOC. IER is an error return flag which should be set to 0 if successful, a positive value if a recoverable error occurred (in which case KINSOL will attempt to correct), or a negative value if FKLOCFN failed unrecoverably (in which case the solution process is halted).
SUBROUTINE FKCOMMFN (NLOC, ULOC, IER) DIMENSION ULOC(*)This routine is to perform the inter-processor communication necessary for the FKLOCFN routine. Each call to FKCOMMFN is preceded by a call to the system function routine FKFUN with the same argument ULOC. IER is an error return flag which should be set to 0 if successful, a positive value if a recoverable error occurred (in which case KINSOL will attempt to correct), or a negative value if FKCOMMFN failed recoverably (in which case the solution process is halted).
The subroutine FKCOMMFN must be supplied even if it is not needed and must return IER=0.
Optionally, the user can supply a routine FKINJTIMES for the evaluation of Jacobian-vector products, as described above in step 4 in §6.3.