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