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 paramter (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,0);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 = (MatMFFD)mat->data; 179 MatMFFD_DS *hctx; 180 181 PetscFunctionBegin; 182 if (!ctx) SETERRQ(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: Requires 2 norms and 1 inner product, but they are computed together 228 so only one parallel collective operation is needed. See MATMFFD_WP for a method 229 (with GMRES) that requires NO collective operations. 230 231 Formula used: 232 F'(u)*a = [F(u+h*a) - F(u)]/h where 233 h = error_rel*u'a/||a||^2 if |u'a| > umin*||a||_{1} 234 = error_rel*umin*sign(u'a)*||a||_{1}/||a||^2 otherwise 235 where 236 error_rel = square root of relative error in function evaluation 237 umin = minimum iterate parameter 238 239 .seealso: MATMFFD, MatCreateMFFD(), MatCreateSNESMF(), MATMFFD_WP, MatMFFDDSSetUmin() 240 241 M*/ 242 PETSC_EXTERN PetscErrorCode MatCreateMFFD_DS(MatMFFD ctx) 243 { 244 MatMFFD_DS *hctx; 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 /* allocate my own private data structure */ 249 ierr = PetscNewLog(ctx,&hctx);CHKERRQ(ierr); 250 ctx->hctx = (void*)hctx; 251 /* set a default for my parameter */ 252 hctx->umin = 1.e-6; 253 254 /* set the functions I am providing */ 255 ctx->ops->compute = MatMFFDCompute_DS; 256 ctx->ops->destroy = MatMFFDDestroy_DS; 257 ctx->ops->view = MatMFFDView_DS; 258 ctx->ops->setfromoptions = MatMFFDSetFromOptions_DS; 259 260 ierr = PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDDSSetUmin_C",MatMFFDDSSetUmin_DS);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 265 266 267 268 269 270