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