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;
8: PetscErrorCode ierr;
11: lmvm->k = -1;
12: lmvm->prev_set = PETSC_FALSE;
13: lmvm->shift = 0.0;
14: if (destructive && lmvm->allocated) {
15: MatLMVMClearJ0(B);
16: B->rmap->n = B->rmap->N = B->cmap->n = B->cmap->N = 0;
17: VecDestroyVecs(lmvm->m, &lmvm->S);
18: VecDestroyVecs(lmvm->m, &lmvm->Y);
19: VecDestroy(&lmvm->Xprev);
20: VecDestroy(&lmvm->Fprev);
21: lmvm->nupdates = 0;
22: lmvm->nrejects = 0;
23: lmvm->m_old = 0;
24: lmvm->allocated = PETSC_FALSE;
25: B->preallocated = PETSC_FALSE;
26: B->assembled = PETSC_FALSE;
27: }
28: ++lmvm->nresets;
29: return(0);
30: }
32: /*------------------------------------------------------------*/
34: PetscErrorCode MatAllocate_LMVM(Mat B, Vec X, Vec F)
35: {
36: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
37: PetscErrorCode ierr;
38: PetscBool same, allocate = PETSC_FALSE;
39: PetscInt m, n, M, N;
40: VecType type;
43: if (lmvm->allocated) {
44: VecCheckMatCompatible(B, X, 2, F, 3);
45: VecGetType(X, &type);
46: PetscObjectTypeCompare((PetscObject)lmvm->Xprev, type, &same);
47: if (!same) {
48: /* Given X vector has a different type than allocated X-type data structures.
49: We need to destroy all of this and duplicate again out of the given vector. */
50: allocate = PETSC_TRUE;
51: MatLMVMReset(B, PETSC_TRUE);
52: }
53: } else {
54: allocate = PETSC_TRUE;
55: }
56: if (allocate) {
57: VecGetLocalSize(X, &n);
58: VecGetSize(X, &N);
59: VecGetLocalSize(F, &m);
60: VecGetSize(F, &M);
61: B->rmap->n = m;
62: B->cmap->n = n;
63: B->rmap->N = M > -1 ? M : B->rmap->N;
64: B->cmap->N = N > -1 ? N : B->cmap->N;
65: VecDuplicate(X, &lmvm->Xprev);
66: VecDuplicate(F, &lmvm->Fprev);
67: if (lmvm->m > 0) {
68: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
69: VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);
70: }
71: lmvm->m_old = lmvm->m;
72: lmvm->allocated = PETSC_TRUE;
73: B->preallocated = PETSC_TRUE;
74: B->assembled = PETSC_TRUE;
75: }
76: return(0);
77: }
79: /*------------------------------------------------------------*/
81: PetscErrorCode MatUpdateKernel_LMVM(Mat B, Vec S, Vec Y)
82: {
83: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
84: PetscErrorCode ierr;
85: PetscInt i;
86: Vec Stmp, Ytmp;
89: if (lmvm->k == lmvm->m-1) {
90: /* We hit the memory limit, so shift all the vectors back one spot
91: and shift the oldest to the front to receive the latest update. */
92: Stmp = lmvm->S[0];
93: Ytmp = lmvm->Y[0];
94: for (i = 0; i < lmvm->k; ++i) {
95: lmvm->S[i] = lmvm->S[i+1];
96: lmvm->Y[i] = lmvm->Y[i+1];
97: }
98: lmvm->S[lmvm->k] = Stmp;
99: lmvm->Y[lmvm->k] = Ytmp;
100: } else {
101: ++lmvm->k;
102: }
103: /* Put the precomputed update into the last vector */
104: VecCopy(S, lmvm->S[lmvm->k]);
105: VecCopy(Y, lmvm->Y[lmvm->k]);
106: ++lmvm->nupdates;
107: return(0);
108: }
110: /*------------------------------------------------------------*/
112: PetscErrorCode MatUpdate_LMVM(Mat B, Vec X, Vec F)
113: {
114: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
115: PetscErrorCode ierr;
118: if (!lmvm->m) return(0);
119: if (lmvm->prev_set) {
120: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
121: VecAXPBY(lmvm->Xprev, 1.0, -1.0, X);
122: VecAXPBY(lmvm->Fprev, 1.0, -1.0, F);
123: /* Update S and Y */
124: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
125: }
127: /* Save the solution and function to be used in the next update */
128: VecCopy(X, lmvm->Xprev);
129: VecCopy(F, lmvm->Fprev);
130: lmvm->prev_set = PETSC_TRUE;
131: return(0);
132: }
134: /*------------------------------------------------------------*/
136: static PetscErrorCode MatMultAdd_LMVM(Mat B, Vec X, Vec Y, Vec Z)
137: {
138: PetscErrorCode ierr;
141: MatMult(B, X, Z);
142: VecAXPY(Z, 1.0, Y);
143: return(0);
144: }
146: /*------------------------------------------------------------*/
148: static PetscErrorCode MatMult_LMVM(Mat B, Vec X, Vec Y)
149: {
150: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
151: PetscErrorCode ierr;
154: VecCheckSameSize(X, 2, Y, 3);
155: VecCheckMatCompatible(B, X, 2, Y, 3);
156: if (!lmvm->allocated) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
157: (*lmvm->ops->mult)(B, X, Y);
158: if (lmvm->shift != 0.0) {
159: VecAXPY(Y, lmvm->shift, X);
160: }
161: return(0);
162: }
164: /*------------------------------------------------------------*/
166: static PetscErrorCode MatCopy_LMVM(Mat B, Mat M, MatStructure str)
167: {
168: Mat_LMVM *bctx = (Mat_LMVM*)B->data;
169: Mat_LMVM *mctx;
170: PetscErrorCode ierr;
171: PetscInt i;
172: PetscBool allocatedM;
175: if (str == DIFFERENT_NONZERO_PATTERN) {
176: MatLMVMReset(M, PETSC_TRUE);
177: MatLMVMAllocate(M, bctx->Xprev, bctx->Fprev);
178: } else {
179: MatLMVMIsAllocated(M, &allocatedM);
180: if (!allocatedM) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_WRONGSTATE, "Target matrix must be allocated first");
181: MatCheckSameSize(B, 1, M, 2);
182: }
184: mctx = (Mat_LMVM*)M->data;
185: if (bctx->user_pc) {
186: MatLMVMSetJ0PC(M, bctx->J0pc);
187: } else if (bctx->user_ksp) {
188: MatLMVMSetJ0KSP(M, bctx->J0ksp);
189: } else if (bctx->J0) {
190: MatLMVMSetJ0(M, bctx->J0);
191: } else if (bctx->user_scale) {
192: if (bctx->J0diag) {
193: MatLMVMSetJ0Diag(M, bctx->J0diag);
194: } else {
195: MatLMVMSetJ0Scale(M, bctx->J0scalar);
196: }
197: }
198: mctx->nupdates = bctx->nupdates;
199: mctx->nrejects = bctx->nrejects;
200: mctx->k = bctx->k;
201: for (i=0; i<=bctx->k; ++i) {
202: VecCopy(bctx->S[i], mctx->S[i]);
203: VecCopy(bctx->Y[i], mctx->Y[i]);
204: VecCopy(bctx->Xprev, mctx->Xprev);
205: VecCopy(bctx->Fprev, mctx->Fprev);
206: }
207: if (bctx->ops->copy) {
208: (*bctx->ops->copy)(B, M, str);
209: }
210: return(0);
211: }
213: /*------------------------------------------------------------*/
215: static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
216: {
217: Mat_LMVM *bctx = (Mat_LMVM*)B->data;
218: Mat_LMVM *mctx;
219: PetscErrorCode ierr;
220: MatType lmvmType;
221: Mat A;
224: MatGetType(B, &lmvmType);
225: MatCreate(PetscObjectComm((PetscObject)B), mat);
226: MatSetType(*mat, lmvmType);
228: A = *mat;
229: mctx = (Mat_LMVM*)A->data;
230: mctx->m = bctx->m;
231: mctx->ksp_max_it = bctx->ksp_max_it;
232: mctx->ksp_rtol = bctx->ksp_rtol;
233: mctx->ksp_atol = bctx->ksp_atol;
234: mctx->shift = bctx->shift;
235: KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it);
237: MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev);
238: if (op == MAT_COPY_VALUES) {
239: MatCopy(B, *mat, SAME_NONZERO_PATTERN);
240: }
241: return(0);
242: }
244: /*------------------------------------------------------------*/
246: static PetscErrorCode MatShift_LMVM(Mat B, PetscScalar a)
247: {
248: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
251: if (!lmvm->allocated) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
252: lmvm->shift += PetscRealPart(a);
253: return(0);
254: }
256: /*------------------------------------------------------------*/
258: static PetscErrorCode MatGetVecs_LMVM(Mat B, Vec *L, Vec *R)
259: {
260: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
261: PetscErrorCode ierr;
264: if (!lmvm->allocated) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ORDER, "LMVM matrix must be allocated first");
265: VecDuplicate(lmvm->Xprev, L);
266: VecDuplicate(lmvm->Fprev, R);
267: return(0);
268: }
270: /*------------------------------------------------------------*/
272: PetscErrorCode MatView_LMVM(Mat B, PetscViewer pv)
273: {
274: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
275: PetscErrorCode ierr;
276: PetscBool isascii;
277: MatType type;
280: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
281: if (isascii) {
282: MatGetType(B, &type);
283: PetscViewerASCIIPrintf(pv,"Max. storage: %D\n",lmvm->m);
284: PetscViewerASCIIPrintf(pv,"Used storage: %D\n",lmvm->k+1);
285: PetscViewerASCIIPrintf(pv,"Number of updates: %D\n",lmvm->nupdates);
286: PetscViewerASCIIPrintf(pv,"Number of rejects: %D\n",lmvm->nrejects);
287: PetscViewerASCIIPrintf(pv,"Number of resets: %D\n",lmvm->nresets);
288: if (lmvm->J0) {
289: PetscViewerASCIIPrintf(pv,"J0 Matrix:\n");
290: PetscViewerPushFormat(pv, PETSC_VIEWER_ASCII_INFO);
291: MatView(lmvm->J0, pv);
292: PetscViewerPopFormat(pv);
293: }
294: }
295: return(0);
296: }
298: /*------------------------------------------------------------*/
300: PetscErrorCode MatSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject, Mat B)
301: {
302: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
303: PetscErrorCode ierr;
306: PetscOptionsHead(PetscOptionsObject,"Limited-memory Variable Metric matrix for approximating Jacobians");
307: PetscOptionsInt("-mat_lmvm_hist_size","number of past updates kept in memory for the approximation","",lmvm->m,&lmvm->m,NULL);
308: 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);
309: PetscOptionsReal("-mat_lmvm_eps","(developer) machine zero definition","",lmvm->eps,&lmvm->eps,NULL);
310: PetscOptionsTail();
311: KSPSetFromOptions(lmvm->J0ksp);
312: return(0);
313: }
315: /*------------------------------------------------------------*/
317: PetscErrorCode MatSetUp_LMVM(Mat B)
318: {
319: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
320: PetscErrorCode ierr;
321: PetscInt m, n, M, N;
322: PetscMPIInt size;
323: MPI_Comm comm = PetscObjectComm((PetscObject)B);
326: MatGetSize(B, &M, &N);
327: if (M == 0 && N == 0) SETERRQ(comm, PETSC_ERR_ORDER, "MatSetSizes() must be called before MatSetUp()");
328: if (!lmvm->allocated) {
329: MPI_Comm_size(comm, &size);
330: if (size == 1) {
331: VecCreateSeq(comm, N, &lmvm->Xprev);
332: VecCreateSeq(comm, M, &lmvm->Fprev);
333: } else {
334: MatGetLocalSize(B, &m, &n);
335: VecCreateMPI(comm, n, N, &lmvm->Xprev);
336: VecCreateMPI(comm, m, M, &lmvm->Fprev);
337: }
338: if (lmvm->m > 0) {
339: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lmvm->S);
340: VecDuplicateVecs(lmvm->Fprev, lmvm->m, &lmvm->Y);
341: }
342: lmvm->m_old = lmvm->m;
343: lmvm->allocated = PETSC_TRUE;
344: B->preallocated = PETSC_TRUE;
345: B->assembled = PETSC_TRUE;
346: }
347: return(0);
348: }
350: /*------------------------------------------------------------*/
352: PetscErrorCode MatDestroy_LMVM(Mat B)
353: {
354: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
355: PetscErrorCode ierr;
358: if (lmvm->allocated) {
359: VecDestroyVecs(lmvm->m, &lmvm->S);
360: VecDestroyVecs(lmvm->m, &lmvm->Y);
361: VecDestroy(&lmvm->Xprev);
362: VecDestroy(&lmvm->Fprev);
363: }
364: KSPDestroy(&lmvm->J0ksp);
365: MatLMVMClearJ0(B);
366: PetscFree(B->data);
367: return(0);
368: }
370: /*------------------------------------------------------------*/
372: PetscErrorCode MatCreate_LMVM(Mat B)
373: {
374: Mat_LMVM *lmvm;
375: PetscErrorCode ierr;
378: PetscNewLog(B, &lmvm);
379: B->data = (void*)lmvm;
381: lmvm->m_old = 0;
382: lmvm->m = 5;
383: lmvm->k = -1;
384: lmvm->nupdates = 0;
385: lmvm->nrejects = 0;
386: lmvm->nresets = 0;
388: lmvm->ksp_max_it = 20;
389: lmvm->ksp_rtol = 0.0;
390: lmvm->ksp_atol = 0.0;
392: lmvm->shift = 0.0;
394: lmvm->eps = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0);
395: lmvm->allocated = PETSC_FALSE;
396: lmvm->prev_set = PETSC_FALSE;
397: lmvm->user_scale = PETSC_FALSE;
398: lmvm->user_pc = PETSC_FALSE;
399: lmvm->user_ksp = PETSC_FALSE;
400: lmvm->square = PETSC_FALSE;
402: B->ops->destroy = MatDestroy_LMVM;
403: B->ops->setfromoptions = MatSetFromOptions_LMVM;
404: B->ops->view = MatView_LMVM;
405: B->ops->setup = MatSetUp_LMVM;
406: B->ops->getvecs = MatGetVecs_LMVM;
407: B->ops->shift = MatShift_LMVM;
408: B->ops->duplicate = MatDuplicate_LMVM;
409: B->ops->mult = MatMult_LMVM;
410: B->ops->multadd = MatMultAdd_LMVM;
411: B->ops->copy = MatCopy_LMVM;
413: lmvm->ops->update = MatUpdate_LMVM;
414: lmvm->ops->allocate = MatAllocate_LMVM;
415: lmvm->ops->reset = MatReset_LMVM;
417: KSPCreate(PetscObjectComm((PetscObject)B), &lmvm->J0ksp);
418: PetscObjectIncrementTabLevel((PetscObject)lmvm->J0ksp, (PetscObject)B, 1);
419: KSPSetOptionsPrefix(lmvm->J0ksp, "mat_lmvm_");
420: KSPSetType(lmvm->J0ksp, KSPGMRES);
421: KSPSetTolerances(lmvm->J0ksp, lmvm->ksp_rtol, lmvm->ksp_atol, PETSC_DEFAULT, lmvm->ksp_max_it);
422: return(0);
423: }