1 #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
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";
11
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;
39
KSPGuessReset_POD(KSPGuess guess)40 static PetscErrorCode KSPGuessReset_POD(KSPGuess guess)
41 {
42 KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
43 PetscLayout Alay = NULL, vlay = NULL;
44 PetscBool cong;
45
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 }
65
KSPGuessSetUp_POD(KSPGuess guess)66 static PetscErrorCode KSPGuessSetUp_POD(KSPGuess guess)
67 {
68 KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
69
70 PetscFunctionBegin;
71 if (!pod->corr) {
72 PetscScalar sdummy;
73 PetscReal rdummy = 0;
74 PetscBLASInt bN, lierr, idummy = 0;
75
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;
97
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;
106
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 }
116
KSPGuessDestroy_POD(KSPGuess guess)117 static PetscErrorCode KSPGuessDestroy_POD(KSPGuess guess)
118 {
119 KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
120
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 }
136
137 static PetscErrorCode KSPGuessUpdate_POD(KSPGuess, Vec, Vec);
138
KSPGuessFormGuess_POD(KSPGuess guess,Vec b,Vec x)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
146 PetscFunctionBegin;
147 PetscCall(PetscCitationsRegister(citation, &cited));
148 if (pod->ndots_iallreduce) { /* complete communication and project the linear system */
149 PetscCall(KSPGuessUpdate_POD(guess, NULL, NULL));
150 }
151 if (!pod->nen) PetscFunctionReturn(PETSC_SUCCESS);
152 /* b_low = S * V^T * X^T * b */
153 PetscCall(VecGetLocalVectorRead(b, pod->bsnap[pod->curr]));
154 PetscCall(VecMDot(pod->bsnap[pod->curr], pod->n, pod->xsnap, pod->swork));
155 PetscCall(VecRestoreLocalVectorRead(b, pod->bsnap[pod->curr]));
156 PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + pod->n, pod->n, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
157 PetscCall(PetscBLASIntCast(pod->n, &bN));
158 PetscCall(PetscBLASIntCast(pod->nen, &bNen));
159 PetscCallBLAS("BLASgemv", BLASgemv_("T", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->swork + pod->n, &ione, &zero, pod->swork, &ione));
160 if (pod->monitor) {
161 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " KSPGuessPOD alphas = "));
162 for (i = 0; i < pod->nen; i++) {
163 #if defined(PETSC_USE_COMPLEX)
164 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g + %g i", (double)PetscRealPart(pod->swork[i]), (double)PetscImaginaryPart(pod->swork[i])));
165 #else
166 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g ", (double)pod->swork[i]));
167 #endif
168 }
169 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "\n"));
170 }
171 /* A_low x_low = b_low */
172 if (!pod->Aspd) { /* A is spd -> LOW = Identity */
173 KSP pksp = guess->ksp;
174 PetscBool tsolve, symm, set;
175
176 if (pod->monitor) {
177 PetscMPIInt rank;
178 Mat L;
179
180 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)guess), &rank));
181 if (rank == 0) {
182 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " L = "));
183 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pod->nen, pod->nen, pod->low, &L));
184 PetscCall(MatView(L, NULL));
185 PetscCall(MatDestroy(&L));
186 }
187 }
188 PetscCall(MatIsSymmetricKnown(guess->A, &set, &symm));
189 tsolve = (set && symm) ? PETSC_FALSE : pksp->transpose_solve;
190 PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&bNen, &bNen, pod->low, &bNen, pod->iwork, &lierr));
191 PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRF Lapack routine %" PetscBLASInt_FMT, lierr);
192 PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_(tsolve ? "T" : "N", &bNen, &ione, pod->low, &bNen, pod->iwork, pod->swork, &bNen, &lierr));
193 PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in GETRS Lapack routine %" PetscBLASInt_FMT, lierr);
194 }
195 /* x = X * V * S * x_low */
196 PetscCallBLAS("BLASgemv", BLASgemv_("N", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->swork, &ione, &zero, pod->swork + pod->n, &ione));
197 if (pod->monitor) {
198 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " KSPGuessPOD sol = "));
199 for (i = 0; i < pod->nen; i++) {
200 #if defined(PETSC_USE_COMPLEX)
201 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g + %g i", (double)PetscRealPart(pod->swork[i + pod->n]), (double)PetscImaginaryPart(pod->swork[i + pod->n])));
202 #else
203 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "%g ", (double)pod->swork[i + pod->n]));
204 #endif
205 }
206 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "\n"));
207 }
208 PetscCall(VecGetLocalVector(x, pod->bsnap[pod->curr]));
209 PetscCall(VecSet(pod->bsnap[pod->curr], 0));
210 PetscCall(VecMAXPY(pod->bsnap[pod->curr], pod->n, pod->swork + pod->n, pod->xsnap));
211 PetscCall(VecRestoreLocalVector(x, pod->bsnap[pod->curr]));
212 PetscFunctionReturn(PETSC_SUCCESS);
213 }
214
KSPGuessUpdate_POD(KSPGuess guess,Vec b,Vec x)215 static PetscErrorCode KSPGuessUpdate_POD(KSPGuess guess, Vec b, Vec x)
216 {
217 KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
218 PetscScalar one = 1, zero = 0;
219 PetscReal toten, parten, reps = 0; /* dlamch? */
220 PetscBLASInt bN, lierr, idummy = 0;
221 PetscInt i;
222 PetscMPIInt podn;
223
224 PetscFunctionBegin;
225 if (pod->ndots_iallreduce) goto complete_request;
226 pod->n = pod->n < pod->maxn ? pod->n + 1 : pod->maxn;
227 PetscCall(PetscMPIIntCast(pod->n, &podn));
228 PetscCall(VecCopy(x, pod->xsnap[pod->curr]));
229 PetscCall(KSP_MatMult(guess->ksp, guess->A, x, pod->work[0]));
230 PetscCall(VecCopy(pod->work[0], pod->bsnap[pod->curr]));
231 if (pod->Aspd) {
232 PetscCall(VecMDot(pod->xsnap[pod->curr], pod->n, pod->bsnap, pod->swork));
233 #if !defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
234 PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 3 * pod->n, podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
235 #else
236 PetscCallMPI(MPI_Iallreduce(pod->swork, pod->dots_iallreduce, podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess), &pod->req_iallreduce));
237 pod->ndots_iallreduce = 1;
238 #endif
239 } else {
240 PetscInt off;
241 PetscBool set, herm;
242
243 #if defined(PETSC_USE_COMPLEX)
244 PetscCall(MatIsHermitianKnown(guess->A, &set, &herm));
245 #else
246 PetscCall(MatIsSymmetricKnown(guess->A, &set, &herm));
247 #endif
248 off = (guess->ksp->transpose_solve && (!set || !herm)) ? 2 * pod->n : pod->n;
249
250 /* TODO: we may want to use a user-defined dot for the correlation matrix */
251 PetscCall(VecMDot(pod->xsnap[pod->curr], pod->n, pod->xsnap, pod->swork));
252 PetscCall(VecMDot(pod->bsnap[pod->curr], pod->n, pod->xsnap, pod->swork + off));
253 if (!set || !herm) {
254 off = (off == pod->n) ? 2 * pod->n : pod->n;
255 PetscCall(VecMDot(pod->xsnap[pod->curr], pod->n, pod->bsnap, pod->swork + off));
256 #if !defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
257 PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 3 * pod->n, 3 * podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
258 #else
259 PetscCallMPI(MPI_Iallreduce(pod->swork, pod->dots_iallreduce, 3 * podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess), &pod->req_iallreduce));
260 pod->ndots_iallreduce = 3;
261 #endif
262 } else {
263 #if !defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
264 PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 3 * pod->n, 2 * podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
265 for (i = 0; i < pod->n; i++) pod->swork[5 * pod->n + i] = pod->swork[4 * pod->n + i];
266 #else
267 PetscCallMPI(MPI_Iallreduce(pod->swork, pod->dots_iallreduce, 2 * podn, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess), &pod->req_iallreduce));
268 pod->ndots_iallreduce = 2;
269 #endif
270 }
271 }
272 if (pod->ndots_iallreduce) PetscFunctionReturn(PETSC_SUCCESS);
273
274 complete_request:
275 if (pod->ndots_iallreduce) {
276 PetscCallMPI(MPI_Wait(&pod->req_iallreduce, MPI_STATUS_IGNORE));
277 switch (pod->ndots_iallreduce) {
278 case 3:
279 for (i = 0; i < pod->n; i++) pod->swork[3 * pod->n + i] = pod->dots_iallreduce[i];
280 for (i = 0; i < pod->n; i++) pod->swork[4 * pod->n + i] = pod->dots_iallreduce[pod->n + i];
281 for (i = 0; i < pod->n; i++) pod->swork[5 * pod->n + i] = pod->dots_iallreduce[2 * pod->n + i];
282 break;
283 case 2:
284 for (i = 0; i < pod->n; i++) pod->swork[3 * pod->n + i] = pod->dots_iallreduce[i];
285 for (i = 0; i < pod->n; i++) pod->swork[4 * pod->n + i] = pod->dots_iallreduce[pod->n + i];
286 for (i = 0; i < pod->n; i++) pod->swork[5 * pod->n + i] = pod->dots_iallreduce[pod->n + i];
287 break;
288 case 1:
289 for (i = 0; i < pod->n; i++) pod->swork[3 * pod->n + i] = pod->dots_iallreduce[i];
290 break;
291 default:
292 SETERRQ(PetscObjectComm((PetscObject)guess), PETSC_ERR_PLIB, "Invalid number of outstanding dots operations: %" PetscInt_FMT, pod->ndots_iallreduce);
293 }
294 }
295 pod->ndots_iallreduce = 0;
296
297 /* correlation matrix and Y^H A Y (Galerkin) */
298 for (i = 0; i < pod->n; i++) {
299 pod->corr[pod->curr * pod->maxn + i] = pod->swork[3 * pod->n + i];
300 pod->corr[i * pod->maxn + pod->curr] = PetscConj(pod->swork[3 * pod->n + i]);
301 if (!pod->Aspd) {
302 pod->yhay[pod->curr * pod->maxn + i] = pod->swork[4 * pod->n + i];
303 pod->yhay[i * pod->maxn + pod->curr] = PetscConj(pod->swork[5 * pod->n + i]);
304 }
305 }
306 /* syevx changes the input matrix */
307 for (i = 0; i < pod->n; i++) {
308 PetscInt j;
309 for (j = i; j < pod->n; j++) pod->swork[i * pod->n + j] = pod->corr[i * pod->maxn + j];
310 }
311 PetscCall(PetscBLASIntCast(pod->n, &bN));
312 #if !defined(PETSC_USE_COMPLEX)
313 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));
314 #else
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->rwork, pod->iwork, pod->iwork + 5 * bN, &lierr));
316 #endif
317 PetscCheck(lierr >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in SYEV Lapack routine: illegal argument %" PetscBLASInt_FMT, -lierr);
318 PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in SYEV Lapack routine: %" PetscBLASInt_FMT " eigenvectors failed to converge", lierr);
319
320 /* dimension of lower dimensional system */
321 pod->st = -1;
322 for (i = 0, toten = 0; i < pod->n; i++) {
323 pod->eigs[i] = PetscMax(pod->eigs[i], 0.0);
324 toten += pod->eigs[i];
325 if (!pod->eigs[i]) pod->st = i;
326 }
327 pod->nen = 0;
328 for (i = pod->n - 1, parten = 0; i > pod->st && toten > 0; i--) {
329 pod->nen++;
330 parten += pod->eigs[i];
331 if (parten + toten * pod->tol >= toten) break;
332 }
333 pod->st = pod->n - pod->nen;
334
335 /* Compute eigv = V * S */
336 for (i = pod->st; i < pod->n; i++) {
337 const PetscReal v = 1.0 / PetscSqrtReal(pod->eigs[i]);
338 const PetscInt st = pod->n * i;
339 PetscInt j;
340
341 for (j = 0; j < pod->n; j++) pod->eigv[st + j] *= v;
342 }
343
344 /* compute S * V^T * X^T * A * X * V * S if needed */
345 if (pod->nen && !pod->Aspd) {
346 PetscBLASInt bNen, bMaxN;
347 PetscInt st = pod->st * pod->n;
348 PetscCall(PetscBLASIntCast(pod->nen, &bNen));
349 PetscCall(PetscBLASIntCast(pod->maxn, &bMaxN));
350 PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &bNen, &bN, &bN, &one, pod->eigv + st, &bN, pod->yhay, &bMaxN, &zero, pod->swork, &bNen));
351 PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &bNen, &bNen, &bN, &one, pod->swork, &bNen, pod->eigv + st, &bN, &zero, pod->low, &bNen));
352 }
353
354 if (pod->monitor) {
355 PetscMPIInt rank;
356 Mat C;
357
358 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)guess), &rank));
359 if (rank == 0) {
360 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " C = "));
361 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pod->n, pod->n, pod->corr, &C));
362 PetscCall(MatDenseSetLDA(C, pod->maxn));
363 PetscCall(MatView(C, NULL));
364 PetscCall(MatDestroy(&C));
365 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " YHAY = "));
366 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, pod->n, pod->n, pod->yhay, &C));
367 PetscCall(MatDenseSetLDA(C, pod->maxn));
368 PetscCall(MatView(C, NULL));
369 PetscCall(MatDestroy(&C));
370 }
371 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " KSPGuessPOD: basis %" PetscBLASInt_FMT ", energy fractions = ", pod->nen));
372 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));
373 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), "\n"));
374 if (PetscDefined(USE_DEBUG)) {
375 for (i = 0; i < pod->n; i++) {
376 Vec v;
377 PetscInt j;
378 PetscBLASInt bNen, ione = 1;
379
380 PetscCall(VecDuplicate(pod->xsnap[i], &v));
381 PetscCall(VecCopy(pod->xsnap[i], v));
382 PetscCall(PetscBLASIntCast(pod->nen, &bNen));
383 PetscCallBLAS("BLASgemv", BLASgemv_("T", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->corr + pod->maxn * i, &ione, &zero, pod->swork, &ione));
384 PetscCallBLAS("BLASgemv", BLASgemv_("N", &bN, &bNen, &one, pod->eigv + pod->st * pod->n, &bN, pod->swork, &ione, &zero, pod->swork + pod->n, &ione));
385 for (j = 0; j < pod->n; j++) pod->swork[j] = -pod->swork[pod->n + j];
386 PetscCall(VecMAXPY(v, pod->n, pod->swork, pod->xsnap));
387 PetscCall(VecDot(v, v, pod->swork));
388 PetscCallMPI(MPIU_Allreduce(pod->swork, pod->swork + 1, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)guess)));
389 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)guess), " Error projection %" PetscInt_FMT ": %g (expected lower than %g)\n", i, (double)PetscRealPart(pod->swork[1]), (double)(toten - parten)));
390 PetscCall(VecDestroy(&v));
391 }
392 }
393 }
394 /* new tip */
395 pod->curr = (pod->curr + 1) % pod->maxn;
396 PetscFunctionReturn(PETSC_SUCCESS);
397 }
398
KSPGuessSetFromOptions_POD(KSPGuess guess)399 static PetscErrorCode KSPGuessSetFromOptions_POD(KSPGuess guess)
400 {
401 KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
402
403 PetscFunctionBegin;
404 PetscOptionsBegin(PetscObjectComm((PetscObject)guess), ((PetscObject)guess)->prefix, "POD initial guess options", "KSPGuess");
405 PetscCall(PetscOptionsInt("-ksp_guess_pod_size", "Number of snapshots", NULL, pod->maxn, &pod->maxn, NULL));
406 PetscCall(PetscOptionsBool("-ksp_guess_pod_monitor", "Monitor initial guess generator", NULL, pod->monitor, &pod->monitor, NULL));
407 PetscCall(PetscOptionsReal("-ksp_guess_pod_tol", "Tolerance to retain eigenvectors", "KSPGuessSetTolerance", pod->tol, &pod->tol, NULL));
408 PetscCall(PetscOptionsBool("-ksp_guess_pod_Ainner", "Use the operator as inner product (must be SPD)", NULL, pod->Aspd, &pod->Aspd, NULL));
409 PetscOptionsEnd();
410 PetscFunctionReturn(PETSC_SUCCESS);
411 }
412
KSPGuessSetTolerance_POD(KSPGuess guess,PetscReal tol)413 static PetscErrorCode KSPGuessSetTolerance_POD(KSPGuess guess, PetscReal tol)
414 {
415 KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
416
417 PetscFunctionBegin;
418 pod->tol = tol;
419 PetscFunctionReturn(PETSC_SUCCESS);
420 }
421
KSPGuessView_POD(KSPGuess guess,PetscViewer viewer)422 static PetscErrorCode KSPGuessView_POD(KSPGuess guess, PetscViewer viewer)
423 {
424 KSPGuessPOD *pod = (KSPGuessPOD *)guess->data;
425 PetscBool isascii;
426
427 PetscFunctionBegin;
428 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
429 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "Max size %" PetscInt_FMT ", tolerance %g, Ainner %d\n", pod->maxn, (double)pod->tol, pod->Aspd));
430 PetscFunctionReturn(PETSC_SUCCESS);
431 }
432
433 /*MC
434 KSPGUESSPOD - Implements a proper orthogonal decomposition based Galerkin scheme for repeated linear system solves.
435
436 Options Database Keys:
437 + -ksp_guess_pod_size <size> - Number of snapshots
438 . -ksp_guess_pod_monitor <true or false> - Monitor initial guess generator
439 . -ksp_guess_pod_tol <tol> - Tolerance to retain eigenvectors
440 - -ksp_guess_pod_Ainner <true> - Use the operator as inner product (must be SPD)
441
442 Level: intermediate
443
444 Note:
445 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`.
446
447 .seealso: [](ch_ksp), `KSPGuess`, `KSPGuessType`, `KSPGuessCreate()`, `KSPSetGuess()`, `KSPGetGuess()`
448 M*/
KSPGuessCreate_POD(KSPGuess guess)449 PetscErrorCode KSPGuessCreate_POD(KSPGuess guess)
450 {
451 KSPGuessPOD *pod;
452
453 PetscFunctionBegin;
454 PetscCall(PetscNew(&pod));
455 pod->maxn = 10;
456 pod->tol = PETSC_MACHINE_EPSILON;
457 guess->data = pod;
458
459 guess->ops->setfromoptions = KSPGuessSetFromOptions_POD;
460 guess->ops->destroy = KSPGuessDestroy_POD;
461 guess->ops->settolerance = KSPGuessSetTolerance_POD;
462 guess->ops->setup = KSPGuessSetUp_POD;
463 guess->ops->view = KSPGuessView_POD;
464 guess->ops->reset = KSPGuessReset_POD;
465 guess->ops->update = KSPGuessUpdate_POD;
466 guess->ops->formguess = KSPGuessFormGuess_POD;
467 PetscFunctionReturn(PETSC_SUCCESS);
468 }
469