00001 /* 00002 * ----------------------------------------------------------------- 00003 * $Revision: 1.2 $ 00004 * $Date: 2006/11/29 00:05:07 $ 00005 * ----------------------------------------------------------------- 00006 * Programmer(s): Scott D. Cohen and Alan C. Hindmarsh @ LLNL 00007 * ----------------------------------------------------------------- 00008 * Copyright (c) 2002, 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 * This header file contains declarations intended for use by 00014 * generic iterative solvers of Ax = b. The enumeration gives 00015 * symbolic names for the type of preconditioning to be used. 00016 * The function type declarations give the prototypes for the 00017 * functions to be called within an iterative linear solver, that 00018 * are responsible for 00019 * multiplying A by a given vector v (ATimesFn), and 00020 * solving the preconditioner equation Pz = r (PSolveFn). 00021 * ----------------------------------------------------------------- 00022 */ 00023 00024 #ifndef _ITERATIVE_H 00025 #define _ITERATIVE_H 00026 00027 #ifdef __cplusplus /* wrapper to enable C++ usage */ 00028 extern "C" { 00029 #endif 00030 00031 #include <sundials/sundials_nvector.h> 00032 00033 00034 /* 00035 * ----------------------------------------------------------------- 00036 * enum : types of preconditioning 00037 * ----------------------------------------------------------------- 00038 * PREC_NONE : The iterative linear solver should not use 00039 * preconditioning. 00040 * 00041 * PREC_LEFT : The iterative linear solver uses preconditioning on 00042 * the left only. 00043 * 00044 * PREC_RIGHT : The iterative linear solver uses preconditioning on 00045 * the right only. 00046 * 00047 * PREC_BOTH : The iterative linear solver uses preconditioning on 00048 * both the left and the right. 00049 * ----------------------------------------------------------------- 00050 */ 00051 00052 enum { PREC_NONE, PREC_LEFT, PREC_RIGHT, PREC_BOTH }; 00053 00054 /* 00055 * ----------------------------------------------------------------- 00056 * enum : types of Gram-Schmidt routines 00057 * ----------------------------------------------------------------- 00058 * MODIFIED_GS : The iterative solver uses the modified 00059 * Gram-Schmidt routine ModifiedGS listed in this 00060 * file. 00061 * 00062 * CLASSICAL_GS : The iterative solver uses the classical 00063 * Gram-Schmidt routine ClassicalGS listed in this 00064 * file. 00065 * ----------------------------------------------------------------- 00066 */ 00067 00068 enum { MODIFIED_GS = 1, CLASSICAL_GS = 2 }; 00069 00070 /* 00071 * ----------------------------------------------------------------- 00072 * Type: ATimesFn 00073 * ----------------------------------------------------------------- 00074 * An ATimesFn multiplies Av and stores the result in z. The 00075 * caller is responsible for allocating memory for the z vector. 00076 * The parameter A_data is a pointer to any information about A 00077 * which the function needs in order to do its job. The vector v 00078 * is unchanged. An ATimesFn returns 0 if successful and a 00079 * non-zero value if unsuccessful. 00080 * ----------------------------------------------------------------- 00081 */ 00082 00083 typedef int (*ATimesFn)(void *A_data, N_Vector v, N_Vector z); 00084 00085 /* 00086 * ----------------------------------------------------------------- 00087 * Type: PSolveFn 00088 * ----------------------------------------------------------------- 00089 * A PSolveFn solves the preconditioner equation Pz = r for the 00090 * vector z. The caller is responsible for allocating memory for 00091 * the z vector. The parameter P_data is a pointer to any 00092 * information about P which the function needs in order to do 00093 * its job. The parameter lr is input, and indicates whether P 00094 * is to be taken as the left preconditioner or the right 00095 * preconditioner: lr = 1 for left and lr = 2 for right. 00096 * If preconditioning is on one side only, lr can be ignored. 00097 * The vector r is unchanged. 00098 * A PSolveFn returns 0 if successful and a non-zero value if 00099 * unsuccessful. On a failure, a negative return value indicates 00100 * an unrecoverable condition, while a positive value indicates 00101 * a recoverable one, in which the calling routine may reattempt 00102 * the solution after updating preconditioner data. 00103 * ----------------------------------------------------------------- 00104 */ 00105 00106 typedef int (*PSolveFn)(void *P_data, N_Vector r, N_Vector z, int lr); 00107 00108 /* 00109 * ----------------------------------------------------------------- 00110 * Function: ModifiedGS 00111 * ----------------------------------------------------------------- 00112 * ModifiedGS performs a modified Gram-Schmidt orthogonalization 00113 * of the N_Vector v[k] against the p unit N_Vectors at 00114 * v[k-1], v[k-2], ..., v[k-p]. 00115 * 00116 * v is an array of (k+1) N_Vectors v[i], i=0, 1, ..., k. 00117 * v[k-1], v[k-2], ..., v[k-p] are assumed to have L2-norm 00118 * equal to 1. 00119 * 00120 * h is the output k by k Hessenberg matrix of inner products. 00121 * This matrix must be allocated row-wise so that the (i,j)th 00122 * entry is h[i][j]. The inner products (v[i],v[k]), 00123 * i=i0, i0+1, ..., k-1, are stored at h[i][k-1]. Here 00124 * i0=MAX(0,k-p). 00125 * 00126 * k is the index of the vector in the v array that needs to be 00127 * orthogonalized against previous vectors in the v array. 00128 * 00129 * p is the number of previous vectors in the v array against 00130 * which v[k] is to be orthogonalized. 00131 * 00132 * new_vk_norm is a pointer to memory allocated by the caller to 00133 * hold the Euclidean norm of the orthogonalized vector v[k]. 00134 * 00135 * If (k-p) < 0, then ModifiedGS uses p=k. The orthogonalized 00136 * v[k] is NOT normalized and is stored over the old v[k]. Once 00137 * the orthogonalization has been performed, the Euclidean norm 00138 * of v[k] is stored in (*new_vk_norm). 00139 * 00140 * ModifiedGS returns 0 to indicate success. It cannot fail. 00141 * ----------------------------------------------------------------- 00142 */ 00143 00144 SUNDIALS_EXPORT int ModifiedGS(N_Vector *v, realtype **h, int k, int p, 00145 realtype *new_vk_norm); 00146 00147 /* 00148 * ----------------------------------------------------------------- 00149 * Function: ClassicalGS 00150 * ----------------------------------------------------------------- 00151 * ClassicalGS performs a classical Gram-Schmidt 00152 * orthogonalization of the N_Vector v[k] against the p unit 00153 * N_Vectors at v[k-1], v[k-2], ..., v[k-p]. The parameters v, h, 00154 * k, p, and new_vk_norm are as described in the documentation 00155 * for ModifiedGS. 00156 * 00157 * temp is an N_Vector which can be used as workspace by the 00158 * ClassicalGS routine. 00159 * 00160 * s is a length k array of realtype which can be used as 00161 * workspace by the ClassicalGS routine. 00162 * 00163 * ClassicalGS returns 0 to indicate success. It cannot fail. 00164 * ----------------------------------------------------------------- 00165 */ 00166 00167 SUNDIALS_EXPORT int ClassicalGS(N_Vector *v, realtype **h, int k, int p, 00168 realtype *new_vk_norm, N_Vector temp, realtype *s); 00169 00170 /* 00171 * ----------------------------------------------------------------- 00172 * Function: QRfact 00173 * ----------------------------------------------------------------- 00174 * QRfact performs a QR factorization of the Hessenberg matrix H. 00175 * 00176 * n is the problem size; the matrix H is (n+1) by n. 00177 * 00178 * h is the (n+1) by n Hessenberg matrix H to be factored. It is 00179 * stored row-wise. 00180 * 00181 * q is an array of length 2*n containing the Givens rotations 00182 * computed by this function. A Givens rotation has the form: 00183 * | c -s | 00184 * | s c |. 00185 * The components of the Givens rotations are stored in q as 00186 * (c, s, c, s, ..., c, s). 00187 * 00188 * job is a control flag. If job==0, then a new QR factorization 00189 * is performed. If job!=0, then it is assumed that the first 00190 * n-1 columns of h have already been factored and only the last 00191 * column needs to be updated. 00192 * 00193 * QRfact returns 0 if successful. If a zero is encountered on 00194 * the diagonal of the triangular factor R, then QRfact returns 00195 * the equation number of the zero entry, where the equations are 00196 * numbered from 1, not 0. If QRsol is subsequently called in 00197 * this situation, it will return an error because it could not 00198 * divide by the zero diagonal entry. 00199 * ----------------------------------------------------------------- 00200 */ 00201 00202 SUNDIALS_EXPORT int QRfact(int n, realtype **h, realtype *q, int job); 00203 00204 /* 00205 * ----------------------------------------------------------------- 00206 * Function: QRsol 00207 * ----------------------------------------------------------------- 00208 * QRsol solves the linear least squares problem 00209 * 00210 * min (b - H*x, b - H*x), x in R^n, 00211 * 00212 * where H is a Hessenberg matrix, and b is in R^(n+1). 00213 * It uses the QR factors of H computed by QRfact. 00214 * 00215 * n is the problem size; the matrix H is (n+1) by n. 00216 * 00217 * h is a matrix (computed by QRfact) containing the upper 00218 * triangular factor R of the original Hessenberg matrix H. 00219 * 00220 * q is an array of length 2*n (computed by QRfact) containing 00221 * the Givens rotations used to factor H. 00222 * 00223 * b is the (n+1)-vector appearing in the least squares problem 00224 * above. 00225 * 00226 * On return, b contains the solution x of the least squares 00227 * problem, if QRsol was successful. 00228 * 00229 * QRsol returns a 0 if successful. Otherwise, a zero was 00230 * encountered on the diagonal of the triangular factor R. 00231 * In this case, QRsol returns the equation number (numbered 00232 * from 1, not 0) of the zero entry. 00233 * ----------------------------------------------------------------- 00234 */ 00235 00236 SUNDIALS_EXPORT int QRsol(int n, realtype **h, realtype *q, realtype *b); 00237 00238 #ifdef __cplusplus 00239 } 00240 #endif 00241 00242 #endif
1.5.6