00001 /* 00002 * ----------------------------------------------------------------- 00003 * $Revision: 1.3 $ 00004 * $Date: 2006/11/29 00:05:06 $ 00005 * ----------------------------------------------------------------- 00006 * Programmer: Radu Serban @ LLNL 00007 * ----------------------------------------------------------------- 00008 * Copyright (c) 2006, The Regents of the University of California. 00009 * Produced at the Lawrence Livermore National Laboratory. 00010 * All rights reserved. 00011 * For details, see the LICENSE file. 00012 * ----------------------------------------------------------------- 00013 * Common header file for the direct linear solvers in IDA. 00014 * ----------------------------------------------------------------- 00015 */ 00016 00017 #ifndef _IDADIRECT_H 00018 #define _IDADIRECT_H 00019 00020 #ifdef __cplusplus /* wrapper to enable C++ usage */ 00021 extern "C" { 00022 #endif 00023 00024 #include <sundials/sundials_direct.h> 00025 #include <sundials/sundials_nvector.h> 00026 00027 /* 00028 * ================================================================= 00029 * I D A D I R E C T C O N S T A N T S 00030 * ================================================================= 00031 */ 00032 00033 /* 00034 * ----------------------------------------------------------------- 00035 * IDADIRECT return values 00036 * ----------------------------------------------------------------- 00037 */ 00038 00039 #define IDADIRECT_SUCCESS 0 00040 #define IDADIRECT_MEM_NULL -1 00041 #define IDADIRECT_LMEM_NULL -2 00042 #define IDADIRECT_ILL_INPUT -3 00043 #define IDADIRECT_MEM_FAIL -4 00044 00045 /* Additional last_flag values */ 00046 00047 #define IDADIRECT_JACFUNC_UNRECVR -5 00048 #define IDADIRECT_JACFUNC_RECVR -6 00049 00050 /* 00051 * ================================================================= 00052 * F U N C T I O N T Y P E S 00053 * ================================================================= 00054 */ 00055 00056 /* 00057 * ----------------------------------------------------------------- 00058 * Types : IDADlsDenseJacFn 00059 * ----------------------------------------------------------------- 00060 * 00061 * A dense Jacobian approximation function djac must be of type 00062 * IDADlsDenseJacFn. 00063 * Its parameters are: 00064 * 00065 * N is the problem size, and length of all vector arguments. 00066 * 00067 * t is the current value of the independent variable t. 00068 * 00069 * y is the current value of the dependent variable vector, 00070 * namely the predicted value of y(t). 00071 * 00072 * yp is the current value of the derivative vector y', 00073 * namely the predicted value of y'(t). 00074 * 00075 * f is the residual vector F(tt,yy,yp). 00076 * 00077 * c_j is the scalar in the system Jacobian, proportional to 00078 * the inverse of the step size h. 00079 * 00080 * jac_data is a pointer to user Jacobian data - the same as the 00081 * jdata parameter passed to IDADlsSetJacFn. 00082 * 00083 * Jac is the dense matrix (of type DlsMat) to be loaded by 00084 * an IDADlsDenseJacFn routine with an approximation to the 00085 * system Jacobian matrix 00086 * J = dF/dy' + gamma*dF/dy 00087 * at the given point (t,y,y'), where the ODE system is 00088 * given by F(t,y,y') = 0. 00089 * Note that Jac is NOT preset to zero! 00090 * 00091 * tmp1, tmp2, tmp3 are pointers to memory allocated for 00092 * N_Vectors which can be used by an IDADlsDenseJacFn routine 00093 * as temporary storage or work space. 00094 * 00095 * A IDADlsDenseJacFn should return 00096 * 0 if successful, 00097 * a positive int if a recoverable error occurred, or 00098 * a negative int if a nonrecoverable error occurred. 00099 * In the case of a recoverable error return, the integrator will 00100 * attempt to recover by reducing the stepsize (which changes cj). 00101 * 00102 * ----------------------------------------------------------------- 00103 * 00104 * NOTE: The following are two efficient ways to load a dense Jac: 00105 * (1) (with macros - no explicit data structure references) 00106 * for (j=0; j < Neq; j++) { 00107 * col_j = LAPACK_DENSE_COL(Jac,j); 00108 * for (i=0; i < Neq; i++) { 00109 * generate J_ij = the (i,j)th Jacobian element 00110 * col_j[i] = J_ij; 00111 * } 00112 * } 00113 * (2) (without macros - explicit data structure references) 00114 * for (j=0; j < Neq; j++) { 00115 * col_j = (Jac->data)[j]; 00116 * for (i=0; i < Neq; i++) { 00117 * generate J_ij = the (i,j)th Jacobian element 00118 * col_j[i] = J_ij; 00119 * } 00120 * } 00121 * A third way, using the LAPACK_DENSE_ELEM(A,i,j) macro, is much less 00122 * efficient in general. It is only appropriate for use in small 00123 * problems in which efficiency of access is NOT a major concern. 00124 * 00125 * NOTE: If the user's Jacobian routine needs other quantities, 00126 * they are accessible as follows: hcur (the current stepsize) 00127 * and ewt (the error weight vector) are accessible through 00128 * IDAGetCurrentStep and IDAGetErrWeights, respectively 00129 * (see ida.h). The unit roundoff is available as 00130 * UNIT_ROUNDOFF defined in sundials_types.h. 00131 * 00132 * ----------------------------------------------------------------- 00133 */ 00134 00135 00136 typedef int (*IDADlsDenseJacFn)(int N, realtype t, realtype c_j, 00137 N_Vector y, N_Vector yp, N_Vector r, 00138 DlsMat Jac, void *jac_data, 00139 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); 00140 00141 /* 00142 * ----------------------------------------------------------------- 00143 * Types : IDADlsBandJacFn 00144 * ----------------------------------------------------------------- 00145 * A banded Jacobian approximation function bjac must have the 00146 * prototype given below. Its parameters are: 00147 * 00148 * Neq is the problem size, and length of all vector arguments. 00149 * 00150 * mupper is the upper bandwidth of the banded Jacobian matrix. 00151 * 00152 * mlower is the lower bandwidth of the banded Jacobian matrix. 00153 * 00154 * tt is the current value of the independent variable t. 00155 * 00156 * yy is the current value of the dependent variable vector, 00157 * namely the predicted value of y(t). 00158 * 00159 * yp is the current value of the derivative vector y', 00160 * namely the predicted value of y'(t). 00161 * 00162 * rr is the residual vector F(tt,yy,yp). 00163 * 00164 * c_j is the scalar in the system Jacobian, proportional to 1/hh. 00165 * 00166 * jac_data is a pointer to user Jacobian data - the same as the 00167 * jdata parameter passed to IDADlsSetJacFn. 00168 * 00169 * Jac is the band matrix (of type BandMat) to be loaded by 00170 * an IDADlsBandJacFn routine with an approximation to the 00171 * system Jacobian matrix 00172 * J = dF/dy + cj*dF/dy' 00173 * at the given point (t,y,y'), where the DAE system is 00174 * given by F(t,y,y') = 0. Jac is preset to zero, so only 00175 * the nonzero elements need to be loaded. See note below. 00176 * 00177 * tmp1, tmp2, tmp3 are pointers to memory allocated for 00178 * N_Vectors which can be used by an IDADlsBandJacFn routine 00179 * as temporary storage or work space. 00180 * 00181 * An IDADlsBandJacFn function should return 00182 * 0 if successful, 00183 * a positive int if a recoverable error occurred, or 00184 * a negative int if a nonrecoverable error occurred. 00185 * In the case of a recoverable error return, the integrator will 00186 * attempt to recover by reducing the stepsize (which changes cj). 00187 * 00188 * ----------------------------------------------------------------- 00189 * 00190 * NOTE: The following are two efficient ways to load Jac: 00191 * 00192 * (1) (with macros - no explicit data structure references) 00193 * for (j=0; j < Neq; j++) { 00194 * col_j = BAND_COL(Jac,j); 00195 * for (i=j-mupper; i <= j+mlower; i++) { 00196 * generate J_ij = the (i,j)th Jacobian element 00197 * BAND_COL_ELEM(col_j,i,j) = J_ij; 00198 * } 00199 * } 00200 * 00201 * (2) (with BAND_COL macro, but without BAND_COL_ELEM macro) 00202 * for (j=0; j < Neq; j++) { 00203 * col_j = BAND_COL(Jac,j); 00204 * for (k=-mupper; k <= mlower; k++) { 00205 * generate J_ij = the (i,j)th Jacobian element, i=j+k 00206 * col_j[k] = J_ij; 00207 * } 00208 * } 00209 * 00210 * A third way, using the BAND_ELEM(A,i,j) macro, is much less 00211 * efficient in general. It is only appropriate for use in small 00212 * problems in which efficiency of access is NOT a major concern. 00213 * 00214 * NOTE: If the user's Jacobian routine needs other quantities, 00215 * they are accessible as follows: hcur (the current stepsize) 00216 * and ewt (the error weight vector) are accessible through 00217 * IDAGetCurrentStep and IDAGetErrWeights, respectively (see 00218 * ida.h). The unit roundoff is available as 00219 * UNIT_ROUNDOFF defined in sundials_types.h 00220 * 00221 * ----------------------------------------------------------------- 00222 */ 00223 00224 typedef int (*IDADlsBandJacFn)(int N, int mupper, int mlower, 00225 realtype t, realtype c_j, 00226 N_Vector y, N_Vector yp, N_Vector r, 00227 DlsMat Jac, void *jac_data, 00228 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); 00229 00230 /* 00231 * ================================================================= 00232 * E X P O R T E D F U N C T I O N S 00233 * ================================================================= 00234 */ 00235 00236 /* 00237 * ----------------------------------------------------------------- 00238 * Optional inputs to the IDADIRECT linear solver 00239 * ----------------------------------------------------------------- 00240 * 00241 * IDADlsSetJacFn specifies the Jacobian approximation routine to be 00242 * used. When using dense Jacobians, a user-supplied jac routine must 00243 * be of type IDADlsDenseJacFn. When using banded Jacobians, a 00244 * user-supplied jac routine must be of type IDADlsBandJacFn. 00245 * By default, a difference quotient approximation, supplied with this 00246 * solver is used. 00247 * IDADlsSetJacFn also specifies a pointer to user data which is 00248 * passed to the user's jac routine every time it is called. 00249 * 00250 * The return value of IDADlsSetJacFn is one of: 00251 * IDADIRECT_SUCCESS if successful 00252 * IDADIRECT_MEM_NULL if the IDA memory was NULL 00253 * IDADIRECT_LMEM_NULL if the linear solver memory was NULL 00254 * ----------------------------------------------------------------- 00255 */ 00256 00257 SUNDIALS_EXPORT int IDADlsSetJacFn(void *ida_mem, void *jac, void *jac_data); 00258 00259 /* 00260 * ----------------------------------------------------------------- 00261 * Optional outputs from the IDADIRECT linear solver 00262 * ----------------------------------------------------------------- 00263 * 00264 * IDADlsGetWorkSpace returns the real and integer workspace used 00265 * by the direct linear solver. 00266 * IDADlsGetNumJacEvals returns the number of calls made to the 00267 * Jacobian evaluation routine jac. 00268 * IDADlsGetNumResEvals returns the number of calls to the user 00269 * f routine due to finite difference Jacobian 00270 * evaluation. 00271 * IDADlsGetLastFlag returns the last error flag set by any of 00272 * the IDADIRECT interface functions. 00273 * 00274 * The return value of IDADlsGet* is one of: 00275 * IDADIRECT_SUCCESS if successful 00276 * IDADIRECT_MEM_NULL if the IDA memory was NULL 00277 * IDADIRECT_LMEM_NULL if the linear solver memory was NULL 00278 * ----------------------------------------------------------------- 00279 */ 00280 00281 SUNDIALS_EXPORT int IDADlsGetWorkSpace(void *ida_mem, long int *lenrwLS, long int *leniwLS); 00282 SUNDIALS_EXPORT int IDADlsGetNumJacEvals(void *ida_mem, long int *njevals); 00283 SUNDIALS_EXPORT int IDADlsGetNumResEvals(void *ida_mem, long int *nfevalsLS); 00284 SUNDIALS_EXPORT int IDADlsGetLastFlag(void *ida_mem, int *flag); 00285 00286 /* 00287 * ----------------------------------------------------------------- 00288 * The following function returns the name of the constant 00289 * associated with a IDADIRECT return flag 00290 * ----------------------------------------------------------------- 00291 */ 00292 00293 SUNDIALS_EXPORT char *IDADlsGetReturnFlagName(int flag); 00294 00295 00296 #ifdef __cplusplus 00297 } 00298 #endif 00299 00300 #endif