1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: snesmfj.c,v 1.54 1997/07/09 20:59:37 balay Exp bsmith $"; 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 PLogEventBarrierBegin(VEC_NormBarrier,0,0,0,0,comm); 126 MPI_Allreduce(ovalues,values,3,MPI_DOUBLE,MPI_SUM,comm ); 127 PLogEventBarrierEnd(VEC_NormBarrier,0,0,0,0,comm); 128 dot = values[0]; sum = values[1]; norm = sqrt(values[2]); 129 PLogEventEnd(VEC_Norm,a,0,0,0); 130 #else 131 { 132 Scalar cvalues[3],covalues[3]; 133 134 PLogEventBegin(VEC_Dot,U,a,0,0); 135 ierr = VecDot_Seq(U,a,covalues); CHKERRQ(ierr); 136 PLogEventEnd(VEC_Dot,U,a,0,0); 137 PLogEventBegin(VEC_Norm,a,0,0,0); 138 ierr = VecNorm_Seq(a,NORM_1,ovalues+1); CHKERRQ(ierr); 139 ierr = VecNorm_Seq(a,NORM_2,ovalues+2); CHKERRQ(ierr); 140 covalues[1] = ovalues[1]; 141 covalues[2] = ovalues[2]*ovalues[2]; 142 PLogEventBarrierBegin(VEC_NormBarrier,0,0,0,0,comm); 143 MPI_Allreduce(covalues,cvalues,6,MPI_DOUBLE,MPI_SUM,comm ); 144 PLogEventBarrierBegin(VEC_NormBarrier,0,0,0,0,comm); 145 dot = cvalues[0]; sum = PetscReal(cvalues[1]); norm = sqrt(PetscReal(cvalues[2])); 146 PLogEventEnd(VEC_Norm,a,0,0,0); 147 } 148 #endif 149 150 151 /* Safeguard for step sizes too small */ 152 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 153 #if defined(PETSC_COMPLEX) 154 else if (abs(dot) < umin*sum && real(dot) >= 0.0) dot = umin*sum; 155 else if (abs(dot) < 0.0 && real(dot) > -umin*sum) dot = -umin*sum; 156 #else 157 else if (dot < umin*sum && dot >= 0.0) dot = umin*sum; 158 else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum; 159 #endif 160 h = ctx->error_rel*dot/(norm*norm); 161 162 /* Evaluate function at F(u + ha) */ 163 ierr = VecWAXPY(&h,a,U,w); CHKERRQ(ierr); 164 ierr = eval_fct(snes,w,y); CHKERRQ(ierr); 165 ierr = VecAXPY(&mone,F,y); CHKERRQ(ierr); 166 h = 1.0/h; 167 ierr = VecScale(&h,y); CHKERRQ(ierr); 168 if (ctx->sp) {ierr = PCNullSpaceRemove(ctx->sp,y); CHKERRQ(ierr);} 169 170 PLogEventEnd(MAT_MatrixFreeMult,a,y,0,0); 171 return 0; 172 } 173 174 #undef __FUNC__ 175 #define __FUNC__ "SNESDefaultMatrixFreeMatCreate" 176 /*@C 177 SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 178 context for use with a SNES solver. This matrix can be used as 179 the Jacobian argument for the routine SNESSetJacobian(). 180 181 Input Parameters: 182 . snes - the SNES context 183 . x - vector where SNES solution is to be stored. 184 185 Output Parameter: 186 . J - the matrix-free matrix 187 188 Notes: 189 The matrix-free matrix context merely contains the function pointers 190 and work space for performing finite difference approximations of 191 Jacobian-vector products, J(u)*a, via 192 193 $ J(u)*a = [J(u+h*a) - J(u)]/h where 194 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 195 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 196 $ where 197 $ error_rel = square root of relative error in 198 $ function evaluation 199 $ umin = minimum iterate parameter 200 201 The user can set these parameters via SNESSetMatrixFreeParameters(). 202 See the nonlinear solvers chapter of the users manual for details. 203 204 The user should call MatDestroy() when finished with the matrix-free 205 matrix context. 206 207 Options Database Keys: 208 $ -snes_mf_err <error_rel> 209 $ -snes_mf_unim <umin> 210 211 .keywords: SNES, default, matrix-free, create, matrix 212 213 .seealso: MatDestroy(), SNESSetMatrixFreeParameters() 214 @*/ 215 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 216 { 217 MPI_Comm comm; 218 MFCtx_Private *mfctx; 219 int n, nloc, ierr, flg; 220 char p[64]; 221 222 mfctx = (MFCtx_Private *) PetscMalloc(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 223 PLogObjectMemory(snes,sizeof(MFCtx_Private)); 224 mfctx->sp = 0; 225 mfctx->snes = snes; 226 mfctx->error_rel = 1.e-8; /* assumes double precision */ 227 mfctx->umin = 1.e-6; 228 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_err",&mfctx->error_rel,&flg); CHKERRQ(ierr); 229 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_umin",&mfctx->umin,&flg); CHKERRQ(ierr); 230 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 231 PetscStrcpy(p,"-"); 232 if (snes->prefix) PetscStrcat(p,snes->prefix); 233 if (flg) { 234 PetscPrintf(snes->comm," %ssnes_mf_err <err>: set sqrt rel error in function (default %g)\n",p,mfctx->error_rel); 235 PetscPrintf(snes->comm," %ssnes_mf_umin <umin> see users manual (default %g)\n",p,mfctx->umin); 236 } 237 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 238 ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr); 239 ierr = VecGetSize(x,&n); CHKERRQ(ierr); 240 ierr = VecGetLocalSize(x,&nloc); CHKERRQ(ierr); 241 ierr = MatCreateShell(comm,nloc,n,n,n,mfctx,J); CHKERRQ(ierr); 242 ierr = MatShellSetOperation(*J,MATOP_MULT,(void*)SNESMatrixFreeMult_Private);CHKERRQ(ierr); 243 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void *)SNESMatrixFreeDestroy_Private);CHKERRQ(ierr); 244 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void *)SNESMatrixFreeView_Private); CHKERRQ(ierr); 245 PLogObjectParent(*J,mfctx->w); 246 PLogObjectParent(snes,*J); 247 return 0; 248 } 249 250 #undef __FUNC__ 251 #define __FUNC__ "SNESSetMatrixFreeParameters" 252 /*@ 253 SNESSetMatrixFreeParameters - Sets the parameters for the approximation of 254 matrix-vector products using finite differences. 255 256 $ J(u)*a = [J(u+h*a) - J(u)]/h where 257 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 258 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 259 $ 260 Input Parameters: 261 . snes - the SNES context 262 . error_rel - relative error (should be set to the square root of 263 the relative error in the function evaluations) 264 . umin - minimum allowable u-value 265 266 .keywords: SNES, matrix-free, parameters 267 268 .seealso: SNESDefaultMatrixFreeMatCreate() 269 @*/ 270 int SNESSetMatrixFreeParameters(SNES snes,double error,double umin) 271 { 272 MFCtx_Private *ctx; 273 int ierr; 274 Mat mat; 275 276 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 277 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 278 if (ctx) { 279 if (error != PETSC_DEFAULT) ctx->error_rel = error; 280 if (umin != PETSC_DEFAULT) ctx->umin = umin; 281 } 282 return 0; 283 } 284 285 #undef __FUNC__ 286 #define __FUNC__ "SNESDefaultMatrixFreeMatAddNullSpace" 287 /*@ 288 SNESDefaultMatrixFreeMatAddNullSpace - Provides a null space that 289 an operator is supposed to have. Since roundoff will create a 290 small component in the null space, if you know the null space 291 you may have it automatically removed. 292 293 Input Parameters: 294 . J - the matrix-free matrix context 295 . has_cnst - PETSC_TRUE or PETSC_FALSE, indicating if null space has constants 296 . n - number of vectors (excluding constant vector) in null space 297 . vecs - the vectors that span the null space (excluding the constant vector); 298 . these vectors must be orthonormal 299 300 .keywords: SNES, matrix-free, null space 301 @*/ 302 int SNESDefaultMatrixFreeMatAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs) 303 { 304 int ierr; 305 MFCtx_Private *ctx; 306 MPI_Comm comm; 307 308 PetscObjectGetComm((PetscObject)J,&comm); 309 310 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 311 /* no context indicates that it is not the "matrix free" matrix type */ 312 if (!ctx) return 0; 313 ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr); 314 return 0; 315 } 316 317 318 319 320