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