Actual source code: lmvmimpl.c
1: #include <../src/ksp/ksp/utils/lmvm/lmvm.h>
3: /*------------------------------------------------------------*/
5: PetscErrorCode MatReset_LMVM(Mat B, PetscBool destructive)
6: {
7: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
9: lmvm->k = -1;
10: lmvm->prev_set = PETSC_FALSE;
11: lmvm->shift = 0.0;
12: if (destructive && lmvm->allocated) {
13: MatLMVMClearJ0(B);
14: B->rmap->n = B->rmap->N = B->cmap->n = B->cmap->N = 0;
15: VecDestroyVecs(lmvm->m, &lmvm->S);
16: VecDestroyVecs(lmvm->m, &lmvm->Y);
17: VecDestroy(&lmvm->Xprev);
18: 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: return 0;
28: }
30: /*------------------------------------------------------------*/
32: PetscErrorCode MatAllocate_LMVM(Mat B, Vec X, Vec F)
33: {
34: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
35: PetscBool same, allocate = PETSC_FALSE;
36: PetscInt m, n, M, N;
37: VecType type;
39: if (lmvm->allocated) {
40: VecCheckMatCompatible(B, X, 2, F, 3);
41: VecGetType(X, &type);
42: PetscObjectTypeCompare((PetscObject)lmvm->Xprev, type, &same);
43: if (!same) {
44: /* Given X vector has a different type than allocated X-type data structures.
45: We need to destroy all of this and duplicate again out of the given vector. */
46: allocate = PETSC_TRUE;
47: MatLMVMReset(B, PETSC_TRUE);
48: }
49: } else {
50: allocate = PETSC_TRUE;
51: }
52: if (allocate) {
53: VecGetLocalSize(X, &n);
54: VecGetSize(X, &N);
55: VecGetLocalSize(F, &m);
56: VecGetSize(F, &M);
57: B->rmap->n = m;
58: B->cmap->n = n;
59: B->rmap->N = M > -1 ? M : B->rmap->N;
60: B->cmap->N = N > -1 ? N : B->cmap->N;
61: VecDuplicate(X, &lmvm->Xprev);
62: VecDuplicate(F, &lmvm->Fprev);
63: if (lmvm->m > 0) {
64: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
65: VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);
66: }
67: lmvm->m_old = lmvm->m;
68: lmvm->allocated = PETSC_TRUE;
69: B->preallocated = PETSC_TRUE;
70: B->assembled = PETSC_TRUE;
71: }
72: return 0;
73: }
75: /*------------------------------------------------------------*/
77: PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
78: {
79: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
80: PetscInt i;
81: Vec Stmp, Ytmp;
83: if (lmvm->k == lmvm->m-1) {
84: /* We hit the memory limit, so shift all the vectors back one spot
85: and shift the oldest to the front to receive the latest update. */
86: Stmp = lmvm->S[0];
87: Ytmp = lmvm->Y[0];
88: for (i = 0; i < lmvm->k; ++i) {
89: lmvm->S[i] = lmvm->S[i+1];
90: lmvm->Y[i] = lmvm->Y[i+1];
91: }
92: lmvm->S[lmvm->k] = Stmp;
93: lmvm->Y[lmvm->k] = Ytmp;
94: } else {
95: ++lmvm->k;
96: }
97: /* Put the precomputed update into the last vector */
98: VecCopy(S, lmvm->S[lmvm->k]);
99: VecCopy(Y, lmvm->Y[lmvm->k]);
100: ++lmvm->nupdates;
101: return 0;
102: }
104: /*------------------------------------------------------------*/
106: PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
107: {
108: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
110: if (!lmvm->m) return 0;
111: if (lmvm->prev_set) {
112: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
113: VecAXPBY(lmvm->Xprev, 1.0, -1.0, X);
114: VecAXPBY(lmvm->Fprev, 1.0, -1.0, F);
115: /* Update S and Y */
116: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
117: }
119: /* Save the solution and function to be used in the next update */
120: VecCopy(X, lmvm->Xprev);
121: VecCopy(F, lmvm->Fprev);
122: lmvm->prev_set = PETSC_TRUE;
123: return 0;
124: }
126: /*------------------------------------------------------------*/
128: static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
129: {
130: MatMult(B, X, Z);
131: VecAXPY(Z, 1.0, Y);
132: return 0;
133: }
135: /*------------------------------------------------------------*/
137: static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
138: {
139: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
141: VecCheckSameSize(X, 2, Y, 3);
142: VecCheckMatCompatible(B, X, 2, Y, 3);
144: (*lmvm->ops->mult)(B, X, Y);
145: if (lmvm->shift != 0.0) {
146: VecAXPY(Y, lmvm->shift, X);
147: }
148: return 0;
149: }
151: /*------------------------------------------------------------*/
153: static PetscErrorCode MatCopy_LMVM(Mat B, Mat M, MatStructure str)
154: {
155: Mat_LMVM *bctx = (Mat_LMVM*)B->data;
156: Mat_LMVM *mctx;
157: PetscInt i;
158: PetscBool allocatedM;
160: if (str == DIFFERENT_NONZERO_PATTERN) {
161: MatLMVMReset(M, PETSC_TRUE);
162: MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev);
163: } else {
164: MatLMVMIsAllocated(M, &allocatedM);
166: MatCheckSameSize(B, 1, M, 2);
167: }
169: mctx = (Mat_LMVM*)M->data;
170: if (bctx->user_pc) {
171: MatLMVMSetJ0PC(M, bctx->J0pc);
172: } else if (bctx->user_ksp) {
173: MatLMVMSetJ0KSP(M, bctx->J0ksp);
174: } else if (bctx->J0) {
175: MatLMVMSetJ0(M, bctx->J0);
176: } else if (bctx->user_scale) {
177: if (bctx->J0diag) {
178: MatLMVMSetJ0Diag(M, bctx->J0diag);
179: } else {
180: MatLMVMSetJ0Scale(M, bctx->J0scalar);
181: }
182: }
183: mctx->nupdates = bctx->nupdates;
184: mctx->nrejects = bctx->nrejects;
185: mctx->k = bctx->k;
186: for (i=0; i<=bctx->k; ++i) {
187: VecCopy(bctx->S[i], mctx->S[i]);
188: VecCopy(bctx->Y[i], mctx->Y[i]);
189: VecCopy(bctx->Xprev, mctx->Xprev);
190: VecCopy(bctx->Fprev, mctx->Fprev);
191: }
192: if (bctx->ops->copy) {
193: (*bctx->ops->copy)(B, M, str);
194: }
195: return 0;
196: }
198: /*------------------------------------------------------------*/
200: static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
201: {
202: Mat_LMVM *bctx = (Mat_LMVM*)B->data;
203: Mat_LMVM *mctx;
204: MatType lmvmType;
205: Mat A;
207: MatGetType(B, &lmvmType);
208: MatCreate(PetscObjectComm((PetscObject)B), mat);
209: MatSetType(*mat, lmvmType);
211: A = *mat;
212: mctx = (Mat_LMVM*)A->data;
213: mctx->m = bctx->m;
214: mctx->ksp_max_it = bctx->ksp_max_it;
215: mctx->ksp_rtol = bctx->ksp_rtol;
216: mctx->ksp_atol = bctx->ksp_atol;
217: mctx->shift = bctx->shift;
218: KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it);
220: MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev);
221: if (op == MAT_COPY_VALUES) {
222: MatCopy(B, *mat, SAME_NONZERO_PATTERN);
223: }
224: return 0;
225: }
227: /*------------------------------------------------------------*/
229: static PetscErrorCode MatShift_LMVM(Mat B, PetscScalar a)
230: {
231: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
234: lmvm->shift += PetscRealPart(a);
235: return 0;
236: }
238: /*------------------------------------------------------------*/
240: static PetscErrorCode MatGetVecs_LMVM(Mat B, Vec *L, Vec *R)
241: {
242: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
245: VecDuplicate(lmvm->Xprev, L);
246: VecDuplicate(lmvm->Fprev, R);
247: return 0;
248: }
250: /*------------------------------------------------------------*/
252: PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
253: {
254: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
255: PetscBool isascii;
256: MatType type;
258: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
259: if (isascii) {
260: MatGetType(B, &type);
261: PetscViewerASCIIPrintf(pv,"Max. storage: %D\n",lmvm->m);
262: PetscViewerASCIIPrintf(pv,"Used storage: %D\n",lmvm->k+1);
263: PetscViewerASCIIPrintf(pv,"Number of updates: %D\n",lmvm->nupdates);
264: PetscViewerASCIIPrintf(pv,"Number of rejects: %D\n",lmvm->nrejects);
265: PetscViewerASCIIPrintf(pv,"Number of resets: %D\n",lmvm->nresets);
266: if (lmvm->J0) {
267: PetscViewerASCIIPrintf(pv,"J0 Matrix:\n");
268: PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO);
269: MatView(lmvm->J0, pv);
270: PetscViewerPopFormat(pv);
271: }
272: }
273: return 0;
274: }
276: /*------------------------------------------------------------*/
278: PetscErrorCode MatSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject, Mat B)
279: {
280: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
282: PetscOptionsHead(PetscOptionsObject,"Limited-memory Variable Metric matrix for approximating Jacobians");
283: PetscOptionsInt("-mat_lmvm_hist_size","number of past updates kept in memory for the approximation","",lmvm->m,&lmvm->m,NULL);
284: 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);
285: PetscOptionsReal("-mat_lmvm_eps","(developer) machine zero definition","",lmvm->eps,&lmvm->eps,NULL);
286: PetscOptionsTail();
287: KSPSetFromOptions(lmvm->J0ksp);
288: return 0;
289: }
291: /*------------------------------------------------------------*/
293: PetscErrorCode MatSetUp_LMVM(Mat B)
294: {
295: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
296: PetscInt m, n, M, N;
297: PetscMPIInt size;
298: MPI_Comm comm = PetscObjectComm((PetscObject)B);
300: MatGetSize(B, &M, &N);
302: if (!lmvm->allocated) {
303: MPI_Comm_size(comm, &size);
304: if (size == 1) {
305: VecCreateSeq(comm, N, &lmvm->Xprev);
306: VecCreateSeq(comm, M, &lmvm->Fprev);
307: } else {
308: MatGetLocalSize(B, &m, &n);
309: VecCreateMPI(comm, n, N, &lmvm->Xprev);
310: VecCreateMPI(comm, m, M, &lmvm->Fprev);
311: }
312: if (lmvm->m > 0) {
313: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
314: VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);
315: }
316: lmvm->m_old = lmvm->m;
317: lmvm->allocated = PETSC_TRUE;
318: B->preallocated = PETSC_TRUE;
319: B->assembled = PETSC_TRUE;
320: }
321: return 0;
322: }
324: /*------------------------------------------------------------*/
326: PetscErrorCode MatDestroy_LMVM(Mat B)
327: {
328: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
330: if (lmvm->allocated) {
331: VecDestroyVecs(lmvm->m, &lmvm->S);
332: VecDestroyVecs(lmvm->m, &lmvm->Y);
333: VecDestroy(&lmvm->Xprev);
334: VecDestroy(&lmvm->Fprev);
335: }
336: KSPDestroy(&lmvm->J0ksp);
337: MatLMVMClearJ0(B);
338: PetscFree(B->data);
339: return 0;
340: }
342: /*------------------------------------------------------------*/
344: PetscErrorCode MatCreate_LMVM(Mat B)
345: {
346: Mat_LMVM *lmvm;
348: PetscNewLog(B, &lmvm);
349: B->data = (void*)lmvm;
351: lmvm->m_old = 0;
352: lmvm->m = 5;
353: lmvm->k = -1;
354: lmvm->nupdates = 0;
355: lmvm->nrejects = 0;
356: lmvm->nresets = 0;
358: lmvm->ksp_max_it = 20;
359: lmvm->ksp_rtol = 0.0;
360: lmvm->ksp_atol = 0.0;
362: lmvm->shift = 0.0;
364: lmvm->eps = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0);
365: lmvm->allocated = PETSC_FALSE;
366: lmvm->prev_set = PETSC_FALSE;
367: lmvm->user_scale = PETSC_FALSE;
368: lmvm->user_pc = PETSC_FALSE;
369: lmvm->user_ksp = PETSC_FALSE;
370: lmvm->square = PETSC_FALSE;
372: B->ops->destroy = MatDestroy_LMVM;
373: B->ops->setfromoptions = MatSetFromOptions_LMVM;
374: B->ops->view = MatView_LMVM;
375: B->ops->setup = MatSetUp_LMVM;
376: B->ops->getvecs = MatGetVecs_LMVM;
377: B->ops->shift = MatShift_LMVM;
378: B->ops->duplicate = MatDuplicate_LMVM;
379: B->ops->mult = MatMult_LMVM;
380: B->ops->multadd = MatMultAdd_LMVM;
381: B->ops->copy = MatCopy_LMVM;
383: lmvm->ops->update = MatUpdate_LMVM;
384: lmvm->ops->allocate = MatAllocate_LMVM;
385: lmvm->ops->reset = MatReset_LMVM;
387: KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp);
388: PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1);
389: KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_");
390: KSPSetType(lmvm->J0ksp, KSPGMRES);
391: KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it);
392: return 0;
393: }