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