1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: snesmfj.c,v 1.53 1997/07/02 22:27:24 bsmith Exp balay $"; 3 #endif 4 5 #include "src/snes/snesimpl.h" /*I "snes.h" I*/ 6 #include "pinclude/pviewer.h" 7 #include <math.h> 8 9 typedef struct { /* default context for matrix-free SNES */ 10 SNES snes; /* SNES context */ 11 Vec w; /* work vector */ 12 PCNullSpace sp; /* null space context */ 13 double error_rel; /* square root of relative error in computing function */ 14 double umin; /* minimum allowable u'a value relative to |u|_1 */ 15 } MFCtx_Private; 16 17 #undef __FUNC__ 18 #define __FUNC__ "SNESMatrixFreeDestroy_Private" /* ADIC Ignore */ 19 int SNESMatrixFreeDestroy_Private(PetscObject obj) 20 { 21 int ierr; 22 Mat mat = (Mat) obj; 23 MFCtx_Private *ctx; 24 25 ierr = MatShellGetContext(mat,(void **)&ctx); 26 ierr = VecDestroy(ctx->w); CHKERRQ(ierr); 27 if (ctx->sp) {ierr = PCNullSpaceDestroy(ctx->sp); CHKERRQ(ierr);} 28 PetscFree(ctx); 29 return 0; 30 } 31 32 #undef __FUNC__ 33 #define __FUNC__ "SNESMatrixFreeView_Private" /* ADIC Ignore */ 34 /* 35 SNESMatrixFreeView_Private - Views matrix-free parameters. 36 */ 37 int SNESMatrixFreeView_Private(Mat J,Viewer viewer) 38 { 39 int ierr; 40 MFCtx_Private *ctx; 41 MPI_Comm comm; 42 FILE *fd; 43 ViewerType vtype; 44 45 PetscObjectGetComm((PetscObject)J,&comm); 46 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 47 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 48 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 49 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 50 PetscFPrintf(comm,fd," SNES matrix-free approximation:\n"); 51 PetscFPrintf(comm,fd," err=%g (relative error in function evaluation)\n",ctx->error_rel); 52 PetscFPrintf(comm,fd," umin=%g (minimum iterate parameter)\n",ctx->umin); 53 } 54 return 0; 55 } 56 57 extern int VecDot_Seq(Vec,Vec,Scalar *); 58 extern int VecNorm_Seq(Vec,NormType,double *); 59 60 #undef __FUNC__ 61 #define __FUNC__ "SNESMatrixFreeMult_Private" 62 /* 63 SNESMatrixFreeMult_Private - Default matrix-free form for Jacobian-vector 64 product, y = F'(u)*a: 65 y = ( F(u + ha) - F(u) ) /h, 66 where F = nonlinear function, as set by SNESSetFunction() 67 u = current iterate 68 h = difference interval 69 */ 70 int SNESMatrixFreeMult_Private(Mat mat,Vec a,Vec y) 71 { 72 MFCtx_Private *ctx; 73 SNES snes; 74 double ovalues[3],values[3],norm, sum, umin; 75 Scalar h, dot, mone = -1.0; 76 Vec w,U,F; 77 int ierr, (*eval_fct)(SNES,Vec,Vec); 78 MPI_Comm comm; 79 80 PLogEventBegin(MAT_MatrixFreeMult,a,y,0,0); 81 82 PetscObjectGetComm((PetscObject)mat,&comm); 83 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 84 snes = ctx->snes; 85 w = ctx->w; 86 umin = ctx->umin; 87 88 /* We log matrix-free matrix-vector products separately, so that we can 89 separate the performance monitoring from the cases that use conventional 90 storage. We may eventually modify event logging to associate events 91 with particular objects, hence alleviating the more general problem. */ 92 93 ierr = SNESGetSolution(snes,&U); CHKERRQ(ierr); 94 if (snes->method_class == SNES_NONLINEAR_EQUATIONS) { 95 eval_fct = SNESComputeFunction; 96 ierr = SNESGetFunction(snes,&F); CHKERRQ(ierr); 97 } 98 else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) { 99 eval_fct = SNESComputeGradient; 100 ierr = SNESGetGradient(snes,&F); CHKERRQ(ierr); 101 } 102 else SETERRQ(1,0,"Invalid method class"); 103 104 /* Determine a "good" step size, h */ 105 106 /* 107 ierr = VecDot(U,a,&dot); CHKERRQ(ierr); 108 ierr = VecNorm(a,NORM_1,&sum); CHKERRQ(ierr); 109 ierr = VecNorm(a,NORM_2,&norm); CHKERRQ(ierr); 110 */ 111 112 /* 113 Call the Seq Vector routines and then do a single reduction 114 to reduce the number of communications required 115 */ 116 117 #if !defined(PETSC_COMPLEX) 118 PLogEventBegin(VEC_Dot,U,a,0,0); 119 ierr = VecDot_Seq(U,a,ovalues); CHKERRQ(ierr); 120 PLogEventEnd(VEC_Dot,U,a,0,0); 121 PLogEventBegin(VEC_Norm,a,0,0,0); 122 ierr = VecNorm_Seq(a,NORM_1,ovalues+1); CHKERRQ(ierr); 123 ierr = VecNorm_Seq(a,NORM_2,ovalues+2); CHKERRQ(ierr); 124 ovalues[2] = ovalues[2]*ovalues[2]; 125 MPI_Allreduce(ovalues,values,3,MPI_DOUBLE,MPI_SUM,comm ); 126 dot = values[0]; sum = values[1]; norm = sqrt(values[2]); 127 PLogEventEnd(VEC_Norm,a,0,0,0); 128 #else 129 { 130 Scalar cvalues[3],covalues[3]; 131 132 PLogEventBegin(VEC_Dot,U,a,0,0); 133 ierr = VecDot_Seq(U,a,covalues); CHKERRQ(ierr); 134 PLogEventEnd(VEC_Dot,U,a,0,0); 135 PLogEventBegin(VEC_Norm,a,0,0,0); 136 ierr = VecNorm_Seq(a,NORM_1,ovalues+1); CHKERRQ(ierr); 137 ierr = VecNorm_Seq(a,NORM_2,ovalues+2); CHKERRQ(ierr); 138 covalues[1] = ovalues[1]; 139 covalues[2] = ovalues[2]*ovalues[2]; 140 MPI_Allreduce(covalues,cvalues,6,MPI_DOUBLE,MPI_SUM,comm ); 141 dot = cvalues[0]; sum = PetscReal(cvalues[1]); norm = sqrt(PetscReal(cvalues[2])); 142 PLogEventEnd(VEC_Norm,a,0,0,0); 143 } 144 #endif 145 146 147 /* Safeguard for step sizes too small */ 148 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 149 #if defined(PETSC_COMPLEX) 150 else if (abs(dot) < umin*sum && real(dot) >= 0.0) dot = umin*sum; 151 else if (abs(dot) < 0.0 && real(dot) > -umin*sum) dot = -umin*sum; 152 #else 153 else if (dot < umin*sum && dot >= 0.0) dot = umin*sum; 154 else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum; 155 #endif 156 h = ctx->error_rel*dot/(norm*norm); 157 158 /* Evaluate function at F(u + ha) */ 159 ierr = VecWAXPY(&h,a,U,w); CHKERRQ(ierr); 160 ierr = eval_fct(snes,w,y); CHKERRQ(ierr); 161 ierr = VecAXPY(&mone,F,y); CHKERRQ(ierr); 162 h = 1.0/h; 163 ierr = VecScale(&h,y); CHKERRQ(ierr); 164 if (ctx->sp) {ierr = PCNullSpaceRemove(ctx->sp,y); CHKERRQ(ierr);} 165 166 PLogEventEnd(MAT_MatrixFreeMult,a,y,0,0); 167 return 0; 168 } 169 170 #undef __FUNC__ 171 #define __FUNC__ "SNESDefaultMatrixFreeMatCreate" 172 /*@C 173 SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 174 context for use with a SNES solver. This matrix can be used as 175 the Jacobian argument for the routine SNESSetJacobian(). 176 177 Input Parameters: 178 . snes - the SNES context 179 . x - vector where SNES solution is to be stored. 180 181 Output Parameter: 182 . J - the matrix-free matrix 183 184 Notes: 185 The matrix-free matrix context merely contains the function pointers 186 and work space for performing finite difference approximations of 187 Jacobian-vector products, J(u)*a, via 188 189 $ J(u)*a = [J(u+h*a) - J(u)]/h where 190 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 191 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 192 $ where 193 $ error_rel = square root of relative error in 194 $ function evaluation 195 $ umin = minimum iterate parameter 196 197 The user can set these parameters via SNESSetMatrixFreeParameters(). 198 See the nonlinear solvers chapter of the users manual for details. 199 200 The user should call MatDestroy() when finished with the matrix-free 201 matrix context. 202 203 Options Database Keys: 204 $ -snes_mf_err <error_rel> 205 $ -snes_mf_unim <umin> 206 207 .keywords: SNES, default, matrix-free, create, matrix 208 209 .seealso: MatDestroy(), SNESSetMatrixFreeParameters() 210 @*/ 211 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 212 { 213 MPI_Comm comm; 214 MFCtx_Private *mfctx; 215 int n, nloc, ierr, flg; 216 char p[64]; 217 218 mfctx = (MFCtx_Private *) PetscMalloc(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 219 PLogObjectMemory(snes,sizeof(MFCtx_Private)); 220 mfctx->sp = 0; 221 mfctx->snes = snes; 222 mfctx->error_rel = 1.e-8; /* assumes double precision */ 223 mfctx->umin = 1.e-6; 224 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_err",&mfctx->error_rel,&flg); CHKERRQ(ierr); 225 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_umin",&mfctx->umin,&flg); CHKERRQ(ierr); 226 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 227 PetscStrcpy(p,"-"); 228 if (snes->prefix) PetscStrcat(p,snes->prefix); 229 if (flg) { 230 PetscPrintf(snes->comm," %ssnes_mf_err <err>: set sqrt rel error in function (default %g)\n",p,mfctx->error_rel); 231 PetscPrintf(snes->comm," %ssnes_mf_umin <umin> see users manual (default %g)\n",p,mfctx->umin); 232 } 233 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 234 ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr); 235 ierr = VecGetSize(x,&n); CHKERRQ(ierr); 236 ierr = VecGetLocalSize(x,&nloc); CHKERRQ(ierr); 237 ierr = MatCreateShell(comm,nloc,n,n,n,mfctx,J); CHKERRQ(ierr); 238 ierr = MatShellSetOperation(*J,MATOP_MULT,(void*)SNESMatrixFreeMult_Private);CHKERRQ(ierr); 239 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void *)SNESMatrixFreeDestroy_Private);CHKERRQ(ierr); 240 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void *)SNESMatrixFreeView_Private); CHKERRQ(ierr); 241 PLogObjectParent(*J,mfctx->w); 242 PLogObjectParent(snes,*J); 243 return 0; 244 } 245 246 #undef __FUNC__ 247 #define __FUNC__ "SNESSetMatrixFreeParameters" 248 /*@ 249 SNESSetMatrixFreeParameters - Sets the parameters for the approximation of 250 matrix-vector products using finite differences. 251 252 $ J(u)*a = [J(u+h*a) - J(u)]/h where 253 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 254 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 255 $ 256 Input Parameters: 257 . snes - the SNES context 258 . error_rel - relative error (should be set to the square root of 259 the relative error in the function evaluations) 260 . umin - minimum allowable u-value 261 262 .keywords: SNES, matrix-free, parameters 263 264 .seealso: SNESDefaultMatrixFreeMatCreate() 265 @*/ 266 int SNESSetMatrixFreeParameters(SNES snes,double error,double umin) 267 { 268 MFCtx_Private *ctx; 269 int ierr; 270 Mat mat; 271 272 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 273 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 274 if (ctx) { 275 if (error != PETSC_DEFAULT) ctx->error_rel = error; 276 if (umin != PETSC_DEFAULT) ctx->umin = umin; 277 } 278 return 0; 279 } 280 281 #undef __FUNC__ 282 #define __FUNC__ "SNESDefaultMatrixFreeMatAddNullSpace" 283 /*@ 284 SNESDefaultMatrixFreeMatAddNullSpace - Provides a null space that 285 an operator is supposed to have. Since roundoff will create a 286 small component in the null space, if you know the null space 287 you may have it automatically removed. 288 289 Input Parameters: 290 . J - the matrix-free matrix context 291 . has_cnst - PETSC_TRUE or PETSC_FALSE, indicating if null space has constants 292 . n - number of vectors (excluding constant vector) in null space 293 . vecs - the vectors that span the null space (excluding the constant vector); 294 . these vectors must be orthonormal 295 296 .keywords: SNES, matrix-free, null space 297 @*/ 298 int SNESDefaultMatrixFreeMatAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs) 299 { 300 int ierr; 301 MFCtx_Private *ctx; 302 MPI_Comm comm; 303 304 PetscObjectGetComm((PetscObject)J,&comm); 305 306 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 307 /* no context indicates that it is not the "matrix free" matrix type */ 308 if (!ctx) return 0; 309 ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr); 310 return 0; 311 } 312 313 314 315 316