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