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