xref: /petsc/src/ksp/ksp/impls/gmres/gmreig.c (revision 750b007cd8d816cecd9de99077bb0a703b4cf61a) 
1 
2 #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>
3 #include <petscblaslapack.h>
4 
5 PetscErrorCode KSPComputeExtremeSingularValues_GMRES(KSP ksp, PetscReal *emax, PetscReal *emin) {
6   KSP_GMRES   *gmres = (KSP_GMRES *)ksp->data;
7   PetscInt     n = gmres->it + 1, i, N = gmres->max_k + 2;
8   PetscBLASInt bn, bN, lwork, idummy, lierr;
9   PetscScalar *R = gmres->Rsvd, *work = R + N * N, sdummy = 0;
10   PetscReal   *realpart = gmres->Dsvd;
11 
12   PetscFunctionBegin;
13   PetscCall(PetscBLASIntCast(n, &bn));
14   PetscCall(PetscBLASIntCast(N, &bN));
15   PetscCall(PetscBLASIntCast(5 * N, &lwork));
16   PetscCall(PetscBLASIntCast(N, &idummy));
17   if (n <= 0) {
18     *emax = *emin = 1.0;
19     PetscFunctionReturn(0);
20   }
21   /* copy R matrix to work space */
22   PetscCall(PetscArraycpy(R, gmres->hh_origin, (gmres->max_k + 2) * (gmres->max_k + 1)));
23 
24   /* zero below diagonal garbage */
25   for (i = 0; i < n; i++) R[i * N + i + 1] = 0.0;
26 
27   /* compute Singular Values */
28   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
29 #if !defined(PETSC_USE_COMPLEX)
30   PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("N", "N", &bn, &bn, R, &bN, realpart, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &lierr));
31 #else
32   PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("N", "N", &bn, &bn, R, &bN, realpart, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, realpart + N, &lierr));
33 #endif
34   PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in SVD Lapack routine %d", (int)lierr);
35   PetscCall(PetscFPTrapPop());
36 
37   *emin = realpart[n - 1];
38   *emax = realpart[0];
39   PetscFunctionReturn(0);
40 }
41 
42 PetscErrorCode KSPComputeEigenvalues_GMRES(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig) {
43 #if !defined(PETSC_USE_COMPLEX)
44   KSP_GMRES   *gmres = (KSP_GMRES *)ksp->data;
45   PetscInt     n = gmres->it + 1, N = gmres->max_k + 1, i, *perm;
46   PetscBLASInt bn, bN, lwork, idummy, lierr = -1;
47   PetscScalar *R = gmres->Rsvd, *work = R + N * N;
48   PetscScalar *realpart = gmres->Dsvd, *imagpart = realpart + N, sdummy = 0;
49 
50   PetscFunctionBegin;
51   PetscCall(PetscBLASIntCast(n, &bn));
52   PetscCall(PetscBLASIntCast(N, &bN));
53   PetscCall(PetscBLASIntCast(5 * N, &lwork));
54   PetscCall(PetscBLASIntCast(N, &idummy));
55   PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
56   *neig = n;
57 
58   if (!n) PetscFunctionReturn(0);
59 
60   /* copy R matrix to work space */
61   PetscCall(PetscArraycpy(R, gmres->hes_origin, N * N));
62 
63   /* compute eigenvalues */
64   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
65   PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "N", &bn, R, &bN, realpart, imagpart, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &lierr));
66   PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int)lierr);
67   PetscCall(PetscFPTrapPop());
68   PetscCall(PetscMalloc1(n, &perm));
69   for (i = 0; i < n; i++) perm[i] = i;
70   PetscCall(PetscSortRealWithPermutation(n, realpart, perm));
71   for (i = 0; i < n; i++) {
72     r[i] = realpart[perm[i]];
73     c[i] = imagpart[perm[i]];
74   }
75   PetscCall(PetscFree(perm));
76 #else
77   KSP_GMRES   *gmres = (KSP_GMRES *)ksp->data;
78   PetscInt     n = gmres->it + 1, N = gmres->max_k + 1, i, *perm;
79   PetscScalar *R = gmres->Rsvd, *work = R + N * N, *eigs = work + 5 * N, sdummy;
80   PetscBLASInt bn, bN, lwork, idummy, lierr = -1;
81 
82   PetscFunctionBegin;
83   PetscCall(PetscBLASIntCast(n, &bn));
84   PetscCall(PetscBLASIntCast(N, &bN));
85   PetscCall(PetscBLASIntCast(5 * N, &lwork));
86   PetscCall(PetscBLASIntCast(N, &idummy));
87   PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
88   *neig = n;
89 
90   if (!n) PetscFunctionReturn(0);
91 
92   /* copy R matrix to work space */
93   PetscCall(PetscArraycpy(R, gmres->hes_origin, N * N));
94 
95   /* compute eigenvalues */
96   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
97   PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "N", &bn, R, &bN, eigs, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, gmres->Dsvd, &lierr));
98   PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine");
99   PetscCall(PetscFPTrapPop());
100   PetscCall(PetscMalloc1(n, &perm));
101   for (i = 0; i < n; i++) perm[i] = i;
102   for (i = 0; i < n; i++) r[i] = PetscRealPart(eigs[i]);
103   PetscCall(PetscSortRealWithPermutation(n, r, perm));
104   for (i = 0; i < n; i++) {
105     r[i] = PetscRealPart(eigs[perm[i]]);
106     c[i] = PetscImaginaryPart(eigs[perm[i]]);
107   }
108   PetscCall(PetscFree(perm));
109 #endif
110   PetscFunctionReturn(0);
111 }
112 
113 PetscErrorCode KSPComputeRitz_GMRES(KSP ksp, PetscBool ritz, PetscBool small, PetscInt *nrit, Vec S[], PetscReal *tetar, PetscReal *tetai) {
114   KSP_GMRES   *gmres = (KSP_GMRES *)ksp->data;
115   PetscInt     NbrRitz, nb = 0, n;
116   PetscInt     i, j, *perm;
117   PetscScalar *H, *Q, *Ht; /* H Hessenberg matrix; Q matrix of eigenvectors of H */
118   PetscScalar *wr, *wi;    /* Real and imaginary part of the Ritz values */
119   PetscScalar *SR, *work;
120   PetscReal   *modul;
121   PetscBLASInt bn, bN, lwork, idummy;
122   PetscScalar *t, sdummy = 0;
123   Mat          A;
124 
125   PetscFunctionBegin;
126   /* Express sizes in PetscBLASInt for LAPACK routines*/
127   PetscCall(PetscBLASIntCast(gmres->fullcycle ? gmres->max_k : gmres->it + 1, &bn)); /* size of the Hessenberg matrix */
128   PetscCall(PetscBLASIntCast(gmres->max_k + 1, &bN));                                /* LDA of the Hessenberg matrix */
129   PetscCall(PetscBLASIntCast(gmres->max_k + 1, &idummy));
130   PetscCall(PetscBLASIntCast(5 * (gmres->max_k + 1) * (gmres->max_k + 1), &lwork));
131 
132   /* NbrRitz: number of (Harmonic) Ritz pairs to extract */
133   NbrRitz = PetscMin(*nrit, bn);
134   PetscCall(KSPGetOperators(ksp, &A, NULL));
135   PetscCall(MatGetSize(A, &n, NULL));
136   NbrRitz = PetscMin(NbrRitz, n);
137 
138   PetscCall(PetscMalloc4(bN * bN, &H, bn * bn, &Q, bn, &wr, bn, &wi));
139 
140   /* copy H matrix to work space */
141   PetscCall(PetscArraycpy(H, gmres->fullcycle ? gmres->hes_ritz : gmres->hes_origin, bN * bN));
142 
143   /* Modify H to compute Harmonic Ritz pairs H = H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
144   if (!ritz) {
145     /* Transpose the Hessenberg matrix => Ht */
146     PetscCall(PetscMalloc1(bn * bn, &Ht));
147     for (i = 0; i < bn; i++) {
148       for (j = 0; j < bn; j++) Ht[i * bn + j] = PetscConj(H[j * bN + i]);
149     }
150     /* Solve the system H^T*t = h^2_{m+1,m}e_m */
151     PetscCall(PetscCalloc1(bn, &t));
152     /* t = h^2_{m+1,m}e_m */
153     if (gmres->fullcycle) t[bn - 1] = PetscSqr(gmres->hes_ritz[(bn - 1) * bN + bn]);
154     else t[bn - 1] = PetscSqr(gmres->hes_origin[(bn - 1) * bN + bn]);
155 
156     /* Call the LAPACK routine dgesv to compute t = H^{-T}*t */
157     {
158       PetscBLASInt  info;
159       PetscBLASInt  nrhs = 1;
160       PetscBLASInt *ipiv;
161       PetscCall(PetscMalloc1(bn, &ipiv));
162       PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&bn, &nrhs, Ht, &bn, ipiv, t, &bn, &info));
163       PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Error while calling the Lapack routine DGESV");
164       PetscCall(PetscFree(ipiv));
165       PetscCall(PetscFree(Ht));
166     }
167     /* Form H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
168     for (i = 0; i < bn; i++) H[(bn - 1) * bn + i] += t[i];
169     PetscCall(PetscFree(t));
170   }
171 
172   /*
173     Compute (Harmonic) Ritz pairs;
174     For a real Ritz eigenvector at wr(j)  Q(:,j) columns contain the real right eigenvector
175     For a complex Ritz pair of eigenvectors at wr(j), wi(j), wr(j+1), and wi(j+1), Q(:,j) + i Q(:,j+1) and Q(:,j) - i Q(:,j+1) are the two eigenvectors
176   */
177   {
178     PetscBLASInt info;
179 #if defined(PETSC_USE_COMPLEX)
180     PetscReal *rwork = NULL;
181 #endif
182     PetscCall(PetscMalloc1(lwork, &work));
183     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
184 #if !defined(PETSC_USE_COMPLEX)
185     PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "V", &bn, H, &bN, wr, wi, &sdummy, &idummy, Q, &bn, work, &lwork, &info));
186 #else
187     PetscCall(PetscMalloc1(2 * n, &rwork));
188     PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "V", &bn, H, &bN, wr, &sdummy, &idummy, Q, &bn, work, &lwork, rwork, &info));
189     PetscCall(PetscFree(rwork));
190 #endif
191     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine");
192     PetscCall(PetscFPTrapPop());
193     PetscCall(PetscFree(work));
194   }
195   /* sort the (Harmonic) Ritz values */
196   PetscCall(PetscMalloc2(bn, &modul, bn, &perm));
197 #if defined(PETSC_USE_COMPLEX)
198   for (i = 0; i < bn; i++) modul[i] = PetscAbsScalar(wr[i]);
199 #else
200   for (i = 0; i < bn; i++) modul[i] = PetscSqrtReal(wr[i] * wr[i] + wi[i] * wi[i]);
201 #endif
202   for (i = 0; i < bn; i++) perm[i] = i;
203   PetscCall(PetscSortRealWithPermutation(bn, modul, perm));
204 
205 #if defined(PETSC_USE_COMPLEX)
206   /* sort extracted (Harmonic) Ritz pairs */
207   nb = NbrRitz;
208   PetscCall(PetscMalloc1(nb * bn, &SR));
209   for (i = 0; i < nb; i++) {
210     if (small) {
211       tetar[i] = PetscRealPart(wr[perm[i]]);
212       tetai[i] = PetscImaginaryPart(wr[perm[i]]);
213       PetscCall(PetscArraycpy(&SR[i * bn], &(Q[perm[i] * bn]), bn));
214     } else {
215       tetar[i] = PetscRealPart(wr[perm[bn - nb + i]]);
216       tetai[i] = PetscImaginaryPart(wr[perm[bn - nb + i]]);
217       PetscCall(PetscArraycpy(&SR[i * bn], &(Q[perm[bn - nb + i] * bn]), bn)); /* permute columns of Q */
218     }
219   }
220 #else
221   /* count the number of extracted (Harmonic) Ritz pairs (with complex conjugates) */
222   if (small) {
223     while (nb < NbrRitz) {
224       if (!wi[perm[nb]]) nb += 1;
225       else {
226         if (nb < NbrRitz - 1) nb += 2;
227         else break;
228       }
229     }
230     PetscCall(PetscMalloc1(nb * bn, &SR));
231     for (i = 0; i < nb; i++) {
232       tetar[i] = wr[perm[i]];
233       tetai[i] = wi[perm[i]];
234       PetscCall(PetscArraycpy(&SR[i * bn], &(Q[perm[i] * bn]), bn));
235     }
236   } else {
237     while (nb < NbrRitz) {
238       if (wi[perm[bn - nb - 1]] == 0) nb += 1;
239       else {
240         if (nb < NbrRitz - 1) nb += 2;
241         else break;
242       }
243     }
244     PetscCall(PetscMalloc1(nb * bn, &SR)); /* bn rows, nb columns */
245     for (i = 0; i < nb; i++) {
246       tetar[i] = wr[perm[bn - nb + i]];
247       tetai[i] = wi[perm[bn - nb + i]];
248       PetscCall(PetscArraycpy(&SR[i * bn], &(Q[perm[bn - nb + i] * bn]), bn)); /* permute columns of Q */
249     }
250   }
251 #endif
252   PetscCall(PetscFree2(modul, perm));
253   PetscCall(PetscFree4(H, Q, wr, wi));
254 
255   /* Form the (Harmonic) Ritz vectors S = SR*V, columns of VV correspond to the basis of the Krylov subspace */
256   for (j = 0; j < nb; j++) {
257     PetscCall(VecZeroEntries(S[j]));
258     PetscCall(VecMAXPY(S[j], bn, &SR[j * bn], gmres->fullcycle ? gmres->vecb : &VEC_VV(0)));
259   }
260 
261   PetscCall(PetscFree(SR));
262   *nrit = nb;
263   PetscFunctionReturn(0);
264 }
265