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