Actual source code: lmvmimpl.c
1: #include <../src/ksp/ksp/utils/lmvm/lmvm.h>
3: PetscErrorCode MatReset_LMVM(Mat B, PetscBool destructive)
4: {
5: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
7: PetscFunctionBegin;
8: lmvm->k = -1;
9: lmvm->prev_set = PETSC_FALSE;
10: lmvm->shift = 0.0;
11: if (destructive && lmvm->allocated) {
12: PetscCall(MatLMVMClearJ0(B));
13: PetscCall(VecDestroyVecs(lmvm->m, &lmvm->S));
14: PetscCall(VecDestroyVecs(lmvm->m, &lmvm->Y));
15: PetscCall(VecDestroy(&lmvm->Xprev));
16: PetscCall(VecDestroy(&lmvm->Fprev));
17: lmvm->nupdates = 0;
18: lmvm->nrejects = 0;
19: lmvm->m_old = 0;
20: lmvm->allocated = PETSC_FALSE;
21: B->preallocated = PETSC_FALSE;
22: B->assembled = PETSC_FALSE;
23: }
24: ++lmvm->nresets;
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: PetscErrorCode MatAllocate_LMVM(Mat B, Vec X, Vec F)
29: {
30: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
31: PetscBool same, allocate = PETSC_FALSE;
32: VecType vtype;
34: PetscFunctionBegin;
35: if (lmvm->allocated) {
36: VecCheckMatCompatible(B, X, 2, F, 3);
37: PetscCall(VecGetType(X, &vtype));
38: PetscCall(PetscObjectTypeCompare((PetscObject)lmvm->Xprev, vtype, &same));
39: if (!same) {
40: /* Given X vector has a different type than allocated X-type data structures.
41: We need to destroy all of this and duplicate again out of the given vector. */
42: allocate = PETSC_TRUE;
43: PetscCall(MatLMVMReset(B, PETSC_TRUE));
44: }
45: } else allocate = PETSC_TRUE;
46: if (allocate) {
47: PetscCall(VecGetType(X, &vtype));
48: PetscCall(MatSetVecType(B, vtype));
49: PetscCall(PetscLayoutReference(F->map, &B->rmap));
50: PetscCall(PetscLayoutReference(X->map, &B->cmap));
51: PetscCall(VecDuplicate(X, &lmvm->Xprev));
52: PetscCall(VecDuplicate(F, &lmvm->Fprev));
53: if (lmvm->m > 0) {
54: PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S));
55: PetscCall(VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y));
56: }
57: lmvm->m_old = lmvm->m;
58: lmvm->allocated = PETSC_TRUE;
59: B->preallocated = PETSC_TRUE;
60: B->assembled = PETSC_TRUE;
61: }
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
66: {
67: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
68: PetscInt i;
69: Vec Stmp, Ytmp;
71: PetscFunctionBegin;
72: if (lmvm->k == lmvm->m - 1) {
73: /* We hit the memory limit, so shift all the vectors back one spot
74: and shift the oldest to the front to receive the latest update. */
75: Stmp = lmvm->S[0];
76: Ytmp = lmvm->Y[0];
77: for (i = 0; i < lmvm->k; ++i) {
78: lmvm->S[i] = lmvm->S[i + 1];
79: lmvm->Y[i] = lmvm->Y[i + 1];
80: }
81: lmvm->S[lmvm->k] = Stmp;
82: lmvm->Y[lmvm->k] = Ytmp;
83: } else {
84: ++lmvm->k;
85: }
86: /* Put the precomputed update into the last vector */
87: PetscCall(VecCopy(S, lmvm->S[lmvm->k]));
88: PetscCall(VecCopy(Y, lmvm->Y[lmvm->k]));
89: ++lmvm->nupdates;
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
94: {
95: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
97: PetscFunctionBegin;
98: if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
99: if (lmvm->prev_set) {
100: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
101: PetscCall(VecAXPBY(lmvm->Xprev, 1.0, -1.0, X));
102: PetscCall(VecAXPBY(lmvm->Fprev, 1.0, -1.0, F));
103: /* Update S and Y */
104: PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
105: }
107: /* Save the solution and function to be used in the next update */
108: PetscCall(VecCopy(X, lmvm->Xprev));
109: PetscCall(VecCopy(F, lmvm->Fprev));
110: lmvm->prev_set = PETSC_TRUE;
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
115: {
116: PetscFunctionBegin;
117: PetscCall(MatMult(B, X, Z));
118: PetscCall(VecAXPY(Z, 1.0, Y));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
123: {
124: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
126: PetscFunctionBegin;
127: VecCheckSameSize(X, 2, Y, 3);
128: VecCheckMatCompatible(B, X, 2, Y, 3);
129: PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
130: PetscCall((*lmvm->ops->mult)(B, X, Y));
131: PetscCall(VecAXPY(Y, lmvm->shift, X));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: static PetscErrorCode MatCopy_LMVM(Mat B, Mat M, MatStructure str)
136: {
137: Mat_LMVM *bctx = (Mat_LMVM *)B->data;
138: Mat_LMVM *mctx;
139: PetscInt i;
140: PetscBool allocatedM;
142: PetscFunctionBegin;
143: if (str == DIFFERENT_NONZERO_PATTERN) {
144: PetscCall(MatLMVMReset(M, PETSC_TRUE));
145: PetscCall(MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev));
146: } else {
147: PetscCall(MatLMVMIsAllocated(M, &allocatedM));
148: PetscCheck(allocatedM, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Target matrix must be allocated first");
149: MatCheckSameSize(B, 1, M, 2);
150: }
152: mctx = (Mat_LMVM *)M->data;
153: if (bctx->user_pc) {
154: PetscCall(MatLMVMSetJ0PC(M, bctx->J0pc));
155: } else if (bctx->user_ksp) {
156: PetscCall(MatLMVMSetJ0KSP(M, bctx->J0ksp));
157: } else if (bctx->J0) {
158: PetscCall(MatLMVMSetJ0(M, bctx->J0));
159: } else if (bctx->user_scale) {
160: if (bctx->J0diag) {
161: PetscCall(MatLMVMSetJ0Diag(M, bctx->J0diag));
162: } else {
163: PetscCall(MatLMVMSetJ0Scale(M, bctx->J0scalar));
164: }
165: }
166: mctx->nupdates = bctx->nupdates;
167: mctx->nrejects = bctx->nrejects;
168: mctx->k = bctx->k;
169: for (i = 0; i <= bctx->k; ++i) {
170: PetscCall(VecCopy(bctx->S[i], mctx->S[i]));
171: PetscCall(VecCopy(bctx->Y[i], mctx->Y[i]));
172: PetscCall(VecCopy(bctx->Xprev, mctx->Xprev));
173: PetscCall(VecCopy(bctx->Fprev, mctx->Fprev));
174: }
175: if (bctx->ops->copy) PetscCall((*bctx->ops->copy)(B, M, str));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
180: {
181: Mat_LMVM *bctx = (Mat_LMVM *)B->data;
182: Mat_LMVM *mctx;
183: MatType lmvmType;
184: Mat A;
186: PetscFunctionBegin;
187: PetscCall(MatGetType(B, &lmvmType));
188: PetscCall(MatCreate(PetscObjectComm((PetscObject)B), mat));
189: PetscCall(MatSetType(*mat, lmvmType));
191: A = *mat;
192: mctx = (Mat_LMVM *)A->data;
193: mctx->m = bctx->m;
194: mctx->ksp_max_it = bctx->ksp_max_it;
195: mctx->ksp_rtol = bctx->ksp_rtol;
196: mctx->ksp_atol = bctx->ksp_atol;
197: mctx->shift = bctx->shift;
198: PetscCall(KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it));
200: PetscCall(MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev));
201: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(B, *mat, SAME_NONZERO_PATTERN));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: static PetscErrorCode MatShift_LMVM(Mat B, PetscScalar a)
206: {
207: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
209: PetscFunctionBegin;
210: PetscCheck(lmvm->allocated, PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
211: lmvm->shift += PetscRealPart(a);
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
216: {
217: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
218: PetscBool isascii;
219: MatType type;
221: PetscFunctionBegin;
222: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
223: if (isascii) {
224: PetscCall(MatGetType(B, &type));
225: PetscCall(PetscViewerASCIIPrintf(pv, "Max. storage: %" PetscInt_FMT "\n", lmvm->m));
226: PetscCall(PetscViewerASCIIPrintf(pv, "Used storage: %" PetscInt_FMT "\n", lmvm->k + 1));
227: PetscCall(PetscViewerASCIIPrintf(pv, "Number of updates: %" PetscInt_FMT "\n", lmvm->nupdates));
228: PetscCall(PetscViewerASCIIPrintf(pv, "Number of rejects: %" PetscInt_FMT "\n", lmvm->nrejects));
229: PetscCall(PetscViewerASCIIPrintf(pv, "Number of resets: %" PetscInt_FMT "\n", lmvm->nresets));
230: if (lmvm->J0) {
231: PetscCall(PetscViewerASCIIPrintf(pv, "J0 Matrix:\n"));
232: PetscCall(PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO));
233: PetscCall(MatView(lmvm->J0, pv));
234: PetscCall(PetscViewerPopFormat(pv));
235: }
236: }
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: PetscErrorCode MatSetFromOptions_LMVM(Mat B, PetscOptionItems *PetscOptionsObject)
241: {
242: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
244: PetscFunctionBegin;
245: PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory Variable Metric matrix for approximating Jacobians");
246: PetscCall(PetscOptionsInt("-mat_lmvm_hist_size", "number of past updates kept in memory for the approximation", "", lmvm->m, &lmvm->m, NULL));
247: 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));
248: PetscCall(PetscOptionsReal("-mat_lmvm_eps", "(developer) machine zero definition", "", lmvm->eps, &lmvm->eps, NULL));
249: PetscOptionsHeadEnd();
250: PetscCall(KSPSetFromOptions(lmvm->J0ksp));
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: PetscErrorCode MatSetUp_LMVM(Mat B)
255: {
256: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
258: PetscFunctionBegin;
259: if (!lmvm->allocated) {
260: PetscCall(PetscLayoutSetUp(B->rmap));
261: PetscCall(PetscLayoutSetUp(B->cmap));
262: PetscCall(MatCreateVecs(B, &lmvm->Xprev, &lmvm->Fprev));
263: if (lmvm->m > 0) {
264: PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S));
265: PetscCall(VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y));
266: }
267: lmvm->m_old = lmvm->m;
268: lmvm->allocated = PETSC_TRUE;
269: B->preallocated = PETSC_TRUE;
270: B->assembled = PETSC_TRUE;
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: PetscErrorCode MatDestroy_LMVM(Mat B)
276: {
277: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
279: PetscFunctionBegin;
280: if (lmvm->allocated) {
281: PetscCall(VecDestroyVecs(lmvm->m, &lmvm->S));
282: PetscCall(VecDestroyVecs(lmvm->m, &lmvm->Y));
283: PetscCall(VecDestroy(&lmvm->Xprev));
284: PetscCall(VecDestroy(&lmvm->Fprev));
285: }
286: PetscCall(KSPDestroy(&lmvm->J0ksp));
287: PetscCall(MatLMVMClearJ0(B));
288: PetscCall(PetscFree(B->data));
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: PetscErrorCode MatCreate_LMVM(Mat B)
293: {
294: Mat_LMVM *lmvm;
296: PetscFunctionBegin;
297: PetscCall(PetscNew(&lmvm));
298: B->data = (void *)lmvm;
300: lmvm->m_old = 0;
301: lmvm->m = 5;
302: lmvm->k = -1;
303: lmvm->nupdates = 0;
304: lmvm->nrejects = 0;
305: lmvm->nresets = 0;
307: lmvm->ksp_max_it = 20;
308: lmvm->ksp_rtol = 0.0;
309: lmvm->ksp_atol = 0.0;
311: lmvm->shift = 0.0;
313: lmvm->eps = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0);
314: lmvm->allocated = PETSC_FALSE;
315: lmvm->prev_set = PETSC_FALSE;
316: lmvm->user_scale = PETSC_FALSE;
317: lmvm->user_pc = PETSC_FALSE;
318: lmvm->user_ksp = PETSC_FALSE;
319: lmvm->square = PETSC_FALSE;
321: B->ops->destroy = MatDestroy_LMVM;
322: B->ops->setfromoptions = MatSetFromOptions_LMVM;
323: B->ops->view = MatView_LMVM;
324: B->ops->setup = MatSetUp_LMVM;
325: B->ops->shift = MatShift_LMVM;
326: B->ops->duplicate = MatDuplicate_LMVM;
327: B->ops->mult = MatMult_LMVM;
328: B->ops->multadd = MatMultAdd_LMVM;
329: B->ops->copy = MatCopy_LMVM;
331: lmvm->ops->update = MatUpdate_LMVM;
332: lmvm->ops->allocate = MatAllocate_LMVM;
333: lmvm->ops->reset = MatReset_LMVM;
335: PetscCall(KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp));
336: PetscCall(PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1));
337: PetscCall(KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_"));
338: PetscCall(KSPSetType(lmvm->J0ksp, KSPGMRES));
339: PetscCall(KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }