Actual source code: lmvmimpl.c
1: #include <petscdevice.h>
2: #include <../src/ksp/ksp/utils/lmvm/lmvm.h>
3: #include <petsc/private/deviceimpl.h>
5: PetscErrorCode MatReset_LMVM(Mat B, PetscBool destructive)
6: {
7: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
9: PetscFunctionBegin;
10: lmvm->k = -1;
11: lmvm->prev_set = PETSC_FALSE;
12: lmvm->shift = 0.0;
13: if (destructive && lmvm->allocated) {
14: PetscCall(MatLMVMClearJ0(B));
15: PetscCall(VecDestroyVecs(lmvm->m, &lmvm->S));
16: PetscCall(VecDestroyVecs(lmvm->m, &lmvm->Y));
17: PetscCall(VecDestroy(&lmvm->Xprev));
18: PetscCall(VecDestroy(&lmvm->Fprev));
19: lmvm->nupdates = 0;
20: lmvm->nrejects = 0;
21: lmvm->m_old = 0;
22: lmvm->allocated = PETSC_FALSE;
23: B->preallocated = PETSC_FALSE;
24: B->assembled = PETSC_FALSE;
25: }
26: ++lmvm->nresets;
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: PetscErrorCode MatAllocate_LMVM(Mat B, Vec X, Vec F)
31: {
32: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
33: PetscBool same, allocate = PETSC_FALSE;
34: VecType vtype;
36: PetscFunctionBegin;
37: if (lmvm->allocated) {
38: VecCheckMatCompatible(B, X, 2, F, 3);
39: PetscCall(VecGetType(X, &vtype));
40: PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->Xprev, vtype, &same));
41: if (!same) {
42: /* Given X vector has a different type than allocated X-type data structures.
43: We need to destroy all of this and duplicate again out of the given vector. */
44: allocate = PETSC_TRUE;
45: PetscCall(MatLMVMReset(B, PETSC_TRUE));
46: }
47: } else allocate = PETSC_TRUE;
48: if (allocate) {
49: PetscCall(VecGetType(X, &vtype));
50: PetscCall(MatSetVecType(B, vtype));
51: PetscCall(PetscLayoutReference(F->map, &B->rmap));
52: PetscCall(PetscLayoutReference(X->map, &B->cmap));
53: PetscCall(VecDuplicate(X, &lmvm->Xprev));
54: PetscCall(VecDuplicate(F, &lmvm->Fprev));
55: if (lmvm->m > 0) {
56: PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S));
57: PetscCall(VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y));
58: }
59: lmvm->m_old = lmvm->m;
60: lmvm->allocated = PETSC_TRUE;
61: B->preallocated = PETSC_TRUE;
62: B->assembled = PETSC_TRUE;
63: }
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
68: {
69: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
70: PetscInt i;
71: Vec Stmp, Ytmp;
73: PetscFunctionBegin;
74: if (lmvm->k == lmvm->m - 1) {
75: /* We hit the memory limit, so shift all the vectors back one spot
76: and shift the oldest to the front to receive the latest update. */
77: Stmp = lmvm->S[0];
78: Ytmp = lmvm->Y[0];
79: for (i = 0; i < lmvm->k; ++i) {
80: lmvm->S[i] = lmvm->S[i + 1];
81: lmvm->Y[i] = lmvm->Y[i + 1];
82: }
83: lmvm->S[lmvm->k] = Stmp;
84: lmvm->Y[lmvm->k] = Ytmp;
85: } else {
86: ++lmvm->k;
87: }
88: /* Put the precomputed update into the last vector */
89: PetscCall(VecCopy(S, lmvm->S[lmvm->k]));
90: PetscCall(VecCopy(Y, lmvm->Y[lmvm->k]));
91: ++lmvm->nupdates;
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
96: {
97: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
99: PetscFunctionBegin;
100: if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
101: if (lmvm->prev_set) {
102: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
103: PetscCall(VecAXPBY(lmvm->Xprev, 1.0, -1.0, X));
104: PetscCall(VecAXPBY(lmvm->Fprev, 1.0, -1.0, F));
105: /* Update S and Y */
106: PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
107: }
109: /* Save the solution and function to be used in the next update */
110: PetscCall(VecCopy(X, lmvm->Xprev));
111: PetscCall(VecCopy(F, lmvm->Fprev));
112: lmvm->prev_set = PETSC_TRUE;
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
117: {
118: PetscFunctionBegin;
119: PetscCall(MatMult(B, X, Z));
120: PetscCall(VecAXPY(Z, 1.0, Y));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
125: {
126: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
128: PetscFunctionBegin;
129: VecCheckSameSize(X, 2, Y, 3);
130: VecCheckMatCompatible(B, X, 2, Y, 3);
131: PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
132: PetscCall((*lmvm->ops->mult)(B, X, Y));
133: if (lmvm->shift) PetscCall(VecAXPY(Y, lmvm->shift, X));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode MatSolve_LMVM(Mat B, Vec F, Vec dX)
138: {
139: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
141: PetscFunctionBegin;
142: VecCheckSameSize(F, 2, dX, 3);
143: VecCheckMatCompatible(B, F, 2, dX, 3);
144: PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
145: PetscCheck(*lmvm->ops->solve, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "LMVM matrix does not have a solution or inversion implementation");
146: PetscCall((*lmvm->ops->solve)(B, F, dX));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode MatCopy_LMVM(Mat B, Mat M, MatStructure str)
151: {
152: Mat_LMVM *bctx = (Mat_LMVM *)B->data;
153: Mat_LMVM *mctx;
154: PetscInt i;
155: PetscBool allocatedM;
157: PetscFunctionBegin;
158: if (str == DIFFERENT_NONZERO_PATTERN) {
159: PetscCall(MatLMVMReset(M, PETSC_TRUE));
160: PetscCall(MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev));
161: } else {
162: PetscCall(MatLMVMIsAllocated(M, &allocatedM));
163: PetscCheck(allocatedM, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Target matrix must be allocated first");
164: MatCheckSameSize(B, 1, M, 2);
165: }
167: mctx = (Mat_LMVM *)M->data;
168: if (bctx->user_pc) {
169: PetscCall(MatLMVMSetJ0PC(M, bctx->J0pc));
170: } else if (bctx->user_ksp) {
171: PetscCall(MatLMVMSetJ0KSP(M, bctx->J0ksp));
172: } else if (bctx->J0) {
173: PetscCall(MatLMVMSetJ0(M, bctx->J0));
174: } else if (bctx->user_scale) {
175: if (bctx->J0diag) {
176: PetscCall(MatLMVMSetJ0Diag(M, bctx->J0diag));
177: } else {
178: PetscCall(MatLMVMSetJ0Scale(M, bctx->J0scalar));
179: }
180: }
181: mctx->nupdates = bctx->nupdates;
182: mctx->nrejects = bctx->nrejects;
183: mctx->k = bctx->k;
184: for (i = 0; i <= bctx->k; ++i) {
185: if (bctx->S) PetscCall(VecCopy(bctx->S[i], mctx->S[i]));
186: if (bctx->Y) PetscCall(VecCopy(bctx->Y[i], mctx->Y[i]));
187: PetscCall(VecCopy(bctx->Xprev, mctx->Xprev));
188: PetscCall(VecCopy(bctx->Fprev, mctx->Fprev));
189: }
190: if (bctx->ops->copy) PetscCall((*bctx->ops->copy)(B, M, str));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
195: {
196: Mat_LMVM *bctx = (Mat_LMVM *)B->data;
197: Mat_LMVM *mctx;
198: MatType lmvmType;
199: Mat A;
201: PetscFunctionBegin;
202: PetscCall(MatGetType(B, &lmvmType));
203: PetscCall(MatCreate(PetscObjectComm((PetscObject)B), mat));
204: PetscCall(MatSetType(*mat, lmvmType));
206: A = *mat;
207: mctx = (Mat_LMVM *)A->data;
208: mctx->m = bctx->m;
209: mctx->ksp_max_it = bctx->ksp_max_it;
210: mctx->ksp_rtol = bctx->ksp_rtol;
211: mctx->ksp_atol = bctx->ksp_atol;
212: mctx->shift = bctx->shift;
213: PetscCall(KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it));
215: PetscCall(MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev));
216: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(B, *mat, SAME_NONZERO_PATTERN));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: static PetscErrorCode MatShift_LMVM(Mat B, PetscScalar a)
221: {
222: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
224: PetscFunctionBegin;
225: PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
226: lmvm->shift += PetscRealPart(a);
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
231: {
232: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
233: PetscBool isascii;
234: MatType type;
236: PetscFunctionBegin;
237: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
238: if (isascii) {
239: PetscCall(MatGetType(B, &type));
240: PetscCall(PetscViewerASCIIPrintf(pv, "Max. storage: %" PetscInt_FMT "\n", lmvm->m));
241: PetscCall(PetscViewerASCIIPrintf(pv, "Used storage: %" PetscInt_FMT "\n", lmvm->k + 1));
242: PetscCall(PetscViewerASCIIPrintf(pv, "Number of updates: %" PetscInt_FMT "\n", lmvm->nupdates));
243: PetscCall(PetscViewerASCIIPrintf(pv, "Number of rejects: %" PetscInt_FMT "\n", lmvm->nrejects));
244: PetscCall(PetscViewerASCIIPrintf(pv, "Number of resets: %" PetscInt_FMT "\n", lmvm->nresets));
245: if (lmvm->J0) {
246: PetscCall(PetscViewerASCIIPrintf(pv, "J0 Matrix:\n"));
247: PetscCall(PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO));
248: PetscCall(MatView(lmvm->J0, pv));
249: PetscCall(PetscViewerPopFormat(pv));
250: }
251: }
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: PetscErrorCode MatSetFromOptions_LMVM(Mat B, PetscOptionItems *PetscOptionsObject)
256: {
257: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
258: PetscInt m_new = lmvm->m;
260: PetscFunctionBegin;
261: PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory Variable Metric matrix for approximating Jacobians");
262: PetscCall(PetscOptionsInt("-mat_lmvm_hist_size", "number of past updates kept in memory for the approximation", "", m_new, &m_new, NULL));
263: PetscCall(PetscOptionsInt("-mat_lmvm_ksp_its", "(developer) fixed number of KSP iterations to take when inverting J0", "", lmvm->ksp_max_it, &lmvm->ksp_max_it, NULL));
264: PetscCall(PetscOptionsReal("-mat_lmvm_eps", "(developer) machine zero definition", "", lmvm->eps, &lmvm->eps, NULL));
265: PetscOptionsHeadEnd();
266: if (m_new != lmvm->m) PetscCall(MatLMVMSetHistorySize(B, m_new));
267: PetscCall(KSPSetFromOptions(lmvm->J0ksp));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: PetscErrorCode MatSetUp_LMVM(Mat B)
272: {
273: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
275: PetscFunctionBegin;
276: if (!lmvm->allocated) {
277: PetscCall(PetscLayoutSetUp(B->rmap));
278: PetscCall(PetscLayoutSetUp(B->cmap));
279: PetscCall(MatCreateVecs(B, &lmvm->Xprev, &lmvm->Fprev));
280: if (lmvm->m > 0) {
281: PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S));
282: PetscCall(VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y));
283: }
284: lmvm->m_old = lmvm->m;
285: lmvm->allocated = PETSC_TRUE;
286: B->preallocated = PETSC_TRUE;
287: B->assembled = PETSC_TRUE;
288: }
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: PetscErrorCode MatDestroy_LMVM(Mat B)
293: {
294: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
296: PetscFunctionBegin;
297: if (lmvm->allocated) {
298: PetscCall(VecDestroyVecs(lmvm->m, &lmvm->S));
299: PetscCall(VecDestroyVecs(lmvm->m, &lmvm->Y));
300: PetscCall(VecDestroy(&lmvm->Xprev));
301: PetscCall(VecDestroy(&lmvm->Fprev));
302: }
303: PetscCall(KSPDestroy(&lmvm->J0ksp));
304: PetscCall(MatLMVMClearJ0(B));
305: PetscCall(PetscFree(B->data));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*MC
310: MATLMVM - MATLMVM = "lmvm" - A matrix type used for Limited-Memory Variable Metric (LMVM) matrices.
312: Level: intermediate
314: Developer notes:
315: Improve this manual page as well as many others in the MATLMVM family.
317: .seealso: [](sec_matlmvm), `Mat`
318: M*/
319: PetscErrorCode MatCreate_LMVM(Mat B)
320: {
321: Mat_LMVM *lmvm;
323: PetscFunctionBegin;
324: PetscCall(PetscNew(&lmvm));
325: B->data = (void *)lmvm;
327: lmvm->m_old = 0;
328: lmvm->m = 5;
329: lmvm->k = -1;
330: lmvm->nupdates = 0;
331: lmvm->nrejects = 0;
332: lmvm->nresets = 0;
334: lmvm->ksp_max_it = 20;
335: lmvm->ksp_rtol = 0.0;
336: lmvm->ksp_atol = 0.0;
338: lmvm->shift = 0.0;
340: lmvm->eps = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0);
341: lmvm->allocated = PETSC_FALSE;
342: lmvm->prev_set = PETSC_FALSE;
343: lmvm->user_scale = PETSC_FALSE;
344: lmvm->user_pc = PETSC_FALSE;
345: lmvm->user_ksp = PETSC_FALSE;
346: lmvm->square = PETSC_FALSE;
348: B->ops->destroy = MatDestroy_LMVM;
349: B->ops->setfromoptions = MatSetFromOptions_LMVM;
350: B->ops->view = MatView_LMVM;
351: B->ops->setup = MatSetUp_LMVM;
352: B->ops->shift = MatShift_LMVM;
353: B->ops->duplicate = MatDuplicate_LMVM;
354: B->ops->mult = MatMult_LMVM;
355: B->ops->multadd = MatMultAdd_LMVM;
356: B->ops->solve = MatSolve_LMVM;
357: B->ops->copy = MatCopy_LMVM;
359: lmvm->ops->update = MatUpdate_LMVM;
360: lmvm->ops->allocate = MatAllocate_LMVM;
361: lmvm->ops->reset = MatReset_LMVM;
363: PetscCall(KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp));
364: PetscCall(PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1));
365: PetscCall(KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_"));
366: PetscCall(KSPSetType(lmvm->J0ksp, KSPGMRES));
367: PetscCall(KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }