1 2 /* 3 Implements the DS PETSc approach for computing the h 4 parameter used with the finite difference based matrix-free 5 Jacobian-vector products. 6 7 To make your own: clone this file and modify for your needs. 8 9 Mandatory functions: 10 ------------------- 11 MatMFFDCompute_ - for a given point and direction computes h 12 13 MatCreateMFFD _ - fills in the MatMFFD data structure 14 for this particular implementation 15 16 17 Optional functions: 18 ------------------- 19 MatMFFDView_ - prints information about the parameters being used. 20 This is called when SNESView() or -snes_view is used. 21 22 MatMFFDSetFromOptions_ - checks the options database for options that 23 apply to this method. 24 25 MatMFFDDestroy_ - frees any space allocated by the routines above 26 27 */ 28 29 /* 30 This include file defines the data structure MatMFFD that 31 includes information about the computation of h. It is shared by 32 all implementations that people provide 33 */ 34 #include <petsc/private/matimpl.h> 35 #include <../src/mat/impls/mffd/mffdimpl.h> /*I "petscmat.h" I*/ 36 37 /* 38 The method has one parameter that is used to 39 "cutoff" very small values. This is stored in a data structure 40 that is only visible to this file. If your method has no parameters 41 it can omit this, if it has several simply reorganize the data structure. 42 The data structure is "hung-off" the MatMFFD data structure in 43 the void *hctx; field. 44 */ 45 typedef struct { 46 PetscReal umin; /* minimum allowable u'a value relative to |u|_1 */ 47 } MatMFFD_DS; 48 49 /* 50 MatMFFDCompute_DS - Standard PETSc code for computing the 51 differencing parameter (h) for use with matrix-free finite differences. 52 53 Input Parameters: 54 + ctx - the matrix free context 55 . U - the location at which you want the Jacobian 56 - a - the direction you want the derivative 57 58 59 Output Parameter: 60 . h - the scale computed 61 62 */ 63 static PetscErrorCode MatMFFDCompute_DS(MatMFFD ctx,Vec U,Vec a,PetscScalar *h,PetscBool *zeroa) 64 { 65 MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx; 66 PetscReal nrm,sum,umin = hctx->umin; 67 PetscScalar dot; 68 PetscErrorCode ierr; 69 70 PetscFunctionBegin; 71 if (!(ctx->count % ctx->recomputeperiod)) { 72 /* 73 This algorithm requires 2 norms and 1 inner product. Rather than 74 use directly the VecNorm() and VecDot() routines (and thus have 75 three separate collective operations, we use the VecxxxBegin/End() routines 76 */ 77 ierr = VecDotBegin(U,a,&dot);CHKERRQ(ierr); 78 ierr = VecNormBegin(a,NORM_1,&sum);CHKERRQ(ierr); 79 ierr = VecNormBegin(a,NORM_2,&nrm);CHKERRQ(ierr); 80 ierr = VecDotEnd(U,a,&dot);CHKERRQ(ierr); 81 ierr = VecNormEnd(a,NORM_1,&sum);CHKERRQ(ierr); 82 ierr = VecNormEnd(a,NORM_2,&nrm);CHKERRQ(ierr); 83 84 if (nrm == 0.0) { 85 *zeroa = PETSC_TRUE; 86 PetscFunctionReturn(0); 87 } 88 *zeroa = PETSC_FALSE; 89 90 /* 91 Safeguard for step sizes that are "too small" 92 */ 93 if (PetscAbsScalar(dot) < umin*sum && PetscRealPart(dot) >= 0.0) dot = umin*sum; 94 else if (PetscAbsScalar(dot) < 0.0 && PetscRealPart(dot) > -umin*sum) dot = -umin*sum; 95 *h = ctx->error_rel*dot/(nrm*nrm); 96 if (PetscIsInfOrNanScalar(*h)) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Differencing parameter is not a number sum = %g dot = %g norm = %g",(double)sum,(double)PetscRealPart(dot),(double)nrm); 97 } else { 98 *h = ctx->currenth; 99 } 100 ctx->count++; 101 PetscFunctionReturn(0); 102 } 103 104 /* 105 MatMFFDView_DS - Prints information about this particular 106 method for computing h. Note that this does not print the general 107 information about the matrix-free method, as such info is printed 108 by the calling routine. 109 110 Input Parameters: 111 + ctx - the matrix free context 112 - viewer - the PETSc viewer 113 */ 114 static PetscErrorCode MatMFFDView_DS(MatMFFD ctx,PetscViewer viewer) 115 { 116 MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx; 117 PetscErrorCode ierr; 118 PetscBool iascii; 119 120 PetscFunctionBegin; 121 /* 122 Currently this only handles the ascii file viewers, others 123 could be added, but for this type of object other viewers 124 make less sense 125 */ 126 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 127 if (iascii) { 128 ierr = PetscViewerASCIIPrintf(viewer," umin=%g (minimum iterate parameter)\n",(double)hctx->umin);CHKERRQ(ierr); 129 } 130 PetscFunctionReturn(0); 131 } 132 133 /* 134 MatMFFDSetFromOptions_DS - Looks in the options database for 135 any options appropriate for this method. 136 137 Input Parameter: 138 . ctx - the matrix free context 139 140 */ 141 static PetscErrorCode MatMFFDSetFromOptions_DS(PetscOptionItems *PetscOptionsObject,MatMFFD ctx) 142 { 143 PetscErrorCode ierr; 144 MatMFFD_DS *hctx = (MatMFFD_DS*)ctx->hctx; 145 146 PetscFunctionBegin; 147 ierr = PetscOptionsHead(PetscOptionsObject,"Finite difference matrix free parameters");CHKERRQ(ierr); 148 ierr = PetscOptionsReal("-mat_mffd_umin","umin","MatMFFDDSSetUmin",hctx->umin,&hctx->umin,NULL);CHKERRQ(ierr); 149 ierr = PetscOptionsTail();CHKERRQ(ierr); 150 PetscFunctionReturn(0); 151 } 152 153 /* 154 MatMFFDDestroy_DS - Frees the space allocated by 155 MatCreateMFFD_DS(). 156 157 Input Parameter: 158 . ctx - the matrix free context 159 160 Notes: 161 Does not free the ctx, that is handled by the calling routine 162 */ 163 static PetscErrorCode MatMFFDDestroy_DS(MatMFFD ctx) 164 { 165 PetscErrorCode ierr; 166 167 PetscFunctionBegin; 168 ierr = PetscFree(ctx->hctx);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 /* 173 The following two routines use the PetscObjectCompose() and PetscObjectQuery() 174 mechanism to allow the user to change the Umin parameter used in this method. 175 */ 176 PetscErrorCode MatMFFDDSSetUmin_DS(Mat mat,PetscReal umin) 177 { 178 MatMFFD ctx=NULL; 179 MatMFFD_DS *hctx; 180 PetscErrorCode ierr; 181 182 PetscFunctionBegin; 183 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 184 if (!ctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatMFFDDSSetUmin() attached to non-shell matrix"); 185 hctx = (MatMFFD_DS*)ctx->hctx; 186 hctx->umin = umin; 187 PetscFunctionReturn(0); 188 } 189 190 /*@ 191 MatMFFDDSSetUmin - Sets the "umin" parameter used by the 192 PETSc routine for computing the differencing parameter, h, which is used 193 for matrix-free Jacobian-vector products. 194 195 Input Parameters: 196 + A - the matrix created with MatCreateSNESMF() 197 - umin - the parameter 198 199 Level: advanced 200 201 Notes: 202 See the manual page for MatCreateSNESMF() for a complete description of the 203 algorithm used to compute h. 204 205 .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF() 206 207 @*/ 208 PetscErrorCode MatMFFDDSSetUmin(Mat A,PetscReal umin) 209 { 210 PetscErrorCode ierr; 211 212 PetscFunctionBegin; 213 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 214 ierr = PetscTryMethod(A,"MatMFFDDSSetUmin_C",(Mat,PetscReal),(A,umin));CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 /*MC 219 MATMFFD_DS - the code for compute the "h" used in the finite difference 220 matrix-free matrix vector product. This code 221 implements the strategy in Dennis and Schnabel, "Numerical Methods for Unconstrained 222 Optimization and Nonlinear Equations". 223 224 Options Database Keys: 225 . -mat_mffd_umin <umin> see MatMFFDDSSetUmin() 226 227 Level: intermediate 228 229 Notes: 230 Requires 2 norms and 1 inner product, but they are computed together 231 so only one parallel collective operation is needed. See MATMFFD_WP for a method 232 (with GMRES) that requires NO collective operations. 233 234 Formula used: 235 F'(u)*a = [F(u+h*a) - F(u)]/h where 236 h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 237 = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 238 where 239 error_rel = square root of relative error in function evaluation 240 umin = minimum iterate parameter 241 242 .seealso: MATMFFD, MatCreateMFFD(), MatCreateSNESMF(), MATMFFD_WP, MatMFFDDSSetUmin() 243 244 M*/ 245 PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx) 246 { 247 MatMFFD_DS *hctx; 248 PetscErrorCode ierr; 249 250 PetscFunctionBegin; 251 /* allocate my own private data structure */ 252 ierr = PetscNewLog(ctx,&hctx);CHKERRQ(ierr); 253 ctx->hctx = (void*)hctx; 254 /* set a default for my parameter */ 255 hctx->umin = 1.e-6; 256 257 /* set the functions I am providing */ 258 ctx->ops->compute = MatMFFDCompute_DS; 259 ctx->ops->destroy = MatMFFDDestroy_DS; 260 ctx->ops->view = MatMFFDView_DS; 261 ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS; 262 263 ierr = PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDDSSetUmin_C",MatMFFDDSSetUmin_DS);CHKERRQ(ierr); 264 PetscFunctionReturn(0); 265 } 266 267 268 269 270 271 272 273