Actual source code: pod.c
petsc-3.9.4 2018-09-11
1: #include <petsc/private/kspimpl.h>
2: #include <petsc/private/matimpl.h>
3: #include <petscblaslapack.h>
4: static PetscBool cited = PETSC_FALSE;
5: static const char citation[] =
6: "@phdthesis{zampini2010non,\n"
7: " title={Non-overlapping Domain Decomposition Methods for Cardiac Reaction-Diffusion Models and Applications},\n"
8: " author={Zampini, S},\n"
9: " year={2010},\n"
10: " school={PhD thesis, Universita degli Studi di Milano}\n"
11: "}\n";
13: typedef struct {
14: PetscInt maxn; /* maximum number of snapshots */
15: PetscInt n; /* number of active snapshots */
16: PetscInt curr; /* current tip of snapshots set */
17: Vec *xsnap; /* snapshots */
18: Vec *bsnap; /* rhs snapshots */
19: PetscScalar *dots_iallreduce;
20: MPI_Request req_iallreduce;
21: PetscInt ndots_iallreduce; /* if we have iallreduce we can hide the VecMDot communications */
22: PetscReal tol; /* relative tolerance to retain eigenvalues */
23: PetscBool Aspd; /* if true, uses the SPD operator as inner product */
24: PetscScalar *corr; /* correlation matrix */
25: PetscReal *eigs; /* eigenvalues */
26: PetscScalar *eigv; /* eigenvectors */
27: PetscBLASInt nen; /* dimension of lower dimensional system */
28: PetscInt st; /* first eigenvector of correlation matrix to be retained */
29: PetscBLASInt *iwork; /* integer work vector */
30: PetscScalar *yhay; /* Y^H * A * Y */
31: PetscScalar *low; /* lower dimensional linear system */
32: #if defined(PETSC_USE_COMPLEX)
33: PetscReal *rwork;
34: #endif
35: PetscBLASInt lwork;
36: PetscScalar *swork;
37: PetscBool monitor;
38: } KSPGuessPOD;
40: static PetscErrorCode KSPGuessReset_POD(KSPGuess guess)
41: {
42: KSPGuessPOD *pod = (KSPGuessPOD*)guess->data;
44: PetscLayout Alay = NULL,vlay = NULL;
45: PetscBool cong;
48: pod->nen = 0;
49: pod->n = 0;
50: pod->curr = 0;
51: /* need to wait for completion of outstanding requests */
52: if (pod->ndots_iallreduce) {
53: MPI_Wait(&pod->req_iallreduce,MPI_STATUS_IGNORE);
54: }
55: pod->ndots_iallreduce = 0;
56: /* destroy vectors if the size of the linear system has changed */
57: if (guess->A) {
58: MatGetLayouts(guess->A,&Alay,NULL);
59: }
60: if (pod->xsnap) {
61: VecGetLayout(pod->xsnap[0],&vlay);
62: }
63: cong = PETSC_FALSE;
64: if (vlay && Alay) {
65: PetscLayoutCompare(Alay,vlay,&cong);
66: }
67: if (!cong) {
68: VecDestroyVecs(pod->maxn,&pod->xsnap);
69: VecDestroyVecs(pod->maxn,&pod->bsnap);
70: }
71: return(0);
72: }
74: static PetscErrorCode KSPGuessSetUp_POD(KSPGuess guess)
75: {
76: KSPGuessPOD *pod = (KSPGuessPOD*)guess->data;
80: if (!pod->corr) {
81: PetscScalar sdummy;
82: PetscReal rdummy = 0;
83: PetscBLASInt bN,lierr,idummy;
85: PetscCalloc6(pod->maxn*pod->maxn,&pod->corr,pod->maxn,&pod->eigs,pod->maxn*pod->maxn,&pod->eigv,
86: 6*pod->maxn,&pod->iwork,pod->maxn*pod->maxn,&pod->yhay,pod->maxn*pod->maxn,&pod->low);
87: #if defined(PETSC_USE_COMPLEX)
88: PetscMalloc1(7*pod->maxn,&pod->rwork);
89: #endif
90: #if defined(PETSC_HAVE_MPI_IALLREDUCE)
91: PetscMalloc1(3*pod->maxn,&pod->dots_iallreduce);
92: #endif
93: pod->lwork = -1;
94: PetscBLASIntCast(pod->maxn,&bN);
95: #if !defined(PETSC_USE_COMPLEX)
96: PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","L",&bN,pod->corr,&bN,&rdummy,&rdummy,&idummy,&idummy,
97: &rdummy,&idummy,pod->eigs,pod->eigv,&bN,&sdummy,&pod->lwork,pod->iwork,pod->iwork+5*bN,&lierr));
98: #else
99: PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","L",&bN,pod->corr,&bN,&rdummy,&rdummy,&idummy,&idummy,
100: &rdummy,&idummy,pod->eigs,pod->eigv,&bN,&sdummy,&pod->lwork,pod->rwork,pod->iwork,pod->iwork+5*bN,&lierr));
101: #endif
102: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to SYEV Lapack routine %d",(int)lierr);
103: PetscBLASIntCast((PetscInt)PetscRealPart(sdummy),&pod->lwork);
104: PetscMalloc1(pod->lwork+PetscMax(bN*bN,6*bN),&pod->swork);
105: }
106: /* work vectors are sequential, we explicitly use MPI_Allreduce */
107: if (!pod->xsnap) {
108: VecType type;
109: Vec *v,vseq;
110: PetscInt n;
112: KSPCreateVecs(guess->ksp,1,&v,0,NULL);
113: VecCreate(PETSC_COMM_SELF,&vseq);
114: VecGetLocalSize(v[0],&n);
115: VecSetSizes(vseq,n,n);
116: VecGetType(v[0],&type);
117: VecSetType(vseq,type);
118: VecDestroyVecs(1,&v);
119: VecDuplicateVecs(vseq,pod->maxn,&pod->xsnap);
120: VecDestroy(&vseq);
121: PetscLogObjectParents(guess,pod->maxn,pod->xsnap);
122: }
123: if (!pod->bsnap) {
124: VecDuplicateVecs(pod->xsnap[0],pod->maxn,&pod->bsnap);
125: PetscLogObjectParents(guess,pod->maxn,pod->bsnap);
126: }
127: return(0);
128: }
130: static PetscErrorCode KSPGuessDestroy_POD(KSPGuess guess)
131: {
132: KSPGuessPOD *pod = (KSPGuessPOD*)guess->data;
133: PetscErrorCode ierr;
136: PetscFree6(pod->corr,pod->eigs,pod->eigv,pod->iwork,
137: pod->yhay,pod->low);
138: #if defined(PETSC_USE_COMPLEX)
139: PetscFree(pod->rwork);
140: #endif
141: /* need to wait for completion before destroying dots_iallreduce */
142: if (pod->ndots_iallreduce) {
143: MPI_Wait(&pod->req_iallreduce,MPI_STATUS_IGNORE);
144: }
145: PetscFree(pod->dots_iallreduce);
146: PetscFree(pod->swork);
147: VecDestroyVecs(pod->maxn,&pod->bsnap);
148: VecDestroyVecs(pod->maxn,&pod->xsnap);
149: PetscFree(pod);
150: return(0);
151: }
153: static PetscErrorCode KSPGuessUpdate_POD(KSPGuess,Vec,Vec);
155: static PetscErrorCode KSPGuessFormGuess_POD(KSPGuess guess,Vec b,Vec x)
156: {
157: KSPGuessPOD *pod = (KSPGuessPOD*)guess->data;
159: PetscScalar one = 1, zero = 0, *array;
160: PetscBLASInt bN,ione = 1,bNen,lierr;
161: PetscInt i;
164: PetscCitationsRegister(citation,&cited);
165: if (pod->ndots_iallreduce) { /* complete communication and project the linear system */
166: KSPGuessUpdate_POD(guess,NULL,NULL);
167: }
168: if (!pod->nen) return(0);
169: /* b_low = S * V^T * X^T * b */
170: VecGetArrayRead(b,(const PetscScalar**)&array);
171: VecPlaceArray(pod->bsnap[pod->curr],array);
172: VecRestoreArrayRead(b,(const PetscScalar**)&array);
173: VecMDot(pod->bsnap[pod->curr],pod->n,pod->xsnap,pod->swork);
174: VecResetArray(pod->bsnap[pod->curr]);
175: MPIU_Allreduce(pod->swork,pod->swork + pod->n,pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
176: PetscBLASIntCast(pod->n,&bN);
177: PetscBLASIntCast(pod->nen,&bNen);
178: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&bN,&bNen,&one,pod->eigv+pod->st*pod->n,&bN,pod->swork+pod->n,&ione,&zero,pod->swork,&ione));
179: if (pod->monitor) {
180: PetscPrintf(PetscObjectComm((PetscObject)guess)," KSPGuessPOD alphas = ");
181: for (i=0; i<pod->nen; i++) {
182: #if defined(PETSC_USE_COMPLEX)
183: PetscPrintf(PetscObjectComm((PetscObject)guess),"%g + %g i",(double)PetscRealPart(pod->swork[i]),(double)PetscImaginaryPart(pod->swork[i]));
184: #else
185: PetscPrintf(PetscObjectComm((PetscObject)guess),"%g ",(double)pod->swork[i]);
186: #endif
187: }
188: PetscPrintf(PetscObjectComm((PetscObject)guess),"\n");
189: }
190: /* A_low x_low = b_low */
191: if (!pod->Aspd) { /* A is spd -> LOW = Identity */
192: KSP pksp = guess->ksp;
193: PetscBool tsolve,symm;
195: if (pod->monitor) {
196: PetscMPIInt rank;
197: Mat L;
199: MPI_Comm_rank(PetscObjectComm((PetscObject)guess),&rank);
200: MatCreateSeqDense(PETSC_COMM_SELF,pod->nen,pod->nen,pod->low,&L);
201: if (!rank) {
202: MatView(L,NULL);
203: }
204: MatDestroy(&L);
205: }
206: MatGetOption(guess->A,MAT_SYMMETRIC,&symm);
207: tsolve = symm ? PETSC_FALSE : pksp->transpose_solve;
208: PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&bNen,&bNen,pod->low,&bNen,pod->iwork,&lierr));
209: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)lierr);
210: PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_(tsolve ? "T" : "N",&bNen,&ione,pod->low,&bNen,pod->iwork,pod->swork,&bNen,&lierr));
211: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRS Lapack routine %d",(int)lierr);
212: }
213: /* x = X * V * S * x_low */
214: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&bN,&bNen,&one,pod->eigv+pod->st*pod->n,&bN,pod->swork,&ione,&zero,pod->swork+pod->n,&ione));
215: if (pod->monitor) {
216: PetscPrintf(PetscObjectComm((PetscObject)guess)," KSPGuessPOD sol = ");
217: for (i=0; i<pod->nen; i++) {
218: #if defined(PETSC_USE_COMPLEX)
219: PetscPrintf(PetscObjectComm((PetscObject)guess),"%g + %g i",(double)PetscRealPart(pod->swork[i+pod->n]),(double)PetscImaginaryPart(pod->swork[i+pod->n]));
220: #else
221: PetscPrintf(PetscObjectComm((PetscObject)guess),"%g ",(double)pod->swork[i+pod->n]);
222: #endif
223: }
224: PetscPrintf(PetscObjectComm((PetscObject)guess),"\n");
225: }
226: VecGetArray(x,&array);
227: VecPlaceArray(pod->bsnap[pod->curr],array);
228: VecRestoreArray(x,&array);
229: VecSet(pod->bsnap[pod->curr],0);
230: VecMAXPY(pod->bsnap[pod->curr],pod->n,pod->swork+pod->n,pod->xsnap);
231: VecResetArray(pod->bsnap[pod->curr]);
232: PetscObjectStateIncrease((PetscObject)x);
233: return(0);
234: }
236: static PetscErrorCode KSPGuessUpdate_POD(KSPGuess guess, Vec b, Vec x)
237: {
238: KSPGuessPOD *pod = (KSPGuessPOD*)guess->data;
239: PetscScalar one = 1, zero = 0,*array;
240: PetscReal toten, parten, reps = 0; /* dlamch? */
241: PetscBLASInt bN,lierr,idummy;
242: PetscInt i;
246: if (pod->ndots_iallreduce) goto complete_request;
247: pod->n = pod->n < pod->maxn ? pod->n+1 : pod->maxn;
248: VecCopy(x,pod->xsnap[pod->curr]);
249: VecGetArray(pod->bsnap[pod->curr],&array);
250: VecPlaceArray(b,array);
251: VecRestoreArray(pod->bsnap[pod->curr],&array);
252: KSP_MatMult(guess->ksp,guess->A,x,b);
253: VecResetArray(b);
254: PetscObjectStateIncrease((PetscObject)pod->bsnap[pod->curr]);
255: if (pod->Aspd) {
256: VecMDot(pod->xsnap[pod->curr],pod->n,pod->bsnap,pod->swork);
257: #if !defined(PETSC_HAVE_MPI_IALLREDUCE)
258: MPIU_Allreduce(pod->swork,pod->swork + 3*pod->n,pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
259: #else
260: MPI_Iallreduce(pod->swork,pod->dots_iallreduce,pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess),&pod->req_iallreduce);
261: pod->ndots_iallreduce = 1;
262: #endif
263: } else {
264: PetscInt off;
265: PetscBool herm;
267: #if defined(PETSC_USE_COMPLEX)
268: MatGetOption(guess->A,MAT_HERMITIAN,&herm);
269: #else
270: MatGetOption(guess->A,MAT_SYMMETRIC,&herm);
271: #endif
272: off = (guess->ksp->transpose_solve && !herm) ? 2*pod->n : pod->n;
274: /* TODO: we may want to use a user-defined dot for the correlation matrix */
275: VecMDot(pod->xsnap[pod->curr],pod->n,pod->xsnap,pod->swork);
276: VecMDot(pod->bsnap[pod->curr],pod->n,pod->xsnap,pod->swork + off);
277: if (!herm) {
278: off = (off == pod->n) ? 2*pod->n : pod->n;
279: VecMDot(pod->xsnap[pod->curr],pod->n,pod->bsnap,pod->swork + off);
280: #if !defined(PETSC_HAVE_MPI_IALLREDUCE)
281: MPIU_Allreduce(pod->swork,pod->swork + 3*pod->n,3*pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
282: #else
283: MPI_Iallreduce(pod->swork,pod->dots_iallreduce,3*pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess),&pod->req_iallreduce);
284: pod->ndots_iallreduce = 3;
285: #endif
286: } else {
287: #if !defined(PETSC_HAVE_MPI_IALLREDUCE)
288: MPIU_Allreduce(pod->swork,pod->swork + 3*pod->n,2*pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
289: for (i=0;i<pod->n;i++) pod->swork[5*pod->n + i] = pod->swork[4*pod->n + i];
290: #else
291: MPI_Iallreduce(pod->swork,pod->dots_iallreduce,2*pod->n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess),&pod->req_iallreduce);
292: pod->ndots_iallreduce = 2;
293: #endif
294: }
295: }
296: if (pod->ndots_iallreduce) return(0);
298: complete_request:
299: if (pod->ndots_iallreduce) {
300: MPI_Wait(&pod->req_iallreduce,MPI_STATUS_IGNORE);
301: switch (pod->ndots_iallreduce) {
302: case 3:
303: for (i=0;i<pod->n;i++) pod->swork[3*pod->n + i] = pod->dots_iallreduce[ i];
304: for (i=0;i<pod->n;i++) pod->swork[4*pod->n + i] = pod->dots_iallreduce[ pod->n+i];
305: for (i=0;i<pod->n;i++) pod->swork[5*pod->n + i] = pod->dots_iallreduce[2*pod->n+i];
306: break;
307: case 2:
308: for (i=0;i<pod->n;i++) pod->swork[3*pod->n + i] = pod->dots_iallreduce[ i];
309: for (i=0;i<pod->n;i++) pod->swork[4*pod->n + i] = pod->dots_iallreduce[pod->n+i];
310: for (i=0;i<pod->n;i++) pod->swork[5*pod->n + i] = pod->dots_iallreduce[pod->n+i];
311: break;
312: case 1:
313: for (i=0;i<pod->n;i++) pod->swork[3*pod->n + i] = pod->dots_iallreduce[i];
314: break;
315: default:
316: SETERRQ1(PetscObjectComm((PetscObject)guess),PETSC_ERR_PLIB,"Invalid number of outstanding dots operations: %D",pod->ndots_iallreduce);
317: break;
318: }
319: }
320: pod->ndots_iallreduce = 0;
322: /* correlation matrix and Y^H A Y (Galerkin) */
323: for (i=0;i<pod->n;i++) {
324: pod->corr[pod->curr*pod->maxn+i] = pod->swork[3*pod->n + i];
325: pod->corr[i*pod->maxn+pod->curr] = PetscConj(pod->swork[3*pod->n + i]);
326: if (!pod->Aspd) {
327: pod->yhay[pod->curr*pod->maxn+i] = pod->swork[4*pod->n + i];
328: pod->yhay[i*pod->maxn+pod->curr] = PetscConj(pod->swork[5*pod->n + i]);
329: }
330: }
331: /* syevx change the input matrix */
332: for (i=0;i<pod->n;i++) {
333: PetscInt j;
334: for (j=i;j<pod->n;j++) pod->swork[i*pod->n+j] = pod->corr[i*pod->maxn+j];
335: }
336: PetscBLASIntCast(pod->n,&bN);
337: #if !defined(PETSC_USE_COMPLEX)
338: PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","L",&bN,pod->swork,&bN,
339: &reps,&reps,&idummy,&idummy,
340: &reps,&idummy,pod->eigs,pod->eigv,&bN,
341: pod->swork+bN*bN,&pod->lwork,pod->iwork,pod->iwork+5*bN,&lierr));
342: #else
343: PetscStackCallBLAS("LAPACKsyevx",LAPACKsyevx_("V","A","L",&bN,pod->swork,&bN,
344: &reps,&reps,&idummy,&idummy,
345: &reps,&idummy,pod->eigs,pod->eigv,&bN,
346: pod->swork+bN*bN,&pod->lwork,pod->rwork,pod->iwork,pod->iwork+5*bN,&lierr));
347: #endif
348: if (lierr<0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYEV Lapack routine: illegal argument %d",-(int)lierr);
349: else if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYEV Lapack routine: %d eigenvectors failed to converge",(int)lierr);
351: /* dimension of lower dimensional system */
352: pod->st = -1;
353: for (i=0,toten=0;i<pod->n;i++) {
354: pod->eigs[i] = PetscMax(pod->eigs[i],0.0);
355: toten += pod->eigs[i];
356: if (!pod->eigs[i]) pod->st = i;
357: }
358: pod->nen = 0;
359: for (i=pod->n-1,parten=0;i>pod->st && toten > 0;i--) {
360: pod->nen++;
361: parten += pod->eigs[i];
362: if (parten + toten*pod->tol >= toten) break;
363: }
364: pod->st = pod->n - pod->nen;
366: /* Compute eigv = V * S */
367: for (i=pod->st;i<pod->n;i++) {
368: const PetscReal v = 1.0/PetscSqrtReal(pod->eigs[i]);
369: const PetscInt st = pod->n*i;
370: PetscInt j;
372: for (j=0;j<pod->n;j++) pod->eigv[st+j] *= v;
373: }
375: /* compute S * V^T * X^T * A * X * V * S if needed */
376: if (pod->nen && !pod->Aspd) {
377: PetscBLASInt bNen,bMaxN;
378: PetscInt st = pod->st*pod->n;
379: PetscBLASIntCast(pod->nen,&bNen);
380: PetscBLASIntCast(pod->maxn,&bMaxN);
381: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&bNen,&bN,&bN,&one,pod->eigv+st,&bN,pod->yhay,&bMaxN,&zero,pod->swork,&bNen));
382: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&bNen,&bNen,&bN,&one,pod->swork,&bNen,pod->eigv+st,&bN,&zero,pod->low,&bNen));
383: }
385: if (pod->monitor) {
386: PetscPrintf(PetscObjectComm((PetscObject)guess)," KSPGuessPOD: basis %D, energy fractions = ",pod->nen);
387: for (i=pod->n-1;i>=0;i--) {
388: PetscPrintf(PetscObjectComm((PetscObject)guess),"%1.6e (%d) ",pod->eigs[i]/toten,i >= pod->st ? 1 : 0);
389: }
390: PetscPrintf(PetscObjectComm((PetscObject)guess),"\n");
391: #if defined(PETSC_USE_DEBUG)
392: for (i=0;i<pod->n;i++) {
393: Vec v;
394: PetscInt j;
395: PetscBLASInt bNen,ione = 1;
397: VecDuplicate(pod->xsnap[i],&v);
398: VecCopy(pod->xsnap[i],v);
399: PetscBLASIntCast(pod->nen,&bNen);
400: PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&bN,&bNen,&one,pod->eigv+pod->st*pod->n,&bN,pod->corr+pod->maxn*i,&ione,&zero,pod->swork,&ione));
401: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&bN,&bNen,&one,pod->eigv+pod->st*pod->n,&bN,pod->swork,&ione,&zero,pod->swork+pod->n,&ione));
402: for (j=0;j<pod->n;j++) pod->swork[j] = -pod->swork[pod->n+j];
403: VecMAXPY(v,pod->n,pod->swork,pod->xsnap);
404: VecDot(v,v,pod->swork);
405: MPIU_Allreduce(pod->swork,pod->swork + 1,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)guess));
406: PetscPrintf(PetscObjectComm((PetscObject)guess)," Error projection %D: %g (expected lower than %g)\n",i,(double)PetscRealPart(pod->swork[1]),(double)(toten-parten));
407: VecDestroy(&v);
408: }
409: #endif
410: }
411: /* new tip */
412: pod->curr = (pod->curr+1)%pod->maxn;
413: return(0);
414: }
416: static PetscErrorCode KSPGuessSetFromOptions_POD(KSPGuess guess)
417: {
418: KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
422: PetscOptionsBegin(PetscObjectComm((PetscObject)guess),((PetscObject)guess)->prefix,"POD initial guess options","KSPGuess");
423: PetscOptionsInt("-ksp_guess_pod_size","Number of snapshots",NULL,pod->maxn,&pod->maxn,NULL);
424: PetscOptionsBool("-ksp_guess_pod_monitor","Monitor initial guess generator",NULL,pod->monitor,&pod->monitor,NULL);
425: PetscOptionsReal("-ksp_guess_pod_tol","Tolerance to retain eigenvectors",NULL,pod->tol,&pod->tol,NULL);
426: PetscOptionsBool("-ksp_guess_pod_Ainner","Use the operator as inner product (must be SPD)",NULL,pod->Aspd,&pod->Aspd,NULL);
427: PetscOptionsEnd();
428: return(0);
429: }
431: static PetscErrorCode KSPGuessView_POD(KSPGuess guess,PetscViewer viewer)
432: {
433: KSPGuessPOD *pod = (KSPGuessPOD*)guess->data;
434: PetscBool isascii;
438: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
439: if (isascii) {
440: PetscViewerASCIIPrintf(viewer,"Max size %D, tolerance %g, Ainner %d\n",pod->maxn,pod->tol,pod->Aspd);
441: }
442: return(0);
443: }
445: /*
446: KSPGUESSPOD - Implements a proper orthogonal decomposition based Galerkin scheme for repeated linear system solves.
448: The initial guess is obtained by solving a small and dense linear system, obtained by Galerkin projection on a lower dimensional space generated by the previous solutions.
449: The number of solutions to be retained and the energy tolerance to construct the lower dimensional basis can be specified at command line by -ksp_guess_pod_tol <real> and -ksp_guess_pod_size <int>.
451: References:
452: . 1. - http://www.math.uni-konstanz.de/numerik/personen/volkwein/teaching/POD-Book.pdf
454: Level: intermediate
456: .seealso: KSPGuess, KSPGuessType, KSPGuessCreate(), KSPSetGuess(), KSPGetGuess()
457: @*/
458: PetscErrorCode KSPGuessCreate_POD(KSPGuess guess)
459: {
460: KSPGuessPOD *pod;
464: PetscNewLog(guess,&pod);
465: pod->maxn = 10;
466: pod->tol = PETSC_MACHINE_EPSILON;
467: guess->data = pod;
469: guess->ops->setfromoptions = KSPGuessSetFromOptions_POD;
470: guess->ops->destroy = KSPGuessDestroy_POD;
471: guess->ops->setup = KSPGuessSetUp_POD;
472: guess->ops->view = KSPGuessView_POD;
473: guess->ops->reset = KSPGuessReset_POD;
474: guess->ops->update = KSPGuessUpdate_POD;
475: guess->ops->formguess = KSPGuessFormGuess_POD;
476: return(0);
477: }