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