1 #include "symbrdn.h" /*I "petscksp.h" I*/
2 #include <petscblaslapack.h>
3
4 /* Compute the forward symmetric Broyden scaling parameter phi corresponding the already known value of the "bad"
5 symmetric Broyden scaling parameter psi */
PhiFromPsi(PetscScalar psi,PetscScalar yts,PetscScalar stBs,PetscScalar ytHy)6 static inline PetscScalar PhiFromPsi(PetscScalar psi, PetscScalar yts, PetscScalar stBs, PetscScalar ytHy)
7 {
8 PetscScalar numer = (1.0 - psi) * PetscRealPart(PetscConj(yts) * yts);
9 PetscScalar phi = numer / (numer + psi * stBs * ytHy);
10 return phi;
11 }
12
13 /* The symmetric Broyden update can be written as
14
15 [ | ] [ a_k | b_k ] [ s_k^T B_k ]
16 B_{k+1} = B_k + [ B_k s_k | y_k ] [-----+-----] [-----------]
17 [ | ] [ b_k | c_k ] [ y_k^T ]
18
19 We can unroll this as
20
21 [ | ] [ a_i | b_i ] [ s_i^T B_i ]
22 B_{k+1} = B_0 + \sum_i [ B_i s_i | y_i ] [-----+-----] [-----------]
23 [ | ] [ b_i | c_i ] [ y_i^T ]
24
25 The a_i, b_i, and c_i values are stored in M00, M01, and M11, and are computed in
26 SymBroydenRecursiveBasisUpdate() below
27 */
SymBroydenKernel_Recursive_Inner(Mat B,MatLMVMMode mode,PetscInt oldest,PetscInt next,Vec X,Vec B0X)28 static PetscErrorCode SymBroydenKernel_Recursive_Inner(Mat B, MatLMVMMode mode, PetscInt oldest, PetscInt next, Vec X, Vec B0X)
29 {
30 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
31 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
32 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
33 LMBasis BkS = lsb->basis[LMVMModeMap(SYMBROYDEN_BASIS_BKS, mode)];
34 LMProducts M00 = lsb->products[LMVMModeMap(SYMBROYDEN_PRODUCTS_M00, mode)];
35 LMProducts M01 = lsb->products[LMVMModeMap(SYMBROYDEN_PRODUCTS_M01, mode)];
36 LMProducts M11 = lsb->products[LMVMModeMap(SYMBROYDEN_PRODUCTS_M11, mode)];
37 LMBasis Y;
38 Vec StBkX, YtX, U, V;
39
40 PetscFunctionBegin;
41 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
42 PetscCall(MatLMVMGetWorkRow(B, &StBkX));
43 PetscCall(MatLMVMGetWorkRow(B, &YtX));
44 PetscCall(MatLMVMGetWorkRow(B, &U));
45 PetscCall(MatLMVMGetWorkRow(B, &V));
46 PetscCall(LMBasisGEMVH(BkS, oldest, next, 1.0, X, 0.0, StBkX));
47 PetscCall(LMBasisGEMVH(Y, oldest, next, 1.0, X, 0.0, YtX));
48 PetscCall(LMProductsMult(M00, oldest, next, 1.0, StBkX, 0.0, U, PETSC_FALSE));
49 PetscCall(LMProductsMult(M01, oldest, next, 1.0, YtX, 1.0, U, PETSC_FALSE));
50 PetscCall(LMProductsMult(M01, oldest, next, 1.0, StBkX, 0.0, V, PETSC_FALSE));
51 PetscCall(LMProductsMult(M11, oldest, next, 1.0, YtX, 1.0, V, PETSC_FALSE));
52 PetscCall(LMBasisGEMV(BkS, oldest, next, 1.0, U, 1.0, B0X));
53 PetscCall(LMBasisGEMV(Y, oldest, next, 1.0, V, 1.0, B0X));
54 PetscCall(MatLMVMRestoreWorkRow(B, &V));
55 PetscCall(MatLMVMRestoreWorkRow(B, &U));
56 PetscCall(MatLMVMRestoreWorkRow(B, &YtX));
57 PetscCall(MatLMVMRestoreWorkRow(B, &StBkX));
58 PetscFunctionReturn(PETSC_SUCCESS);
59 }
60
MatLMVMSymBroydenGetConvexFactor(Mat B,SymBroydenProductsType Phi_t,LMProducts * Phi)61 static PetscErrorCode MatLMVMSymBroydenGetConvexFactor(Mat B, SymBroydenProductsType Phi_t, LMProducts *Phi)
62 {
63 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
64 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
65 PetscScalar phi = Phi_t == SYMBROYDEN_PRODUCTS_PHI ? lsb->phi_scalar : lsb->psi_scalar;
66
67 PetscFunctionBegin;
68 if (!lsb->products[Phi_t]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lsb->products[Phi_t]));
69 *Phi = lsb->products[Phi_t];
70 if (phi != PETSC_DETERMINE) {
71 PetscInt oldest, next, start;
72
73 PetscCall(MatLMVMGetRange(B, &oldest, &next));
74 start = PetscMax((*Phi)->k, oldest);
75 (*Phi)->k = start;
76 for (PetscInt i = start; i < next; i++) PetscCall(LMProductsInsertNextDiagonalValue(*Phi, i, phi));
77 }
78 PetscFunctionReturn(PETSC_SUCCESS);
79 }
80
SymBroydenRecursiveBasisUpdate(Mat B,MatLMVMMode mode,PetscBool update_phi_from_psi)81 static PetscErrorCode SymBroydenRecursiveBasisUpdate(Mat B, MatLMVMMode mode, PetscBool update_phi_from_psi)
82 {
83 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
84 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
85 MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
86 MatLMVMBasisType B0S_t = LMVMModeMap(LMBASIS_B0S, mode);
87 SymBroydenProductsType Phi_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_PHI, mode);
88 SymBroydenProductsType Psi_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_PSI, mode);
89 SymBroydenProductsType StBkS_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_STBKS, mode);
90 SymBroydenProductsType YtHkY_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_YTHKY, mode);
91 SymBroydenProductsType M00_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_M00, mode);
92 SymBroydenProductsType M01_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_M01, mode);
93 SymBroydenProductsType M11_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_M11, mode);
94 SymBroydenProductsType M_t[3] = {M00_t, M01_t, M11_t};
95 LMProducts M[3];
96 LMProducts Phi, Psi = NULL;
97 SymBroydenBasisType BkS_t = LMVMModeMap(SYMBROYDEN_BASIS_BKS, mode);
98 LMBasis BkS;
99 LMProducts StBkS, YtHkY = NULL;
100 PetscInt oldest, start, next;
101 PetscInt products_oldest;
102 LMBasis S;
103
104 PetscFunctionBegin;
105 PetscCall(MatLMVMGetRange(B, &oldest, &next));
106 for (PetscInt i = 0; i < 3; i++) {
107 if (lsb->products[M_t[i]] && lsb->products[M_t[i]]->block_type != LMBLOCK_DIAGONAL) PetscCall(LMProductsDestroy(&lsb->products[M_t[i]]));
108 if (!lsb->products[M_t[i]]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lsb->products[M_t[i]]));
109 M[i] = lsb->products[M_t[i]];
110 }
111 if (!lsb->basis[BkS_t]) PetscCall(LMBasisCreate(MatLMVMBasisSizeOf(B0S_t) == LMBASIS_S ? lmvm->Xprev : lmvm->Fprev, lmvm->m, &lsb->basis[BkS_t]));
112 BkS = lsb->basis[BkS_t];
113 PetscCall(MatLMVMSymBroydenGetConvexFactor(B, Phi_t, &Phi));
114 if (!lsb->products[StBkS_t]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lsb->products[StBkS_t]));
115 StBkS = lsb->products[StBkS_t];
116 PetscCall(LMProductsPrepare(StBkS, lmvm->J0, oldest, next));
117 products_oldest = PetscMax(0, StBkS->k - lmvm->m);
118 if (oldest > products_oldest) {
119 // recursion is starting from a different starting index, it must be recomputed
120 StBkS->k = oldest;
121 }
122 BkS->k = start = StBkS->k;
123 for (PetscInt i = 0; i < 3; i++) M[i]->k = start;
124 if (start == next) PetscFunctionReturn(PETSC_SUCCESS);
125
126 if (update_phi_from_psi) {
127 Phi->k = start;
128 // we have to first make sure that the inverse data is up to date
129 PetscCall(SymBroydenRecursiveBasisUpdate(B, (MatLMVMMode)(mode ^ 1), PETSC_FALSE));
130 PetscCall(MatLMVMSymBroydenGetConvexFactor(B, Psi_t, &Psi));
131 YtHkY = lsb->products[YtHkY_t];
132 }
133 PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
134
135 for (PetscInt j = start; j < next; j++) {
136 Vec p_j, s_j, B0s_j;
137 PetscScalar alpha, sjtbjsj;
138 PetscScalar m_00, m_01, m_11, yjtsj;
139 PetscScalar phi_j;
140
141 PetscCall(LMBasisGetWorkVec(BkS, &p_j));
142 // p_j starts as B_0 * s_j
143 PetscCall(MatLMVMBasisGetVecRead(B, B0S_t, j, &B0s_j, &alpha));
144 PetscCall(VecAXPBY(p_j, alpha, 0.0, B0s_j));
145 PetscCall(MatLMVMBasisRestoreVecRead(B, B0S_t, j, &B0s_j, &alpha));
146
147 // Use the matmult kernel to compute p_j = B_j * p_j
148 PetscCall(LMBasisGetVecRead(S, j, &s_j));
149 if (j > oldest) PetscCall(SymBroydenKernel_Recursive_Inner(B, mode, oldest, j, s_j, p_j));
150 PetscCall(VecDot(p_j, s_j, &sjtbjsj));
151 PetscCall(LMBasisRestoreVecRead(S, j, &s_j));
152 PetscCall(LMProductsInsertNextDiagonalValue(StBkS, j, sjtbjsj));
153 PetscCall(LMBasisSetNextVec(BkS, p_j));
154 PetscCall(LMBasisRestoreWorkVec(BkS, &p_j));
155
156 PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_Y, LMBASIS_S, j, &yjtsj));
157 if (update_phi_from_psi) {
158 PetscScalar psi_j;
159 PetscScalar yjthjyj;
160
161 PetscCall(LMProductsGetDiagonalValue(YtHkY, j, &yjthjyj));
162 PetscCall(LMProductsGetDiagonalValue(Psi, j, &psi_j));
163
164 phi_j = PhiFromPsi(psi_j, yjtsj, sjtbjsj, yjthjyj);
165 PetscCall(LMProductsInsertNextDiagonalValue(Phi, j, phi_j));
166 } else PetscCall(LMProductsGetDiagonalValue(Phi, j, &phi_j));
167
168 /* The symmetric Broyden update can be represented as
169
170 [ | ][ (phi - 1) / s_j^T p_j | -phi / y_j^T s_j ][ p_j^T ]
171 [ p_j | y_j ][-----------------------+-----------------------------------------------][-------]
172 [ | ][ -phi / y_j^T s_j | (y_j^T s_j + phi * s_j^T p_j) / (y_j^T s_j)^2 ][ y_j^T ]
173
174 We store diagonal vectors with these values
175 */
176
177 m_00 = (phi_j - 1.0) / sjtbjsj;
178 m_01 = -phi_j / yjtsj;
179 m_11 = (yjtsj + phi_j * sjtbjsj) / (yjtsj * yjtsj);
180 PetscCall(LMProductsInsertNextDiagonalValue(M[0], j, m_00));
181 PetscCall(LMProductsInsertNextDiagonalValue(M[1], j, m_01));
182 PetscCall(LMProductsInsertNextDiagonalValue(M[2], j, m_11));
183 }
184 PetscFunctionReturn(PETSC_SUCCESS);
185 }
186
SymBroydenKernel_Recursive(Mat B,MatLMVMMode mode,Vec X,Vec Y,PetscBool update_phi_from_psi)187 PETSC_INTERN PetscErrorCode SymBroydenKernel_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec Y, PetscBool update_phi_from_psi)
188 {
189 PetscInt oldest, next;
190
191 PetscFunctionBegin;
192 PetscCall(MatLMVMApplyJ0Mode(mode)(B, X, Y));
193 PetscCall(MatLMVMGetRange(B, &oldest, &next));
194 if (next > oldest) {
195 PetscCall(SymBroydenRecursiveBasisUpdate(B, mode, update_phi_from_psi));
196 PetscCall(SymBroydenKernel_Recursive_Inner(B, mode, oldest, next, X, Y));
197 }
198 PetscFunctionReturn(PETSC_SUCCESS);
199 }
200
MatMult_LMVMSymBrdn_Recursive(Mat B,Vec X,Vec Y)201 static PetscErrorCode MatMult_LMVMSymBrdn_Recursive(Mat B, Vec X, Vec Y)
202 {
203 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
204 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
205
206 PetscFunctionBegin;
207 PetscCheck(lsb->phi_scalar != PETSC_DETERMINE, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Phi must first be set using MatLMVMSymBroydenSetPhi()");
208 if (lsb->phi_scalar == 0.0) {
209 PetscCall(BFGSKernel_Recursive(B, MATLMVM_MODE_PRIMAL, X, Y));
210 } else if (lsb->phi_scalar == 1.0) {
211 PetscCall(DFPKernel_Recursive(B, MATLMVM_MODE_PRIMAL, X, Y));
212 } else {
213 PetscCall(SymBroydenKernel_Recursive(B, MATLMVM_MODE_PRIMAL, X, Y, PETSC_FALSE));
214 }
215 PetscFunctionReturn(PETSC_SUCCESS);
216 }
217
MatSolve_LMVMSymBrdn_Recursive(Mat B,Vec X,Vec Y)218 static PetscErrorCode MatSolve_LMVMSymBrdn_Recursive(Mat B, Vec X, Vec Y)
219 {
220 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
221 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
222
223 PetscFunctionBegin;
224 PetscCheck(lsb->phi_scalar != PETSC_DETERMINE, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Phi must first be set using MatLMVMSymBroydenSetPhi()");
225 if (lsb->phi_scalar == 0.0) {
226 PetscCall(DFPKernel_Recursive(B, MATLMVM_MODE_DUAL, X, Y));
227 } else if (lsb->phi_scalar == 1.0) {
228 PetscCall(BFGSKernel_Recursive(B, MATLMVM_MODE_DUAL, X, Y));
229 } else {
230 PetscCall(SymBroydenKernel_Recursive(B, MATLMVM_MODE_DUAL, X, Y, PETSC_TRUE));
231 }
232 PetscFunctionReturn(PETSC_SUCCESS);
233 }
234
235 PetscBool ErwayMarciaCite = PETSC_FALSE;
236 const char ErwayMarciaCitation[] = "@article{Erway2015,"
237 " title = {On Efficiently Computing the Eigenvalues of Limited-Memory Quasi-Newton Matrices},"
238 " volume = {36},"
239 " ISSN = {1095-7162},"
240 " url = {http://dx.doi.org/10.1137/140997737},"
241 " DOI = {10.1137/140997737},"
242 " number = {3},"
243 " journal = {SIAM Journal on Matrix Analysis and Applications},"
244 " publisher = {Society for Industrial & Applied Mathematics (SIAM)},"
245 " author = {Erway, Jennifer B. and Marcia, Roummel F.},"
246 " year = {2015},"
247 " month = jan,"
248 " pages = {1338-1359}"
249 "}\n";
250
251 // TODO: on device (e.g. cuBLAS) implementation?
SymBroydenCompactDenseUpdateArrays(PetscBLASInt m,PetscBLASInt oldest,PetscBLASInt next,PetscScalar M00[],PetscBLASInt lda00,PetscScalar M01[],PetscBLASInt lda01,PetscScalar M11[],PetscBLASInt lda11,const PetscScalar StB0S[],PetscBLASInt ldasbs,const PetscScalar YtS[],PetscBLASInt ldays,const PetscScalar Phi[],PetscScalar p0[],PetscScalar p1[],const PetscScalar Psi[],const PetscScalar YtHkY[],PetscScalar StBkS[])252 static PetscErrorCode SymBroydenCompactDenseUpdateArrays(PetscBLASInt m, PetscBLASInt oldest, PetscBLASInt next, PetscScalar M00[], PetscBLASInt lda00, PetscScalar M01[], PetscBLASInt lda01, PetscScalar M11[], PetscBLASInt lda11, const PetscScalar StB0S[], PetscBLASInt ldasbs, const PetscScalar YtS[], PetscBLASInt ldays, const PetscScalar Phi[], PetscScalar p0[], PetscScalar p1[], const PetscScalar Psi[], const PetscScalar YtHkY[], PetscScalar StBkS[])
253 {
254 PetscBLASInt i;
255 PetscScalar alpha, beta, delta;
256 PetscBLASInt ione = 1;
257 PetscScalar sone = 1.0;
258 PetscScalar szero = 0.0;
259 PetscScalar sBis;
260 PetscScalar yts;
261 PetscScalar phi;
262
263 PetscFunctionBegin;
264 if (next <= oldest || m <= 0) PetscFunctionReturn(PETSC_SUCCESS);
265
266 PetscCall(PetscCitationsRegister(ErwayMarciaCitation, &ErwayMarciaCite));
267
268 PetscCall(PetscArrayzero(M00, m * lda00));
269 PetscCall(PetscArrayzero(M01, m * lda01));
270 PetscCall(PetscArrayzero(M11, m * lda11));
271
272 i = oldest % m;
273
274 // Base case entries
275 sBis = StB0S[i + i * ldasbs];
276 if (StBkS) StBkS[i] = sBis;
277 yts = YtS[i + i * ldays];
278 if (Psi) {
279 phi = PhiFromPsi(Psi[i], yts, sBis, YtHkY[i]);
280 } else {
281 phi = Phi[i];
282 }
283 alpha = PetscRealPart(-(1.0 - phi) / sBis);
284 beta = -phi / yts;
285 delta = (1.0 + phi * sBis / yts) / yts;
286 M00[i + i * lda00] = alpha;
287 M01[i + i * lda01] = beta;
288 M11[i + i * lda11] = delta;
289 for (PetscBLASInt i_ = oldest + 1; i_ < next; i_++) {
290 const PetscScalar *q0, *q1;
291 PetscScalar qp;
292
293 i = i_ % m;
294 q0 = &StB0S[0 + i * ldasbs];
295 q1 = &YtS[0 + i * ldays];
296
297 // p_i <- M_{i-1} q_i
298
299 PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &m, &sone, M00, &lda00, q0, &ione, &szero, p0, &ione));
300 PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &m, &sone, M01, &lda01, q1, &ione, &sone, p0, &ione));
301 PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &m, &sone, M01, &lda01, q0, &ione, &szero, p1, &ione));
302 PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &m, &sone, M11, &lda11, q1, &ione, &sone, p1, &ione));
303
304 // q'p
305 PetscCallBLAS("BLASdot", qp = BLASdot_(&m, q0, &ione, p0, &ione));
306 PetscCallBLAS("BLASdot", qp += BLASdot_(&m, q1, &ione, p1, &ione));
307
308 sBis = StB0S[i + i * ldasbs] + qp;
309 if (StBkS) StBkS[i] = sBis;
310 yts = YtS[i + i * ldays];
311 if (Psi) {
312 phi = PhiFromPsi(Psi[i], yts, sBis, YtHkY[i]);
313 } else {
314 phi = Phi[i];
315 }
316
317 alpha = PetscRealPart(-(1.0 - phi) / sBis);
318 beta = -phi / yts;
319 delta = (1.0 + phi * sBis / yts) / yts;
320
321 PetscCallBLAS("LAPACKgerc", LAPACKgerc_(&m, &m, &alpha, p0, &ione, p0, &ione, M00, &lda00));
322 for (PetscInt j = 0; j < m; j++) M00[j + i * lda00] = alpha * p0[j];
323 for (PetscInt j = 0; j < m; j++) M00[i + j * lda00] = PetscConj(alpha * p0[j]);
324 M00[i + i * lda00] = alpha;
325
326 PetscCallBLAS("LAPACKgerc", LAPACKgerc_(&m, &m, &alpha, p0, &ione, p1, &ione, M01, &lda01));
327 for (PetscBLASInt j = 0; j < m; j++) M01[j + i * lda01] = beta * p0[j];
328 for (PetscBLASInt j = 0; j < m; j++) M01[i + j * lda01] = PetscConj(alpha * p1[j]);
329 M01[i + i * lda01] = beta;
330
331 PetscCallBLAS("LAPACKgerc", LAPACKgerc_(&m, &m, &alpha, p1, &ione, p1, &ione, M11, &lda11));
332 for (PetscInt j = 0; j < m; j++) M11[j + i * lda11] = beta * p1[j];
333 for (PetscInt j = 0; j < m; j++) M11[i + j * lda11] = PetscConj(beta * p1[j]);
334 M11[i + i * lda11] = delta;
335 }
336 PetscFunctionReturn(PETSC_SUCCESS);
337 }
338
SymBroydenCompactProductsUpdate(Mat B,MatLMVMMode mode,PetscBool update_phi_from_psi)339 static PetscErrorCode SymBroydenCompactProductsUpdate(Mat B, MatLMVMMode mode, PetscBool update_phi_from_psi)
340 {
341 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
342 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
343 PetscInt oldest, next;
344 MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
345 MatLMVMBasisType B0S_t = LMVMModeMap(LMBASIS_B0S, mode);
346 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
347 SymBroydenProductsType Phi_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_PHI, mode);
348 SymBroydenProductsType Psi_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_PSI, mode);
349 SymBroydenProductsType StBkS_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_STBKS, mode);
350 SymBroydenProductsType YtHkY_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_YTHKY, mode);
351 SymBroydenProductsType M00_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_M00, mode);
352 SymBroydenProductsType M01_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_M01, mode);
353 SymBroydenProductsType M11_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_M11, mode);
354 SymBroydenProductsType M_t[3] = {M00_t, M01_t, M11_t};
355 Mat M_local[3];
356 LMProducts M[3], Phi, Psi, YtS, StB0S, StBkS, YtHkY;
357 PetscInt products_oldest, k;
358 PetscBool local_is_nonempty;
359
360 PetscFunctionBegin;
361 PetscCall(MatLMVMGetRange(B, &oldest, &next));
362 for (PetscInt i = 0; i < 3; i++) {
363 if (lsb->products[M_t[i]] && lsb->products[M_t[i]]->block_type != LMBLOCK_FULL) PetscCall(LMProductsDestroy(&lsb->products[M_t[i]]));
364 if (!lsb->products[M_t[i]]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_FULL, &lsb->products[M_t[i]]));
365 M[i] = lsb->products[M_t[i]];
366 }
367 PetscCall(MatLMVMSymBroydenGetConvexFactor(B, Phi_t, &Phi));
368 PetscCall(MatLMVMSymBroydenGetConvexFactor(B, Psi_t, &Psi));
369 if (!lsb->products[StBkS_t]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lsb->products[StBkS_t]));
370 StBkS = lsb->products[StBkS_t];
371 if (!lsb->products[YtHkY_t]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lsb->products[YtHkY_t]));
372 YtHkY = lsb->products[YtHkY_t];
373 PetscCall(LMProductsPrepare(M[0], lmvm->J0, oldest, next));
374 PetscCall(LMProductsGetLocalMatrix(M[0], &M_local[0], &k, &local_is_nonempty));
375 products_oldest = PetscMax(0, k - lmvm->m);
376 if (products_oldest < oldest) k = oldest;
377 if (k < next) {
378 Mat StB0S_local, YtS_local;
379 Vec Psi_local = NULL, Phi_local = NULL, YtHkY_local = NULL, StBkS_local = NULL;
380
381 PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, S_t, LMBLOCK_UPPER_TRIANGLE, &YtS));
382 PetscCall(MatLMVMGetUpdatedProducts(B, S_t, B0S_t, LMBLOCK_UPPER_TRIANGLE, &StB0S));
383 PetscCall(LMProductsGetLocalMatrix(StB0S, &StB0S_local, NULL, NULL));
384 PetscCall(LMProductsGetLocalMatrix(YtS, &YtS_local, NULL, NULL));
385 PetscCall(LMProductsGetLocalMatrix(M[1], &M_local[1], NULL, NULL));
386 PetscCall(LMProductsGetLocalMatrix(M[2], &M_local[2], NULL, NULL));
387 if (update_phi_from_psi) {
388 PetscCall(LMProductsGetLocalDiagonal(Psi, &Psi_local));
389 PetscCall(LMProductsGetLocalDiagonal(YtHkY, &YtHkY_local));
390 } else {
391 PetscCall(LMProductsGetLocalDiagonal(Phi, &Phi_local));
392 PetscCall(LMProductsGetLocalDiagonal(StBkS, &StBkS_local));
393 }
394 if (local_is_nonempty) {
395 PetscInt lda;
396 PetscBLASInt M_lda[3], StB0S_lda, YtS_lda, m_blas, oldest_blas, next_blas;
397 const PetscScalar *StB0S_;
398 const PetscScalar *YtS_;
399 PetscScalar *M_[3];
400 const PetscScalar *Phi_ = NULL;
401 const PetscScalar *Psi_ = NULL;
402 const PetscScalar *YtHkY_ = NULL;
403 PetscScalar *StBkS_ = NULL;
404 PetscScalar *p0, *p1;
405
406 for (PetscInt i = 0; i < 3; i++) {
407 PetscCall(MatDenseGetLDA(M_local[i], &lda));
408 PetscCall(PetscBLASIntCast(lda, &M_lda[i]));
409 PetscCall(MatDenseGetArrayWrite(M_local[i], &M_[i]));
410 }
411
412 PetscCall(MatDenseGetArrayRead(StB0S_local, &StB0S_));
413 PetscCall(MatDenseGetLDA(StB0S_local, &lda));
414 PetscCall(PetscBLASIntCast(lda, &StB0S_lda));
415
416 PetscCall(MatDenseGetArrayRead(YtS_local, &YtS_));
417 PetscCall(MatDenseGetLDA(YtS_local, &lda));
418 PetscCall(PetscBLASIntCast(lda, &YtS_lda));
419
420 PetscCall(PetscBLASIntCast(lmvm->m, &m_blas));
421 PetscCall(PetscBLASIntCast(oldest, &oldest_blas));
422 PetscCall(PetscBLASIntCast(next, &next_blas));
423
424 if (update_phi_from_psi) {
425 PetscCall(VecGetArrayRead(Psi_local, &Psi_));
426 PetscCall(VecGetArrayRead(YtHkY_local, &YtHkY_));
427 } else {
428 PetscCall(VecGetArrayRead(Phi_local, &Phi_));
429 PetscCall(VecGetArrayWrite(StBkS_local, &StBkS_));
430 }
431
432 PetscCall(PetscMalloc2(lmvm->m, &p0, lmvm->m, &p1));
433 PetscCall(SymBroydenCompactDenseUpdateArrays(m_blas, oldest_blas, next_blas, M_[0], M_lda[0], M_[1], M_lda[1], M_[2], M_lda[2], StB0S_, StB0S_lda, YtS_, YtS_lda, Phi_, p0, p1, Psi_, YtHkY_, StBkS_));
434 PetscCall(PetscFree2(p0, p1));
435
436 if (update_phi_from_psi) {
437 PetscCall(VecRestoreArrayRead(YtHkY_local, &YtHkY_));
438 PetscCall(VecRestoreArrayRead(Psi_local, &Psi_));
439 } else {
440 PetscCall(VecRestoreArrayWrite(StBkS_local, &StBkS_));
441 PetscCall(VecRestoreArrayRead(Phi_local, &Phi_));
442 }
443
444 for (PetscInt i = 0; i < 3; i++) PetscCall(MatDenseRestoreArrayWrite(M_local[i], &M_[i]));
445 PetscCall(MatDenseRestoreArrayRead(YtS_local, &YtS_));
446 PetscCall(MatDenseRestoreArrayRead(StB0S_local, &StB0S_));
447 }
448 if (update_phi_from_psi) {
449 PetscCall(LMProductsRestoreLocalDiagonal(YtHkY, &YtHkY_local));
450 PetscCall(LMProductsRestoreLocalDiagonal(Psi, &Psi_local));
451 } else {
452 PetscCall(LMProductsRestoreLocalDiagonal(StBkS, &StBkS_local));
453 PetscCall(LMProductsRestoreLocalDiagonal(Phi, &Phi_local));
454 }
455 PetscCall(LMProductsRestoreLocalMatrix(M[2], &M_local[2], &next));
456 PetscCall(LMProductsRestoreLocalMatrix(M[1], &M_local[1], &next));
457 PetscCall(LMProductsRestoreLocalMatrix(YtS, &YtS_local, NULL));
458 PetscCall(LMProductsRestoreLocalMatrix(StB0S, &StB0S_local, NULL));
459 }
460 PetscCall(LMProductsRestoreLocalMatrix(M[0], &M_local[0], &next));
461 PetscFunctionReturn(PETSC_SUCCESS);
462 }
463
SymBroydenCompactDenseKernelUseB0S(Mat B,MatLMVMMode mode,Vec X,PetscBool * use_B0S)464 PETSC_INTERN PetscErrorCode SymBroydenCompactDenseKernelUseB0S(Mat B, MatLMVMMode mode, Vec X, PetscBool *use_B0S)
465 {
466 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
467 PetscBool is_scalar;
468 PetscScalar J0_scale;
469 LMBasis B0S;
470 Mat J0 = lmvm->J0;
471 PetscObjectId id;
472 PetscObjectState state;
473
474 /*
475 The one time where we would want to compute S^T B_0 X as (B_0 S)^T X instead of S^T (B_0 X)
476 is if (B_0 S)^T X is cached.
477 */
478 PetscFunctionBegin;
479 *use_B0S = PETSC_FALSE;
480 PetscCall(MatLMVMGetJ0Scalar(B, &is_scalar, &J0_scale));
481 B0S = lmvm->basis[is_scalar ? LMVMModeMap(LMBASIS_S, mode) : LMVMModeMap(LMBASIS_B0S, mode)];
482 if ((B0S->k < lmvm->k) || (B0S->cached_product == NULL)) PetscFunctionReturn(PETSC_SUCCESS);
483 if (!is_scalar) {
484 PetscCall(PetscObjectGetId((PetscObject)J0, &id));
485 PetscCall(PetscObjectStateGet((PetscObject)J0, &state));
486 if (id != B0S->operator_id || state != B0S->operator_state) PetscFunctionReturn(PETSC_SUCCESS);
487 }
488 PetscCall(PetscObjectGetId((PetscObject)X, &id));
489 PetscCall(PetscObjectStateGet((PetscObject)X, &state));
490 if (id == B0S->cached_vec_id && state == B0S->cached_vec_state) *use_B0S = PETSC_TRUE;
491 PetscFunctionReturn(PETSC_SUCCESS);
492 }
493
SymBroydenKernel_CompactDense(Mat B,MatLMVMMode mode,Vec X,Vec BX,PetscBool update_phi_from_psi)494 PETSC_INTERN PetscErrorCode SymBroydenKernel_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BX, PetscBool update_phi_from_psi)
495 {
496 PetscInt oldest, next;
497
498 PetscFunctionBegin;
499 PetscCall(MatLMVMApplyJ0Mode(mode)(B, X, BX));
500 PetscCall(MatLMVMGetRange(B, &oldest, &next));
501 if (next > oldest) {
502 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
503 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
504 Vec StB0X, YtX, u, v;
505 MatLMVMBasisType S_t = LMVMModeMap(LMBASIS_S, mode);
506 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode);
507 MatLMVMBasisType B0S_t = LMVMModeMap(LMBASIS_B0S, mode);
508 SymBroydenProductsType M00_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_M00, mode);
509 SymBroydenProductsType M01_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_M01, mode);
510 SymBroydenProductsType M11_t = LMVMModeMap(SYMBROYDEN_PRODUCTS_M11, mode);
511 LMProducts M00, M01, M11;
512 LMBasis S, Y;
513 PetscBool use_B0S;
514
515 if (update_phi_from_psi) PetscCall(SymBroydenCompactProductsUpdate(B, (MatLMVMMode)(mode ^ 1), PETSC_FALSE));
516 PetscCall(SymBroydenCompactProductsUpdate(B, mode, update_phi_from_psi));
517 M00 = lsb->products[M00_t];
518 M01 = lsb->products[M01_t];
519 M11 = lsb->products[M11_t];
520
521 PetscCall(MatLMVMGetUpdatedBasis(B, S_t, &S, NULL, NULL));
522 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL));
523
524 PetscCall(MatLMVMGetWorkRow(B, &StB0X));
525 PetscCall(MatLMVMGetWorkRow(B, &YtX));
526 PetscCall(MatLMVMGetWorkRow(B, &u));
527 PetscCall(MatLMVMGetWorkRow(B, &v));
528
529 PetscCall(SymBroydenCompactDenseKernelUseB0S(B, mode, X, &use_B0S));
530 if (use_B0S) PetscCall(MatLMVMBasisGEMVH(B, B0S_t, oldest, next, 1.0, X, 0.0, StB0X));
531 else PetscCall(LMBasisGEMVH(S, oldest, next, 1.0, BX, 0.0, StB0X));
532
533 PetscCall(LMBasisGEMVH(Y, oldest, next, 1.0, X, 0.0, YtX));
534
535 PetscCall(LMProductsMult(M00, oldest, next, 1.0, StB0X, 0.0, u, PETSC_FALSE));
536 PetscCall(LMProductsMult(M01, oldest, next, 1.0, YtX, 1.0, u, PETSC_FALSE));
537 PetscCall(LMProductsMult(M01, oldest, next, 1.0, StB0X, 0.0, v, PETSC_TRUE));
538 PetscCall(LMProductsMult(M11, oldest, next, 1.0, YtX, 1.0, v, PETSC_FALSE));
539
540 PetscCall(LMBasisGEMV(Y, oldest, next, 1.0, v, 1.0, BX));
541 PetscCall(MatLMVMBasisGEMV(B, B0S_t, oldest, next, 1.0, u, 1.0, BX));
542
543 PetscCall(MatLMVMRestoreWorkRow(B, &v));
544 PetscCall(MatLMVMRestoreWorkRow(B, &u));
545 PetscCall(MatLMVMRestoreWorkRow(B, &YtX));
546 PetscCall(MatLMVMRestoreWorkRow(B, &StB0X));
547 }
548 PetscFunctionReturn(PETSC_SUCCESS);
549 }
550
MatMult_LMVMSymBrdn_CompactDense(Mat B,Vec X,Vec BX)551 static PetscErrorCode MatMult_LMVMSymBrdn_CompactDense(Mat B, Vec X, Vec BX)
552 {
553 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
554 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
555
556 PetscFunctionBegin;
557 PetscCheck(lsb->phi_scalar != PETSC_DETERMINE, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Phi must first be set using MatLMVMSymBroydenSetPhi()");
558 if (lsb->phi_scalar == 0.0) {
559 PetscCall(BFGSKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, X, BX));
560 } else if (lsb->phi_scalar == 1.0) {
561 PetscCall(DFPKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, X, BX));
562 } else {
563 PetscCall(SymBroydenKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, X, BX, PETSC_FALSE));
564 }
565 PetscFunctionReturn(PETSC_SUCCESS);
566 }
567
MatSolve_LMVMSymBrdn_CompactDense(Mat B,Vec X,Vec HX)568 static PetscErrorCode MatSolve_LMVMSymBrdn_CompactDense(Mat B, Vec X, Vec HX)
569 {
570 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
571 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
572
573 PetscFunctionBegin;
574 PetscCheck(lsb->phi_scalar != PETSC_DETERMINE, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Phi must first be set using MatLMVMSymBroydenSetPhi()");
575 if (lsb->phi_scalar == 0.0) {
576 PetscCall(DFPKernel_CompactDense(B, MATLMVM_MODE_DUAL, X, HX));
577 } else if (lsb->phi_scalar == 1.0) {
578 PetscCall(BFGSKernel_CompactDense(B, MATLMVM_MODE_DUAL, X, HX));
579 } else {
580 PetscCall(SymBroydenKernel_CompactDense(B, MATLMVM_MODE_DUAL, X, HX, PETSC_TRUE));
581 }
582 PetscFunctionReturn(PETSC_SUCCESS);
583 }
584
MatUpdate_LMVMSymBrdn(Mat B,Vec X,Vec F)585 static PetscErrorCode MatUpdate_LMVMSymBrdn(Mat B, Vec X, Vec F)
586 {
587 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
588 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
589 PetscReal curvtol;
590 PetscScalar curvature, ststmp;
591 PetscInt oldest, next;
592 PetscBool cache_StFprev = (lmvm->mult_alg != MAT_LMVM_MULT_RECURSIVE) ? lmvm->cache_gradient_products : PETSC_FALSE;
593 PetscBool cache_YtH0Fprev = cache_StFprev;
594 LMBasis S = NULL, H0Y = NULL;
595 PetscScalar H0_alpha = 1.0;
596 MatLMVMBasisType H0Y_t = LMBASIS_H0Y;
597
598 PetscFunctionBegin;
599 if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
600 // BFGS using the dense algorithm does not need to cache YtH0F products
601 if (lsb->phi_scalar == 0.0 && lmvm->mult_alg == MAT_LMVM_MULT_DENSE) cache_YtH0Fprev = PETSC_FALSE;
602 // Caching is no use if we are diagonally updating
603 if (lsb->rescale->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) cache_YtH0Fprev = PETSC_FALSE;
604 PetscCall(MatLMVMGetRange(B, &oldest, &next));
605 if (lmvm->prev_set) {
606 LMBasis Y = NULL;
607 Vec Fprev_old = NULL;
608
609 if (cache_StFprev || cache_YtH0Fprev) {
610 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &Y, NULL, NULL));
611 PetscCall(LMBasisGetWorkVec(Y, &Fprev_old));
612 PetscCall(VecCopy(lmvm->Fprev, Fprev_old));
613 }
614
615 /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
616 PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
617 PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
618
619 /* Test if the updates can be accepted */
620 {
621 Vec sy[2] = {lmvm->Xprev, lmvm->Fprev};
622 PetscScalar stsy[2];
623
624 PetscCall(VecMDot(lmvm->Xprev, 2, sy, stsy));
625 ststmp = stsy[0];
626 curvature = stsy[1];
627 }
628 curvtol = lmvm->eps * PetscRealPart(ststmp);
629
630 if (PetscRealPart(curvature) > curvtol) { /* Update is good, accept it */
631 LMProducts StY = NULL;
632 LMProducts YtH0Y = NULL;
633 Vec StFprev_old = NULL;
634 Vec YtH0Fprev_old = NULL;
635 PetscInt oldest_new, next_new;
636
637 lsb->watchdog = 0;
638 if (cache_StFprev) {
639 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_S, &S, NULL, NULL));
640 if (!lsb->StFprev) PetscCall(LMBasisCreateRow(S, &lsb->StFprev));
641 PetscCall(MatLMVMGetUpdatedProducts(B, LMBASIS_S, LMBASIS_Y, LMBLOCK_UPPER_TRIANGLE, &StY));
642 PetscCall(LMProductsGetNextColumn(StY, &StFprev_old));
643 PetscCall(VecCopy(lsb->StFprev, StFprev_old));
644 }
645 if (cache_YtH0Fprev) {
646 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_H0Y, &H0Y, &H0Y_t, &H0_alpha));
647 if (!lsb->YtH0Fprev) PetscCall(LMBasisCreateRow(H0Y, &lsb->YtH0Fprev));
648 PetscCall(MatLMVMGetUpdatedProducts(B, LMBASIS_Y, H0Y_t, LMBLOCK_UPPER_TRIANGLE, &YtH0Y));
649 PetscCall(LMProductsGetNextColumn(YtH0Y, &YtH0Fprev_old));
650 if (lsb->YtH0Fprev == H0Y->cached_product) {
651 PetscCall(VecCopy(lsb->YtH0Fprev, YtH0Fprev_old));
652 } else {
653 if (next > oldest) {
654 // need to recalculate
655 PetscCall(LMBasisGEMVH(H0Y, oldest, next, 1.0, Fprev_old, 0.0, YtH0Fprev_old));
656 } else {
657 PetscCall(VecZeroEntries(YtH0Fprev_old));
658 }
659 }
660 }
661
662 PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
663 PetscCall(MatLMVMGetRange(B, &oldest_new, &next_new));
664 if (cache_StFprev) {
665 // compute the one new s_i^T F_old value
666 PetscCall(LMBasisGEMVH(S, next, next_new, 1.0, Fprev_old, 0.0, StFprev_old));
667 PetscCall(LMBasisGEMVH(S, oldest_new, next_new, 1.0, F, 0.0, lsb->StFprev));
668 PetscCall(LMBasisSetCachedProduct(S, F, lsb->StFprev));
669 PetscCall(VecAXPBY(StFprev_old, 1.0, -1.0, lsb->StFprev));
670 PetscCall(LMProductsRestoreNextColumn(StY, &StFprev_old));
671 }
672 if (cache_YtH0Fprev) {
673 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_H0Y, &H0Y, &H0Y_t, &H0_alpha));
674 // compute the one new (H_0 y_i)^T F_old value
675 PetscCall(LMBasisGEMVH(H0Y, next, next_new, 1.0, Fprev_old, 0.0, YtH0Fprev_old));
676 PetscCall(LMBasisGEMVH(H0Y, oldest_new, next_new, 1.0, F, 0.0, lsb->YtH0Fprev));
677 PetscCall(LMBasisSetCachedProduct(H0Y, F, lsb->YtH0Fprev));
678 PetscCall(VecAXPBY(YtH0Fprev_old, 1.0, -1.0, lsb->YtH0Fprev));
679 PetscCall(LMProductsRestoreNextColumn(YtH0Y, &YtH0Fprev_old));
680 }
681
682 PetscCall(MatLMVMProductsInsertDiagonalValue(B, LMBASIS_Y, LMBASIS_S, next, PetscRealPart(curvature)));
683 PetscCall(MatLMVMProductsInsertDiagonalValue(B, LMBASIS_S, LMBASIS_Y, next, PetscRealPart(curvature)));
684 PetscCall(MatLMVMProductsInsertDiagonalValue(B, LMBASIS_S, LMBASIS_S, next, ststmp));
685 PetscCall(SymBroydenRescaleUpdate(B, lsb->rescale));
686 } else {
687 /* Update is bad, skip it */
688 PetscCall(PetscInfo(B, "Rejecting update: curvature %g, ||s||^2 %g\n", (double)PetscRealPart(curvature), (double)PetscRealPart(ststmp)));
689 lmvm->nrejects++;
690 lsb->watchdog++;
691 if (cache_StFprev) {
692 // we still need to update the cached product
693 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_S, &S, NULL, NULL));
694 PetscCall(LMBasisGEMVH(S, oldest, next, 1.0, F, 0.0, lsb->StFprev));
695 PetscCall(LMBasisSetCachedProduct(S, F, lsb->StFprev));
696 }
697 if (cache_YtH0Fprev) {
698 // we still need to update the cached product
699 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_H0Y, &H0Y, &H0Y_t, &H0_alpha));
700 PetscCall(LMBasisGEMVH(H0Y, oldest, next, 1.0, F, 0.0, lsb->YtH0Fprev));
701 PetscCall(LMBasisSetCachedProduct(H0Y, F, lsb->StFprev));
702 }
703 }
704 if (cache_StFprev || cache_YtH0Fprev) PetscCall(LMBasisRestoreWorkVec(Y, &Fprev_old));
705 } else {
706 if (cache_StFprev) {
707 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_S, &S, NULL, NULL));
708 if (!lsb->StFprev) PetscCall(LMBasisCreateRow(S, &lsb->StFprev));
709 }
710 if (cache_YtH0Fprev) {
711 MatLMVMBasisType H0Y_t;
712 PetscScalar H0_alpha;
713
714 PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_H0Y, &H0Y, &H0Y_t, &H0_alpha));
715 if (!lsb->YtH0Fprev) PetscCall(LMBasisCreateRow(H0Y, &lsb->YtH0Fprev));
716 }
717 }
718
719 if (lsb->watchdog > lsb->max_seq_rejects) PetscCall(MatLMVMReset(B, PETSC_FALSE));
720
721 /* Save the solution and function to be used in the next update */
722 PetscCall(VecCopy(X, lmvm->Xprev));
723 PetscCall(VecCopy(F, lmvm->Fprev));
724 lmvm->prev_set = PETSC_TRUE;
725 PetscFunctionReturn(PETSC_SUCCESS);
726 }
727
MatCopy_LMVMSymBrdn(Mat B,Mat M,MatStructure str)728 static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
729 {
730 Mat_LMVM *bdata = (Mat_LMVM *)B->data;
731 Mat_SymBrdn *blsb = (Mat_SymBrdn *)bdata->ctx;
732 Mat_LMVM *mdata = (Mat_LMVM *)M->data;
733 Mat_SymBrdn *mlsb = (Mat_SymBrdn *)mdata->ctx;
734
735 PetscFunctionBegin;
736 mlsb->phi_scalar = blsb->phi_scalar;
737 mlsb->psi_scalar = blsb->psi_scalar;
738 mlsb->watchdog = blsb->watchdog;
739 mlsb->max_seq_rejects = blsb->max_seq_rejects;
740 PetscCall(SymBroydenRescaleCopy(blsb->rescale, mlsb->rescale));
741 PetscFunctionReturn(PETSC_SUCCESS);
742 }
743
MatReset_LMVMSymBrdn_Internal(Mat B,MatLMVMResetMode mode)744 static PetscErrorCode MatReset_LMVMSymBrdn_Internal(Mat B, MatLMVMResetMode mode)
745 {
746 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
747 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
748
749 PetscFunctionBegin;
750 if (MatLMVMResetClearsBases(mode)) {
751 for (PetscInt i = 0; i < SYMBROYDEN_BASIS_COUNT; i++) PetscCall(LMBasisDestroy(&lsb->basis[i]));
752 for (PetscInt i = 0; i < SYMBROYDEN_PRODUCTS_COUNT; i++) PetscCall(LMProductsDestroy(&lsb->products[i]));
753 PetscCall(VecDestroy(&lsb->StFprev));
754 PetscCall(VecDestroy(&lsb->YtH0Fprev));
755 } else {
756 for (PetscInt i = 0; i < SYMBROYDEN_BASIS_COUNT; i++) PetscCall(LMBasisReset(lsb->basis[i]));
757 for (PetscInt i = 0; i < SYMBROYDEN_PRODUCTS_COUNT; i++) PetscCall(LMProductsReset(lsb->products[i]));
758 if (lsb->StFprev) PetscCall(VecZeroEntries(lsb->StFprev));
759 if (lsb->YtH0Fprev) PetscCall(VecZeroEntries(lsb->YtH0Fprev));
760 }
761 PetscFunctionReturn(PETSC_SUCCESS);
762 }
763
MatReset_LMVMSymBrdn(Mat B,MatLMVMResetMode mode)764 static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, MatLMVMResetMode mode)
765 {
766 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
767 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
768
769 PetscFunctionBegin;
770 lsb->watchdog = 0;
771 PetscCall(SymBroydenRescaleReset(B, lsb->rescale, mode));
772 PetscCall(MatReset_LMVMSymBrdn_Internal(B, mode));
773 PetscFunctionReturn(PETSC_SUCCESS);
774 }
775
MatDestroy_LMVMSymBrdn(Mat B)776 static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
777 {
778 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
779 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
780
781 PetscFunctionBegin;
782 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenGetPhi_C", NULL));
783 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetPhi_C", NULL));
784 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBadBroydenGetPsi_C", NULL));
785 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBadBroydenSetPsi_C", NULL));
786 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetDelta_C", NULL));
787 PetscCall(SymBroydenRescaleDestroy(&lsb->rescale));
788 PetscCall(MatReset_LMVMSymBrdn_Internal(B, MAT_LMVM_RESET_ALL));
789 PetscCall(PetscFree(lmvm->ctx));
790 PetscCall(MatDestroy_LMVM(B));
791 PetscFunctionReturn(PETSC_SUCCESS);
792 }
793
MatSetUp_LMVMSymBrdn(Mat B)794 static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
795 {
796 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
797 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
798
799 PetscFunctionBegin;
800 PetscCall(MatSetUp_LMVM(B));
801 PetscCall(SymBroydenRescaleInitializeJ0(B, lsb->rescale));
802 PetscFunctionReturn(PETSC_SUCCESS);
803 }
804
MatView_LMVMSymBrdn(Mat B,PetscViewer pv)805 static PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
806 {
807 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
808 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
809 PetscBool isascii;
810
811 PetscFunctionBegin;
812 PetscCall(MatView_LMVM(B, pv));
813 PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
814 if (isascii) {
815 PetscBool is_other;
816
817 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &is_other, MATLMVMBFGS, MATLMVMDFP, ""));
818 if (!is_other) {
819 if (lsb->phi_scalar != PETSC_DETERMINE) PetscCall(PetscViewerASCIIPrintf(pv, "Convex factor phi = %g\n", (double)lsb->phi_scalar));
820 if (lsb->psi_scalar != PETSC_DETERMINE) PetscCall(PetscViewerASCIIPrintf(pv, "Dual convex factor psi = %g\n", (double)lsb->psi_scalar));
821 }
822 }
823 PetscCall(SymBroydenRescaleView(lsb->rescale, pv));
824 PetscFunctionReturn(PETSC_SUCCESS);
825 }
826
MatLMVMSetMultAlgorithm_SymBrdn(Mat B)827 static PetscErrorCode MatLMVMSetMultAlgorithm_SymBrdn(Mat B)
828 {
829 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
830
831 PetscFunctionBegin;
832 switch (lmvm->mult_alg) {
833 case MAT_LMVM_MULT_RECURSIVE:
834 lmvm->ops->mult = MatMult_LMVMSymBrdn_Recursive;
835 lmvm->ops->solve = MatSolve_LMVMSymBrdn_Recursive;
836 break;
837 case MAT_LMVM_MULT_DENSE:
838 case MAT_LMVM_MULT_COMPACT_DENSE:
839 lmvm->ops->mult = MatMult_LMVMSymBrdn_CompactDense;
840 lmvm->ops->solve = MatSolve_LMVMSymBrdn_CompactDense;
841 break;
842 }
843 lmvm->ops->multht = lmvm->ops->mult;
844 lmvm->ops->solveht = lmvm->ops->solve;
845 PetscFunctionReturn(PETSC_SUCCESS);
846 }
847
MatSetFromOptions_LMVMSymBrdn(Mat B,PetscOptionItems PetscOptionsObject)848 static PetscErrorCode MatSetFromOptions_LMVMSymBrdn(Mat B, PetscOptionItems PetscOptionsObject)
849 {
850 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
851 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
852
853 PetscFunctionBegin;
854 PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
855 PetscOptionsHeadBegin(PetscOptionsObject, "Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
856 PetscCall(PetscOptionsReal("-mat_lmvm_phi", "convex ratio between BFGS and DFP components of the update", "", lsb->phi_scalar, &lsb->phi_scalar, NULL));
857 PetscCheck(lsb->phi_scalar >= 0.0 && lsb->phi_scalar <= 1.0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the update formula cannot be outside the range of [0, 1]");
858 PetscCall(SymBroydenRescaleSetFromOptions(B, lsb->rescale, PetscOptionsObject));
859 PetscOptionsHeadEnd();
860 PetscFunctionReturn(PETSC_SUCCESS);
861 }
862
863 /*@
864 MatLMVMSymBroydenGetPhi - Get the phi parameter for a Broyden class quasi-Newton update matrix
865
866 Input Parameter:
867 . B - The matrix
868
869 Output Parameter:
870 . phi - a number defining an update that is an affine combination of the BFGS update (phi = 0) and DFP update (phi = 1)
871
872 Level: advanced
873
874 Note:
875 If `B` does not have a constant value of `phi` for all iterations this will
876 return `phi` = `PETSC_DETERMINE` = -1, a negative value that `phi` cannot
877 attain for a valid general Broyden update.
878 This is the case if `B` is a `MATLMVMSYMBADBROYDEN`, where `phi`'s dual value
879 `psi` is constant and `phi` changes from iteration to iteration.
880
881 .seealso: [](ch_ksp),
882 `MATLMVMSYMBROYDEN`, `MATLMVMSYMBADBROYDEN`,
883 `MATLMVMDFP`, `MATLMVMBFGS`,
884 `MatLMVMSymBroydenSetPhi()`,
885 `MatLMVMSymBadBroydenGetPsi()`, `MatLMVMSymBadBroydenSetPsi()`
886 @*/
MatLMVMSymBroydenGetPhi(Mat B,PetscReal * phi)887 PetscErrorCode MatLMVMSymBroydenGetPhi(Mat B, PetscReal *phi)
888 {
889 PetscFunctionBegin;
890 *phi = PETSC_DETERMINE;
891 PetscUseMethod(B, "MatLMVMSymBroydenGetPhi_C", (Mat, PetscReal *), (B, phi));
892 PetscFunctionReturn(PETSC_SUCCESS);
893 }
894
MatLMVMSymBroydenGetPhi_SymBrdn(Mat B,PetscReal * phi)895 static PetscErrorCode MatLMVMSymBroydenGetPhi_SymBrdn(Mat B, PetscReal *phi)
896 {
897 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
898 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
899
900 PetscFunctionBegin;
901 *phi = lsb->phi_scalar;
902 PetscFunctionReturn(PETSC_SUCCESS);
903 }
904
905 /*@
906 MatLMVMSymBroydenSetPhi - Get the phi parameter for a Broyden class quasi-Newton update matrix
907
908 Input Parameters:
909 + B - The matrix
910 - phi - a number defining an update that is a convex combination of the BFGS update (phi = 0) and DFP update (phi = 1)
911
912 Level: advanced
913
914 Note:
915 If `B` cannot have a constant value of `phi` for all iterations this will be ignored.
916 This is the case if `B` is a `MATLMVMSYMBADBROYDEN`, where `phi`'s dual value
917 `psi` is constant and `phi` changes from iteration to iteration.
918
919 .seealso: [](ch_ksp),
920 `MATLMVMSYMBROYDEN`, `MATLMVMSYMBADBROYDEN`,
921 `MATLMVMDFP`, `MATLMVMBFGS`,
922 `MatLMVMSymBroydenGetPhi()`,
923 `MatLMVMSymBadBroydenGetPsi()`, `MatLMVMSymBadBroydenSetPsi()`
924 @*/
MatLMVMSymBroydenSetPhi(Mat B,PetscReal phi)925 PetscErrorCode MatLMVMSymBroydenSetPhi(Mat B, PetscReal phi)
926 {
927 PetscFunctionBegin;
928 PetscTryMethod(B, "MatLMVMSymBroydenSetPhi_C", (Mat, PetscReal), (B, phi));
929 PetscFunctionReturn(PETSC_SUCCESS);
930 }
931
MatLMVMSymBroydenSetPhi_SymBrdn(Mat B,PetscReal phi)932 static PetscErrorCode MatLMVMSymBroydenSetPhi_SymBrdn(Mat B, PetscReal phi)
933 {
934 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
935 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
936
937 PetscFunctionBegin;
938 lsb->phi_scalar = phi;
939 PetscFunctionReturn(PETSC_SUCCESS);
940 }
941
942 /*@
943 MatLMVMSymBadBroydenGetPsi - Get the psi parameter for a Broyden class quasi-Newton update matrix
944
945 Input Parameter:
946 . B - The matrix
947
948 Output Parameter:
949 . psi - a number defining an update that is an affine combination of the BFGS update (psi = 1) and DFP update (psi = 0)
950
951 Level: advanced
952
953 Note:
954 If B does not have a constant value of `psi` for all iterations this will
955 return `psi` = `PETSC_DETERMINE` = -1, a negative value that `psi` cannot
956 attain for a valid general Broyden update.
957 This is the case if `B` is a `MATLMVMSYMBROYDEN`, where `psi`'s dual value
958 `phi` is constant and `psi` changes from iteration to iteration.
959
960 .seealso: [](ch_ksp),
961 `MATLMVMSYMBROYDEN`, `MATLMVMSYMBADBROYDEN`,
962 `MATLMVMDFP`, `MATLMVMBFGS`,
963 `MatLMVMSymBadBroydenSetPsi()`,
964 `MatLMVMSymBroydenGetPhi()`, `MatLMVMSymBroydenSetPhi()`
965 @*/
MatLMVMSymBadBroydenGetPsi(Mat B,PetscReal * psi)966 PetscErrorCode MatLMVMSymBadBroydenGetPsi(Mat B, PetscReal *psi)
967 {
968 PetscFunctionBegin;
969 *psi = PETSC_DETERMINE;
970 PetscTryMethod(B, "MatLMVMSymBadBroydenGetPsi_C", (Mat, PetscReal *), (B, psi));
971 PetscFunctionReturn(PETSC_SUCCESS);
972 }
973
MatLMVMSymBadBroydenGetPsi_SymBrdn(Mat B,PetscReal * psi)974 static PetscErrorCode MatLMVMSymBadBroydenGetPsi_SymBrdn(Mat B, PetscReal *psi)
975 {
976 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
977 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
978
979 PetscFunctionBegin;
980 *psi = lsb->psi_scalar;
981 PetscFunctionReturn(PETSC_SUCCESS);
982 }
983
984 /*@
985 MatLMVMSymBadBroydenSetPsi - Get the psi parameter for a Broyden class quasi-Newton update matrix
986
987 Input Parameters:
988 + B - The matrix
989 - psi - a number defining an update that is a convex combination of the BFGS update (psi = 1) and DFP update (psi = 0)
990
991 Level: developer
992
993 Note:
994 If `B` cannot have a constant value of `psi` for all iterations this will
995 be ignored.
996 This is the case if `B` is a `MATLMVMSYMBROYDEN`, where `psi`'s dual value
997 `phi` is constant and `psi` changes from iteration to iteration.
998
999 .seealso: [](ch_ksp),
1000 `MATLMVMSYMBROYDEN`, `MATLMVMSYMBADBROYDEN`,
1001 `MATLMVMDFP`, `MATLMVMBFGS`,
1002 `MatLMVMSymBadBroydenGetPsi()`,
1003 `MatLMVMSymBroydenGetPhi()`, `MatLMVMSymBroydenSetPhi()`
1004 @*/
MatLMVMSymBadBroydenSetPsi(Mat B,PetscReal psi)1005 PetscErrorCode MatLMVMSymBadBroydenSetPsi(Mat B, PetscReal psi)
1006 {
1007 PetscFunctionBegin;
1008 PetscTryMethod(B, "MatLMVMSymBadBroydenSetPsi_C", (Mat, PetscReal), (B, psi));
1009 PetscFunctionReturn(PETSC_SUCCESS);
1010 }
1011
MatLMVMSymBadBroydenSetPsi_SymBrdn(Mat B,PetscReal psi)1012 static PetscErrorCode MatLMVMSymBadBroydenSetPsi_SymBrdn(Mat B, PetscReal psi)
1013 {
1014 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1015 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
1016
1017 PetscFunctionBegin;
1018 lsb->psi_scalar = psi;
1019 PetscFunctionReturn(PETSC_SUCCESS);
1020 }
1021
MatLMVMSymBroydenSetDelta_SymBrdn(Mat B,PetscScalar delta)1022 static PetscErrorCode MatLMVMSymBroydenSetDelta_SymBrdn(Mat B, PetscScalar delta)
1023 {
1024 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1025 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
1026
1027 PetscFunctionBegin;
1028 PetscCall(SymBroydenRescaleSetDelta(B, lsb->rescale, PetscAbsReal(PetscRealPart(delta))));
1029 PetscFunctionReturn(PETSC_SUCCESS);
1030 }
1031
MatCreate_LMVMSymBrdn(Mat B)1032 PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
1033 {
1034 Mat_LMVM *lmvm;
1035 Mat_SymBrdn *lsb;
1036
1037 PetscFunctionBegin;
1038 PetscCall(MatCreate_LMVM(B));
1039 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN));
1040 PetscCall(MatSetOption(B, MAT_HERMITIAN, PETSC_TRUE));
1041 PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE)); // TODO: change to HPD when available
1042 PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
1043 B->ops->view = MatView_LMVMSymBrdn;
1044 B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
1045 B->ops->setup = MatSetUp_LMVMSymBrdn;
1046 B->ops->destroy = MatDestroy_LMVMSymBrdn;
1047
1048 lmvm = (Mat_LMVM *)B->data;
1049 lmvm->ops->reset = MatReset_LMVMSymBrdn;
1050 lmvm->ops->update = MatUpdate_LMVMSymBrdn;
1051 lmvm->ops->copy = MatCopy_LMVMSymBrdn;
1052 lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_SymBrdn;
1053 lmvm->cache_gradient_products = PETSC_TRUE;
1054 PetscCall(MatLMVMSetMultAlgorithm_SymBrdn(B));
1055
1056 PetscCall(PetscNew(&lsb));
1057 lmvm->ctx = (void *)lsb;
1058 lsb->phi_scalar = 0.125;
1059 lsb->psi_scalar = PETSC_DETERMINE;
1060 lsb->watchdog = 0;
1061 lsb->max_seq_rejects = lmvm->m / 2;
1062
1063 PetscCall(SymBroydenRescaleCreate(&lsb->rescale));
1064 lsb->rescale->theta = lsb->phi_scalar;
1065 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenGetPhi_C", MatLMVMSymBroydenGetPhi_SymBrdn));
1066 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetPhi_C", MatLMVMSymBroydenSetPhi_SymBrdn));
1067 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBadBroydenGetPsi_C", MatLMVMSymBadBroydenGetPsi_SymBrdn));
1068 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBadBroydenSetPsi_C", MatLMVMSymBadBroydenSetPsi_SymBrdn));
1069 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatLMVMSymBroydenSetDelta_C", MatLMVMSymBroydenSetDelta_SymBrdn));
1070 PetscFunctionReturn(PETSC_SUCCESS);
1071 }
1072
1073 /*@
1074 MatLMVMSymBroydenSetDelta - Sets the starting value for the diagonal scaling vector computed
1075 in the SymBrdn approximations (also works for BFGS and DFP).
1076
1077 Input Parameters:
1078 + B - `MATLMVM` matrix
1079 - delta - initial value for diagonal scaling
1080
1081 Level: intermediate
1082
1083 .seealso: [](ch_ksp), `MATLMVMSYMBROYDEN`
1084 @*/
MatLMVMSymBroydenSetDelta(Mat B,PetscScalar delta)1085 PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
1086 {
1087 PetscFunctionBegin;
1088 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
1089 PetscTryMethod(B, "MatLMVMSymBroydenSetDelta_C", (Mat, PetscScalar), (B, delta));
1090 PetscFunctionReturn(PETSC_SUCCESS);
1091 }
1092
1093 /*@
1094 MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.
1095
1096 Input Parameters:
1097 + B - the `MATLMVM` matrix
1098 - stype - scale type, see `MatLMVMSymBroydenScaleType`
1099
1100 Options Database Key:
1101 . -mat_lmvm_scale_type <none,scalar,diagonal> - set the scaling type
1102
1103 Level: intermediate
1104
1105 MatLMVMSymBrdnScaleTypes\:
1106 + `MAT_LMVM_SYMBROYDEN_SCALE_NONE` - use whatever initial Hessian is already there (will be the identity if the user does nothing)
1107 . `MAT_LMVM_SYMBROYDEN_SCALE_SCALAR` - use the Shanno scalar as the initial Hessian
1108 . `MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL` - use a diagonalized BFGS update as the initial Hessian
1109 . `MAT_LMVM_SYMBROYDEN_SCALE_USER` - same as `MAT_LMVM_SYMBROYDEN_NONE`
1110 - `MAT_LMVM_SYMBROYDEN_SCALE_DECIDE` - let PETSc decide
1111
1112 .seealso: [](ch_ksp), `MATLMVMSYMBROYDEN`, `MatCreateLMVMSymBroyden()`, `MatLMVMSymBroydenScaleType`
1113 @*/
MatLMVMSymBroydenSetScaleType(Mat B,MatLMVMSymBroydenScaleType stype)1114 PetscErrorCode MatLMVMSymBroydenSetScaleType(Mat B, MatLMVMSymBroydenScaleType stype)
1115 {
1116 Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
1117 Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
1118
1119 PetscFunctionBegin;
1120 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
1121 PetscCall(SymBroydenRescaleSetType(lsb->rescale, stype));
1122 PetscFunctionReturn(PETSC_SUCCESS);
1123 }
1124
1125 /*@
1126 MatCreateLMVMSymBroyden - Creates a limited-memory Symmetric Broyden-type matrix used
1127 for approximating Jacobians.
1128
1129 Collective
1130
1131 Input Parameters:
1132 + comm - MPI communicator, set to `PETSC_COMM_SELF`
1133 . n - number of local rows for storage vectors
1134 - N - global size of the storage vectors
1135
1136 Output Parameter:
1137 . B - the matrix
1138
1139 Options Database Keys:
1140 + -mat_lmvm_hist_size - the number of history vectors to keep
1141 . -mat_lmvm_phi - convex ratio between BFGS and DFP components of the update
1142 . -mat_lmvm_scale_type - type of scaling applied to J0 (none, scalar, diagonal)
1143 . -mat_lmvm_mult_algorithm - the algorithm to use for multiplication (recursive, dense, compact_dense)
1144 . -mat_lmvm_cache_J0_products - whether products between the base Jacobian J0 and history vectors should be cached or recomputed
1145 . -mat_lmvm_eps - (developer) numerical zero tolerance for testing when an update should be skipped
1146 . -mat_lmvm_debug - (developer) perform internal debugging checks
1147 . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
1148 . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
1149 . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
1150 . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
1151 - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
1152
1153 Level: intermediate
1154
1155 Notes:
1156 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()` paradigm instead of this
1157 routine directly.
1158
1159 L-SymBrdn is a convex combination of L-DFP and L-BFGS such that $B = (1 - \phi)B_{\text{BFGS}} + \phi B_{\text{DFP}}$.
1160 The combination factor $\phi$ is restricted to the range $[0, 1]$, where the L-SymBrdn matrix is guaranteed to be
1161 symmetric positive-definite.
1162
1163 To use the L-SymBrdn matrix with other vector types, the matrix must be created using `MatCreate()` and `MatSetType()`,
1164 followed by `MatLMVMAllocate()`. This ensures that the internal storage and work vectors are duplicated from the
1165 correct type of vector.
1166
1167 .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMSYMBROYDEN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
1168 `MatCreateLMVMBFGS()`, `MatCreateLMVMBroyden()`, `MatCreateLMVMBadBroyden()`
1169 @*/
MatCreateLMVMSymBroyden(MPI_Comm comm,PetscInt n,PetscInt N,Mat * B)1170 PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
1171 {
1172 PetscFunctionBegin;
1173 PetscCall(KSPInitializePackage());
1174 PetscCall(MatCreate(comm, B));
1175 PetscCall(MatSetSizes(*B, n, n, N, N));
1176 PetscCall(MatSetType(*B, MATLMVMSYMBROYDEN));
1177 PetscCall(MatSetUp(*B));
1178 PetscFunctionReturn(PETSC_SUCCESS);
1179 }
1180