1 2 #ifndef lint 3 static char vcid[] = "$Id: snesmfj.c,v 1.39 1996/12/16 20:46:01 balay Exp balay $"; 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 __FUNCTION__ 19 #define __FUNCTION__ "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 __FUNCTION__ 34 #define __FUNCTION__ "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 __FUNCTION__ 59 #define __FUNCTION__ "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,"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 __FUNCTION__ 128 #define __FUNCTION__ "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 __FUNCTION__ 186 #define __FUNCTION__ "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 198 . umin - minimum allowable u-value 199 200 .keywords: SNES, matrix-free, parameters 201 @*/ 202 int SNESSetMatrixFreeParameters(SNES snes,double error,double umin) 203 { 204 MFCtx_Private *ctx; 205 int ierr; 206 Mat mat; 207 208 ierr = SNESGetJacobian(snes,&mat,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); 209 ierr = MatShellGetContext(mat,(void **)&ctx); CHKERRQ(ierr); 210 if (ctx) { 211 if (error != PETSC_DEFAULT) ctx->error_rel = error; 212 if (umin != PETSC_DEFAULT) ctx->umin = umin; 213 } 214 return 0; 215 } 216 217 #undef __FUNCTION__ 218 #define __FUNCTION__ "SNESDefaultMatrixFreeMatAddNullSpace" 219 /*@ 220 SNESDefaultMatrixFreeMatAddNullSpace - Provides a null space that 221 an operator is supposed to have. Since roundoff will create a 222 small component in the null space, if you know the null space 223 you may have it automatically removed. 224 225 Input Parameters: 226 . J - the matrix-free matrix context 227 . has_cnst - PETSC_TRUE or PETSC_FALSE, indicating if null space has constants 228 . n - number of vectors (excluding constant vector) in null space 229 . vecs - the vectors that span the null space (excluding the constant vector); 230 . these vectors must be orthonormal 231 232 .keywords: SNES, matrix-free, null space 233 @*/ 234 int SNESDefaultMatrixFreeMatAddNullSpace(Mat J,int has_cnst,int n,Vec *vecs) 235 { 236 int ierr; 237 MFCtx_Private *ctx; 238 MPI_Comm comm; 239 240 PetscObjectGetComm((PetscObject)J,&comm); 241 242 ierr = MatShellGetContext(J,(void **)&ctx); CHKERRQ(ierr); 243 /* no context indicates that it is not the "matrix free" matrix type */ 244 if (!ctx) return 0; 245 ierr = PCNullSpaceCreate(comm,has_cnst,n,vecs,&ctx->sp); CHKERRQ(ierr); 246 return 0; 247 } 248 249 250 251 252