next up previous contents index
Next: 7. Description of the Up: User Documentation for KINSOL Previous: 5. Using KINSOL for   Contents   Index

Subsections


6. FKINSOL, an Interface Module for FORTRAN Applications

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.


6.1 Important note on portability

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.


6.2 FKINSOL routines

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
In contrast to the case of direct use of KINSOL, the names of all user-supplied routines here are fixed, in order to maximize portability for the resulting mixed-language program.


6.3 Usage of the FKINSOL interface module

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.

  1. Nonlinear system function specification

    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).

  2. NVECTOR module initialization

    [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.

  3. Problem specification

    To set various problem and solution parameters and allocate internal memory, make the following call:

    FKINMALLOC
    Call
       CALL FKINMALLOC (IOUT, ROUT, IER)    
    Description
    This function specifies the optional output arrays, allocates internal memory, and initializes KINSOL.
    Arguments
    IOUT
    is an integer array for integer optional outputs.
    ROUT
    is a real array for real optional outputs.
    Return value
    IER is the return completion flag. Values are 0 for successful return and -1 otherwise. See printed message for details in case of failure.
    Notes
    The user integer data array IOUT must be declared as INTEGER*4 or INTEGER*8 according to the C type long int.

    The optional outputs associated with the main KINSOL integrator are listed in Table 6.2.

  4. Linear solver specification

    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

    The user must make the call:

          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 = $ \partial$F/$ \partial$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 $ \neq$ 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

    The user must make the call:

          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 = $ \partial$F/$ \partial$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 $ \neq$ 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 $ \neq$ 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 $ \neq$ 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 = $ \partial$F/$ \partial$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.

  5. Problem solution

    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).

  6. Memory deallocation

    To free the internal memory created by the call to FKINMALLOC, make the call

          CALL FKINFREE
    

6.4 FKINSOL optional input and output

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.



Table 6.1: Keys for setting FKINSOL optional inputs
Integer optional inputs FKINSETIIN
Key Optional input Default value
PRNT_LEVEL Verbosity level of output 0
MAX_NITER Maximum no. of nonlinear iterations 200
ETA_FORM Form of $ \eta$ coefficient 1 (KIN_ETACHOICE1)
MAX_SETUPS Maximum no. of iterations without prec. setup 10
MAX_SP_SETUPS Maximum no. of iterations without residual check 5
NO_INIT_SETUP No initial preconditioner setup FALSE
NO_MIN_EPS Lower bound on $ \epsilon$ FALSE
NO_RES_MON No residual monitoring FALSE
 
Real optional inputs (FKINSETRIN)
Key Optional input Default value
FNORM_TOL Function-norm stopping tolerance $ \sqrt[3]{{\text{uround}}}$
SSTEP_TOL Scaled-step stopping tolerance uround2/3
MAX_STEP Max. scaled length of Newton step 1000| Duu0|2
RERR_FUNC Relative error for F.D. Jv $ \sqrt{{\text{uround}}}$
ETA_CONST Constant value of $ \eta$ 0.1
ETA_PARAMS Values of $ \gamma$ and $ \alpha$ 0.9 and 2.0
RMON_CONST Constant value of $ \omega$ 0.9
RMON_PARAMS Values of $ \omega_{{min}}^{}$ and $ \omega_{{max}}^{}$ 0.00001 and 0.9



Table 6.2: Description of the FKINSOL optional output arrays IOUT and ROUT
Integer output array IOUT
Index Optional output KINSOL function
KINSOL main solver
1 LENRW KINGetWorkSpace
2 LENIW KINGetWorkSpace
3 NNI KINGetNumNonlinSolvIters
4 NFE KINGetNumFuncEvals
5 NBCF KINGetNumBetaCondFails
6 NBKTRK KINGetNumBacktrackOps
KINDENSE linear solver
7 LENRWLS KINDenseGetWorkSpace
8 LENIWLS KINDenseGetWorkSpace
9 LS_FLAG KINDenseGetLastFlag
10 NFELS KINDenseGetNumFuncEvals
11 NJE KINDenseGetNumJacEvals
KINBAND linear solver
7 LENRWLS KINBandGetWorkSpace
8 LENIWLS KINBandGetWorkSpace
9 LS_FLAG KINBandGetLastFlag
10 NFELS KINBandGetNumFuncEvals
11 NJE KINBandGetNumJacEvals
KINSPGMR, KINSPBCG, KINSPTFQMR linear solvers
7 LENRWLS KINSpilsGetWorkSpace
8 LENIWLS KINSpilsGetWorkSpace
9 LS_FLAG KINSpilsGetLastFlag
10 NFELS KINSpilsGetNumFuncEvals
11 NJTV KINSpilsGetNumJacEvals
12 NPE KINSpilsGetNumPrecEvals
13 NPS KINSpilsGetNumPrecSolves
14 NLI KINSpilsGetNumLinIters
15 NCFL KINSpilsGetNumConvFails
 
Real output array ROUT
Index Optional output KINSOL function
1 FNORM KINGetFuncNorm
2 SSTEP KINGetStepLength

6.5 Usage of the FKINBBD interface to KINBBDPRE

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
As with the rest of the FKINSOL routines, the names of all user-supplied routines here are fixed, in order to maximize portability for the resulting mixed-language program. Additionally, based on flags discussed above in §6.2, the names of the user-supplied routines are mapped to actual values through a series of definitions in the header file fkinbbd.h.

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.

  1. Nonlinear system function specification

  2. NVECTOR module initialization

  3. Problem specification

  4. Linear solver specification

    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 $ \neq$ 0. (see step 4 in §6.3).

  5. Problem solution

  6. KINBBDPRE Optional outputs

    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.

  7. Memory deallocation

    To free the internal memory created by the call to FKINBBDINIT, before calling FKINFREE and FNVFREEP, the user must call

          CALL FKINBBDFREE
    

  8. User-supplied routines

    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.




next up previous contents index
Next: 7. Description of the Up: User Documentation for KINSOL Previous: 5. Using KINSOL for   Contents   Index
Radu Serban 2006-11-06