1 #include <petscdevice.h>
2 #include <../src/ksp/ksp/utils/lmvm/rescale/symbrdnrescale.h> /*I "petscksp.h" I*/
3
MatSolve_DiagBrdn(Mat B,Vec F,Vec dX)4 static PetscErrorCode MatSolve_DiagBrdn(Mat B, Vec F, Vec dX)
5 {
6 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
7
8 PetscFunctionBegin;
9 PetscCall(MatSolve(lmvm->J0, F, dX));
10 PetscFunctionReturn(PETSC_SUCCESS);
11 }
12
MatMult_DiagBrdn(Mat B,Vec X,Vec Z)13 static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
14 {
15 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
16
17 PetscFunctionBegin;
18 PetscCall(MatMult(lmvm->J0, X, Z));
19 PetscFunctionReturn(PETSC_SUCCESS);
20 }
21
MatUpdate_DiagBrdn(Mat B,Vec X,Vec F)22 static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
23 {
24 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
25
26 PetscFunctionBegin;
27 if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
28 if (lmvm->prev_set) {
29 SymBroydenRescale ldb = (SymBroydenRescale)lmvm->ctx;
30 PetscScalar curvature;
31 PetscReal curvtol, ststmp;
32 PetscInt oldest, next;
33
34 PetscCall(MatLMVMGetRange(B, &oldest, &next));
35 /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
36 PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
37 PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
38
39 /* Test if the updates can be accepted */
40 PetscCall(VecDotNorm2(lmvm->Fprev, lmvm->Xprev, &curvature, &ststmp));
41 if (ststmp < lmvm->eps) curvtol = 0.0;
42 else curvtol = lmvm->eps * ststmp;
43
44 /* Test the curvature for the update */
45 if (PetscRealPart(curvature) > curvtol) {
46 /* Update is good so we accept it */
47 PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
48 PetscCall(MatLMVMProductsInsertDiagonalValue(B, LMBASIS_Y, LMBASIS_S, next, PetscRealPart(curvature)));
49 PetscCall(MatLMVMProductsInsertDiagonalValue(B, LMBASIS_S, LMBASIS_S, next, ststmp));
50 PetscCall(SymBroydenRescaleUpdate(B, ldb));
51 } else {
52 /* reset */
53 PetscCall(SymBroydenRescaleReset(B, ldb, MAT_LMVM_RESET_HISTORY));
54 }
55 /* End DiagBrdn update */
56 }
57 /* Save the solution and function to be used in the next update */
58 PetscCall(VecCopy(X, lmvm->Xprev));
59 PetscCall(VecCopy(F, lmvm->Fprev));
60 lmvm->prev_set = PETSC_TRUE;
61 PetscFunctionReturn(PETSC_SUCCESS);
62 }
63
MatCopy_DiagBrdn(Mat B,Mat M,MatStructure str)64 static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
65 {
66 Mat_LMVM *bdata = (Mat_LMVM *)B->data;
67 SymBroydenRescale bctx = (SymBroydenRescale)bdata->ctx;
68 Mat_LMVM *mdata = (Mat_LMVM *)M->data;
69 SymBroydenRescale mctx = (SymBroydenRescale)mdata->ctx;
70
71 PetscFunctionBegin;
72 PetscCall(SymBroydenRescaleCopy(bctx, mctx));
73 PetscFunctionReturn(PETSC_SUCCESS);
74 }
75
MatView_DiagBrdn(Mat B,PetscViewer pv)76 static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
77 {
78 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
79 SymBroydenRescale ldb = (SymBroydenRescale)lmvm->ctx;
80
81 PetscFunctionBegin;
82 PetscCall(MatView_LMVM(B, pv));
83 PetscCall(SymBroydenRescaleView(ldb, pv));
84 PetscFunctionReturn(PETSC_SUCCESS);
85 }
86
MatSetFromOptions_DiagBrdn(Mat B,PetscOptionItems PetscOptionsObject)87 static PetscErrorCode MatSetFromOptions_DiagBrdn(Mat B, PetscOptionItems PetscOptionsObject)
88 {
89 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
90 SymBroydenRescale ldb = (SymBroydenRescale)lmvm->ctx;
91
92 PetscFunctionBegin;
93 PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
94 PetscCall(SymBroydenRescaleSetFromOptions(B, ldb, PetscOptionsObject));
95 PetscFunctionReturn(PETSC_SUCCESS);
96 }
97
MatReset_DiagBrdn(Mat B,MatLMVMResetMode mode)98 static PetscErrorCode MatReset_DiagBrdn(Mat B, MatLMVMResetMode mode)
99 {
100 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
101 SymBroydenRescale ldb = (SymBroydenRescale)lmvm->ctx;
102
103 PetscFunctionBegin;
104 PetscCall(SymBroydenRescaleReset(B, ldb, mode));
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
MatDestroy_DiagBrdn(Mat B)108 static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
109 {
110 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
111
112 PetscFunctionBegin;
113 PetscCall(SymBroydenRescaleDestroy((SymBroydenRescale *)&lmvm->ctx));
114 PetscCall(MatDestroy_LMVM(B));
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
MatSetUp_DiagBrdn(Mat B)118 static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
119 {
120 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
121 SymBroydenRescale ldb = (SymBroydenRescale)lmvm->ctx;
122
123 PetscFunctionBegin;
124 PetscCall(MatSetUp_LMVM(B));
125 PetscCall(SymBroydenRescaleInitializeJ0(B, ldb));
126 PetscFunctionReturn(PETSC_SUCCESS);
127 }
128
MatCreate_LMVMDiagBrdn(Mat B)129 PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
130 {
131 Mat_LMVM *lmvm;
132 SymBroydenRescale ldb;
133
134 PetscFunctionBegin;
135 PetscCall(MatCreate_LMVM(B));
136 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN));
137 PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
138 PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
139 PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
140 B->ops->setup = MatSetUp_DiagBrdn;
141 B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
142 B->ops->destroy = MatDestroy_DiagBrdn;
143 B->ops->view = MatView_DiagBrdn;
144
145 lmvm = (Mat_LMVM *)B->data;
146 lmvm->ops->reset = MatReset_DiagBrdn;
147 lmvm->ops->mult = MatMult_DiagBrdn;
148 lmvm->ops->solve = MatSolve_DiagBrdn;
149 lmvm->ops->update = MatUpdate_DiagBrdn;
150 lmvm->ops->copy = MatCopy_DiagBrdn;
151
152 PetscCall(SymBroydenRescaleCreate(&ldb));
153 lmvm->ctx = (void *)ldb;
154
155 PetscCall(MatLMVMSetHistorySize(B, 1));
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158
159 /*@
160 MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
161 for approximating Hessians.
162
163 Collective
164
165 Input Parameters:
166 + comm - MPI communicator
167 . n - number of local rows for storage vectors
168 - N - global size of the storage vectors
169
170 Output Parameter:
171 . B - the matrix
172
173 Options Database Keys:
174 + -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
175 . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
176 . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
177 . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
178 . -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
179 . -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
180 - -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal
181
182 Level: intermediate
183
184 Notes:
185 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
186 paradigm instead of this routine directly.
187
188 It consists of a convex combination of DFP and BFGS
189 diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
190 To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
191 We also ensure positive definiteness by taking the `VecAbs()` of the final vector.
192
193 There are two ways of approximating the diagonal: using the forward (B) update
194 schemes for BFGS and DFP and then taking the inverse, or directly working with
195 the inverse (H) update schemes for the BFGS and DFP updates, derived using the
196 Sherman-Morrison-Woodbury formula. We have implemented both, controlled by a
197 parameter below.
198
199 In order to use the DiagBrdn matrix with other vector types, i.e. doing matrix-vector products
200 and matrix solves, the matrix must first be created using `MatCreate()` and `MatSetType()`,
201 followed by `MatLMVMAllocate()`. Then it will be available for updating
202 (via `MatLMVMUpdate()`) in one's favored solver implementation.
203
204 .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDIAGBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
205 `MatCreateLMVMBFGS()`, `MatCreateLMVMBroyden()`, `MatCreateLMVMSymBroyden()`
206 @*/
MatCreateLMVMDiagBroyden(MPI_Comm comm,PetscInt n,PetscInt N,Mat * B)207 PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
208 {
209 PetscFunctionBegin;
210 PetscCall(KSPInitializePackage());
211 PetscCall(MatCreate(comm, B));
212 PetscCall(MatSetSizes(*B, n, n, N, N));
213 PetscCall(MatSetType(*B, MATLMVMDIAGBROYDEN));
214 PetscCall(MatSetUp(*B));
215 PetscFunctionReturn(PETSC_SUCCESS);
216 }
217