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