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: lmvm->k = -1;
8: lmvm->prev_set = PETSC_FALSE;
9: lmvm->shift = 0.0;
10: if (destructive && lmvm->allocated) {
11: MatLMVMClearJ0(B);
12: B->rmap->n = B->rmap->N = B->cmap->n = B->cmap->N = 0;
13: VecDestroyVecs(lmvm->m, &lmvm->S);
14: VecDestroyVecs(lmvm->m, &lmvm->Y);
15: VecDestroy(&lmvm->Xprev);
16: 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: return 0;
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: PetscInt m, n, M, N;
33: VecType type;
35: if (lmvm->allocated) {
36: VecCheckMatCompatible(B, X, 2, F, 3);
37: VecGetType(X, &type);
38: PetscObjectTypeCompare((PetscObject)lmvm->Xprev, type, &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: MatLMVMReset(B, PETSC_TRUE);
44: }
45: } else {
46: allocate = PETSC_TRUE;
47: }
48: if (allocate) {
49: VecGetLocalSize(X, &n);
50: VecGetSize(X, &N);
51: VecGetLocalSize(F, &m);
52: VecGetSize(F, &M);
53: B->rmap->n = m;
54: B->cmap->n = n;
55: B->rmap->N = M > -1 ? M : B->rmap->N;
56: B->cmap->N = N > -1 ? N : B->cmap->N;
57: VecDuplicate(X, &lmvm->Xprev);
58: VecDuplicate(F, &lmvm->Fprev);
59: if (lmvm->m > 0) {
60: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
61: VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);
62: }
63: lmvm->m_old = lmvm->m;
64: lmvm->allocated = PETSC_TRUE;
65: B->preallocated = PETSC_TRUE;
66: B->assembled = PETSC_TRUE;
67: }
68: return 0;
69: }
71: PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
72: {
73: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
74: PetscInt i;
75: Vec Stmp, Ytmp;
77: if (lmvm->k == lmvm->m - 1) {
78: /* We hit the memory limit, so shift all the vectors back one spot
79: and shift the oldest to the front to receive the latest update. */
80: Stmp = lmvm->S[0];
81: Ytmp = lmvm->Y[0];
82: for (i = 0; i < lmvm->k; ++i) {
83: lmvm->S[i] = lmvm->S[i + 1];
84: lmvm->Y[i] = lmvm->Y[i + 1];
85: }
86: lmvm->S[lmvm->k] = Stmp;
87: lmvm->Y[lmvm->k] = Ytmp;
88: } else {
89: ++lmvm->k;
90: }
91: /* Put the precomputed update into the last vector */
92: VecCopy(S, lmvm->S[lmvm->k]);
93: VecCopy(Y, lmvm->Y[lmvm->k]);
94: ++lmvm->nupdates;
95: return 0;
96: }
98: PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
99: {
100: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
102: if (!lmvm->m) return 0;
103: if (lmvm->prev_set) {
104: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
105: VecAXPBY(lmvm->Xprev, 1.0, -1.0, X);
106: VecAXPBY(lmvm->Fprev, 1.0, -1.0, F);
107: /* Update S and Y */
108: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
109: }
111: /* Save the solution and function to be used in the next update */
112: VecCopy(X, lmvm->Xprev);
113: VecCopy(F, lmvm->Fprev);
114: lmvm->prev_set = PETSC_TRUE;
115: return 0;
116: }
118: static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
119: {
120: MatMult(B, X, Z);
121: VecAXPY(Z, 1.0, Y);
122: return 0;
123: }
125: static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
126: {
127: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
129: VecCheckSameSize(X, 2, Y, 3);
130: VecCheckMatCompatible(B, X, 2, Y, 3);
132: (*lmvm->ops->mult)(B, X, Y);
133: if (lmvm->shift != 0.0) VecAXPY(Y, lmvm->shift, X);
134: return 0;
135: }
137: static PetscErrorCode MatCopy_LMVM(Mat B, Mat M, MatStructure str)
138: {
139: Mat_LMVM *bctx = (Mat_LMVM *)B->data;
140: Mat_LMVM *mctx;
141: PetscInt i;
142: PetscBool allocatedM;
144: if (str == DIFFERENT_NONZERO_PATTERN) {
145: MatLMVMReset(M, PETSC_TRUE);
146: MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev);
147: } else {
148: MatLMVMIsAllocated(M, &allocatedM);
150: MatCheckSameSize(B, 1, M, 2);
151: }
153: mctx = (Mat_LMVM *)M->data;
154: if (bctx->user_pc) {
155: MatLMVMSetJ0PC(M, bctx->J0pc);
156: } else if (bctx->user_ksp) {
157: MatLMVMSetJ0KSP(M, bctx->J0ksp);
158: } else if (bctx->J0) {
159: MatLMVMSetJ0(M, bctx->J0);
160: } else if (bctx->user_scale) {
161: if (bctx->J0diag) {
162: MatLMVMSetJ0Diag(M, bctx->J0diag);
163: } else {
164: MatLMVMSetJ0Scale(M, bctx->J0scalar);
165: }
166: }
167: mctx->nupdates = bctx->nupdates;
168: mctx->nrejects = bctx->nrejects;
169: mctx->k = bctx->k;
170: for (i = 0; i <= bctx->k; ++i) {
171: VecCopy(bctx->S[i], mctx->S[i]);
172: VecCopy(bctx->Y[i], mctx->Y[i]);
173: VecCopy(bctx->Xprev, mctx->Xprev);
174: VecCopy(bctx->Fprev, mctx->Fprev);
175: }
176: if (bctx->ops->copy) (*bctx->ops->copy)(B, M, str);
177: return 0;
178: }
180: static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
181: {
182: Mat_LMVM *bctx = (Mat_LMVM *)B->data;
183: Mat_LMVM *mctx;
184: MatType lmvmType;
185: Mat A;
187: MatGetType(B, &lmvmType);
188: MatCreate(PetscObjectComm((PetscObject)B), mat);
189: 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: KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it);
200: MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev);
201: if (op == MAT_COPY_VALUES) MatCopy(B, *mat, SAME_NONZERO_PATTERN);
202: return 0;
203: }
205: static PetscErrorCode MatShift_LMVM(Mat B, PetscScalar a)
206: {
207: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
210: lmvm->shift += PetscRealPart(a);
211: return 0;
212: }
214: static PetscErrorCode MatCreateVecs_LMVM(Mat B, Vec *L, Vec *R)
215: {
216: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
219: if (L) VecDuplicate(lmvm->Xprev, L);
220: if (R) VecDuplicate(lmvm->Fprev, R);
221: return 0;
222: }
224: PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
225: {
226: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
227: PetscBool isascii;
228: MatType type;
230: PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii);
231: if (isascii) {
232: MatGetType(B, &type);
233: PetscViewerASCIIPrintf(pv, "Max. storage: %" PetscInt_FMT "\n", lmvm->m);
234: PetscViewerASCIIPrintf(pv, "Used storage: %" PetscInt_FMT "\n", lmvm->k + 1);
235: PetscViewerASCIIPrintf(pv, "Number of updates: %" PetscInt_FMT "\n", lmvm->nupdates);
236: PetscViewerASCIIPrintf(pv, "Number of rejects: %" PetscInt_FMT "\n", lmvm->nrejects);
237: PetscViewerASCIIPrintf(pv, "Number of resets: %" PetscInt_FMT "\n", lmvm->nresets);
238: if (lmvm->J0) {
239: PetscViewerASCIIPrintf(pv, "J0 Matrix:\n");
240: PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO);
241: MatView(lmvm->J0, pv);
242: PetscViewerPopFormat(pv);
243: }
244: }
245: return 0;
246: }
248: PetscErrorCode MatSetFromOptions_LMVM(Mat B, PetscOptionItems *PetscOptionsObject)
249: {
250: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
252: PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory Variable Metric matrix for approximating Jacobians");
253: PetscOptionsInt("-mat_lmvm_hist_size", "number of past updates kept in memory for the approximation", "", lmvm->m, &lmvm->m, NULL);
254: 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);
255: PetscOptionsReal("-mat_lmvm_eps", "(developer) machine zero definition", "", lmvm->eps, &lmvm->eps, NULL);
256: PetscOptionsHeadEnd();
257: KSPSetFromOptions(lmvm->J0ksp);
258: return 0;
259: }
261: PetscErrorCode MatSetUp_LMVM(Mat B)
262: {
263: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
264: PetscInt m, n, M, N;
265: PetscMPIInt size;
266: MPI_Comm comm = PetscObjectComm((PetscObject)B);
268: MatGetSize(B, &M, &N);
270: if (!lmvm->allocated) {
271: MPI_Comm_size(comm, &size);
272: if (size == 1) {
273: VecCreateSeq(comm, N, &lmvm->Xprev);
274: VecCreateSeq(comm, M, &lmvm->Fprev);
275: } else {
276: MatGetLocalSize(B, &m, &n);
277: VecCreateMPI(comm, n, N, &lmvm->Xprev);
278: VecCreateMPI(comm, m, M, &lmvm->Fprev);
279: }
280: if (lmvm->m > 0) {
281: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
282: 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: return 0;
290: }
292: PetscErrorCode MatDestroy_LMVM(Mat B)
293: {
294: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
296: if (lmvm->allocated) {
297: VecDestroyVecs(lmvm->m, &lmvm->S);
298: VecDestroyVecs(lmvm->m, &lmvm->Y);
299: VecDestroy(&lmvm->Xprev);
300: VecDestroy(&lmvm->Fprev);
301: }
302: KSPDestroy(&lmvm->J0ksp);
303: MatLMVMClearJ0(B);
304: PetscFree(B->data);
305: return 0;
306: }
308: PetscErrorCode MatCreate_LMVM(Mat B)
309: {
310: Mat_LMVM *lmvm;
312: PetscNew(&lmvm);
313: B->data = (void *)lmvm;
315: lmvm->m_old = 0;
316: lmvm->m = 5;
317: lmvm->k = -1;
318: lmvm->nupdates = 0;
319: lmvm->nrejects = 0;
320: lmvm->nresets = 0;
322: lmvm->ksp_max_it = 20;
323: lmvm->ksp_rtol = 0.0;
324: lmvm->ksp_atol = 0.0;
326: lmvm->shift = 0.0;
328: lmvm->eps = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0);
329: lmvm->allocated = PETSC_FALSE;
330: lmvm->prev_set = PETSC_FALSE;
331: lmvm->user_scale = PETSC_FALSE;
332: lmvm->user_pc = PETSC_FALSE;
333: lmvm->user_ksp = PETSC_FALSE;
334: lmvm->square = PETSC_FALSE;
336: B->ops->destroy = MatDestroy_LMVM;
337: B->ops->setfromoptions = MatSetFromOptions_LMVM;
338: B->ops->view = MatView_LMVM;
339: B->ops->setup = MatSetUp_LMVM;
340: B->ops->getvecs = MatCreateVecs_LMVM;
341: B->ops->shift = MatShift_LMVM;
342: B->ops->duplicate = MatDuplicate_LMVM;
343: B->ops->mult = MatMult_LMVM;
344: B->ops->multadd = MatMultAdd_LMVM;
345: B->ops->copy = MatCopy_LMVM;
347: lmvm->ops->update = MatUpdate_LMVM;
348: lmvm->ops->allocate = MatAllocate_LMVM;
349: lmvm->ops->reset = MatReset_LMVM;
351: KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp);
352: PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1);
353: KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_");
354: KSPSetType(lmvm->J0ksp, KSPGMRES);
355: KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it);
356: return 0;
357: }