1 #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h> /*I "petscksp.h" I*/
2 #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
3
4 const char *const MatLMVMSymBroydenScaleTypes[] = {"NONE","SCALAR","DIAGONAL","USER","MatLMVMSymBrdnScaleType","MAT_LMVM_SYMBROYDEN_SCALING_",NULL};
5
6 /*------------------------------------------------------------*/
7
8 /*
9 The solution method below is the matrix-free implementation of
10 Equation 8.6a in Dennis and More "Quasi-Newton Methods, Motivation
11 and Theory" (https://epubs.siam.org/doi/abs/10.1137/1019005).
12
13 Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
14 the matrix is updated with a new (S[i], Y[i]) pair. This allows
15 repeated calls of MatSolve without incurring redundant computation.
16
17 dX <- J0^{-1} * F
18
19 for i=0,1,2,...,k
20 # Q[i] = (B_i)^T{-1} Y[i]
21
22 rho = 1.0 / (Y[i]^T S[i])
23 alpha = rho * (S[i]^T F)
24 zeta = 1.0 / (Y[i]^T Q[i])
25 gamma = zeta * (Y[i]^T dX)
26
27 dX <- dX - (gamma * Q[i]) + (alpha * Y[i])
28 W <- (rho * S[i]) - (zeta * Q[i])
29 dX <- dX + (psi[i] * (Y[i]^T Q[i]) * (W^T F) * W)
30 end
31 */
MatSolve_LMVMSymBrdn(Mat B,Vec F,Vec dX)32 static PetscErrorCode MatSolve_LMVMSymBrdn(Mat B, Vec F, Vec dX)
33 {
34 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
35 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
36 PetscErrorCode ierr;
37 PetscInt i, j;
38 PetscReal numer;
39 PetscScalar sjtpi, yjtsi, wtsi, yjtqi, sjtyi, wtyi, ytx, stf, wtf, stp, ytq;
40
41 PetscFunctionBegin;
42 /* Efficient shortcuts for pure BFGS and pure DFP configurations */
43 if (lsb->phi == 0.0) {
44 ierr = MatSolve_LMVMBFGS(B, F, dX);CHKERRQ(ierr);
45 PetscFunctionReturn(0);
46 }
47 if (lsb->phi == 1.0) {
48 ierr = MatSolve_LMVMDFP(B, F, dX);CHKERRQ(ierr);
49 PetscFunctionReturn(0);
50 }
51
52 VecCheckSameSize(F, 2, dX, 3);
53 VecCheckMatCompatible(B, dX, 3, F, 2);
54
55 if (lsb->needP) {
56 /* Start the loop for (P[k] = (B_k) * S[k]) */
57 for (i = 0; i <= lmvm->k; ++i) {
58 ierr = MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);CHKERRQ(ierr);
59 for (j = 0; j <= i-1; ++j) {
60 /* Compute the necessary dot products */
61 ierr = VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);CHKERRQ(ierr);
62 ierr = VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);CHKERRQ(ierr);
63 ierr = VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);CHKERRQ(ierr);
64 ierr = VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);CHKERRQ(ierr);
65 /* Compute the pure BFGS component of the forward product */
66 ierr = VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi)/lsb->stp[j], PetscRealPart(yjtsi)/lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]);CHKERRQ(ierr);
67 /* Tack on the convexly scaled extras to the forward product */
68 if (lsb->phi > 0.0) {
69 ierr = VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);CHKERRQ(ierr);
70 ierr = VecDot(lsb->work, lmvm->S[i], &wtsi);CHKERRQ(ierr);
71 ierr = VecAXPY(lsb->P[i], lsb->phi*lsb->stp[j]*PetscRealPart(wtsi), lsb->work);CHKERRQ(ierr);
72 }
73 }
74 ierr = VecDot(lmvm->S[i], lsb->P[i], &stp);CHKERRQ(ierr);
75 lsb->stp[i] = PetscRealPart(stp);
76 }
77 lsb->needP = PETSC_FALSE;
78 }
79 if (lsb->needQ) {
80 /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
81 for (i = 0; i <= lmvm->k; ++i) {
82 ierr = MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]);CHKERRQ(ierr);
83 for (j = 0; j <= i-1; ++j) {
84 /* Compute the necessary dot products */
85 ierr = VecDotBegin(lmvm->Y[j], lsb->Q[i], &yjtqi);CHKERRQ(ierr);
86 ierr = VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);CHKERRQ(ierr);
87 ierr = VecDotEnd(lmvm->Y[j], lsb->Q[i], &yjtqi);CHKERRQ(ierr);
88 ierr = VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);CHKERRQ(ierr);
89 /* Compute the pure DFP component of the inverse application*/
90 ierr = VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi)/lsb->ytq[j], PetscRealPart(sjtyi)/lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]);CHKERRQ(ierr);
91 /* Tack on the convexly scaled extras to the inverse application*/
92 if (lsb->psi[j] > 0.0) {
93 ierr = VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]);CHKERRQ(ierr);
94 ierr = VecDot(lsb->work, lmvm->Y[i], &wtyi);CHKERRQ(ierr);
95 ierr = VecAXPY(lsb->Q[i], lsb->psi[j]*lsb->ytq[j]*PetscRealPart(wtyi), lsb->work);CHKERRQ(ierr);
96 }
97 }
98 ierr = VecDot(lmvm->Y[i], lsb->Q[i], &ytq);CHKERRQ(ierr);
99 lsb->ytq[i] = PetscRealPart(ytq);
100 if (lsb->phi == 1.0) {
101 lsb->psi[i] = 0.0;
102 } else if (lsb->phi == 0.0) {
103 lsb->psi[i] = 1.0;
104 } else {
105 numer = (1.0 - lsb->phi)*lsb->yts[i]*lsb->yts[i];
106 lsb->psi[i] = numer / (numer + (lsb->phi*lsb->ytq[i]*lsb->stp[i]));
107 }
108 }
109 lsb->needQ = PETSC_FALSE;
110 }
111
112 /* Start the outer iterations for ((B^{-1}) * dX) */
113 ierr = MatSymBrdnApplyJ0Inv(B, F, dX);CHKERRQ(ierr);
114 for (i = 0; i <= lmvm->k; ++i) {
115 /* Compute the necessary dot products -- store yTs and yTp for inner iterations later */
116 ierr = VecDotBegin(lmvm->Y[i], dX, &ytx);CHKERRQ(ierr);
117 ierr = VecDotBegin(lmvm->S[i], F, &stf);CHKERRQ(ierr);
118 ierr = VecDotEnd(lmvm->Y[i], dX, &ytx);CHKERRQ(ierr);
119 ierr = VecDotEnd(lmvm->S[i], F, &stf);CHKERRQ(ierr);
120 /* Compute the pure DFP component */
121 ierr = VecAXPBYPCZ(dX, -PetscRealPart(ytx)/lsb->ytq[i], PetscRealPart(stf)/lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]);CHKERRQ(ierr);
122 /* Tack on the convexly scaled extras */
123 ierr = VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]);CHKERRQ(ierr);
124 ierr = VecDot(lsb->work, F, &wtf);CHKERRQ(ierr);
125 ierr = VecAXPY(dX, lsb->psi[i]*lsb->ytq[i]*PetscRealPart(wtf), lsb->work);CHKERRQ(ierr);
126 }
127
128 PetscFunctionReturn(0);
129 }
130
131 /*------------------------------------------------------------*/
132
133 /*
134 The forward-product below is the matrix-free implementation of
135 Equation 16 in Dennis and Wolkowicz "Sizing and Least Change Secant
136 Methods" (http://www.caam.rice.edu/caam/trs/90/TR90-05.pdf).
137
138 P[i] = (B_i)*S[i] terms are computed ahead of time whenever
139 the matrix is updated with a new (S[i], Y[i]) pair. This allows
140 repeated calls of MatMult inside KSP solvers without unnecessarily
141 recomputing P[i] terms in expensive nested-loops.
142
143 Z <- J0 * X
144
145 for i=0,1,2,...,k
146 # P[i] = (B_k) * S[i]
147
148 rho = 1.0 / (Y[i]^T S[i])
149 alpha = rho * (Y[i]^T F)
150 zeta = 1.0 / (S[i]^T P[i])
151 gamma = zeta * (S[i]^T dX)
152
153 dX <- dX - (gamma * P[i]) + (alpha * S[i])
154 W <- (rho * Y[i]) - (zeta * P[i])
155 dX <- dX + (phi * (S[i]^T P[i]) * (W^T F) * W)
156 end
157 */
MatMult_LMVMSymBrdn(Mat B,Vec X,Vec Z)158 static PetscErrorCode MatMult_LMVMSymBrdn(Mat B, Vec X, Vec Z)
159 {
160 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
161 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
162 PetscErrorCode ierr;
163 PetscInt i, j;
164 PetscScalar sjtpi, yjtsi, wtsi, stz, ytx, wtx, stp;
165
166
167 PetscFunctionBegin;
168 /* Efficient shortcuts for pure BFGS and pure DFP configurations */
169 if (lsb->phi == 0.0) {
170 ierr = MatMult_LMVMBFGS(B, X, Z);CHKERRQ(ierr);
171 PetscFunctionReturn(0);
172 }
173 if (lsb->phi == 1.0) {
174 ierr = MatMult_LMVMDFP(B, X, Z);CHKERRQ(ierr);
175 PetscFunctionReturn(0);
176 }
177
178 VecCheckSameSize(X, 2, Z, 3);
179 VecCheckMatCompatible(B, X, 2, Z, 3);
180
181 if (lsb->needP) {
182 /* Start the loop for (P[k] = (B_k) * S[k]) */
183 for (i = 0; i <= lmvm->k; ++i) {
184 ierr = MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);CHKERRQ(ierr);
185 for (j = 0; j <= i-1; ++j) {
186 /* Compute the necessary dot products */
187 ierr = VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);CHKERRQ(ierr);
188 ierr = VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);CHKERRQ(ierr);
189 ierr = VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);CHKERRQ(ierr);
190 ierr = VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);CHKERRQ(ierr);
191 /* Compute the pure BFGS component of the forward product */
192 ierr = VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi)/lsb->stp[j], PetscRealPart(yjtsi)/lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]);CHKERRQ(ierr);
193 /* Tack on the convexly scaled extras to the forward product */
194 if (lsb->phi > 0.0) {
195 ierr = VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);CHKERRQ(ierr);
196 ierr = VecDot(lsb->work, lmvm->S[i], &wtsi);CHKERRQ(ierr);
197 ierr = VecAXPY(lsb->P[i], lsb->phi*lsb->stp[j]*PetscRealPart(wtsi), lsb->work);CHKERRQ(ierr);
198 }
199 }
200 ierr = VecDot(lmvm->S[i], lsb->P[i], &stp);CHKERRQ(ierr);
201 lsb->stp[i] = PetscRealPart(stp);
202 }
203 lsb->needP = PETSC_FALSE;
204 }
205
206 /* Start the outer iterations for (B * X) */
207 ierr = MatSymBrdnApplyJ0Fwd(B, X, Z);CHKERRQ(ierr);
208 for (i = 0; i <= lmvm->k; ++i) {
209 /* Compute the necessary dot products */
210 ierr = VecDotBegin(lmvm->S[i], Z, &stz);CHKERRQ(ierr);
211 ierr = VecDotBegin(lmvm->Y[i], X, &ytx);CHKERRQ(ierr);
212 ierr = VecDotEnd(lmvm->S[i], Z, &stz);CHKERRQ(ierr);
213 ierr = VecDotEnd(lmvm->Y[i], X, &ytx);CHKERRQ(ierr);
214 /* Compute the pure BFGS component */
215 ierr = VecAXPBYPCZ(Z, -PetscRealPart(stz)/lsb->stp[i], PetscRealPart(ytx)/lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]);CHKERRQ(ierr);
216 /* Tack on the convexly scaled extras */
217 ierr = VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]);CHKERRQ(ierr);
218 ierr = VecDot(lsb->work, X, &wtx);CHKERRQ(ierr);
219 ierr = VecAXPY(Z, lsb->phi*lsb->stp[i]*PetscRealPart(wtx), lsb->work);CHKERRQ(ierr);
220 }
221 PetscFunctionReturn(0);
222 }
223
224 /*------------------------------------------------------------*/
225
MatUpdate_LMVMSymBrdn(Mat B,Vec X,Vec F)226 static PetscErrorCode MatUpdate_LMVMSymBrdn(Mat B, Vec X, Vec F)
227 {
228 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
229 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
230 Mat_LMVM *dbase;
231 Mat_DiagBrdn *dctx;
232 PetscErrorCode ierr;
233 PetscInt old_k, i;
234 PetscReal curvtol;
235 PetscScalar curvature, ytytmp, ststmp;
236
237 PetscFunctionBegin;
238 if (!lmvm->m) PetscFunctionReturn(0);
239 if (lmvm->prev_set) {
240 /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
241 ierr = VecAYPX(lmvm->Xprev, -1.0, X);CHKERRQ(ierr);
242 ierr = VecAYPX(lmvm->Fprev, -1.0, F);CHKERRQ(ierr);
243 /* Test if the updates can be accepted */
244 ierr = VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);CHKERRQ(ierr);
245 ierr = VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);CHKERRQ(ierr);
246 ierr = VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);CHKERRQ(ierr);
247 ierr = VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);CHKERRQ(ierr);
248 if (PetscRealPart(ststmp) < lmvm->eps) {
249 curvtol = 0.0;
250 } else {
251 curvtol = lmvm->eps * PetscRealPart(ststmp);
252 }
253 if (PetscRealPart(curvature) > curvtol) {
254 /* Update is good, accept it */
255 lsb->watchdog = 0;
256 lsb->needP = lsb->needQ = PETSC_TRUE;
257 old_k = lmvm->k;
258 ierr = MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);CHKERRQ(ierr);
259 /* If we hit the memory limit, shift the yts, yty and sts arrays */
260 if (old_k == lmvm->k) {
261 for (i = 0; i <= lmvm->k-1; ++i) {
262 lsb->yts[i] = lsb->yts[i+1];
263 lsb->yty[i] = lsb->yty[i+1];
264 lsb->sts[i] = lsb->sts[i+1];
265 }
266 }
267 /* Update history of useful scalars */
268 ierr = VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp);CHKERRQ(ierr);
269 lsb->yts[lmvm->k] = PetscRealPart(curvature);
270 lsb->yty[lmvm->k] = PetscRealPart(ytytmp);
271 lsb->sts[lmvm->k] = PetscRealPart(ststmp);
272 /* Compute the scalar scale if necessary */
273 if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) {
274 ierr = MatSymBrdnComputeJ0Scalar(B);CHKERRQ(ierr);
275 }
276 } else {
277 /* Update is bad, skip it */
278 ++lmvm->nrejects;
279 ++lsb->watchdog;
280 }
281 } else {
282 switch (lsb->scale_type) {
283 case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
284 dbase = (Mat_LMVM*)lsb->D->data;
285 dctx = (Mat_DiagBrdn*)dbase->ctx;
286 ierr = VecSet(dctx->invD, lsb->delta);CHKERRQ(ierr);
287 break;
288 case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
289 lsb->sigma = lsb->delta;
290 break;
291 case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
292 lsb->sigma = 1.0;
293 break;
294 default:
295 break;
296 }
297 }
298
299 /* Update the scaling */
300 if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
301 ierr = MatLMVMUpdate(lsb->D, X, F);CHKERRQ(ierr);
302 }
303
304 if (lsb->watchdog > lsb->max_seq_rejects) {
305 ierr = MatLMVMReset(B, PETSC_FALSE);CHKERRQ(ierr);
306 if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
307 ierr = MatLMVMReset(lsb->D, PETSC_FALSE);CHKERRQ(ierr);
308 }
309 }
310
311 /* Save the solution and function to be used in the next update */
312 ierr = VecCopy(X, lmvm->Xprev);CHKERRQ(ierr);
313 ierr = VecCopy(F, lmvm->Fprev);CHKERRQ(ierr);
314 lmvm->prev_set = PETSC_TRUE;
315 PetscFunctionReturn(0);
316 }
317
318 /*------------------------------------------------------------*/
319
MatCopy_LMVMSymBrdn(Mat B,Mat M,MatStructure str)320 static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
321 {
322 Mat_LMVM *bdata = (Mat_LMVM*)B->data;
323 Mat_SymBrdn *blsb = (Mat_SymBrdn*)bdata->ctx;
324 Mat_LMVM *mdata = (Mat_LMVM*)M->data;
325 Mat_SymBrdn *mlsb = (Mat_SymBrdn*)mdata->ctx;
326 PetscErrorCode ierr;
327 PetscInt i;
328
329 PetscFunctionBegin;
330 mlsb->phi = blsb->phi;
331 mlsb->needP = blsb->needP;
332 mlsb->needQ = blsb->needQ;
333 for (i=0; i<=bdata->k; ++i) {
334 mlsb->stp[i] = blsb->stp[i];
335 mlsb->ytq[i] = blsb->ytq[i];
336 mlsb->yts[i] = blsb->yts[i];
337 mlsb->psi[i] = blsb->psi[i];
338 ierr = VecCopy(blsb->P[i], mlsb->P[i]);CHKERRQ(ierr);
339 ierr = VecCopy(blsb->Q[i], mlsb->Q[i]);CHKERRQ(ierr);
340 }
341 mlsb->scale_type = blsb->scale_type;
342 mlsb->alpha = blsb->alpha;
343 mlsb->beta = blsb->beta;
344 mlsb->rho = blsb->rho;
345 mlsb->delta = blsb->delta;
346 mlsb->sigma_hist = blsb->sigma_hist;
347 mlsb->watchdog = blsb->watchdog;
348 mlsb->max_seq_rejects = blsb->max_seq_rejects;
349 switch (blsb->scale_type) {
350 case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
351 mlsb->sigma = blsb->sigma;
352 break;
353 case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
354 ierr = MatCopy(blsb->D, mlsb->D, SAME_NONZERO_PATTERN);CHKERRQ(ierr);
355 break;
356 case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
357 mlsb->sigma = 1.0;
358 break;
359 default:
360 break;
361 }
362 PetscFunctionReturn(0);
363 }
364
365 /*------------------------------------------------------------*/
366
MatReset_LMVMSymBrdn(Mat B,PetscBool destructive)367 static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
368 {
369 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
370 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
371 Mat_LMVM *dbase;
372 Mat_DiagBrdn *dctx;
373 PetscErrorCode ierr;
374
375 PetscFunctionBegin;
376 lsb->watchdog = 0;
377 lsb->needP = lsb->needQ = PETSC_TRUE;
378 if (lsb->allocated) {
379 if (destructive) {
380 ierr = VecDestroy(&lsb->work);CHKERRQ(ierr);
381 ierr = PetscFree5(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts);CHKERRQ(ierr);
382 ierr = PetscFree(lsb->psi);CHKERRQ(ierr);
383 ierr = VecDestroyVecs(lmvm->m, &lsb->P);CHKERRQ(ierr);
384 ierr = VecDestroyVecs(lmvm->m, &lsb->Q);CHKERRQ(ierr);
385 switch (lsb->scale_type) {
386 case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
387 ierr = MatLMVMReset(lsb->D, PETSC_TRUE);CHKERRQ(ierr);
388 break;
389 default:
390 break;
391 }
392 lsb->allocated = PETSC_FALSE;
393 } else {
394 ierr = PetscMemzero(lsb->psi, lmvm->m);CHKERRQ(ierr);
395 switch (lsb->scale_type) {
396 case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
397 lsb->sigma = lsb->delta;
398 break;
399 case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
400 ierr = MatLMVMReset(lsb->D, PETSC_FALSE);CHKERRQ(ierr);
401 dbase = (Mat_LMVM*)lsb->D->data;
402 dctx = (Mat_DiagBrdn*)dbase->ctx;
403 ierr = VecSet(dctx->invD, lsb->delta);CHKERRQ(ierr);
404 break;
405 case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
406 lsb->sigma = 1.0;
407 break;
408 default:
409 break;
410 }
411 }
412 }
413 ierr = MatReset_LMVM(B, destructive);CHKERRQ(ierr);
414 PetscFunctionReturn(0);
415 }
416
417 /*------------------------------------------------------------*/
418
MatAllocate_LMVMSymBrdn(Mat B,Vec X,Vec F)419 static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
420 {
421 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
422 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
423 PetscErrorCode ierr;
424
425 PetscFunctionBegin;
426 ierr = MatAllocate_LMVM(B, X, F);CHKERRQ(ierr);
427 if (!lsb->allocated) {
428 ierr = VecDuplicate(X, &lsb->work);CHKERRQ(ierr);
429 if (lmvm->m > 0) {
430 ierr = PetscMalloc5(lmvm->m,&lsb->stp,lmvm->m,&lsb->ytq,lmvm->m,&lsb->yts,lmvm->m,&lsb->yty,lmvm->m,&lsb->sts);CHKERRQ(ierr);
431 ierr = PetscCalloc1(lmvm->m,&lsb->psi);CHKERRQ(ierr);
432 ierr = VecDuplicateVecs(X, lmvm->m, &lsb->P);CHKERRQ(ierr);
433 ierr = VecDuplicateVecs(X, lmvm->m, &lsb->Q);CHKERRQ(ierr);
434 }
435 switch (lsb->scale_type) {
436 case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
437 ierr = MatLMVMAllocate(lsb->D, X, F);CHKERRQ(ierr);
438 break;
439 default:
440 break;
441 }
442 lsb->allocated = PETSC_TRUE;
443 }
444 PetscFunctionReturn(0);
445 }
446
447 /*------------------------------------------------------------*/
448
MatDestroy_LMVMSymBrdn(Mat B)449 static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
450 {
451 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
452 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
453 PetscErrorCode ierr;
454
455 PetscFunctionBegin;
456 if (lsb->allocated) {
457 ierr = VecDestroy(&lsb->work);CHKERRQ(ierr);
458 ierr = PetscFree5(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts);CHKERRQ(ierr);
459 ierr = PetscFree(lsb->psi);CHKERRQ(ierr);
460 ierr = VecDestroyVecs(lmvm->m, &lsb->P);CHKERRQ(ierr);
461 ierr = VecDestroyVecs(lmvm->m, &lsb->Q);CHKERRQ(ierr);
462 lsb->allocated = PETSC_FALSE;
463 }
464 ierr = MatDestroy(&lsb->D);CHKERRQ(ierr);
465 ierr = PetscFree(lmvm->ctx);CHKERRQ(ierr);
466 ierr = MatDestroy_LMVM(B);CHKERRQ(ierr);
467 PetscFunctionReturn(0);
468 }
469
470 /*------------------------------------------------------------*/
471
MatSetUp_LMVMSymBrdn(Mat B)472 static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
473 {
474 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
475 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
476 PetscErrorCode ierr;
477 PetscInt n, N;
478
479 PetscFunctionBegin;
480 ierr = MatSetUp_LMVM(B);CHKERRQ(ierr);
481 if (!lsb->allocated) {
482 ierr = VecDuplicate(lmvm->Xprev, &lsb->work);CHKERRQ(ierr);
483 if (lmvm->m > 0) {
484 ierr = PetscMalloc5(lmvm->m,&lsb->stp,lmvm->m,&lsb->ytq,lmvm->m,&lsb->yts,lmvm->m,&lsb->yty,lmvm->m,&lsb->sts);CHKERRQ(ierr);
485 ierr = PetscCalloc1(lmvm->m,&lsb->psi);CHKERRQ(ierr);
486 ierr = VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->P);CHKERRQ(ierr);
487 ierr = VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->Q);CHKERRQ(ierr);
488 }
489 switch (lsb->scale_type) {
490 case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
491 ierr = MatGetLocalSize(B, &n, &n);CHKERRQ(ierr);
492 ierr = MatGetSize(B, &N, &N);CHKERRQ(ierr);
493 ierr = MatSetSizes(lsb->D, n, n, N, N);CHKERRQ(ierr);
494 ierr = MatSetUp(lsb->D);CHKERRQ(ierr);
495 break;
496 default:
497 break;
498 }
499 lsb->allocated = PETSC_TRUE;
500 }
501 PetscFunctionReturn(0);
502 }
503
504 /*------------------------------------------------------------*/
505
MatView_LMVMSymBrdn(Mat B,PetscViewer pv)506 PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
507 {
508 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
509 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
510 PetscErrorCode ierr;
511 PetscBool isascii;
512
513 PetscFunctionBegin;
514 ierr = PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
515 if (isascii) {
516 ierr = PetscViewerASCIIPrintf(pv,"Scale type: %s\n",MatLMVMSymBroydenScaleTypes[lsb->scale_type]);CHKERRQ(ierr);
517 ierr = PetscViewerASCIIPrintf(pv,"Scale history: %d\n",lsb->sigma_hist);CHKERRQ(ierr);
518 ierr = PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)lsb->alpha, (double)lsb->beta, (double)lsb->rho);CHKERRQ(ierr);
519 ierr = PetscViewerASCIIPrintf(pv,"Convex factors: phi=%g, theta=%g\n",(double)lsb->phi, (double)lsb->theta);CHKERRQ(ierr);
520 }
521 ierr = MatView_LMVM(B, pv);CHKERRQ(ierr);
522 if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
523 ierr = MatView(lsb->D, pv);CHKERRQ(ierr);
524 }
525 PetscFunctionReturn(0);
526 }
527
528 /*------------------------------------------------------------*/
529
MatSetFromOptions_LMVMSymBrdn(PetscOptionItems * PetscOptionsObject,Mat B)530 PetscErrorCode MatSetFromOptions_LMVMSymBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
531 {
532 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
533 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
534 PetscErrorCode ierr;
535
536 PetscFunctionBegin;
537 ierr = MatSetFromOptions_LMVM(PetscOptionsObject, B);CHKERRQ(ierr);
538 ierr = PetscOptionsHead(PetscOptionsObject,"Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");CHKERRQ(ierr);
539 ierr = PetscOptionsReal("-mat_lmvm_phi","(developer) convex ratio between BFGS and DFP components of the update","",lsb->phi,&lsb->phi,NULL);CHKERRQ(ierr);
540 if ((lsb->phi < 0.0) || (lsb->phi > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the update formula cannot be outside the range of [0, 1]");
541 ierr = MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);CHKERRQ(ierr);
542 ierr = PetscOptionsTail();CHKERRQ(ierr);
543 PetscFunctionReturn(0);
544 }
545
MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionItems * PetscOptionsObject,Mat B)546 PetscErrorCode MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionItems *PetscOptionsObject, Mat B)
547 {
548 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
549 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
550 Mat_LMVM *dbase;
551 Mat_DiagBrdn *dctx;
552 MatLMVMSymBroydenScaleType stype = lsb->scale_type;
553 PetscBool flg;
554 PetscErrorCode ierr;
555
556 PetscFunctionBegin;
557 ierr = PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",lsb->beta,&lsb->beta,NULL);CHKERRQ(ierr);
558 ierr = PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",lsb->theta,&lsb->theta,NULL);CHKERRQ(ierr);
559 if ((lsb->theta < 0.0) || (lsb->theta > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
560 ierr = PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",lsb->rho,&lsb->rho,NULL);CHKERRQ(ierr);
561 if ((lsb->rho < 0.0) || (lsb->rho > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "update limiter in the J0 scaling cannot be outside the range of [0, 1]");
562 ierr = PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",lsb->alpha,&lsb->alpha,NULL);CHKERRQ(ierr);
563 if ((lsb->alpha < 0.0) || (lsb->alpha > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
564 ierr = PetscOptionsBoundedInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",lsb->sigma_hist,&lsb->sigma_hist,NULL,1);CHKERRQ(ierr);
565 ierr = PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0","MatLMVMSymBrdnScaleType",MatLMVMSymBroydenScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);CHKERRQ(ierr);
566 if (flg) ierr = MatLMVMSymBroydenSetScaleType(B, stype);CHKERRQ(ierr);
567 if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
568 ierr = MatSetFromOptions(lsb->D);CHKERRQ(ierr);
569 dbase = (Mat_LMVM*)lsb->D->data;
570 dctx = (Mat_DiagBrdn*)dbase->ctx;
571 dctx->delta_min = lsb->delta_min;
572 dctx->delta_max = lsb->delta_max;
573 dctx->theta = lsb->theta;
574 dctx->rho = lsb->rho;
575 dctx->alpha = lsb->alpha;
576 dctx->beta = lsb->beta;
577 dctx->sigma_hist = lsb->sigma_hist;
578 }
579 PetscFunctionReturn(0);
580 }
581
582 /*------------------------------------------------------------*/
583
MatCreate_LMVMSymBrdn(Mat B)584 PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
585 {
586 Mat_LMVM *lmvm;
587 Mat_SymBrdn *lsb;
588 PetscErrorCode ierr;
589
590 PetscFunctionBegin;
591 ierr = MatCreate_LMVM(B);CHKERRQ(ierr);
592 ierr = PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN);CHKERRQ(ierr);
593 ierr = MatSetOption(B, MAT_SPD, PETSC_TRUE);CHKERRQ(ierr);
594 B->ops->view = MatView_LMVMSymBrdn;
595 B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
596 B->ops->setup = MatSetUp_LMVMSymBrdn;
597 B->ops->destroy = MatDestroy_LMVMSymBrdn;
598 B->ops->solve = MatSolve_LMVMSymBrdn;
599
600 lmvm = (Mat_LMVM*)B->data;
601 lmvm->square = PETSC_TRUE;
602 lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
603 lmvm->ops->reset = MatReset_LMVMSymBrdn;
604 lmvm->ops->update = MatUpdate_LMVMSymBrdn;
605 lmvm->ops->mult = MatMult_LMVMSymBrdn;
606 lmvm->ops->copy = MatCopy_LMVMSymBrdn;
607
608 ierr = PetscNewLog(B, &lsb);CHKERRQ(ierr);
609 lmvm->ctx = (void*)lsb;
610 lsb->allocated = PETSC_FALSE;
611 lsb->needP = lsb->needQ = PETSC_TRUE;
612 lsb->phi = 0.125;
613 lsb->theta = 0.125;
614 lsb->alpha = 1.0;
615 lsb->rho = 1.0;
616 lsb->beta = 0.5;
617 lsb->sigma = 1.0;
618 lsb->delta = 1.0;
619 lsb->delta_min = 1e-7;
620 lsb->delta_max = 100.0;
621 lsb->sigma_hist = 1;
622 lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL;
623 lsb->watchdog = 0;
624 lsb->max_seq_rejects = lmvm->m/2;
625
626 ierr = MatCreate(PetscObjectComm((PetscObject)B), &lsb->D);CHKERRQ(ierr);
627 ierr = MatSetType(lsb->D, MATLMVMDIAGBROYDEN);CHKERRQ(ierr);
628 ierr = MatSetOptionsPrefix(lsb->D, "J0_");CHKERRQ(ierr);
629 PetscFunctionReturn(0);
630 }
631
632 /*------------------------------------------------------------*/
633
634 /*@
635 MatLMVMSymBroydenSetDelta - Sets the starting value for the diagonal scaling vector computed
636 in the SymBrdn approximations (also works for BFGS and DFP).
637
638 Input Parameters:
639 + B - LMVM matrix
640 - delta - initial value for diagonal scaling
641
642 Level: intermediate
643 @*/
644
MatLMVMSymBroydenSetDelta(Mat B,PetscScalar delta)645 PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
646 {
647 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
648 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
649 PetscErrorCode ierr;
650 PetscBool is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn;
651
652 PetscFunctionBegin;
653 ierr = PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs);CHKERRQ(ierr);
654 ierr = PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp);CHKERRQ(ierr);
655 ierr = PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn);CHKERRQ(ierr);
656 ierr = PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn);CHKERRQ(ierr);
657 if (!is_bfgs && !is_dfp && !is_symbrdn && !is_symbadbrdn) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "diagonal scaling is only available for DFP, BFGS and SymBrdn matrices");
658 lsb->delta = PetscAbsReal(PetscRealPart(delta));
659 lsb->delta = PetscMin(lsb->delta, lsb->delta_max);CHKERRQ(ierr);
660 lsb->delta = PetscMax(lsb->delta, lsb->delta_min);CHKERRQ(ierr);
661 PetscFunctionReturn(0);
662 }
663
664 /*------------------------------------------------------------*/
665
666 /*@
667 MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.
668
669 Input Parameters:
670 + snes - the iterative context
671 - rtype - restart type
672
673 Options Database:
674 . -mat_lmvm_scale_type <none,scalar,diagonal> - set the scaling type
675
676 Level: intermediate
677
678 MatLMVMSymBrdnScaleTypes:
679 + MAT_LMVM_SYMBROYDEN_SCALE_NONE - initial Hessian is the identity matrix
680 . MAT_LMVM_SYMBROYDEN_SCALE_SCALAR - use the Shanno scalar as the initial Hessian
681 - MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL - use a diagonalized BFGS update as the initial Hessian
682
683 .seealso: MATLMVMSYMBROYDEN, MatCreateLMVMSymBroyden()
684 @*/
MatLMVMSymBroydenSetScaleType(Mat B,MatLMVMSymBroydenScaleType stype)685 PetscErrorCode MatLMVMSymBroydenSetScaleType(Mat B, MatLMVMSymBroydenScaleType stype)
686 {
687 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
688 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
689
690 PetscFunctionBegin;
691 PetscValidHeaderSpecific(B,MAT_CLASSID,1);
692 lsb->scale_type = stype;
693 PetscFunctionReturn(0);
694 }
695
696 /*------------------------------------------------------------*/
697
698 /*@
699 MatCreateLMVMSymBroyden - Creates a limited-memory Symmetric Broyden-type matrix used
700 for approximating Jacobians. L-SymBrdn is a convex combination of L-DFP and
701 L-BFGS such that SymBrdn = (1 - phi)*BFGS + phi*DFP. The combination factor
702 phi is restricted to the range [0, 1], where the L-SymBrdn matrix is guaranteed
703 to be symmetric positive-definite.
704
705 The provided local and global sizes must match the solution and function vectors
706 used with MatLMVMUpdate() and MatSolve(). The resulting L-SymBrdn matrix will have
707 storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
708 parallel. To use the L-SymBrdn matrix with other vector types, the matrix must be
709 created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
710 This ensures that the internal storage and work vectors are duplicated from the
711 correct type of vector.
712
713 Collective
714
715 Input Parameters:
716 + comm - MPI communicator, set to PETSC_COMM_SELF
717 . n - number of local rows for storage vectors
718 - N - global size of the storage vectors
719
720 Output Parameter:
721 . B - the matrix
722
723 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
724 paradigm instead of this routine directly.
725
726 Options Database Keys:
727 + -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
728 . -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
729 . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
730 . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
731 . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
732 . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
733 . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
734 - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
735
736 Level: intermediate
737
738 .seealso: MatCreate(), MATLMVM, MATLMVMSYMBROYDEN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
739 MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn()
740 @*/
MatCreateLMVMSymBroyden(MPI_Comm comm,PetscInt n,PetscInt N,Mat * B)741 PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
742 {
743 PetscErrorCode ierr;
744
745 PetscFunctionBegin;
746 ierr = MatCreate(comm, B);CHKERRQ(ierr);
747 ierr = MatSetSizes(*B, n, n, N, N);CHKERRQ(ierr);
748 ierr = MatSetType(*B, MATLMVMSYMBROYDEN);CHKERRQ(ierr);
749 ierr = MatSetUp(*B);CHKERRQ(ierr);
750 PetscFunctionReturn(0);
751 }
752
753 /*------------------------------------------------------------*/
754
MatSymBrdnApplyJ0Fwd(Mat B,Vec X,Vec Z)755 PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
756 {
757 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
758 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
759 PetscErrorCode ierr;
760
761 PetscFunctionBegin;
762 if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
763 lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
764 ierr = MatLMVMApplyJ0Fwd(B, X, Z);CHKERRQ(ierr);
765 } else {
766 switch (lsb->scale_type) {
767 case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
768 ierr = VecCopy(X, Z);CHKERRQ(ierr);
769 ierr = VecScale(Z, 1.0/lsb->sigma);CHKERRQ(ierr);
770 break;
771 case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
772 ierr = MatMult(lsb->D, X, Z);CHKERRQ(ierr);
773 break;
774 case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
775 default:
776 ierr = VecCopy(X, Z);CHKERRQ(ierr);
777 break;
778 }
779 }
780 PetscFunctionReturn(0);
781 }
782
783 /*------------------------------------------------------------*/
784
MatSymBrdnApplyJ0Inv(Mat B,Vec F,Vec dX)785 PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
786 {
787 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
788 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
789 PetscErrorCode ierr;
790
791 PetscFunctionBegin;
792 if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
793 lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
794 ierr = MatLMVMApplyJ0Inv(B, F, dX);CHKERRQ(ierr);
795 } else {
796 switch (lsb->scale_type) {
797 case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
798 ierr = VecCopy(F, dX);CHKERRQ(ierr);
799 ierr = VecScale(dX, lsb->sigma);CHKERRQ(ierr);
800 break;
801 case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
802 ierr = MatSolve(lsb->D, F, dX);CHKERRQ(ierr);
803 break;
804 case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
805 default:
806 ierr = VecCopy(F, dX);CHKERRQ(ierr);
807 break;
808 }
809 }
810 PetscFunctionReturn(0);
811 }
812
813 /*------------------------------------------------------------*/
814
MatSymBrdnComputeJ0Scalar(Mat B)815 PetscErrorCode MatSymBrdnComputeJ0Scalar(Mat B)
816 {
817 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
818 Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
819 PetscInt i, start;
820 PetscReal a, b, c, sig1, sig2, signew;
821
822 PetscFunctionBegin;
823 if (lsb->sigma_hist == 0) {
824 signew = 1.0;
825 } else {
826 start = PetscMax(0, lmvm->k-lsb->sigma_hist+1);
827 signew = 0.0;
828 if (lsb->alpha == 1.0) {
829 for (i = start; i <= lmvm->k; ++i) {
830 signew += lsb->yts[i]/lsb->yty[i];
831 }
832 } else if (lsb->alpha == 0.5) {
833 for (i = start; i <= lmvm->k; ++i) {
834 signew += lsb->sts[i]/lsb->yty[i];
835 }
836 signew = PetscSqrtReal(signew);
837 } else if (lsb->alpha == 0.0) {
838 for (i = start; i <= lmvm->k; ++i) {
839 signew += lsb->sts[i]/lsb->yts[i];
840 }
841 } else {
842 /* compute coefficients of the quadratic */
843 a = b = c = 0.0;
844 for (i = start; i <= lmvm->k; ++i) {
845 a += lsb->yty[i];
846 b += lsb->yts[i];
847 c += lsb->sts[i];
848 }
849 a *= lsb->alpha;
850 b *= -(2.0*lsb->alpha - 1.0);
851 c *= lsb->alpha - 1.0;
852 /* use quadratic formula to find roots */
853 sig1 = (-b + PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
854 sig2 = (-b - PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
855 /* accept the positive root as the scalar */
856 if (sig1 > 0.0) {
857 signew = sig1;
858 } else if (sig2 > 0.0) {
859 signew = sig2;
860 } else {
861 SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
862 }
863 }
864 }
865 lsb->sigma = lsb->rho*signew + (1.0 - lsb->rho)*lsb->sigma;
866 PetscFunctionReturn(0);
867 }
868