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