00001 /* 00002 * ----------------------------------------------------------------- 00003 * $Revision: 1.2 $ 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 CPODES. 00014 * ----------------------------------------------------------------- 00015 */ 00016 00017 #ifndef _CPDIRECT_H 00018 #define _CPDIRECT_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 * C P D I R E C T C O N S T A N T S 00030 * ================================================================= 00031 */ 00032 00033 /* 00034 * ----------------------------------------------------------------- 00035 * CPDIRECT input constants 00036 * ----------------------------------------------------------------- 00037 * fact_type: is the type of constraint Jacobian factorization used 00038 * for the projection. For independent constraints (i.e. Jacobian 00039 * with full row rank) any of the following three options can be 00040 * used (although their computational cost varies, depending on 00041 * the number of states and constraints) 00042 * CPDIRECT_LU: use LU decomposition of G^T. 00043 * CPDIRECT_QR: use QR decomposition of G^T 00044 * CPDIRECT_SC: use Schur complement 00045 * If it is known (or suspected) that some constraints are redundant, 00046 * the following option should be used: 00047 * CPDIRECT_QRP: use QR with column pivoting on G^T. 00048 * ----------------------------------------------------------------- 00049 */ 00050 00051 /* fact_type */ 00052 00053 #define CPDIRECT_LU 1 00054 #define CPDIRECT_QR 2 00055 #define CPDIRECT_SC 3 00056 #define CPDIRECT_QRP 4 00057 00058 /* 00059 * ----------------------------------------------------------------- 00060 * CPDIRECT return values 00061 * ----------------------------------------------------------------- 00062 */ 00063 00064 #define CPDIRECT_SUCCESS 0 00065 #define CPDIRECT_MEM_NULL -1 00066 #define CPDIRECT_LMEM_NULL -2 00067 #define CPDIRECT_ILL_INPUT -3 00068 #define CPDIRECT_MEM_FAIL -4 00069 00070 /* Additional last_flag values */ 00071 00072 #define CPDIRECT_JACFUNC_UNRECVR -5 00073 #define CPDIRECT_JACFUNC_RECVR -6 00074 00075 /* 00076 * ================================================================= 00077 * F U N C T I O N T Y P E S 00078 * ================================================================= 00079 */ 00080 00081 /* 00082 * ----------------------------------------------------------------- 00083 * Types: CPDlsDenseJacExplFn and CPDlsDenseJacImplFn 00084 * ----------------------------------------------------------------- 00085 * 00086 * If the ODE is given in explicit form, a dense Jacobian 00087 * approximation function Jac must be of type CPDlsDenseJacFn. 00088 * Its parameters are: 00089 * 00090 * N is the problem size. 00091 * 00092 * Jac is the dense matrix (of type DlsMat) that will be loaded 00093 * by a CPDlsDenseJacFn with an approximation to the Jacobian 00094 * matrix J = (df_i/dy_j) at the point (t,y). 00095 * 00096 * t is the current value of the independent variable. 00097 * 00098 * y is the current value of the dependent variable vector, 00099 * namely the predicted value of y(t). 00100 * 00101 * fy is the vector f(t,y). 00102 * 00103 * jac_data is a pointer to user data - the same as the jac_data 00104 * parameter passed to CPDlsSetJacFn. 00105 * 00106 * tmp1, tmp2, and tmp3 are pointers to memory allocated for 00107 * vectors of length N which can be used by a CPDlsDenseJacFn 00108 * as temporary storage or work space. 00109 * 00110 * A CPDlsDenseJacFn should return 0 if successful, a positive 00111 * value if a recoverable error occurred, and a negative value if 00112 * an unrecoverable error occurred. 00113 * 00114 * ----------------------------------------------------------------- 00115 * 00116 * If the ODE is given in implicit form, a dense Jacobian 00117 * approximation function djac must be of type CPDlsDenseJacImplFn. 00118 * Its parameters are: 00119 * 00120 * N is the problem size, and length of all vector arguments. 00121 * 00122 * t is the current value of the independent variable t. 00123 * 00124 * y is the current value of the dependent variable vector, 00125 * namely the predicted value of y(t). 00126 * 00127 * yp is the current value of the derivative vector y', 00128 * namely the predicted value of y'(t). 00129 * 00130 * f is the residual vector F(tt,yy,yp). 00131 * 00132 * gm is the scalar in the system Jacobian, proportional to 00133 * the step size h. 00134 * 00135 * jac_data is a pointer to user Jacobian data - the same as the 00136 * jdata parameter passed to CPDenseSetJacFn. 00137 * 00138 * Jac is the dense matrix (of type DenseMat) to be loaded by 00139 * an CPDenseJacImplFn routine with an approximation to the 00140 * system Jacobian matrix 00141 * J = dF/dy' + gamma*dF/dy 00142 * at the given point (t,y,y'), where the ODE system is 00143 * given by F(t,y,y') = 0. Jac is preset to zero, so only 00144 * the nonzero elements need to be loaded. See note below. 00145 * 00146 * tmp1, tmp2, tmp3 are pointers to memory allocated for 00147 * N_Vectors which can be used by an CPDenseJacImplFn routine 00148 * as temporary storage or work space. 00149 * 00150 * A CPDenseJacImplFn should return 00151 * 0 if successful, 00152 * a positive int if a recoverable error occurred, or 00153 * a negative int if a nonrecoverable error occurred. 00154 * In the case of a recoverable error return, the integrator will 00155 * attempt to recover by reducing the stepsize (which changes cj). 00156 * 00157 * ----------------------------------------------------------------- 00158 * 00159 * NOTE: The following are two efficient ways to load a dense Jac: 00160 * (1) (with macros - no explicit data structure references) 00161 * for (j=0; j < Neq; j++) { 00162 * col_j = DENSE_COL(Jac,j); 00163 * for (i=0; i < Neq; i++) { 00164 * generate J_ij = the (i,j)th Jacobian element 00165 * col_j[i] = J_ij; 00166 * } 00167 * } 00168 * (2) (without macros - explicit data structure references) 00169 * for (j=0; j < Neq; j++) { 00170 * col_j = (Jac->data)[j]; 00171 * for (i=0; i < Neq; i++) { 00172 * generate J_ij = the (i,j)th Jacobian element 00173 * col_j[i] = J_ij; 00174 * } 00175 * } 00176 * A third way, using the DENSE_ELEM(A,i,j) macro, is much less 00177 * efficient in general. It is only appropriate for use in small 00178 * problems in which efficiency of access is NOT a major concern. 00179 * 00180 * NOTE: If the user's Jacobian routine needs other quantities, 00181 * they are accessible as follows: hcur (the current stepsize) 00182 * and ewt (the error weight vector) are accessible through 00183 * CPodeGetCurrentStep and CPodeGetErrWeights, respectively 00184 * (see cvode.h). The unit roundoff is available as 00185 * UNIT_ROUNDOFF defined in sundials_types.h. 00186 * 00187 * ----------------------------------------------------------------- 00188 */ 00189 00190 00191 typedef int (*CPDlsDenseJacExplFn)(int N, realtype t, 00192 N_Vector y, N_Vector fy, 00193 DlsMat Jac, void *jac_data, 00194 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); 00195 typedef int (*CPDlsDenseJacImplFn)(int N, realtype t, realtype gm, 00196 N_Vector y, N_Vector yp, N_Vector r, 00197 DlsMat Jac, void *jac_data, 00198 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); 00199 00200 /* 00201 * ----------------------------------------------------------------- 00202 * Types: CPDlsBandJacExplFn and CPDlsBandJacImplFn 00203 * ----------------------------------------------------------------- 00204 * 00205 * If the ODE is given in explicit form, a band Jacobian 00206 * approximation function Jac must be of type CPDlsBandJacExplFn. 00207 * Its parameters are: 00208 * 00209 * N is the length of all vector arguments. 00210 * 00211 * mupper is the upper half-bandwidth of the approximate banded 00212 * Jacobian. This parameter is the same as the mupper parameter 00213 * passed by the user to the linear solver initialization function. 00214 * 00215 * mlower is the lower half-bandwidth of the approximate banded 00216 * Jacobian. This parameter is the same as the mlower parameter 00217 * passed by the user to the linear solver initialization function. 00218 * 00219 * t is the current value of the independent variable. 00220 * 00221 * y is the current value of the dependent variable vector, 00222 * namely the predicted value of y(t). 00223 * 00224 * fy is the vector f(t,y). 00225 * 00226 * Jac is the band matrix (of type DlsMat) that will be loaded 00227 * by a CPDlsBandJacFn with an approximation to the Jacobian matrix 00228 * Jac = (df_i/dy_j) at the point (t,y). 00229 * jac_data is a pointer to user data - the same as the jac_data 00230 * parameter passed to CPDlsSetJacFn. 00231 * 00232 * A CPDlsBandJacExplFn should return 0 if successful, a positive 00233 * value if a recoverable error occurred, and a negative value if 00234 * an unrecoverable error occurred. 00235 * 00236 * ----------------------------------------------------------------- 00237 * 00238 * If the ODE is given in explicit form, a band Jacobian 00239 * approximation function Jac must be of type CPDlsBandJacImplFn. 00240 * Its parameters are: 00241 * 00242 * N is the problem size, and length of all vector arguments. 00243 * 00244 * mupper is the upper bandwidth of the banded Jacobian matrix. 00245 * 00246 * mlower is the lower bandwidth of the banded Jacobian matrix. 00247 * 00248 * t is the current value of the independent variable t. 00249 * 00250 * y is the current value of the dependent variable vector, 00251 * namely the predicted value of y(t). 00252 * 00253 * yp is the current value of the derivative vector y', 00254 * namely the predicted value of y'(t). 00255 * 00256 * r is the residual vector F(tt,yy,yp). 00257 * 00258 * gm is the scalar in the system Jacobian, proportional to 00259 * the step size h. 00260 * 00261 * jac_data is a pointer to user Jacobian data - the same as the 00262 * jdata parameter passed to CPBandSetJacFn. 00263 * 00264 * J is the band matrix (of type BandMat) to be loaded by 00265 * with an approximation to the system Jacobian matrix 00266 * J = dF/dy' + gamma*dF/dy 00267 * at the given point (t,y,y'), where the ODE system is 00268 * given by F(t,y,y') = 0. Jac is preset to zero, so only 00269 * the nonzero elements need to be loaded. See note below. 00270 * 00271 * tmp1, tmp2, tmp3 are pointers to memory allocated for 00272 * N_Vectors which can be used by an CPBandJacImplFn routine 00273 * as temporary storage or work space. 00274 * 00275 * A CPDlsBandJacImplFn should return 0 if successful, a positive 00276 * value if a recoverable error occurred, and a negative value if 00277 * an unrecoverable error occurred. 00278 * 00279 * ----------------------------------------------------------------- 00280 * 00281 * NOTE: Three efficient ways to load J are: 00282 * 00283 * (1) (with macros - no explicit data structure references) 00284 * for (j=0; j < n; j++) { 00285 * col_j = BAND_COL(Jac,j); 00286 * for (i=j-mupper; i <= j+mlower; i++) { 00287 * generate J_ij = the (i,j)th Jacobian element 00288 * BAND_COL_ELEM(col_j,i,j) = J_ij; 00289 * } 00290 * } 00291 * 00292 * (2) (with BAND_COL macro, but without BAND_COL_ELEM macro) 00293 * for (j=0; j < n; j++) { 00294 * col_j = BAND_COL(Jac,j); 00295 * for (k=-mupper; k <= mlower; k++) { 00296 * generate J_ij = the (i,j)th Jacobian element, i=j+k 00297 * col_j[k] = J_ij; 00298 * } 00299 * } 00300 * 00301 * (3) (without macros - explicit data structure references) 00302 * offset = Jac->smu; 00303 * for (j=0; j < n; j++) { 00304 * col_j = ((Jac->data)[j])+offset; 00305 * for (k=-mupper; k <= mlower; k++) { 00306 * generate J_ij = the (i,j)th Jacobian element, i=j+k 00307 * col_j[k] = J_ij; 00308 * } 00309 * } 00310 * Caution: Jac->smu is generally NOT the same as mupper. 00311 * 00312 * The BAND_ELEM(A,i,j) macro is appropriate for use in small 00313 * problems in which efficiency of access is NOT a major concern. 00314 * 00315 * NOTE: If the user's Jacobian routine needs other quantities, 00316 * they are accessible as follows: hcur (the current stepsize) 00317 * and ewt (the error weight vector) are accessible through 00318 * CPodeGetCurrentStep and CPodeGetErrWeights, respectively 00319 * (see cvode.h). The unit roundoff is available as 00320 * UNIT_ROUNDOFF defined in sundials_types.h 00321 * 00322 * tmp1, tmp2, and tmp3 are pointers to memory allocated for 00323 * vectors of length N which can be used by a CPDlsBandJacFn 00324 * as temporary storage or work space. 00325 * 00326 * ----------------------------------------------------------------- 00327 */ 00328 00329 typedef int (*CPDlsBandJacExplFn)(int N, int mupper, int mlower, 00330 realtype t, N_Vector y, N_Vector fy, 00331 DlsMat Jac, void *jac_data, 00332 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); 00333 typedef int (*CPDlsBandJacImplFn)(int N, int mupper, int mlower, 00334 realtype t, realtype gm, 00335 N_Vector y, N_Vector yp, N_Vector r, 00336 DlsMat Jac, void *jac_data, 00337 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); 00338 00339 /* 00340 * ----------------------------------------------------------------- 00341 * Type: CPDlsDenseProjJacFn 00342 * ----------------------------------------------------------------- 00343 * 00344 * 00345 * ----------------------------------------------------------------- 00346 */ 00347 00348 typedef int (*CPDlsDenseProjJacFn)(int Nc, int Ny, 00349 realtype t, N_Vector y, N_Vector cy, 00350 DlsMat Jac, void *jac_data, 00351 N_Vector tmp1, N_Vector tmp2); 00352 00353 00354 00355 /* 00356 * ================================================================= 00357 * E X P O R T E D F U N C T I O N S 00358 * ================================================================= 00359 */ 00360 00361 /* 00362 * ----------------------------------------------------------------- 00363 * Optional inputs to a CPDIRECT linear solver 00364 * for implicit integration 00365 * ----------------------------------------------------------------- 00366 * 00367 * CPDlsSetJacFn specifies the Jacobian approximation routine to be 00368 * used. When using dense Jacobians, a user-supplied jac routine must 00369 * be of type CPDlsDenseJacFn. When using banded Jacobians, a 00370 * user-supplied jac routine must be of type CPDlsBandJacFn. 00371 * By default, a difference quotient approximation, supplied with this 00372 * solver is used. 00373 * CPDlsSetJacFn also specifies a pointer to user data which is 00374 * passed to the user's jac routine every time it is called. 00375 * 00376 * The return value of CPDlsSetJacFn is one of: 00377 * CPDIRECT_SUCCESS if successful 00378 * CPDIRECT_MEM_NULL if the CPODES memory was NULL 00379 * CPDIRECT_LMEM_NULL if the linear solver memory was NULL 00380 * ----------------------------------------------------------------- 00381 */ 00382 00383 SUNDIALS_EXPORT int CPDlsSetJacFn(void *cvode_mem, void *jac, void *jac_data); 00384 00385 /* 00386 * ----------------------------------------------------------------- 00387 * Optional outputs from a CPDIRECT linear solver 00388 * for implicit integration 00389 * ----------------------------------------------------------------- 00390 * 00391 * CPDlsGetWorkSpace returns the real and integer workspace used 00392 * by the direct linear solver. 00393 * CPDlsGetNumJacEvals returns the number of calls made to the 00394 * Jacobian evaluation routine jac. 00395 * CPDlsGetNumFctEvals returns the number of calls to the user 00396 * f routine due to finite difference Jacobian 00397 * evaluation. 00398 * CPDlsGetLastFlag returns the last error flag set by any of 00399 * the CPDIRECT interface functions. 00400 * 00401 * The return value of CPDlsGet* is one of: 00402 * CPDIRECT_SUCCESS if successful 00403 * CPDIRECT_MEM_NULL if the CPODES memory was NULL 00404 * CPDIRECT_LMEM_NULL if the linear solver memory was NULL 00405 * ----------------------------------------------------------------- 00406 */ 00407 00408 SUNDIALS_EXPORT int CPDlsGetWorkSpace(void *cvode_mem, long int *lenrwLS, long int *leniwLS); 00409 SUNDIALS_EXPORT int CPDlsGetNumJacEvals(void *cvode_mem, long int *njevals); 00410 SUNDIALS_EXPORT int CPDlsGetNumFctEvals(void *cvode_mem, long int *nfevalsLS); 00411 SUNDIALS_EXPORT int CPDlsGetLastFlag(void *cvode_mem, int *flag); 00412 00413 /* 00414 * ----------------------------------------------------------------- 00415 * The following function returns the name of the constant 00416 * associated with a CPDIRECT return flag 00417 * ----------------------------------------------------------------- 00418 */ 00419 00420 SUNDIALS_EXPORT char *CPDlsGetReturnFlagName(int flag); 00421 00422 /* 00423 * ----------------------------------------------------------------- 00424 * Optional I/O functions for a CPDIRECT linear solver 00425 * for projection 00426 * ----------------------------------------------------------------- 00427 */ 00428 00429 SUNDIALS_EXPORT int CPDlsProjSetJacFn(void *cpode_mem, void *jacP, void *jacP_data); 00430 00431 SUNDIALS_EXPORT int CPDlsProjGetNumJacEvals(void *cpode_mem, long int *njPevals); 00432 SUNDIALS_EXPORT int CPDlsProjGetNumFctEvals(void *cpode_mem, long int *ncevalsLS); 00433 00434 #ifdef __cplusplus 00435 } 00436 #endif 00437 00438 #endif
1.5.6