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