1 #include <../src/ksp/ksp/utils/lmvm/lmvm.h> /*I "petscksp.h" I*/
2
3 /*@
4 MatLMVMUpdate - Adds (X-Xprev) and (F-Fprev) updates to an LMVM matrix.
5 The first time the function is called for an LMVM matrix, no update is
6 applied, but the given X and F vectors are stored for use as Xprev and
7 Fprev in the next update.
8
9 If the user has provided another LMVM matrix in place of J0, the J0
10 matrix is also updated recursively.
11
12 Input Parameters:
13 + B - An LMVM-type matrix
14 . X - Solution vector
15 - F - Function vector
16
17 Level: intermediate
18
19 .seealso: MatLMVMReset(), MatLMVMAllocate()
20 @*/
MatLMVMUpdate(Mat B,Vec X,Vec F)21 PetscErrorCode MatLMVMUpdate(Mat B, Vec X, Vec F)
22 {
23 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
24 PetscErrorCode ierr;
25 PetscBool same;
26
27 PetscFunctionBegin;
28 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
29 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
30 PetscValidHeaderSpecific(F, VEC_CLASSID, 3);
31 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
32 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
33 if (!lmvm->allocated) {
34 ierr = MatLMVMAllocate(B, X, F);CHKERRQ(ierr);
35 } else {
36 VecCheckMatCompatible(B, X, 2, F, 3);
37 }
38 if (lmvm->J0) {
39 /* If the user provided an LMVM-type matrix as J0, then trigger its update as well */
40 ierr = PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);CHKERRQ(ierr);
41 if (same) {
42 ierr = MatLMVMUpdate(lmvm->J0, X, F);CHKERRQ(ierr);
43 }
44 }
45 ierr = (*lmvm->ops->update)(B, X, F);CHKERRQ(ierr);
46 PetscFunctionReturn(0);
47 }
48
49 /*------------------------------------------------------------*/
50
51 /*@
52 MatLMVMClearJ0 - Removes all definitions of J0 and reverts to
53 an identity matrix (scale = 1.0).
54
55 Input Parameters:
56 . B - An LMVM-type matrix
57
58 Level: advanced
59
60 .seealso: MatLMVMSetJ0()
61 @*/
MatLMVMClearJ0(Mat B)62 PetscErrorCode MatLMVMClearJ0(Mat B)
63 {
64 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
65 PetscErrorCode ierr;
66 PetscBool same;
67 MPI_Comm comm = PetscObjectComm((PetscObject)B);
68
69 PetscFunctionBegin;
70 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
71 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
72 if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
73 lmvm->user_pc = PETSC_FALSE;
74 lmvm->user_ksp = PETSC_FALSE;
75 lmvm->user_scale = PETSC_FALSE;
76 lmvm->J0scalar = 1.0;
77 ierr = VecDestroy(&lmvm->J0diag);CHKERRQ(ierr);
78 ierr = MatDestroy(&lmvm->J0);CHKERRQ(ierr);
79 ierr = PCDestroy(&lmvm->J0pc);CHKERRQ(ierr);
80 PetscFunctionReturn(0);
81 }
82
83 /*------------------------------------------------------------*/
84
85 /*@
86 MatLMVMSetJ0Scale - Allows the user to define a scalar value
87 mu such that J0 = mu*I.
88
89 Input Parameters:
90 + B - An LMVM-type matrix
91 - scale - Scalar value mu that defines the initial Jacobian
92
93 Level: advanced
94
95 .seealso: MatLMVMSetDiagScale(), MatLMVMSetJ0()
96 @*/
MatLMVMSetJ0Scale(Mat B,PetscReal scale)97 PetscErrorCode MatLMVMSetJ0Scale(Mat B, PetscReal scale)
98 {
99 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
100 PetscErrorCode ierr;
101 PetscBool same;
102 MPI_Comm comm = PetscObjectComm((PetscObject)B);
103
104 PetscFunctionBegin;
105 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
106 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
107 if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
108 if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Scaling is available only for square LMVM matrices");
109 ierr = MatLMVMClearJ0(B);CHKERRQ(ierr);
110 lmvm->J0scalar = scale;
111 lmvm->user_scale = PETSC_TRUE;
112 PetscFunctionReturn(0);
113 }
114
115 /*------------------------------------------------------------*/
116
117 /*@
118 MatLMVMSetJ0Diag - Allows the user to define a vector
119 V such that J0 = diag(V).
120
121 Input Parameters:
122 + B - An LMVM-type matrix
123 - V - Vector that defines the diagonal of the initial Jacobian
124
125 Level: advanced
126
127 .seealso: MatLMVMSetScale(), MatLMVMSetJ0()
128 @*/
MatLMVMSetJ0Diag(Mat B,Vec V)129 PetscErrorCode MatLMVMSetJ0Diag(Mat B, Vec V)
130 {
131 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
132 PetscErrorCode ierr;
133 PetscBool same;
134 MPI_Comm comm = PetscObjectComm((PetscObject)B);
135
136 PetscFunctionBegin;
137 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
138 PetscValidHeaderSpecific(V, VEC_CLASSID, 2);
139 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
140 if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
141 if (!lmvm->allocated) SETERRQ(comm, PETSC_ERR_ORDER, "Matrix must be allocated before setting diagonal scaling");
142 if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Diagonal scaling is available only for square LMVM matrices");
143 VecCheckSameSize(V, 2, lmvm->Fprev, 3);CHKERRQ(ierr);
144 ierr = MatLMVMClearJ0(B);CHKERRQ(ierr);
145 if (!lmvm->J0diag) {
146 ierr = VecDuplicate(V, &lmvm->J0diag);CHKERRQ(ierr);
147 }
148 ierr = VecCopy(V, lmvm->J0diag);CHKERRQ(ierr);
149 lmvm->user_scale = PETSC_TRUE;
150 PetscFunctionReturn(0);
151 }
152
153 /*------------------------------------------------------------*/
154
155 /*@
156 MatLMVMSetJ0 - Allows the user to define the initial
157 Jacobian matrix from which the LMVM approximation is
158 built up. Inverse of this initial Jacobian is applied
159 using an internal KSP solver, which defaults to GMRES.
160 This internal KSP solver has the "mat_lmvm_" option
161 prefix.
162
163 Note that another LMVM matrix can be used in place of
164 J0, in which case updating the outer LMVM matrix will
165 also trigger the update for the inner LMVM matrix. This
166 is useful in cases where a full-memory diagonal approximation
167 such as MATLMVMDIAGBRDN is used in place of J0.
168
169 Input Parameters:
170 + B - An LMVM-type matrix
171 - J0 - The initial Jacobian matrix
172
173 Level: advanced
174
175 .seealso: MatLMVMSetJ0PC(), MatLMVMSetJ0KSP()
176 @*/
MatLMVMSetJ0(Mat B,Mat J0)177 PetscErrorCode MatLMVMSetJ0(Mat B, Mat J0)
178 {
179 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
180 PetscErrorCode ierr;
181 PetscBool same;
182 MPI_Comm comm = PetscObjectComm((PetscObject)B);
183
184 PetscFunctionBegin;
185 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
186 PetscValidHeaderSpecific(J0, MAT_CLASSID, 2);
187 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
188 if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
189 ierr = MatLMVMClearJ0(B);CHKERRQ(ierr);
190 ierr = MatDestroy(&lmvm->J0);CHKERRQ(ierr);
191 ierr = PetscObjectReference((PetscObject)J0);CHKERRQ(ierr);
192 lmvm->J0 = J0;
193 ierr = PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);CHKERRQ(ierr);
194 if (!same && lmvm->square) {
195 ierr = KSPSetOperators(lmvm->J0ksp, lmvm->J0, lmvm->J0);CHKERRQ(ierr);
196 }
197 PetscFunctionReturn(0);
198 }
199
200 /*------------------------------------------------------------*/
201
202 /*@
203 MatLMVMSetJ0PC - Allows the user to define a PC object that
204 acts as the initial inverse-Jacobian matrix. This PC should
205 already contain all the operators necessary for its application.
206 The LMVM matrix only calls PCApply() without changing any other
207 options.
208
209 Input Parameters:
210 + B - An LMVM-type matrix
211 - J0pc - PC object where PCApply defines an inverse application for J0
212
213 Level: advanced
214
215 .seealso: MatLMVMGetJ0PC()
216 @*/
MatLMVMSetJ0PC(Mat B,PC J0pc)217 PetscErrorCode MatLMVMSetJ0PC(Mat B, PC J0pc)
218 {
219 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
220 PetscErrorCode ierr;
221 PetscBool same;
222 MPI_Comm comm = PetscObjectComm((PetscObject)B);
223
224 PetscFunctionBegin;
225 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
226 PetscValidHeaderSpecific(J0pc, PC_CLASSID, 2);
227 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
228 if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
229 if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
230 ierr = MatLMVMClearJ0(B);CHKERRQ(ierr);
231 ierr = PetscObjectReference((PetscObject)J0pc);CHKERRQ(ierr);
232 lmvm->J0pc = J0pc;
233 lmvm->user_pc = PETSC_TRUE;
234 PetscFunctionReturn(0);
235 }
236
237 /*------------------------------------------------------------*/
238
239 /*@
240 MatLMVMSetJ0KSP - Allows the user to provide a pre-configured
241 KSP solver for the initial inverse-Jacobian approximation.
242 This KSP solver should already contain all the operators
243 necessary to perform the inversion. The LMVM matrix only
244 calls KSPSolve() without changing any other options.
245
246 Input Parameters:
247 + B - An LMVM-type matrix
248 - J0ksp - KSP solver for the initial inverse-Jacobian application
249
250 Level: advanced
251
252 .seealso: MatLMVMGetJ0KSP()
253 @*/
MatLMVMSetJ0KSP(Mat B,KSP J0ksp)254 PetscErrorCode MatLMVMSetJ0KSP(Mat B, KSP J0ksp)
255 {
256 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
257 PetscErrorCode ierr;
258 PetscBool same;
259 MPI_Comm comm = PetscObjectComm((PetscObject)B);
260
261 PetscFunctionBegin;
262 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
263 PetscValidHeaderSpecific(J0ksp, KSP_CLASSID, 2);
264 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
265 if (!same) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
266 if (!lmvm->square) SETERRQ(comm, PETSC_ERR_SUP, "Inverse J0 can be defined only for square LMVM matrices");
267 ierr = MatLMVMClearJ0(B);CHKERRQ(ierr);
268 ierr = KSPDestroy(&lmvm->J0ksp);CHKERRQ(ierr);
269 ierr = PetscObjectReference((PetscObject)J0ksp);CHKERRQ(ierr);
270 lmvm->J0ksp = J0ksp;
271 lmvm->user_ksp = PETSC_TRUE;
272 PetscFunctionReturn(0);
273 }
274
275 /*------------------------------------------------------------*/
276
277 /*@
278 MatLMVMGetJ0 - Returns a pointer to the internal J0 matrix.
279
280 Input Parameters:
281 . B - An LMVM-type matrix
282
283 Output Parameter:
284 . J0 - Mat object for defining the initial Jacobian
285
286 Level: advanced
287
288 .seealso: MatLMVMSetJ0()
289 @*/
MatLMVMGetJ0(Mat B,Mat * J0)290 PetscErrorCode MatLMVMGetJ0(Mat B, Mat *J0)
291 {
292 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
293 PetscErrorCode ierr;
294 PetscBool same;
295
296 PetscFunctionBegin;
297 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
298 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
299 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
300 *J0 = lmvm->J0;
301 PetscFunctionReturn(0);
302 }
303
304 /*------------------------------------------------------------*/
305
306 /*@
307 MatLMVMGetJ0PC - Returns a pointer to the internal PC object
308 associated with the initial Jacobian.
309
310 Input Parameters:
311 . B - An LMVM-type matrix
312
313 Output Parameter:
314 . J0pc - PC object for defining the initial inverse-Jacobian
315
316 Level: advanced
317
318 .seealso: MatLMVMSetJ0PC()
319 @*/
MatLMVMGetJ0PC(Mat B,PC * J0pc)320 PetscErrorCode MatLMVMGetJ0PC(Mat B, PC *J0pc)
321 {
322 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
323 PetscErrorCode ierr;
324 PetscBool same;
325
326 PetscFunctionBegin;
327 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
328 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
329 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
330 if (lmvm->J0pc) {
331 *J0pc = lmvm->J0pc;
332 } else {
333 ierr = KSPGetPC(lmvm->J0ksp, J0pc);CHKERRQ(ierr);
334 }
335 PetscFunctionReturn(0);
336 }
337
338 /*------------------------------------------------------------*/
339
340 /*@
341 MatLMVMGetJ0KSP - Returns a pointer to the internal KSP solver
342 associated with the initial Jacobian.
343
344 Input Parameters:
345 . B - An LMVM-type matrix
346
347 Output Parameter:
348 . J0ksp - KSP solver for defining the initial inverse-Jacobian
349
350 Level: advanced
351
352 .seealso: MatLMVMSetJ0KSP()
353 @*/
MatLMVMGetJ0KSP(Mat B,KSP * J0ksp)354 PetscErrorCode MatLMVMGetJ0KSP(Mat B, KSP *J0ksp)
355 {
356 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
357 PetscErrorCode ierr;
358 PetscBool same;
359
360 PetscFunctionBegin;
361 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
362 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
363 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
364 *J0ksp = lmvm->J0ksp;
365 PetscFunctionReturn(0);
366 }
367
368 /*------------------------------------------------------------*/
369
370 /*@
371 MatLMVMApplyJ0Fwd - Applies an approximation of the forward
372 matrix-vector product with the initial Jacobian.
373
374 Input Parameters:
375 + B - An LMVM-type matrix
376 - X - vector to multiply with J0
377
378 Output Parameter:
379 . Y - resulting vector for the operation
380
381 Level: advanced
382
383 .seealso: MatLMVMSetJ0(), MatLMVMSetJ0Scale(), MatLMVMSetJ0ScaleDiag(),
384 MatLMVMSetJ0PC(), MatLMVMSetJ0KSP(), MatLMVMApplyJ0Inv()
385 @*/
MatLMVMApplyJ0Fwd(Mat B,Vec X,Vec Y)386 PetscErrorCode MatLMVMApplyJ0Fwd(Mat B, Vec X, Vec Y)
387 {
388 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
389 PetscErrorCode ierr;
390 PetscBool same, hasMult;
391 MPI_Comm comm = PetscObjectComm((PetscObject)B);
392 Mat Amat, Pmat;
393
394 PetscFunctionBegin;
395 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
396 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
397 PetscValidHeaderSpecific(Y, VEC_CLASSID, 3);
398 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
399 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
400 if (!lmvm->allocated) SETERRQ(comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
401 VecCheckMatCompatible(B, X, 2, Y, 3);
402 if (lmvm->user_pc || lmvm->user_ksp || lmvm->J0) {
403 /* User may have defined a PC or KSP for J0^{-1} so let's try to use its operators. */
404 if (lmvm->user_pc){
405 ierr = PCGetOperators(lmvm->J0pc, &Amat, &Pmat);CHKERRQ(ierr);
406 } else if (lmvm->user_ksp) {
407 ierr = KSPGetOperators(lmvm->J0ksp, &Amat, &Pmat);CHKERRQ(ierr);
408 } else {
409 Amat = lmvm->J0;
410 }
411 ierr = MatHasOperation(Amat, MATOP_MULT, &hasMult);CHKERRQ(ierr);
412 if (hasMult) {
413 /* product is available, use it */
414 ierr = MatMult(Amat, X, Y);CHKERRQ(ierr);
415 } else {
416 /* there's no product, so treat J0 as identity */
417 ierr = VecCopy(X, Y);CHKERRQ(ierr);
418 }
419 } else if (lmvm->user_scale) {
420 if (lmvm->J0diag) {
421 /* User has defined a diagonal vector for J0 */
422 ierr = VecPointwiseMult(X, lmvm->J0diag, Y);CHKERRQ(ierr);
423 } else {
424 /* User has defined a scalar value for J0 */
425 ierr = VecCopy(X, Y);CHKERRQ(ierr);
426 ierr = VecScale(Y, lmvm->J0scalar);CHKERRQ(ierr);
427 }
428 } else {
429 /* There is no J0 representation so just apply an identity matrix */
430 ierr = VecCopy(X, Y);CHKERRQ(ierr);
431 }
432 PetscFunctionReturn(0);
433 }
434
435 /*------------------------------------------------------------*/
436
437 /*@
438 MatLMVMApplyJ0Inv - Applies some estimation of the initial Jacobian
439 inverse to the given vector. The specific form of the application
440 depends on whether the user provided a scaling factor, a J0 matrix,
441 a J0 PC, or a J0 KSP object. If no form of the initial Jacobian is
442 provided, the function simply does an identity matrix application
443 (vector copy).
444
445 Input Parameters:
446 + B - An LMVM-type matrix
447 - X - vector to "multiply" with J0^{-1}
448
449 Output Parameter:
450 . Y - resulting vector for the operation
451
452 Level: advanced
453
454 .seealso: MatLMVMSetJ0(), MatLMVMSetJ0Scale(), MatLMVMSetJ0ScaleDiag(),
455 MatLMVMSetJ0PC(), MatLMVMSetJ0KSP(), MatLMVMApplyJ0Fwd()
456 @*/
MatLMVMApplyJ0Inv(Mat B,Vec X,Vec Y)457 PetscErrorCode MatLMVMApplyJ0Inv(Mat B, Vec X, Vec Y)
458 {
459 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
460 PetscErrorCode ierr;
461 PetscBool same, hasSolve;
462 MPI_Comm comm = PetscObjectComm((PetscObject)B);
463
464 PetscFunctionBegin;
465 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
466 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
467 PetscValidHeaderSpecific(Y, VEC_CLASSID, 3);
468 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
469 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
470 if (!lmvm->allocated) SETERRQ(comm, PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
471 VecCheckMatCompatible(B, X, 2, Y, 3);
472 /* Invert the initial Jacobian onto q (or apply scaling) */
473 if (lmvm->user_pc) {
474 /* User has defined a J0 inverse so we can directly apply it as a preconditioner */
475 ierr = PCApply(lmvm->J0pc, X, Y);CHKERRQ(ierr);
476 } else if (lmvm->user_ksp) {
477 /* User has defined a J0 or a custom KSP so just perform a solution */
478 ierr = KSPSolve(lmvm->J0ksp, X, Y);CHKERRQ(ierr);
479 } else if (lmvm->J0) {
480 ierr = MatHasOperation(lmvm->J0, MATOP_SOLVE, &hasSolve);CHKERRQ(ierr);
481 if (hasSolve) {
482 ierr = MatSolve(lmvm->J0, X, Y);CHKERRQ(ierr);
483 } else {
484 ierr = KSPSolve(lmvm->J0ksp, X, Y);CHKERRQ(ierr);
485 }
486 } else if (lmvm->user_scale) {
487 if (lmvm->J0diag) {
488 ierr = VecPointwiseDivide(X, Y, lmvm->J0diag);CHKERRQ(ierr);
489 } else {
490 ierr = VecCopy(X, Y);CHKERRQ(ierr);
491 ierr = VecScale(Y, 1.0/lmvm->J0scalar);CHKERRQ(ierr);
492 }
493 } else {
494 /* There is no J0 representation so just apply an identity matrix */
495 ierr = VecCopy(X, Y);CHKERRQ(ierr);
496 }
497 PetscFunctionReturn(0);
498 }
499
500 /*------------------------------------------------------------*/
501
502 /*@
503 MatLMVMIsAllocated - Returns a boolean flag that shows whether
504 the necessary data structures for the underlying matrix is allocated.
505
506 Input Parameters:
507 . B - An LMVM-type matrix
508
509 Output Parameter:
510 . flg - PETSC_TRUE if allocated, PETSC_FALSE otherwise
511
512 Level: intermediate
513
514 .seealso: MatLMVMAllocate(), MatLMVMReset()
515 @*/
516
MatLMVMIsAllocated(Mat B,PetscBool * flg)517 PetscErrorCode MatLMVMIsAllocated(Mat B, PetscBool *flg)
518 {
519 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
520 PetscErrorCode ierr;
521 PetscBool same;
522
523 PetscFunctionBegin;
524 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
525 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
526 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
527 *flg = PETSC_FALSE;
528 if (lmvm->allocated && B->preallocated && B->assembled) *flg = PETSC_TRUE;
529 PetscFunctionReturn(0);
530 }
531
532 /*------------------------------------------------------------*/
533
534 /*@
535 MatLMVMAllocate - Produces all necessary common memory for
536 LMVM approximations based on the solution and function vectors
537 provided. If MatSetSizes() and MatSetUp() have not been called
538 before MatLMVMAllocate(), the allocation will read sizes from
539 the provided vectors and update the matrix.
540
541 Input Parameters:
542 + B - An LMVM-type matrix
543 . X - Solution vector
544 - F - Function vector
545
546 Level: intermediate
547
548 .seealso: MatLMVMReset(), MatLMVMUpdate()
549 @*/
MatLMVMAllocate(Mat B,Vec X,Vec F)550 PetscErrorCode MatLMVMAllocate(Mat B, Vec X, Vec F)
551 {
552 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
553 PetscErrorCode ierr;
554 PetscBool same;
555
556 PetscFunctionBegin;
557 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
558 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
559 PetscValidHeaderSpecific(F, VEC_CLASSID, 3);
560 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
561 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
562 ierr = (*lmvm->ops->allocate)(B, X, F);CHKERRQ(ierr);
563 if (lmvm->J0) {
564 ierr = PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);CHKERRQ(ierr);
565 if (same) {
566 ierr = MatLMVMAllocate(lmvm->J0, X, F);CHKERRQ(ierr);
567 }
568 }
569 PetscFunctionReturn(0);
570 }
571
572 /*------------------------------------------------------------*/
573
574 /*@
575 MatLMVMResetShift - Zero the shift factor.
576
577 Input Parameters:
578 . B - An LMVM-type matrix
579
580 Level: intermediate
581
582 .seealso: MatLMVMAllocate(), MatLMVMUpdate()
583 @*/
MatLMVMResetShift(Mat B)584 PetscErrorCode MatLMVMResetShift(Mat B)
585 {
586 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
587 PetscErrorCode ierr;
588 PetscBool same;
589
590 PetscFunctionBegin;
591 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
592 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
593 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
594 lmvm->shift = 0.0;
595 PetscFunctionReturn(0);
596 }
597
598 /*------------------------------------------------------------*/
599
600 /*@
601 MatLMVMReset - Flushes all of the accumulated updates out of
602 the LMVM approximation. In practice, this will not actually
603 destroy the data associated with the updates. It simply resets
604 counters, which leads to existing data being overwritten, and
605 MatSolve() being applied as if there are no updates. A boolean
606 flag is available to force destruction of the update vectors.
607
608 If the user has provided another LMVM matrix as J0, the J0
609 matrix is also reset in this function.
610
611 Input Parameters:
612 + B - An LMVM-type matrix
613 - destructive - flag for enabling destruction of data structures
614
615 Level: intermediate
616
617 .seealso: MatLMVMAllocate(), MatLMVMUpdate()
618 @*/
MatLMVMReset(Mat B,PetscBool destructive)619 PetscErrorCode MatLMVMReset(Mat B, PetscBool destructive)
620 {
621 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
622 PetscErrorCode ierr;
623 PetscBool same;
624
625 PetscFunctionBegin;
626 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
627 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
628 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
629 ierr = (*lmvm->ops->reset)(B, destructive);CHKERRQ(ierr);
630 if (lmvm->J0) {
631 ierr = PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &same);CHKERRQ(ierr);
632 if (same) {
633 ierr = MatLMVMReset(lmvm->J0, destructive);CHKERRQ(ierr);
634 }
635 }
636 PetscFunctionReturn(0);
637 }
638
639 /*------------------------------------------------------------*/
640
641 /*@
642 MatLMVMSetHistorySize - Set the number of past iterates to be
643 stored for the construction of the limited-memory QN update.
644
645 Input Parameters:
646 + B - An LMVM-type matrix
647 - hist_size - number of past iterates (default 5)
648
649 Options Database:
650 . -mat_lmvm_hist_size <m>
651
652 Level: beginner
653
654 .seealso: MatLMVMGetUpdateCount()
655 @*/
656
MatLMVMSetHistorySize(Mat B,PetscInt hist_size)657 PetscErrorCode MatLMVMSetHistorySize(Mat B, PetscInt hist_size)
658 {
659 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
660 PetscErrorCode ierr;
661 PetscBool same;
662 Vec X, F;
663
664 PetscFunctionBegin;
665 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
666 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
667 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
668 if (hist_size > 0) {
669 lmvm->m = hist_size;
670 if (lmvm->allocated && lmvm->m != lmvm->m_old) {
671 ierr = VecDuplicate(lmvm->Xprev, &X);CHKERRQ(ierr);
672 ierr = VecDuplicate(lmvm->Fprev, &F);CHKERRQ(ierr);
673 ierr = MatLMVMReset(B, PETSC_TRUE);CHKERRQ(ierr);
674 ierr = MatLMVMAllocate(B, X, F);CHKERRQ(ierr);
675 ierr = VecDestroy(&X);CHKERRQ(ierr);
676 ierr = VecDestroy(&F);CHKERRQ(ierr);
677 }
678 } else {
679 SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "QN history size must be a positive integer.");
680 }
681
682 PetscFunctionReturn(0);
683 }
684
685 /*------------------------------------------------------------*/
686
687 /*@
688 MatLMVMGetUpdateCount - Returns the number of accepted updates.
689 This number may be greater than the total number of update vectors
690 stored in the matrix. The counters are reset when MatLMVMReset()
691 is called.
692
693 Input Parameters:
694 . B - An LMVM-type matrix
695
696 Output Parameter:
697 . nupdates - number of accepted updates
698
699 Level: intermediate
700
701 .seealso: MatLMVMGetRejectCount(), MatLMVMReset()
702 @*/
MatLMVMGetUpdateCount(Mat B,PetscInt * nupdates)703 PetscErrorCode MatLMVMGetUpdateCount(Mat B, PetscInt *nupdates)
704 {
705 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
706 PetscErrorCode ierr;
707 PetscBool same;
708
709 PetscFunctionBegin;
710 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
711 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
712 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
713 *nupdates = lmvm->nupdates;
714 PetscFunctionReturn(0);
715 }
716
717 /*------------------------------------------------------------*/
718
719 /*@
720 MatLMVMGetRejectCount - Returns the number of rejected updates.
721 The counters are reset when MatLMVMReset() is called.
722
723 Input Parameters:
724 . B - An LMVM-type matrix
725
726 Output Parameter:
727 . nrejects - number of rejected updates
728
729 Level: intermediate
730
731 .seealso: MatLMVMGetRejectCount(), MatLMVMReset()
732 @*/
MatLMVMGetRejectCount(Mat B,PetscInt * nrejects)733 PetscErrorCode MatLMVMGetRejectCount(Mat B, PetscInt *nrejects)
734 {
735 Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
736 PetscErrorCode ierr;
737 PetscBool same;
738
739 PetscFunctionBegin;
740 PetscValidHeaderSpecific(B, MAT_CLASSID, 1);
741 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr);
742 if (!same) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
743 *nrejects = lmvm->nrejects;
744 PetscFunctionReturn(0);
745 }
746