Actual source code: pod.c
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[] = "@phdthesis{zampini2010non,\n"
6: " title={Non-overlapping Domain Decomposition Methods for Cardiac Reaction-Diffusion Models and Applications},\n"
7: " author={Zampini, S},\n"
8: " year={2010},\n"
9: " school={PhD thesis, Universita degli Studi di Milano}\n"
10: "}\n";
12: typedef struct {
13: PetscInt maxn; /* maximum number of snapshots */
14: PetscInt n; /* number of active snapshots */
15: PetscInt curr; /* current tip of snapshots set */
16: Vec *xsnap; /* snapshots */
17: Vec *bsnap; /* rhs snapshots */
18: Vec *work; /* parallel work vectors */
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;
43: PetscLayout Alay = NULL, vlay = NULL;
44: PetscBool cong;
46: PetscFunctionBegin;
47: pod->nen = 0;
48: pod->n = 0;
49: pod->curr = 0;
50: /* need to wait for completion of outstanding requests */
51: if (pod->ndots_iallreduce) PetscCallMPI(MPI_Wait(&pod->req_iallreduce, MPI_STATUS_IGNORE));
52: pod->ndots_iallreduce = 0;
53: /* destroy vectors if the size of the linear system has changed */
54: if (guess->A) PetscCall(MatGetLayouts(guess->A, &Alay, NULL));
55: if (pod->xsnap) PetscCall(VecGetLayout(pod->xsnap[0], &vlay));
56: cong = PETSC_FALSE;
57: if (vlay && Alay) PetscCall(PetscLayoutCompare(Alay, vlay, &cong));
58: if (!cong) {
59: PetscCall(VecDestroyVecs(pod->maxn, &pod->xsnap));
60: PetscCall(VecDestroyVecs(pod->maxn, &pod->bsnap));
61: PetscCall(VecDestroyVecs(1, &pod->work));
62: }
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: static PetscErrorCode KSPGuessSetUp_POD(KSPGuess guess)
67: {
68: KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
70: PetscFunctionBegin;
71: if (!pod->corr) {
72: PetscScalar sdummy;
73: PetscReal rdummy = 0;
74: PetscBLASInt bN, lierr, idummy;
76: PetscCall(PetscCalloc6(pod->maxn * pod->maxn, &pod->corr, pod->maxn, &pod->eigs, pod->maxn * pod->maxn, &pod->eigv, 6 * pod->maxn, &pod->iwork, pod->maxn * pod->maxn, &pod->yhay, pod->maxn * pod->maxn, &pod->low));
77: #if defined(PETSC_USE_COMPLEX)
78: PetscCall(PetscMalloc1(7 * pod->maxn, &pod->rwork));
79: #endif
80: #if defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
81: PetscCall(PetscMalloc1(3 * pod->maxn, &pod->dots_iallreduce));
82: #endif
83: pod->lwork = -1;
84: PetscCall(PetscBLASIntCast(pod->maxn, &bN));
85: #if !defined(PETSC_USE_COMPLEX)
86: PetscCallBLAS("LAPACKsyevx", LAPACKsyevx_("V", "A", "L", &bN, pod->corr, &bN, &rdummy, &rdummy, &idummy, &idummy, &rdummy, &idummy, pod->eigs, pod->eigv, &bN, &sdummy, &pod->lwork, pod->iwork, pod->iwork + 5 * bN, &lierr));
87: #else
88: PetscCallBLAS("LAPACKsyevx", LAPACKsyevx_("V", "A", "L", &bN, pod->corr, &bN, &rdummy, &rdummy, &idummy, &idummy, &rdummy, &idummy, pod->eigs, pod->eigv, &bN, &sdummy, &pod->lwork, pod->rwork, pod->iwork, pod->iwork + 5 * bN, &lierr));
89: #endif
90: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in query to SYEV Lapack routine %" PetscBLASInt_FMT, lierr);
91: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(sdummy), &pod->lwork));
92: PetscCall(PetscMalloc1(pod->lwork + PetscMax(bN * bN, 6 * bN), &pod->swork));
93: }
94: /* work vectors are sequential, we explicitly use MPI_Allreduce */
95: if (!pod->xsnap) {
96: Vec *v, vseq;
98: PetscCall(KSPCreateVecs(guess->ksp, 1, &v, 0, NULL));
99: PetscCall(VecCreateLocalVector(v[0], &vseq));
100: PetscCall(VecDestroyVecs(1, &v));
101: PetscCall(VecDuplicateVecs(vseq, pod->maxn, &pod->xsnap));
102: PetscCall(VecDestroy(&vseq));
103: }
104: if (!pod->bsnap) {
105: Vec *v, vseq;
107: PetscCall(KSPCreateVecs(guess->ksp, 0, NULL, 1, &v));
108: PetscCall(VecCreateLocalVector(v[0], &vseq));
109: PetscCall(VecDestroyVecs(1, &v));
110: PetscCall(VecDuplicateVecs(vseq, pod->maxn, &pod->bsnap));
111: PetscCall(VecDestroy(&vseq));
112: }
113: if (!pod->work) PetscCall(KSPCreateVecs(guess->ksp, 1, &pod->work, 0, NULL));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode KSPGuessDestroy_POD(KSPGuess guess)
118: {
119: KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
121: PetscFunctionBegin;
122: PetscCall(PetscFree6(pod->corr, pod->eigs, pod->eigv, pod->iwork, pod->yhay, pod->low));
123: #if defined(PETSC_USE_COMPLEX)
124: PetscCall(PetscFree(pod->rwork));
125: #endif
126: /* need to wait for completion before destroying dots_iallreduce */
127: if (pod->ndots_iallreduce) PetscCallMPI(MPI_Wait(&pod->req_iallreduce, MPI_STATUS_IGNORE));
128: PetscCall(PetscFree(pod->dots_iallreduce));
129: PetscCall(PetscFree(pod->swork));
130: PetscCall(VecDestroyVecs(pod->maxn, &pod->bsnap));
131: PetscCall(VecDestroyVecs(pod->maxn, &pod->xsnap));
132: PetscCall(VecDestroyVecs(1, &pod->work));
133: PetscCall(PetscFree(pod));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode KSPGuessUpdate_POD(KSPGuess, Vec, Vec);
139: static PetscErrorCode KSPGuessFormGuess_POD(KSPGuess guess, Vec b, Vec x)
140: {
141: KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
142: PetscScalar one = 1, zero = 0;
143: PetscBLASInt bN, ione = 1, bNen, lierr;
144: PetscInt i;
145: PetscMPIInt podn;
147: PetscFunctionBegin;
148: PetscCall(PetscCitationsRegister(citation, &cited));
149: if (pod->ndots_iallreduce) { /* complete communication and project the linear system */
150: PetscCall(KSPGuessUpdate_POD(guess, NULL, NULL));
151: }
152: if (!pod->nen) PetscFunctionReturn(PETSC_SUCCESS);
153: /* b_low = S * V^T * X^T * b */
154: PetscCall(VecGetLocalVectorRead(b, pod->bsnap[pod->curr]));
155: PetscCall(VecMDot(pod->bsnap[pod->curr], pod->n, pod->xsnap, pod->swork));
156: PetscCall(VecRestoreLocalVectorRead(b, pod->bsnap[pod->curr]));
157: PetscCall(PetscMPIIntCast(pod->n, &podn));
158: PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + pod->n, podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
159: PetscCall(PetscBLASIntCast(pod->n, &bN));
160: PetscCall(PetscBLASIntCast(pod->nen, &bNen));
161: PetscCallBLAS("BLASgemv", BLASgemv_("T", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->swork + pod->n, &ione, &zero, pod->swork, &ione));
162: if (pod->monitor) {
163: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " KSPGuessPOD alphas = "));
164: for (i = 0; i < pod->nen; i++) {
165: #if defined(PETSC_USE_COMPLEX)
166: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g + %g i", (double)PetscRealPart(pod->swork[i]), (double)PetscImaginaryPart(pod->swork[i])));
167: #else
168: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g ", (double)pod->swork[i]));
169: #endif
170: }
171: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "\n"));
172: }
173: /* A_low x_low = b_low */
174: if (!pod->Aspd) { /* A is spd -> LOW = Identity */
175: KSP pksp = guess->ksp;
176: PetscBool tsolve, symm, set;
178: if (pod->monitor) {
179: PetscMPIInt rank;
180: Mat L;
182: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)guess), &rank));
183: if (rank == 0) {
184: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " L = "));
185: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pod->nen, pod->nen, pod->low, &L));
186: PetscCall(MatView(L, NULL));
187: PetscCall(MatDestroy(&L));
188: }
189: }
190: PetscCall(MatIsSymmetricKnown(guess->A, &set, &symm));
191: tsolve = (set && symm) ? PETSC_FALSE : pksp->transpose_solve;
192: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bNen, &bNen, pod->low, &bNen, pod->iwork, &lierr));
193: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRF Lapack routine %" PetscBLASInt_FMT, lierr);
194: PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_(tsolve ? "T" : "N", &bNen, &ione, pod->low, &bNen, pod->iwork, pod->swork, &bNen, &lierr));
195: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRS Lapack routine %" PetscBLASInt_FMT, lierr);
196: }
197: /* x = X * V * S * x_low */
198: PetscCallBLAS("BLASgemv", BLASgemv_("N", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->swork, &ione, &zero, pod->swork + pod->n, &ione));
199: if (pod->monitor) {
200: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " KSPGuessPOD sol = "));
201: for (i = 0; i < pod->nen; i++) {
202: #if defined(PETSC_USE_COMPLEX)
203: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g + %g i", (double)PetscRealPart(pod->swork[i + pod->n]), (double)PetscImaginaryPart(pod->swork[i + pod->n])));
204: #else
205: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g ", (double)pod->swork[i + pod->n]));
206: #endif
207: }
208: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "\n"));
209: }
210: PetscCall(VecGetLocalVector(x, pod->bsnap[pod->curr]));
211: PetscCall(VecSet(pod->bsnap[pod->curr], 0));
212: PetscCall(VecMAXPY(pod->bsnap[pod->curr], pod->n, pod->swork + pod->n, pod->xsnap));
213: PetscCall(VecRestoreLocalVector(x, pod->bsnap[pod->curr]));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static PetscErrorCode KSPGuessUpdate_POD(KSPGuess guess, Vec b, Vec x)
218: {
219: KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
220: PetscScalar one = 1, zero = 0;
221: PetscReal toten, parten, reps = 0; /* dlamch? */
222: PetscBLASInt bN, lierr, idummy;
223: PetscInt i;
224: PetscMPIInt podn;
226: PetscFunctionBegin;
227: if (pod->ndots_iallreduce) goto complete_request;
228: pod->n = pod->n < pod->maxn ? pod->n + 1 : pod->maxn;
229: PetscCall(PetscMPIIntCast(pod->n, &podn));
230: PetscCall(VecCopy(x, pod->xsnap[pod->curr]));
231: PetscCall(KSP_MatMult(guess->ksp, guess->A, x, pod->work[0]));
232: PetscCall(VecCopy(pod->work[0], pod->bsnap[pod->curr]));
233: if (pod->Aspd) {
234: PetscCall(VecMDot(pod->xsnap[pod->curr], pod->n, pod->bsnap, pod->swork));
235: #if !defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
236: PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 3 * pod->n, podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
237: #else
238: PetscCallMPI(MPI_Iallreduce(pod->swork, pod->dots_iallreduce, podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess), &pod->req_iallreduce));
239: pod->ndots_iallreduce = 1;
240: #endif
241: } else {
242: PetscInt off;
243: PetscBool set, herm;
245: #if defined(PETSC_USE_COMPLEX)
246: PetscCall(MatIsHermitianKnown(guess->A, &set, &herm));
247: #else
248: PetscCall(MatIsSymmetricKnown(guess->A, &set, &herm));
249: #endif
250: off = (guess->ksp->transpose_solve && (!set || !herm)) ? 2 * pod->n : pod->n;
252: /* TODO: we may want to use a user-defined dot for the correlation matrix */
253: PetscCall(VecMDot(pod->xsnap[pod->curr], pod->n, pod->xsnap, pod->swork));
254: PetscCall(VecMDot(pod->bsnap[pod->curr], pod->n, pod->xsnap, pod->swork + off));
255: if (!set || !herm) {
256: off = (off == pod->n) ? 2 * pod->n : pod->n;
257: PetscCall(VecMDot(pod->xsnap[pod->curr], pod->n, pod->bsnap, pod->swork + off));
258: #if !defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
259: PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 3 * pod->n, 3 * podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
260: #else
261: PetscCallMPI(MPI_Iallreduce(pod->swork, pod->dots_iallreduce, 3 * podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess), &pod->req_iallreduce));
262: pod->ndots_iallreduce = 3;
263: #endif
264: } else {
265: #if !defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
266: PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 3 * pod->n, 2 * podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
267: for (i = 0; i < pod->n; i++) pod->swork[5 * pod->n + i] = pod->swork[4 * pod->n + i];
268: #else
269: PetscCallMPI(MPI_Iallreduce(pod->swork, pod->dots_iallreduce, 2 * podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess), &pod->req_iallreduce));
270: pod->ndots_iallreduce = 2;
271: #endif
272: }
273: }
274: if (pod->ndots_iallreduce) PetscFunctionReturn(PETSC_SUCCESS);
276: complete_request:
277: if (pod->ndots_iallreduce) {
278: PetscCallMPI(MPI_Wait(&pod->req_iallreduce, MPI_STATUS_IGNORE));
279: switch (pod->ndots_iallreduce) {
280: case 3:
281: for (i = 0; i < pod->n; i++) pod->swork[3 * pod->n + i] = pod->dots_iallreduce[i];
282: for (i = 0; i < pod->n; i++) pod->swork[4 * pod->n + i] = pod->dots_iallreduce[pod->n + i];
283: for (i = 0; i < pod->n; i++) pod->swork[5 * pod->n + i] = pod->dots_iallreduce[2 * pod->n + i];
284: break;
285: case 2:
286: for (i = 0; i < pod->n; i++) pod->swork[3 * pod->n + i] = pod->dots_iallreduce[i];
287: for (i = 0; i < pod->n; i++) pod->swork[4 * pod->n + i] = pod->dots_iallreduce[pod->n + i];
288: for (i = 0; i < pod->n; i++) pod->swork[5 * pod->n + i] = pod->dots_iallreduce[pod->n + i];
289: break;
290: case 1:
291: for (i = 0; i < pod->n; i++) pod->swork[3 * pod->n + i] = pod->dots_iallreduce[i];
292: break;
293: default:
294: SETERRQ(PetscObjectComm((PetscObject)guess), PETSC_ERR_PLIB, "Invalid number of outstanding dots operations: %" PetscInt_FMT, pod->ndots_iallreduce);
295: }
296: }
297: pod->ndots_iallreduce = 0;
299: /* correlation matrix and Y^H A Y (Galerkin) */
300: for (i = 0; i < pod->n; i++) {
301: pod->corr[pod->curr * pod->maxn + i] = pod->swork[3 * pod->n + i];
302: pod->corr[i * pod->maxn + pod->curr] = PetscConj(pod->swork[3 * pod->n + i]);
303: if (!pod->Aspd) {
304: pod->yhay[pod->curr * pod->maxn + i] = pod->swork[4 * pod->n + i];
305: pod->yhay[i * pod->maxn + pod->curr] = PetscConj(pod->swork[5 * pod->n + i]);
306: }
307: }
308: /* syevx changes the input matrix */
309: for (i = 0; i < pod->n; i++) {
310: PetscInt j;
311: for (j = i; j < pod->n; j++) pod->swork[i * pod->n + j] = pod->corr[i * pod->maxn + j];
312: }
313: PetscCall(PetscBLASIntCast(pod->n, &bN));
314: #if !defined(PETSC_USE_COMPLEX)
315: PetscCallBLAS("LAPACKsyevx", LAPACKsyevx_("V", "A", "L", &bN, pod->swork, &bN, &reps, &reps, &idummy, &idummy, &reps, &idummy, pod->eigs, pod->eigv, &bN, pod->swork + bN * bN, &pod->lwork, pod->iwork, pod->iwork + 5 * bN, &lierr));
316: #else
317: PetscCallBLAS("LAPACKsyevx", LAPACKsyevx_("V", "A", "L", &bN, pod->swork, &bN, &reps, &reps, &idummy, &idummy, &reps, &idummy, pod->eigs, pod->eigv, &bN, pod->swork + bN * bN, &pod->lwork, pod->rwork, pod->iwork, pod->iwork + 5 * bN, &lierr));
318: #endif
319: PetscCheck(lierr >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in SYEV Lapack routine: illegal argument %" PetscBLASInt_FMT, -lierr);
320: PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in SYEV Lapack routine: %" PetscBLASInt_FMT " eigenvectors failed to converge", lierr);
322: /* dimension of lower dimensional system */
323: pod->st = -1;
324: for (i = 0, toten = 0; i < pod->n; i++) {
325: pod->eigs[i] = PetscMax(pod->eigs[i], 0.0);
326: toten += pod->eigs[i];
327: if (!pod->eigs[i]) pod->st = i;
328: }
329: pod->nen = 0;
330: for (i = pod->n - 1, parten = 0; i > pod->st && toten > 0; i--) {
331: pod->nen++;
332: parten += pod->eigs[i];
333: if (parten + toten * pod->tol >= toten) break;
334: }
335: pod->st = pod->n - pod->nen;
337: /* Compute eigv = V * S */
338: for (i = pod->st; i < pod->n; i++) {
339: const PetscReal v = 1.0 / PetscSqrtReal(pod->eigs[i]);
340: const PetscInt st = pod->n * i;
341: PetscInt j;
343: for (j = 0; j < pod->n; j++) pod->eigv[st + j] *= v;
344: }
346: /* compute S * V^T * X^T * A * X * V * S if needed */
347: if (pod->nen && !pod->Aspd) {
348: PetscBLASInt bNen, bMaxN;
349: PetscInt st = pod->st * pod->n;
350: PetscCall(PetscBLASIntCast(pod->nen, &bNen));
351: PetscCall(PetscBLASIntCast(pod->maxn, &bMaxN));
352: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &bNen, &bN, &bN, &one, pod->eigv + st, &bN, pod->yhay, &bMaxN, &zero, pod->swork, &bNen));
353: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &bNen, &bNen, &bN, &one, pod->swork, &bNen, pod->eigv + st, &bN, &zero, pod->low, &bNen));
354: }
356: if (pod->monitor) {
357: PetscMPIInt rank;
358: Mat C;
360: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)guess), &rank));
361: if (rank == 0) {
362: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " C = "));
363: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pod->n, pod->n, pod->corr, &C));
364: PetscCall(MatDenseSetLDA(C, pod->maxn));
365: PetscCall(MatView(C, NULL));
366: PetscCall(MatDestroy(&C));
367: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " YHAY = "));
368: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pod->n, pod->n, pod->yhay, &C));
369: PetscCall(MatDenseSetLDA(C, pod->maxn));
370: PetscCall(MatView(C, NULL));
371: PetscCall(MatDestroy(&C));
372: }
373: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " KSPGuessPOD: basis %" PetscBLASInt_FMT ", energy fractions = ", pod->nen));
374: for (i = pod->n - 1; i >= 0; i--) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%1.6e (%d) ", (double)(pod->eigs[i] / toten), i >= pod->st ? 1 : 0));
375: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "\n"));
376: if (PetscDefined(USE_DEBUG)) {
377: for (i = 0; i < pod->n; i++) {
378: Vec v;
379: PetscInt j;
380: PetscBLASInt bNen, ione = 1;
382: PetscCall(VecDuplicate(pod->xsnap[i], &v));
383: PetscCall(VecCopy(pod->xsnap[i], v));
384: PetscCall(PetscBLASIntCast(pod->nen, &bNen));
385: PetscCallBLAS("BLASgemv", BLASgemv_("T", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->corr + pod->maxn * i, &ione, &zero, pod->swork, &ione));
386: PetscCallBLAS("BLASgemv", BLASgemv_("N", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->swork, &ione, &zero, pod->swork + pod->n, &ione));
387: for (j = 0; j < pod->n; j++) pod->swork[j] = -pod->swork[pod->n + j];
388: PetscCall(VecMAXPY(v, pod->n, pod->swork, pod->xsnap));
389: PetscCall(VecDot(v, v, pod->swork));
390: PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 1, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
391: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " Error projection %" PetscInt_FMT ": %g (expected lower than %g)\n", i, (double)PetscRealPart(pod->swork[1]), (double)(toten - parten)));
392: PetscCall(VecDestroy(&v));
393: }
394: }
395: }
396: /* new tip */
397: pod->curr = (pod->curr + 1) % pod->maxn;
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: static PetscErrorCode KSPGuessSetFromOptions_POD(KSPGuess guess)
402: {
403: KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
405: PetscFunctionBegin;
406: PetscOptionsBegin(PetscObjectComm((PetscObject)guess), ((PetscObject)guess)->prefix, "POD initial guess options", "KSPGuess");
407: PetscCall(PetscOptionsInt("-ksp_guess_pod_size", "Number of snapshots", NULL, pod->maxn, &pod->maxn, NULL));
408: PetscCall(PetscOptionsBool("-ksp_guess_pod_monitor", "Monitor initial guess generator", NULL, pod->monitor, &pod->monitor, NULL));
409: PetscCall(PetscOptionsReal("-ksp_guess_pod_tol", "Tolerance to retain eigenvectors", "KSPGuessSetTolerance", pod->tol, &pod->tol, NULL));
410: PetscCall(PetscOptionsBool("-ksp_guess_pod_Ainner", "Use the operator as inner product (must be SPD)", NULL, pod->Aspd, &pod->Aspd, NULL));
411: PetscOptionsEnd();
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode KSPGuessSetTolerance_POD(KSPGuess guess, PetscReal tol)
416: {
417: KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
419: PetscFunctionBegin;
420: pod->tol = tol;
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: static PetscErrorCode KSPGuessView_POD(KSPGuess guess, PetscViewer viewer)
425: {
426: KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
427: PetscBool isascii;
429: PetscFunctionBegin;
430: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
431: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "Max size %" PetscInt_FMT ", tolerance %g, Ainner %d\n", pod->maxn, (double)pod->tol, pod->Aspd));
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: /*MC
436: KSPGUESSPOD - Implements a proper orthogonal decomposition based Galerkin scheme for repeated linear system solves.
438: Options Database Keys:
439: + -ksp_guess_pod_size <size> - Number of snapshots
440: . -ksp_guess_pod_monitor <true or false> - Monitor initial guess generator
441: . -ksp_guess_pod_tol <tol> - Tolerance to retain eigenvectors
442: - -ksp_guess_pod_Ainner <true> - Use the operator as inner product (must be SPD)
444: Level: intermediate
446: Note:
447: 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 as presented in {cite}`volkwein2013proper`.
449: .seealso: [](ch_ksp), `KSPGuess`, `KSPGuessType`, `KSPGuessCreate()`, `KSPSetGuess()`, `KSPGetGuess()`
450: M*/
451: PetscErrorCode KSPGuessCreate_POD(KSPGuess guess)
452: {
453: KSPGuessPOD *pod;
455: PetscFunctionBegin;
456: PetscCall(PetscNew(&pod));
457: pod->maxn = 10;
458: pod->tol = PETSC_MACHINE_EPSILON;
459: guess->data = pod;
461: guess->ops->setfromoptions = KSPGuessSetFromOptions_POD;
462: guess->ops->destroy = KSPGuessDestroy_POD;
463: guess->ops->settolerance = KSPGuessSetTolerance_POD;
464: guess->ops->setup = KSPGuessSetUp_POD;
465: guess->ops->view = KSPGuessView_POD;
466: guess->ops->reset = KSPGuessReset_POD;
467: guess->ops->update = KSPGuessUpdate_POD;
468: guess->ops->formguess = KSPGuessFormGuess_POD;
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }