1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.43 1997/01/21 20:12:37 curfman Exp curfman $"; 4 #endif 5 6 #include "draw.h" /*I "draw.h" I*/ 7 #include "src/snes/snesimpl.h" /*I "snes.h" I*/ 8 #include "pinclude/pviewer.h" 9 10 typedef struct { /* default context for matrix-free SNES */ 11 SNES snes; /* SNES context */ 12 Vec w; /* work vector */ 13 PCNullSpace sp; /* null space context */ 14 double error_rel; /* square root of relative error in computing function */ 15 double umin; /* minimum allowable u'a value relative to |u|_1 */ 16 } MFCtx_Private; 17 18 #undef __FUNC__ 19 #define __FUNC__ "SNESMatrixFreeDestroy_Private" 20 int SNESMatrixFreeDestroy_Private(PetscObject obj) 21 { 22 int ierr; 23 Mat mat = (Mat) obj; 24 MFCtx_Private *ctx; 25 26 ierr = MatShellGetContext(mat,(void **)&ctx); 27 ierr = VecDestroy(ctx->w); CHKERRQ(ierr); 28 if (ctx->sp) {ierr = PCNullSpaceDestroy(ctx->sp); CHKERRQ(ierr);} 29 PetscFree(ctx); 30 return 0; 31 } 32 33 #undef __FUNC__ 34 #define __FUNC__ "SNESMatrixFreeView_Private" 35 /* 36 SNESMatrixFreeView_Private - Views matrix-free parameters. 37 */ 38 int SNESMatrixFreeView_Private(Mat J,Viewer viewer) 39 { 40 int ierr; 41 MFCtx_Private *ctx; 42 MPI_Comm comm; 43 FILE *fd; 44 ViewerType vtype; 45 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 } 55 return 0; 56 } 57 58 #undef __FUNC__ 59 #define __FUNC__ "SNESMatrixFreeMult_Private" 60 /* 61 SNESMatrixFreeMult_Private - Default matrix-free form for Jacobian-vector 62 product, y = F'(u)*a: 63 y = ( F(u + ha) - F(u) ) /h, 64 where F = nonlinear function, as set by SNESSetFunction() 65 u = current iterate 66 h = difference interval 67 */ 68 int SNESMatrixFreeMult_Private(Mat mat,Vec a,Vec y) 69 { 70 MFCtx_Private *ctx; 71 SNES snes; 72 double norm, sum, umin; 73 Scalar h, dot, mone = -1.0; 74 Vec w,U,F; 75 int ierr, (*eval_fct)(SNES,Vec,Vec); 76 77 MatShellGetContext(mat,(void **)&ctx); 78 snes = ctx->snes; 79 w = ctx->w; 80 umin = ctx->umin; 81 82 /* We log matrix-free matrix-vector products separately, so that we can 83 separate the performance monitoring from the cases that use conventional 84 storage. We may eventually modify event logging to associate events 85 with particular objects, hence alleviating the more general problem. */ 86 PLogEventBegin(MAT_MatrixFreeMult,a,y,0,0); 87 88 ierr = SNESGetSolution(snes,&U); CHKERRQ(ierr); 89 if (snes->method_class == SNES_NONLINEAR_EQUATIONS) { 90 eval_fct = SNESComputeFunction; 91 ierr = SNESGetFunction(snes,&F); CHKERRQ(ierr); 92 } 93 else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) { 94 eval_fct = SNESComputeGradient; 95 ierr = SNESGetGradient(snes,&F); CHKERRQ(ierr); 96 } 97 else SETERRQ(1,0,"Invalid method class"); 98 99 /* Determine a "good" step size, h */ 100 ierr = VecDot(U,a,&dot); CHKERRQ(ierr); 101 ierr = VecNorm(a,NORM_1,&sum); CHKERRQ(ierr); 102 ierr = VecNorm(a,NORM_2,&norm); CHKERRQ(ierr); 103 104 /* Safeguard for step sizes too small */ 105 if (sum == 0.0) {dot = 1.0; norm = 1.0;} 106 #if defined(PETSC_COMPLEX) 107 else if (abs(dot) < umin*sum && real(dot) >= 0.0) dot = umin*sum; 108 else if (abs(dot) < 0.0 && real(dot) > -umin*sum) dot = -umin*sum; 109 #else 110 else if (dot < umin*sum && dot >= 0.0) dot = umin*sum; 111 else if (dot < 0.0 && dot > -umin*sum) dot = -umin*sum; 112 #endif 113 h = ctx->error_rel*dot/(norm*norm); 114 115 /* Evaluate function at F(u + ha) */ 116 ierr = VecWAXPY(&h,a,U,w); CHKERRQ(ierr); 117 ierr = eval_fct(snes,w,y); CHKERRQ(ierr); 118 ierr = VecAXPY(&mone,F,y); CHKERRQ(ierr); 119 h = 1.0/h; 120 ierr = VecScale(&h,y); CHKERRQ(ierr); 121 if (ctx->sp) {ierr = PCNullSpaceRemove(ctx->sp,y); CHKERRQ(ierr);} 122 123 PLogEventEnd(MAT_MatrixFreeMult,a,y,0,0); 124 return 0; 125 } 126 127 #undef __FUNC__ 128 #define __FUNC__ "SNESDefaultMatrixFreeMatCreate" 129 /*@C 130 SNESDefaultMatrixFreeMatCreate - Creates a matrix-free matrix 131 context for use with a SNES solver. This matrix can be used as 132 the Jacobian argument for the routine SNESSetJacobian(). 133 134 Input Parameters: 135 . snes - the SNES context 136 . x - vector where SNES solution is to be stored. 137 138 Output Parameter: 139 . J - the matrix-free matrix 140 141 Notes: 142 The matrix-free matrix context merely contains the function pointers 143 and work space for performing finite difference approximations of 144 matrix operations such as matrix-vector products. 145 146 The user should call MatDestroy() when finished with the matrix-free 147 matrix context. 148 149 .keywords: SNES, default, matrix-free, create, matrix 150 151 .seealso: MatDestroy() 152 @*/ 153 int SNESDefaultMatrixFreeMatCreate(SNES snes,Vec x, Mat *J) 154 { 155 MPI_Comm comm; 156 MFCtx_Private *mfctx; 157 int n, nloc, ierr, flg; 158 char p[64]; 159 160 mfctx = (MFCtx_Private *) PetscMalloc(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 161 PLogObjectMemory(snes,sizeof(MFCtx_Private)); 162 mfctx->sp = 0; 163 mfctx->snes = snes; 164 mfctx->error_rel = 1.e-8; /* assumes double precision */ 165 mfctx->umin = 1.e-6; 166 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_err",&mfctx->error_rel,&flg); CHKERRQ(ierr); 167 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_umin",&mfctx->umin,&flg); CHKERRQ(ierr); 168 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 169 PetscStrcpy(p,"-"); 170 if (snes->prefix) PetscStrcat(p,snes->prefix); 171 if (flg) { 172 PetscPrintf(snes->comm," %ssnes_mf_err <err>: set sqrt rel error in function (default %g)\n",p,mfctx->error_rel); 173 PetscPrintf(snes->comm," %ssnes_mf_umin <umin> see users manual (default %g)\n",p,mfctx->umin); 174 } 175 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 176 ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr); 177 ierr = VecGetSize(x,&n); CHKERRQ(ierr); 178 ierr = VecGetLocalSize(x,&nloc); CHKERRQ(ierr); 179 ierr = MatCreateShell(comm,nloc,n,n,n,(void*)mfctx,J); CHKERRQ(ierr); 180 ierr = MatShellSetOperation(*J,MATOP_MULT,(void*)SNESMatrixFreeMult_Private); CHKERRQ(ierr); 181 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void *)SNESMatrixFreeDestroy_Private); CHKERRQ(ierr); 182 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void *)SNESMatrixFreeView_Private); CHKERRQ(ierr); 183 PLogObjectParent(*J,mfctx->w); 184 PLogObjectParent(snes,*J); 185 return 0; 186 } 187 188 #undef __FUNC__ 189 #define __FUNC__ "SNESSetMatrixFreeParameters" 190 /*@ 191 SNESSetMatrixFreeParameters - Sets the parameters for the approximation of 192 matrix-vector products using finite differences. 193 194 $ J(u)*a = [J(u+h*a) - J(u)]/h where 195 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 196 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 197 $ 198 Input Parameters: 199 . snes - the SNES context 200 . error_rel - relative error (should be set to the square root of 201 the relative error in the function evaluations) 202 . umin - minimum allowable u-value 203 204 .keywords: SNES, matrix-free, parameters 205 @*/ 206 int SNESSetMatrixFreeParameters(SNES snes,double error,double umin) 207 { 208 MFCtx_Private *ctx; 209 int ierr; 210 Mat mat; 211 212 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 213 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 214 if (ctx) { 215 if (error != PETSC_DEFAULT) ctx->error_rel = error; 216 if (umin != PETSC_DEFAULT) ctx->umin = umin; 217 } 218 return 0; 219 } 220 221 #undef __FUNC__ 222 #define __FUNC__ "SNESDefaultMatrixFreeMatAddNullSpace" 223 /*@ 224 SNESDefaultMatrixFreeMatAddNullSpace - Provides a null space that 225 an operator is supposed to have. Since roundoff will create a 226 small component in the null space, if you know the null space 227 you may have it automatically removed. 228 229 Input Parameters: 230 . J - the matrix-free matrix context 231 . has_cnst - PETSC_TRUE or PETSC_FALSE, indicating if null space has constants 232 . n - number of vectors (excluding constant vector) in null space 233 . vecs - the vectors that span the null space (excluding the constant vector); 234 . these vectors must be orthonormal 235 236 .keywords: SNES, matrix-free, null space 237 @*/ 238 int SNESDefaultMatrixFreeMatAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs) 239 { 240 int ierr; 241 MFCtx_Private *ctx; 242 MPI_Comm comm; 243 244 PetscObjectGetComm((PetscObject)J,&comm); 245 246 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 247 /* no context indicates that it is not the "matrix free" matrix type */ 248 if (!ctx) return 0; 249 ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr); 250 return 0; 251 } 252 253 254 255 256