1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: snesmfj.c,v 1.62 1998/04/03 23:17:44 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 } 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); 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 (abs(dot) < umin*sum && real(dot) >= 0.0) dot = umin*sum; 162 else if (abs(dot) < 0.0 && real(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 ierr = VecAXPY(&mone,F,y); CHKERRQ(ierr); 173 h = 1.0/h; 174 ierr = VecScale(&h,y); CHKERRQ(ierr); 175 if (ctx->sp) {ierr = PCNullSpaceRemove(ctx->sp,y); CHKERRQ(ierr);} 176 177 PLogEventEnd(MAT_MatrixFreeMult,a,y,0,0); 178 PetscFunctionReturn(0); 179 } 180 181 #undef __FUNC__ 182 #define __FUNC__ "SNESDefaultMatrixFreeMatCreate" 183 /*@C 184 SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 185 context for use with a SNES solver. This matrix can be used as 186 the Jacobian argument for the routine SNESSetJacobian(). 187 188 Input Parameters: 189 . snes - the SNES context 190 . x - vector where SNES solution is to be stored. 191 192 Output Parameter: 193 . J - the matrix-free matrix 194 195 Notes: 196 The matrix-free matrix context merely contains the function pointers 197 and work space for performing finite difference approximations of 198 Jacobian-vector products, J(u)*a, via 199 200 $ J(u)*a = [J(u+h*a) - J(u)]/h where 201 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 202 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 203 $ where 204 $ error_rel = square root of relative error in 205 $ function evaluation 206 $ umin = minimum iterate parameter 207 208 The user can set these parameters via SNESSetMatrixFreeParameters(). 209 See the nonlinear solvers chapter of the users manual for details. 210 211 The user should call MatDestroy() when finished with the matrix-free 212 matrix context. 213 214 Options Database Keys: 215 $ -snes_mf_err <error_rel> 216 $ -snes_mf_unim <umin> 217 218 .keywords: SNES, default, matrix-free, create, matrix 219 220 .seealso: MatDestroy(), SNESSetMatrixFreeParameters() 221 @*/ 222 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 223 { 224 MPI_Comm comm; 225 MFCtx_Private *mfctx; 226 int n, nloc, ierr, flg; 227 char p[64]; 228 229 PetscFunctionBegin; 230 mfctx = (MFCtx_Private *) PetscMalloc(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 231 PLogObjectMemory(snes,sizeof(MFCtx_Private)); 232 mfctx->sp = 0; 233 mfctx->snes = snes; 234 mfctx->error_rel = 1.e-8; /* assumes double precision */ 235 mfctx->umin = 1.e-6; 236 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_err",&mfctx->error_rel,&flg); CHKERRQ(ierr); 237 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_umin",&mfctx->umin,&flg); CHKERRQ(ierr); 238 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 239 PetscStrcpy(p,"-"); 240 if (snes->prefix) PetscStrcat(p,snes->prefix); 241 if (flg) { 242 (*PetscHelpPrintf)(snes->comm," %ssnes_mf_err <err>: set sqrt rel error in function (default %g)\n",p,mfctx->error_rel); 243 (*PetscHelpPrintf)(snes->comm," %ssnes_mf_umin <umin> see users manual (default %g)\n",p,mfctx->umin); 244 } 245 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 246 ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr); 247 ierr = VecGetSize(x,&n); CHKERRQ(ierr); 248 ierr = VecGetLocalSize(x,&nloc); CHKERRQ(ierr); 249 ierr = MatCreateShell(comm,nloc,n,n,n,mfctx,J); CHKERRQ(ierr); 250 ierr = MatShellSetOperation(*J,MATOP_MULT,(void*)SNESMatrixFreeMult_Private);CHKERRQ(ierr); 251 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void *)SNESMatrixFreeDestroy_Private);CHKERRQ(ierr); 252 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void *)SNESMatrixFreeView_Private); CHKERRQ(ierr); 253 PLogObjectParent(*J,mfctx->w); 254 PLogObjectParent(snes,*J); 255 PetscFunctionReturn(0); 256 } 257 258 #undef __FUNC__ 259 #define __FUNC__ "SNESSetMatrixFreeParameters" 260 /*@ 261 SNESSetMatrixFreeParameters - Sets the parameters for the approximation of 262 matrix-vector products using finite differences. 263 264 $ J(u)*a = [J(u+h*a) - J(u)]/h where 265 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 266 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 267 $ 268 Input Parameters: 269 . snes - the SNES context 270 . error_rel - relative error (should be set to the square root of 271 the relative error in the function evaluations) 272 . umin - minimum allowable u-value 273 274 .keywords: SNES, matrix-free, parameters 275 276 .seealso: SNESDefaultMatrixFreeMatCreate() 277 @*/ 278 int SNESSetMatrixFreeParameters(SNES snes,double error,double umin) 279 { 280 MFCtx_Private *ctx; 281 int ierr; 282 Mat mat; 283 284 PetscFunctionBegin; 285 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 286 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 287 if (ctx) { 288 if (error != PETSC_DEFAULT) ctx->error_rel = error; 289 if (umin != PETSC_DEFAULT) ctx->umin = umin; 290 } 291 PetscFunctionReturn(0); 292 } 293 294 #undef __FUNC__ 295 #define __FUNC__ "SNESDefaultMatrixFreeMatAddNullSpace" 296 /*@ 297 SNESDefaultMatrixFreeMatAddNullSpace - Provides a null space that 298 an operator is supposed to have. Since roundoff will create a 299 small component in the null space, if you know the null space 300 you may have it automatically removed. 301 302 Input Parameters: 303 . J - the matrix-free matrix context 304 . has_cnst - PETSC_TRUE or PETSC_FALSE, indicating if null space has constants 305 . n - number of vectors (excluding constant vector) in null space 306 . vecs - the vectors that span the null space (excluding the constant vector); 307 . these vectors must be orthonormal 308 309 .keywords: SNES, matrix-free, null space 310 @*/ 311 int SNESDefaultMatrixFreeMatAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs) 312 { 313 int ierr; 314 MFCtx_Private *ctx; 315 MPI_Comm comm; 316 317 PetscFunctionBegin; 318 PetscObjectGetComm((PetscObject)J,&comm); 319 320 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 321 /* no context indicates that it is not the "matrix free" matrix type */ 322 if (!ctx) PetscFunctionReturn(0); 323 ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr); 324 PetscFunctionReturn(0); 325 } 326 327 328 329 330