1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.42 1997/01/06 20:29:45 balay 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 159 mfctx = (MFCtx_Private *) PetscMalloc(sizeof(MFCtx_Private)); CHKPTRQ(mfctx); 160 PLogObjectMemory(snes,sizeof(MFCtx_Private)); 161 mfctx->sp = 0; 162 mfctx->snes = snes; 163 mfctx->error_rel = 1.e-8; /* assumes double precision */ 164 mfctx->umin = 1.e-6; 165 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_err",&mfctx->error_rel,&flg); CHKERRQ(ierr); 166 ierr = OptionsGetDouble(snes->prefix,"-snes_mf_umin",&mfctx->umin,&flg); CHKERRQ(ierr); 167 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 168 if (flg) { 169 PetscPrintf(snes->comm,"-snes_mf_err <err> set sqrt rel tol in function. Default 1.e-8\n"); 170 PetscPrintf(snes->comm,"-snes_mf_umin <umin> see users manual. Default 1.e-8\n"); 171 } 172 ierr = VecDuplicate(x,&mfctx->w); CHKERRQ(ierr); 173 ierr = PetscObjectGetComm((PetscObject)x,&comm); CHKERRQ(ierr); 174 ierr = VecGetSize(x,&n); CHKERRQ(ierr); 175 ierr = VecGetLocalSize(x,&nloc); CHKERRQ(ierr); 176 ierr = MatCreateShell(comm,nloc,n,n,n,(void*)mfctx,J); CHKERRQ(ierr); 177 ierr = MatShellSetOperation(*J,MATOP_MULT,(void*)SNESMatrixFreeMult_Private); CHKERRQ(ierr); 178 ierr = MatShellSetOperation(*J,MATOP_DESTROY,(void *)SNESMatrixFreeDestroy_Private); CHKERRQ(ierr); 179 ierr = MatShellSetOperation(*J,MATOP_VIEW,(void *)SNESMatrixFreeView_Private); CHKERRQ(ierr); 180 PLogObjectParent(*J,mfctx->w); 181 PLogObjectParent(snes,*J); 182 return 0; 183 } 184 185 #undef __FUNC__ 186 #define __FUNC__ "SNESSetMatrixFreeParameters" 187 /*@ 188 SNESSetMatrixFreeParameters - Sets the parameters for the approximation of 189 matrix-vector products using finite differences. 190 191 $ J(u)*a = [J(u+h*a) - J(u)]/h where 192 $ h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 193 $ = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 else 194 $ 195 Input Parameters: 196 . snes - the SNES context 197 . error_rel - relative error (should be set to the square root of 198 the relative error in the function evaluations) 199 . umin - minimum allowable u-value 200 201 .keywords: SNES, matrix-free, parameters 202 @*/ 203 int SNESSetMatrixFreeParameters(SNES snes,double error,double umin) 204 { 205 MFCtx_Private *ctx; 206 int ierr; 207 Mat mat; 208 209 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 210 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 211 if (ctx) { 212 if (error != PETSC_DEFAULT) ctx->error_rel = error; 213 if (umin != PETSC_DEFAULT) ctx->umin = umin; 214 } 215 return 0; 216 } 217 218 #undef __FUNC__ 219 #define __FUNC__ "SNESDefaultMatrixFreeMatAddNullSpace" 220 /*@ 221 SNESDefaultMatrixFreeMatAddNullSpace - Provides a null space that 222 an operator is supposed to have. Since roundoff will create a 223 small component in the null space, if you know the null space 224 you may have it automatically removed. 225 226 Input Parameters: 227 . J - the matrix-free matrix context 228 . has_cnst - PETSC_TRUE or PETSC_FALSE, indicating if null space has constants 229 . n - number of vectors (excluding constant vector) in null space 230 . vecs - the vectors that span the null space (excluding the constant vector); 231 . these vectors must be orthonormal 232 233 .keywords: SNES, matrix-free, null space 234 @*/ 235 int SNESDefaultMatrixFreeMatAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs) 236 { 237 int ierr; 238 MFCtx_Private *ctx; 239 MPI_Comm comm; 240 241 PetscObjectGetComm((PetscObject)J,&comm); 242 243 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 244 /* no context indicates that it is not the "matrix free" matrix type */ 245 if (!ctx) return 0; 246 ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr); 247 return 0; 248 } 249 250 251 252 253