1 #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h> /*I "petscksp.h" I*/
2
3 /*------------------------------------------------------------*/
4
MatSolve_DiagBrdn(Mat B,Vec F,Vec dX)5 static PetscErrorCode MatSolve_DiagBrdn(Mat B, Vec F, Vec dX)
6 {
7 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
8 Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
9 PetscErrorCode ierr;
10
11 PetscFunctionBegin;
12 VecCheckSameSize(F, 2, dX, 3);
13 VecCheckMatCompatible(B, dX, 3, F, 2);
14 ierr = VecPointwiseMult(dX, ldb->invD, F);CHKERRQ(ierr);
15 PetscFunctionReturn(0);
16 }
17
18 /*------------------------------------------------------------*/
19
MatMult_DiagBrdn(Mat B,Vec X,Vec Z)20 static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
21 {
22 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
23 Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
24 PetscErrorCode ierr;
25
26 PetscFunctionBegin;
27 VecCheckSameSize(X, 2, Z, 3);
28 VecCheckMatCompatible(B, X, 2, Z, 3);
29 ierr = VecPointwiseDivide(Z, X, ldb->invD);CHKERRQ(ierr);
30 PetscFunctionReturn(0);
31 }
32
33 /*------------------------------------------------------------*/
34
MatUpdate_DiagBrdn(Mat B,Vec X,Vec F)35 static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
36 {
37 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
38 Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
39 PetscErrorCode ierr;
40 PetscInt old_k, i, start;
41 PetscScalar yty, ststmp, curvature, ytDy, stDs, ytDs;
42 PetscReal curvtol, sigma, yy_sum, ss_sum, ys_sum, denom;
43
44 PetscFunctionBegin;
45 if (!lmvm->m) PetscFunctionReturn(0);
46 if (lmvm->prev_set) {
47 /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
48 ierr = VecAYPX(lmvm->Xprev, -1.0, X);CHKERRQ(ierr);
49 ierr = VecAYPX(lmvm->Fprev, -1.0, F);CHKERRQ(ierr);
50 /* Compute tolerance for accepting the update */
51 ierr = VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);CHKERRQ(ierr);
52 ierr = VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);CHKERRQ(ierr);
53 ierr = VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);CHKERRQ(ierr);
54 ierr = VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);CHKERRQ(ierr);
55 if (PetscRealPart(ststmp) < lmvm->eps){
56 curvtol = 0.0;
57 } else {
58 curvtol = lmvm->eps * PetscRealPart(ststmp);
59 }
60 /* Test the curvature for the update */
61 if (PetscRealPart(curvature) > curvtol) {
62 /* Update is good so we accept it */
63 old_k = lmvm->k;
64 ierr = MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);CHKERRQ(ierr);
65 /* If we hit the memory limit, shift the yty and yts arrays */
66 if (old_k == lmvm->k) {
67 for (i = 0; i <= lmvm->k-1; ++i) {
68 ldb->yty[i] = ldb->yty[i+1];
69 ldb->yts[i] = ldb->yts[i+1];
70 ldb->sts[i] = ldb->sts[i+1];
71 }
72 }
73 /* Accept dot products into the history */
74 ierr = VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);CHKERRQ(ierr);
75 ldb->yty[lmvm->k] = PetscRealPart(yty);
76 ldb->yts[lmvm->k] = PetscRealPart(curvature);
77 ldb->sts[lmvm->k] = PetscRealPart(ststmp);
78 if (ldb->forward) {
79 /* We are doing diagonal scaling of the forward Hessian B */
80 /* BFGS = DFP = inv(D); */
81 ierr = VecCopy(ldb->invD, ldb->invDnew);CHKERRQ(ierr);
82 ierr = VecReciprocal(ldb->invDnew);CHKERRQ(ierr);
83
84 /* V = y*y */
85 ierr = VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]);CHKERRQ(ierr);
86
87 /* W = inv(D)*s */
88 ierr = VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]);CHKERRQ(ierr);
89 ierr = VecDot(ldb->W, lmvm->S[lmvm->k], &stDs);CHKERRQ(ierr);
90
91 /* Safeguard stDs */
92 stDs = PetscMax(PetscRealPart(stDs), ldb->tol);
93
94 if (1.0 != ldb->theta) {
95 /* BFGS portion of the update */
96 /* U = (inv(D)*s)*(inv(D)*s) */
97 ierr = VecPointwiseMult(ldb->U, ldb->W, ldb->W);CHKERRQ(ierr);
98
99 /* Assemble */
100 ierr = VecAXPBY(ldb->BFGS, -1.0/stDs, 0.0, ldb->U);CHKERRQ(ierr);
101 }
102 if (0.0 != ldb->theta) {
103 /* DFP portion of the update */
104 /* U = inv(D)*s*y */
105 ierr = VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]);CHKERRQ(ierr);
106
107 /* Assemble */
108 ierr = VecAXPBY(ldb->DFP, stDs/ldb->yts[lmvm->k], 0.0, ldb->V);CHKERRQ(ierr);
109 ierr = VecAXPY(ldb->DFP, -2.0, ldb->U);CHKERRQ(ierr);
110 }
111
112 if (0.0 == ldb->theta) {
113 ierr = VecAXPY(ldb->invDnew, 1.0, ldb->BFGS);CHKERRQ(ierr);
114 } else if (1.0 == ldb->theta) {
115 ierr = VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->DFP);CHKERRQ(ierr);
116 } else {
117 /* Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
118 ierr = VecAXPBYPCZ(ldb->invDnew, 1.0-ldb->theta, (ldb->theta)/ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP);CHKERRQ(ierr);
119 }
120
121 ierr = VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);CHKERRQ(ierr);
122 /* Obtain inverse and ensure positive definite */
123 ierr = VecReciprocal(ldb->invDnew);CHKERRQ(ierr);
124 ierr = VecAbs(ldb->invDnew);CHKERRQ(ierr);
125
126 } else {
127 /* Inverse Hessian update instead. */
128 ierr = VecCopy(ldb->invD, ldb->invDnew);CHKERRQ(ierr);
129
130 /* V = s*s */
131 ierr = VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]);CHKERRQ(ierr);
132
133 /* W = D*y */
134 ierr = VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]);CHKERRQ(ierr);
135 ierr = VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy);CHKERRQ(ierr);
136
137 /* Safeguard ytDy */
138 ytDy = PetscMax(PetscRealPart(ytDy), ldb->tol);
139
140 if (1.0 != ldb->theta) {
141 /* BFGS portion of the update */
142 /* U = s*Dy */
143 ierr = VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]);CHKERRQ(ierr);
144
145 /* Assemble */
146 ierr = VecAXPBY(ldb->BFGS, ytDy/ldb->yts[lmvm->k], 0.0, ldb->V);CHKERRQ(ierr);
147 ierr = VecAXPY(ldb->BFGS, -2.0, ldb->U);CHKERRQ(ierr);
148 }
149 if (0.0 != ldb->theta) {
150 /* DFP portion of the update */
151
152 /* U = (inv(D)*y)*(inv(D)*y) */
153 ierr = VecPointwiseMult(ldb->U, ldb->W, ldb->W);CHKERRQ(ierr);
154
155 /* Assemble */
156 ierr = VecAXPBY(ldb->DFP, -1.0/ytDy, 0.0, ldb->U);CHKERRQ(ierr);
157 }
158
159 if (0.0 == ldb->theta) {
160 ierr = VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->BFGS);CHKERRQ(ierr);
161 } else if (1.0 == ldb->theta) {
162 ierr = VecAXPY(ldb->invDnew, 1.0, ldb->DFP);CHKERRQ(ierr);
163 } else {
164 /* Broyden update U=(1-theta)*P + theta*Q */
165 ierr = VecAXPBYPCZ(ldb->invDnew, (1.0-ldb->theta)/ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP);CHKERRQ(ierr);
166 }
167 ierr = VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);CHKERRQ(ierr);
168 /* Ensure positive definite */
169 ierr = VecAbs(ldb->invDnew);CHKERRQ(ierr);
170 }
171 if (ldb->sigma_hist > 0){
172 /* Start with re-scaling on the newly computed diagonal */
173 if (0.5 == ldb->beta) {
174 if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
175 ierr = VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew);CHKERRQ(ierr);
176 ierr = VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);CHKERRQ(ierr);
177
178 ierr = VecDotBegin(ldb->V, lmvm->Y[0], &ytDy);CHKERRQ(ierr);
179 ierr = VecDotBegin(ldb->W, lmvm->S[0], &stDs);CHKERRQ(ierr);
180 ierr = VecDotEnd(ldb->V, lmvm->Y[0], &ytDy);CHKERRQ(ierr);
181 ierr = VecDotEnd(ldb->W, lmvm->S[0], &stDs);CHKERRQ(ierr);
182
183 ss_sum = PetscRealPart(stDs);
184 yy_sum = PetscRealPart(ytDy);
185 ys_sum = ldb->yts[0];
186 } else {
187 ierr = VecCopy(ldb->invDnew, ldb->U);CHKERRQ(ierr);
188 ierr = VecReciprocal(ldb->U);CHKERRQ(ierr);
189
190 /* Compute summations for scalar scaling */
191 yy_sum = 0; /* No safeguard required */
192 ys_sum = 0; /* No safeguard required */
193 ss_sum = 0; /* No safeguard required */
194 start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
195 for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
196 ierr = VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U);CHKERRQ(ierr);
197 ierr = VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);CHKERRQ(ierr);
198
199 ierr = VecDotBegin(ldb->W, lmvm->S[i], &stDs);CHKERRQ(ierr);
200 ierr = VecDotBegin(ldb->V, lmvm->Y[i], &ytDy);CHKERRQ(ierr);
201 ierr = VecDotEnd(ldb->W, lmvm->S[i], &stDs);CHKERRQ(ierr);
202 ierr = VecDotEnd(ldb->V, lmvm->Y[i], &ytDy);CHKERRQ(ierr);
203
204 ss_sum += PetscRealPart(stDs);
205 ys_sum += ldb->yts[i];
206 yy_sum += PetscRealPart(ytDy);
207 }
208 }
209 } else if (0.0 == ldb->beta) {
210 if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
211 /* Compute summations for scalar scaling */
212 ierr = VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);CHKERRQ(ierr);
213
214 ierr = VecDotBegin(ldb->W, lmvm->Y[0], &ytDs);CHKERRQ(ierr);
215 ierr = VecDotBegin(ldb->W, ldb->W, &stDs);CHKERRQ(ierr);
216 ierr = VecDotEnd(ldb->W, lmvm->Y[0], &ytDs);CHKERRQ(ierr);
217 ierr = VecDotEnd(ldb->W, ldb->W, &stDs);CHKERRQ(ierr);
218
219 ys_sum = PetscRealPart(ytDs);
220 ss_sum = PetscRealPart(stDs);
221 yy_sum = ldb->yty[0];
222 } else {
223 ierr = VecCopy(ldb->invDnew, ldb->U);CHKERRQ(ierr);
224 ierr = VecReciprocal(ldb->U);CHKERRQ(ierr);
225
226 /* Compute summations for scalar scaling */
227 yy_sum = 0; /* No safeguard required */
228 ys_sum = 0; /* No safeguard required */
229 ss_sum = 0; /* No safeguard required */
230 start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
231 for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
232 ierr = VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);CHKERRQ(ierr);
233
234 ierr = VecDotBegin(ldb->W, lmvm->Y[i], &ytDs);CHKERRQ(ierr);
235 ierr = VecDotBegin(ldb->W, ldb->W, &stDs);CHKERRQ(ierr);
236 ierr = VecDotEnd(ldb->W, lmvm->Y[i], &ytDs);CHKERRQ(ierr);
237 ierr = VecDotEnd(ldb->W, ldb->W, &stDs);CHKERRQ(ierr);
238
239 ss_sum += PetscRealPart(stDs);
240 ys_sum += PetscRealPart(ytDs);
241 yy_sum += ldb->yty[i];
242 }
243 }
244 } else if (1.0 == ldb->beta) {
245 /* Compute summations for scalar scaling */
246 yy_sum = 0; /* No safeguard required */
247 ys_sum = 0; /* No safeguard required */
248 ss_sum = 0; /* No safeguard required */
249 start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
250 for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
251 ierr = VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew);CHKERRQ(ierr);
252
253 ierr = VecDotBegin(ldb->V, lmvm->S[i], &ytDs);CHKERRQ(ierr);
254 ierr = VecDotBegin(ldb->V, ldb->V, &ytDy);CHKERRQ(ierr);
255 ierr = VecDotEnd(ldb->V, lmvm->S[i], &ytDs);CHKERRQ(ierr);
256 ierr = VecDotEnd(ldb->V, ldb->V, &ytDy);CHKERRQ(ierr);
257
258 yy_sum += PetscRealPart(ytDy);
259 ys_sum += PetscRealPart(ytDs);
260 ss_sum += ldb->sts[i];
261 }
262 } else {
263 ierr = VecCopy(ldb->invDnew, ldb->U);CHKERRQ(ierr);
264 ierr = VecPow(ldb->U, ldb->beta-1);CHKERRQ(ierr);
265
266 /* Compute summations for scalar scaling */
267 yy_sum = 0; /* No safeguard required */
268 ys_sum = 0; /* No safeguard required */
269 ss_sum = 0; /* No safeguard required */
270 start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
271 for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
272 ierr = VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]);CHKERRQ(ierr);
273 ierr = VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]);CHKERRQ(ierr);
274
275 ierr = VecDotBegin(ldb->V, ldb->W, &ytDs);CHKERRQ(ierr);
276 ierr = VecDotBegin(ldb->V, ldb->V, &ytDy);CHKERRQ(ierr);
277 ierr = VecDotBegin(ldb->W, ldb->W, &stDs);CHKERRQ(ierr);
278 ierr = VecDotEnd(ldb->V, ldb->W, &ytDs);CHKERRQ(ierr);
279 ierr = VecDotEnd(ldb->V, ldb->V, &ytDy);CHKERRQ(ierr);
280 ierr = VecDotEnd(ldb->W, ldb->W, &stDs);CHKERRQ(ierr);
281
282 yy_sum += PetscRealPart(ytDy);
283 ys_sum += PetscRealPart(ytDs);
284 ss_sum += PetscRealPart(stDs);
285 }
286 }
287
288 if (0.0 == ldb->alpha) {
289 /* Safeguard ys_sum */
290 ys_sum = PetscMax(ldb->tol, ys_sum);
291
292 sigma = ss_sum / ys_sum;
293 } else if (1.0 == ldb->alpha) {
294 /* yy_sum is never 0; if it were, we'd be at the minimum */
295 sigma = ys_sum / yy_sum;
296 } else {
297 denom = 2.0*ldb->alpha*yy_sum;
298
299 /* Safeguard denom */
300 denom = PetscMax(ldb->tol, denom);
301
302 sigma = ((2.0*ldb->alpha-1)*ys_sum + PetscSqrtReal((2.0*ldb->alpha-1)*(2.0*ldb->alpha-1)*ys_sum*ys_sum - 4.0*ldb->alpha*(ldb->alpha-1)*yy_sum*ss_sum)) / denom;
303 }
304 } else {
305 sigma = 1.0;
306 }
307 /* If Q has small values, then Q^(r_beta - 1)
308 can have very large values. Hence, ys_sum
309 and ss_sum can be infinity. In this case,
310 sigma can either be not-a-number or infinity. */
311
312 if (PetscIsInfOrNanScalar(sigma)) {
313 /* sigma is not-a-number; skip rescaling */
314 } else if (0.0 == sigma) {
315 /* sigma is zero; this is a bad case; skip rescaling */
316 } else {
317 /* sigma is positive */
318 ierr = VecScale(ldb->invDnew, sigma);CHKERRQ(ierr);
319 }
320
321 /* Combine the old diagonal and the new diagonal using a convex limiter */
322 if (1.0 == ldb->rho) {
323 ierr = VecCopy(ldb->invDnew, ldb->invD);CHKERRQ(ierr);
324 } else if (ldb->rho) {
325 ierr = VecAXPBY(ldb->invD, 1.0-ldb->rho, ldb->rho, ldb->invDnew);CHKERRQ(ierr);
326 }
327 } else {
328 ierr = MatLMVMReset(B, PETSC_FALSE);CHKERRQ(ierr);
329 }
330 /* End DiagBrdn update */
331
332 }
333 /* Save the solution and function to be used in the next update */
334 ierr = VecCopy(X, lmvm->Xprev);CHKERRQ(ierr);
335 ierr = VecCopy(F, lmvm->Fprev);CHKERRQ(ierr);
336 lmvm->prev_set = PETSC_TRUE;
337 PetscFunctionReturn(0);
338 }
339
340 /*------------------------------------------------------------*/
341
MatCopy_DiagBrdn(Mat B,Mat M,MatStructure str)342 static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
343 {
344 Mat_LMVM *bdata = (Mat_LMVM*)B->data;
345 Mat_DiagBrdn *bctx = (Mat_DiagBrdn*)bdata->ctx;
346 Mat_LMVM *mdata = (Mat_LMVM*)M->data;
347 Mat_DiagBrdn *mctx = (Mat_DiagBrdn*)mdata->ctx;
348 PetscErrorCode ierr;
349 PetscInt i;
350
351 PetscFunctionBegin;
352 mctx->theta = bctx->theta;
353 mctx->alpha = bctx->alpha;
354 mctx->beta = bctx->beta;
355 mctx->rho = bctx->rho;
356 mctx->delta = bctx->delta;
357 mctx->delta_min = bctx->delta_min;
358 mctx->delta_max = bctx->delta_max;
359 mctx->tol = bctx->tol;
360 mctx->sigma = bctx->sigma;
361 mctx->sigma_hist = bctx->sigma_hist;
362 mctx->forward = bctx->forward;
363 ierr = VecCopy(bctx->invD, mctx->invD);CHKERRQ(ierr);
364 for (i=0; i<=bdata->k; ++i) {
365 mctx->yty[i] = bctx->yty[i];
366 mctx->yts[i] = bctx->yts[i];
367 mctx->sts[i] = bctx->sts[i];
368 }
369 PetscFunctionReturn(0);
370 }
371
372 /*------------------------------------------------------------*/
373
MatView_DiagBrdn(Mat B,PetscViewer pv)374 static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
375 {
376 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
377 Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
378 PetscErrorCode ierr;
379 PetscBool isascii;
380
381 PetscFunctionBegin;
382 ierr = PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
383 if (isascii) {
384 ierr = PetscViewerASCIIPrintf(pv,"Scale history: %d\n",ldb->sigma_hist);CHKERRQ(ierr);
385 ierr = PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)ldb->alpha, (double)ldb->beta, (double)ldb->rho);CHKERRQ(ierr);
386 ierr = PetscViewerASCIIPrintf(pv,"Convex factor: theta=%g\n", (double)ldb->theta);CHKERRQ(ierr);
387 }
388 ierr = MatView_LMVM(B, pv);CHKERRQ(ierr);
389 PetscFunctionReturn(0);
390 }
391
392 /*------------------------------------------------------------*/
393
MatSetFromOptions_DiagBrdn(PetscOptionItems * PetscOptionsObject,Mat B)394 static PetscErrorCode MatSetFromOptions_DiagBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
395 {
396 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
397 Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
398 PetscErrorCode ierr;
399
400 PetscFunctionBegin;
401 ierr = MatSetFromOptions_LMVM(PetscOptionsObject, B);CHKERRQ(ierr);
402 ierr = PetscOptionsHead(PetscOptionsObject,"Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");CHKERRQ(ierr);
403 ierr = PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",ldb->theta,&ldb->theta,NULL);CHKERRQ(ierr);
404 ierr = PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",ldb->rho,&ldb->rho,NULL);CHKERRQ(ierr);
405 ierr = PetscOptionsReal("-mat_lmvm_tol","(developer) tolerance for bounding rescaling denominator","",ldb->tol,&ldb->tol,NULL);CHKERRQ(ierr);
406 ierr = PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",ldb->alpha,&ldb->alpha,NULL);CHKERRQ(ierr);
407 ierr = PetscOptionsBool("-mat_lmvm_forward","Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.","",ldb->forward,&ldb->forward,NULL);CHKERRQ(ierr);
408 ierr = PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",ldb->beta,&ldb->beta,NULL);CHKERRQ(ierr);
409 ierr = PetscOptionsInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",ldb->sigma_hist,&ldb->sigma_hist,NULL);CHKERRQ(ierr);
410 ierr = PetscOptionsTail();CHKERRQ(ierr);
411 if ((ldb->theta < 0.0) || (ldb->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]");
412 if ((ldb->alpha < 0.0) || (ldb->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]");
413 if ((ldb->rho < 0.0) || (ldb->rho > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex update limiter in the J0 scaling cannot be outside the range of [0, 1]");
414 if (ldb->sigma_hist < 0) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
415 PetscFunctionReturn(0);
416 }
417
418 /*------------------------------------------------------------*/
419
MatReset_DiagBrdn(Mat B,PetscBool destructive)420 static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
421 {
422 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
423 Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
424 PetscErrorCode ierr;
425
426 PetscFunctionBegin;
427 ierr = VecSet(ldb->invD, ldb->delta);CHKERRQ(ierr);
428 if (destructive && ldb->allocated) {
429 ierr = PetscFree3(ldb->yty, ldb->yts, ldb->sts);CHKERRQ(ierr);
430 ierr = VecDestroy(&ldb->invDnew);CHKERRQ(ierr);
431 ierr = VecDestroy(&ldb->invD);CHKERRQ(ierr);
432 ierr = VecDestroy(&ldb->BFGS);CHKERRQ(ierr);
433 ierr = VecDestroy(&ldb->DFP);CHKERRQ(ierr);
434 ierr = VecDestroy(&ldb->U);CHKERRQ(ierr);
435 ierr = VecDestroy(&ldb->V);CHKERRQ(ierr);
436 ierr = VecDestroy(&ldb->W);CHKERRQ(ierr);
437 ldb->allocated = PETSC_FALSE;
438 }
439 ierr = MatReset_LMVM(B, destructive);CHKERRQ(ierr);
440 PetscFunctionReturn(0);
441 }
442
443 /*------------------------------------------------------------*/
444
MatAllocate_DiagBrdn(Mat B,Vec X,Vec F)445 static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
446 {
447 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
448 Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
449 PetscErrorCode ierr;
450
451 PetscFunctionBegin;
452 ierr = MatAllocate_LMVM(B, X, F);CHKERRQ(ierr);
453 if (!ldb->allocated) {
454 ierr = PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);CHKERRQ(ierr);
455 ierr = VecDuplicate(lmvm->Xprev, &ldb->invDnew);CHKERRQ(ierr);
456 ierr = VecDuplicate(lmvm->Xprev, &ldb->invD);CHKERRQ(ierr);
457 ierr = VecDuplicate(lmvm->Xprev, &ldb->BFGS);CHKERRQ(ierr);
458 ierr = VecDuplicate(lmvm->Xprev, &ldb->DFP);CHKERRQ(ierr);
459 ierr = VecDuplicate(lmvm->Xprev, &ldb->U);CHKERRQ(ierr);
460 ierr = VecDuplicate(lmvm->Xprev, &ldb->V);CHKERRQ(ierr);
461 ierr = VecDuplicate(lmvm->Xprev, &ldb->W);CHKERRQ(ierr);
462 ldb->allocated = PETSC_TRUE;
463 }
464 PetscFunctionReturn(0);
465 }
466
467 /*------------------------------------------------------------*/
468
MatDestroy_DiagBrdn(Mat B)469 static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
470 {
471 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
472 Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
473 PetscErrorCode ierr;
474
475 PetscFunctionBegin;
476 if (ldb->allocated) {
477 ierr = PetscFree3(ldb->yty, ldb->yts, ldb->sts);CHKERRQ(ierr);
478 ierr = VecDestroy(&ldb->invDnew);CHKERRQ(ierr);
479 ierr = VecDestroy(&ldb->invD);CHKERRQ(ierr);
480 ierr = VecDestroy(&ldb->BFGS);CHKERRQ(ierr);
481 ierr = VecDestroy(&ldb->DFP);CHKERRQ(ierr);
482 ierr = VecDestroy(&ldb->U);CHKERRQ(ierr);
483 ierr = VecDestroy(&ldb->V);CHKERRQ(ierr);
484 ierr = VecDestroy(&ldb->W);CHKERRQ(ierr);
485 ldb->allocated = PETSC_FALSE;
486 }
487 ierr = PetscFree(lmvm->ctx);CHKERRQ(ierr);
488 ierr = MatDestroy_LMVM(B);CHKERRQ(ierr);
489 PetscFunctionReturn(0);
490 }
491
492 /*------------------------------------------------------------*/
493
MatSetUp_DiagBrdn(Mat B)494 static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
495 {
496 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
497 Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
498 PetscErrorCode ierr;
499
500 PetscFunctionBegin;
501 ierr = MatSetUp_LMVM(B);CHKERRQ(ierr);
502 if (!ldb->allocated) {
503 ierr = PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);CHKERRQ(ierr);
504 ierr = VecDuplicate(lmvm->Xprev, &ldb->invDnew);CHKERRQ(ierr);
505 ierr = VecDuplicate(lmvm->Xprev, &ldb->invD);CHKERRQ(ierr);
506 ierr = VecDuplicate(lmvm->Xprev, &ldb->BFGS);CHKERRQ(ierr);
507 ierr = VecDuplicate(lmvm->Xprev, &ldb->DFP);CHKERRQ(ierr);
508 ierr = VecDuplicate(lmvm->Xprev, &ldb->U);CHKERRQ(ierr);
509 ierr = VecDuplicate(lmvm->Xprev, &ldb->V);CHKERRQ(ierr);
510 ierr = VecDuplicate(lmvm->Xprev, &ldb->W);CHKERRQ(ierr);
511 ldb->allocated = PETSC_TRUE;
512 }
513 PetscFunctionReturn(0);
514 }
515
516 /*------------------------------------------------------------*/
517
MatCreate_LMVMDiagBrdn(Mat B)518 PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
519 {
520 Mat_LMVM *lmvm;
521 Mat_DiagBrdn *ldb;
522 PetscErrorCode ierr;
523
524 PetscFunctionBegin;
525 ierr = MatCreate_LMVM(B);CHKERRQ(ierr);
526 ierr = PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN);CHKERRQ(ierr);
527 B->ops->setup = MatSetUp_DiagBrdn;
528 B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
529 B->ops->destroy = MatDestroy_DiagBrdn;
530 B->ops->solve = MatSolve_DiagBrdn;
531 B->ops->view = MatView_DiagBrdn;
532
533 lmvm = (Mat_LMVM*)B->data;
534 lmvm->square = PETSC_TRUE;
535 lmvm->m = 1;
536 lmvm->ops->allocate = MatAllocate_DiagBrdn;
537 lmvm->ops->reset = MatReset_DiagBrdn;
538 lmvm->ops->mult = MatMult_DiagBrdn;
539 lmvm->ops->update = MatUpdate_DiagBrdn;
540 lmvm->ops->copy = MatCopy_DiagBrdn;
541
542 ierr = PetscNewLog(B, &ldb);CHKERRQ(ierr);
543 lmvm->ctx = (void*)ldb;
544 ldb->theta = 0.0;
545 ldb->alpha = 1.0;
546 ldb->rho = 1.0;
547 ldb->forward = PETSC_TRUE;
548 ldb->beta = 0.5;
549 ldb->sigma = 1.0;
550 ldb->delta = 1.0;
551 ldb->delta_min = 1e-7;
552 ldb->delta_max = 100.0;
553 ldb->tol = 1e-8;
554 ldb->sigma_hist = 1;
555 ldb->allocated = PETSC_FALSE;
556 PetscFunctionReturn(0);
557 }
558
559 /*------------------------------------------------------------*/
560
561 /*@
562 MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
563 for approximating Hessians. It consists of a convex combination of DFP and BFGS
564 diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
565 To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
566 We also ensure positive definiteness by taking the VecAbs() of the final vector.
567
568 There are two ways of approximating the diagonal: using the forward (B) update
569 schemes for BFGS and DFP and then taking the inverse, or directly working with
570 the inverse (H) update schemes for the BFGS and DFP updates, derived using the
571 Sherman-Morrison-Woodbury formula. We have implemented both, controlled by a
572 parameter below.
573
574 In order to use the DiagBrdn matrix with other vector types, i.e. doing MatMults
575 and MatSolves, the matrix must first be created using MatCreate() and MatSetType(),
576 followed by MatLMVMAllocate(). Then it will be available for updating
577 (via MatLMVMUpdate) in one's favored solver implementation.
578 This allows for MPI compatibility.
579
580 Collective
581
582 Input Parameters:
583 + comm - MPI communicator, set to PETSC_COMM_SELF
584 . n - number of local rows for storage vectors
585 - N - global size of the storage vectors
586
587 Output Parameter:
588 . B - the matrix
589
590 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
591 paradigm instead of this routine directly.
592
593 Options Database Keys:
594 + -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
595 . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
596 . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
597 . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
598 . -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
599 . -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
600 - -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal
601
602 Level: intermediate
603
604 .seealso: MatCreate(), MATLMVM, MATLMVMDIAGBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
605 MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMSymBrdn()
606 @*/
MatCreateLMVMDiagBroyden(MPI_Comm comm,PetscInt n,PetscInt N,Mat * B)607 PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
608 {
609 PetscErrorCode ierr;
610
611 PetscFunctionBegin;
612 ierr = MatCreate(comm, B);CHKERRQ(ierr);
613 ierr = MatSetSizes(*B, n, n, N, N);CHKERRQ(ierr);
614 ierr = MatSetType(*B, MATLMVMDIAGBROYDEN);CHKERRQ(ierr);
615 ierr = MatSetUp(*B);CHKERRQ(ierr);
616 PetscFunctionReturn(0);
617 }
618