Actual source code: gmreig.c


  2: #include <../src/ksp/ksp/impls/gmres/gmresimpl.h>
  3: #include <petscblaslapack.h>

  5: PetscErrorCode KSPComputeExtremeSingularValues_GMRES(KSP ksp, PetscReal *emax, PetscReal *emin)
  6: {
  7:   KSP_GMRES   *gmres = (KSP_GMRES *)ksp->data;
  8:   PetscInt     n = gmres->it + 1, i, N = gmres->max_k + 2;
  9:   PetscBLASInt bn, bN, lwork, idummy, lierr;
 10:   PetscScalar *R = gmres->Rsvd, *work = R + N * N, sdummy = 0;
 11:   PetscReal   *realpart = gmres->Dsvd;

 13:   PetscFunctionBegin;
 14:   PetscCall(PetscBLASIntCast(n, &bn));
 15:   PetscCall(PetscBLASIntCast(N, &bN));
 16:   PetscCall(PetscBLASIntCast(5 * N, &lwork));
 17:   PetscCall(PetscBLASIntCast(N, &idummy));
 18:   if (n <= 0) {
 19:     *emax = *emin = 1.0;
 20:     PetscFunctionReturn(PETSC_SUCCESS);
 21:   }
 22:   /* copy R matrix to work space */
 23:   PetscCall(PetscArraycpy(R, gmres->hh_origin, (gmres->max_k + 2) * (gmres->max_k + 1)));

 25:   /* zero below diagonal garbage */
 26:   for (i = 0; i < n; i++) R[i * N + i + 1] = 0.0;

 28:   /* compute Singular Values */
 29:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 30: #if !defined(PETSC_USE_COMPLEX)
 31:   PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("N", "N", &bn, &bn, R, &bN, realpart, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &lierr));
 32: #else
 33:   PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("N", "N", &bn, &bn, R, &bN, realpart, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, realpart + N, &lierr));
 34: #endif
 35:   PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in SVD Lapack routine %d", (int)lierr);
 36:   PetscCall(PetscFPTrapPop());

 38:   *emin = realpart[n - 1];
 39:   *emax = realpart[0];
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: PetscErrorCode KSPComputeEigenvalues_GMRES(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
 44: {
 45: #if !defined(PETSC_USE_COMPLEX)
 46:   KSP_GMRES   *gmres = (KSP_GMRES *)ksp->data;
 47:   PetscInt     n = gmres->it + 1, N = gmres->max_k + 1, i, *perm;
 48:   PetscBLASInt bn, bN, lwork, idummy, lierr = -1;
 49:   PetscScalar *R = gmres->Rsvd, *work = R + N * N;
 50:   PetscScalar *realpart = gmres->Dsvd, *imagpart = realpart + N, sdummy = 0;

 52:   PetscFunctionBegin;
 53:   PetscCall(PetscBLASIntCast(n, &bn));
 54:   PetscCall(PetscBLASIntCast(N, &bN));
 55:   PetscCall(PetscBLASIntCast(5 * N, &lwork));
 56:   PetscCall(PetscBLASIntCast(N, &idummy));
 57:   PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
 58:   *neig = n;

 60:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);

 62:   /* copy R matrix to work space */
 63:   PetscCall(PetscArraycpy(R, gmres->hes_origin, N * N));

 65:   /* compute eigenvalues */
 66:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 67:   PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "N", &bn, R, &bN, realpart, imagpart, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &lierr));
 68:   PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int)lierr);
 69:   PetscCall(PetscFPTrapPop());
 70:   PetscCall(PetscMalloc1(n, &perm));
 71:   for (i = 0; i < n; i++) perm[i] = i;
 72:   PetscCall(PetscSortRealWithPermutation(n, realpart, perm));
 73:   for (i = 0; i < n; i++) {
 74:     r[i] = realpart[perm[i]];
 75:     c[i] = imagpart[perm[i]];
 76:   }
 77:   PetscCall(PetscFree(perm));
 78: #else
 79:   KSP_GMRES   *gmres = (KSP_GMRES *)ksp->data;
 80:   PetscInt     n = gmres->it + 1, N = gmres->max_k + 1, i, *perm;
 81:   PetscScalar *R = gmres->Rsvd, *work = R + N * N, *eigs = work + 5 * N, sdummy;
 82:   PetscBLASInt bn, bN, lwork, idummy, lierr = -1;

 84:   PetscFunctionBegin;
 85:   PetscCall(PetscBLASIntCast(n, &bn));
 86:   PetscCall(PetscBLASIntCast(N, &bN));
 87:   PetscCall(PetscBLASIntCast(5 * N, &lwork));
 88:   PetscCall(PetscBLASIntCast(N, &idummy));
 89:   PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
 90:   *neig = n;

 92:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);

 94:   /* copy R matrix to work space */
 95:   PetscCall(PetscArraycpy(R, gmres->hes_origin, N * N));

 97:   /* compute eigenvalues */
 98:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 99:   PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "N", &bn, R, &bN, eigs, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, gmres->Dsvd, &lierr));
100:   PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine");
101:   PetscCall(PetscFPTrapPop());
102:   PetscCall(PetscMalloc1(n, &perm));
103:   for (i = 0; i < n; i++) perm[i] = i;
104:   for (i = 0; i < n; i++) r[i] = PetscRealPart(eigs[i]);
105:   PetscCall(PetscSortRealWithPermutation(n, r, perm));
106:   for (i = 0; i < n; i++) {
107:     r[i] = PetscRealPart(eigs[perm[i]]);
108:     c[i] = PetscImaginaryPart(eigs[perm[i]]);
109:   }
110:   PetscCall(PetscFree(perm));
111: #endif
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }

115: PetscErrorCode KSPComputeRitz_GMRES(KSP ksp, PetscBool ritz, PetscBool small, PetscInt *nrit, Vec S[], PetscReal *tetar, PetscReal *tetai)
116: {
117:   KSP_GMRES   *gmres = (KSP_GMRES *)ksp->data;
118:   PetscInt     NbrRitz, nb = 0, n;
119:   PetscInt     i, j, *perm;
120:   PetscScalar *H, *Q, *Ht; /* H Hessenberg matrix; Q matrix of eigenvectors of H */
121:   PetscScalar *wr, *wi;    /* Real and imaginary part of the Ritz values */
122:   PetscScalar *SR, *work;
123:   PetscReal   *modul;
124:   PetscBLASInt bn, bN, lwork, idummy;
125:   PetscScalar *t, sdummy = 0;
126:   Mat          A;

128:   PetscFunctionBegin;
129:   /* Express sizes in PetscBLASInt for LAPACK routines*/
130:   PetscCall(PetscBLASIntCast(gmres->fullcycle ? gmres->max_k : gmres->it + 1, &bn)); /* size of the Hessenberg matrix */
131:   PetscCall(PetscBLASIntCast(gmres->max_k + 1, &bN));                                /* LDA of the Hessenberg matrix */
132:   PetscCall(PetscBLASIntCast(gmres->max_k + 1, &idummy));
133:   PetscCall(PetscBLASIntCast(5 * (gmres->max_k + 1) * (gmres->max_k + 1), &lwork));

135:   /* NbrRitz: number of (Harmonic) Ritz pairs to extract */
136:   NbrRitz = PetscMin(*nrit, bn);
137:   PetscCall(KSPGetOperators(ksp, &A, NULL));
138:   PetscCall(MatGetSize(A, &n, NULL));
139:   NbrRitz = PetscMin(NbrRitz, n);

141:   PetscCall(PetscMalloc4(bN * bN, &H, bn * bn, &Q, bn, &wr, bn, &wi));

143:   /* copy H matrix to work space */
144:   PetscCall(PetscArraycpy(H, gmres->fullcycle ? gmres->hes_ritz : gmres->hes_origin, bN * bN));

146:   /* Modify H to compute Harmonic Ritz pairs H = H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
147:   if (!ritz) {
148:     /* Transpose the Hessenberg matrix => Ht */
149:     PetscCall(PetscMalloc1(bn * bn, &Ht));
150:     for (i = 0; i < bn; i++) {
151:       for (j = 0; j < bn; j++) Ht[i * bn + j] = PetscConj(H[j * bN + i]);
152:     }
153:     /* Solve the system H^T*t = h^2_{m+1,m}e_m */
154:     PetscCall(PetscCalloc1(bn, &t));
155:     /* t = h^2_{m+1,m}e_m */
156:     if (gmres->fullcycle) t[bn - 1] = PetscSqr(gmres->hes_ritz[(bn - 1) * bN + bn]);
157:     else t[bn - 1] = PetscSqr(gmres->hes_origin[(bn - 1) * bN + bn]);

159:     /* Call the LAPACK routine dgesv to compute t = H^{-T}*t */
160:     {
161:       PetscBLASInt  info;
162:       PetscBLASInt  nrhs = 1;
163:       PetscBLASInt *ipiv;
164:       PetscCall(PetscMalloc1(bn, &ipiv));
165:       PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&bn, &nrhs, Ht, &bn, ipiv, t, &bn, &info));
166:       PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Error while calling the Lapack routine DGESV");
167:       PetscCall(PetscFree(ipiv));
168:       PetscCall(PetscFree(Ht));
169:     }
170:     /* Form H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
171:     for (i = 0; i < bn; i++) H[(bn - 1) * bn + i] += t[i];
172:     PetscCall(PetscFree(t));
173:   }

175:   /*
176:     Compute (Harmonic) Ritz pairs;
177:     For a real Ritz eigenvector at wr(j)  Q(:,j) columns contain the real right eigenvector
178:     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
179:   */
180:   {
181:     PetscBLASInt info;
182: #if defined(PETSC_USE_COMPLEX)
183:     PetscReal *rwork = NULL;
184: #endif
185:     PetscCall(PetscMalloc1(lwork, &work));
186:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
187: #if !defined(PETSC_USE_COMPLEX)
188:     PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "V", &bn, H, &bN, wr, wi, &sdummy, &idummy, Q, &bn, work, &lwork, &info));
189: #else
190:     PetscCall(PetscMalloc1(2 * n, &rwork));
191:     PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "V", &bn, H, &bN, wr, &sdummy, &idummy, Q, &bn, work, &lwork, rwork, &info));
192:     PetscCall(PetscFree(rwork));
193: #endif
194:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine");
195:     PetscCall(PetscFPTrapPop());
196:     PetscCall(PetscFree(work));
197:   }
198:   /* sort the (Harmonic) Ritz values */
199:   PetscCall(PetscMalloc2(bn, &modul, bn, &perm));
200: #if defined(PETSC_USE_COMPLEX)
201:   for (i = 0; i < bn; i++) modul[i] = PetscAbsScalar(wr[i]);
202: #else
203:   for (i = 0; i < bn; i++) modul[i] = PetscSqrtReal(wr[i] * wr[i] + wi[i] * wi[i]);
204: #endif
205:   for (i = 0; i < bn; i++) perm[i] = i;
206:   PetscCall(PetscSortRealWithPermutation(bn, modul, perm));

208: #if defined(PETSC_USE_COMPLEX)
209:   /* sort extracted (Harmonic) Ritz pairs */
210:   nb = NbrRitz;
211:   PetscCall(PetscMalloc1(nb * bn, &SR));
212:   for (i = 0; i < nb; i++) {
213:     if (small) {
214:       tetar[i] = PetscRealPart(wr[perm[i]]);
215:       tetai[i] = PetscImaginaryPart(wr[perm[i]]);
216:       PetscCall(PetscArraycpy(&SR[i * bn], &(Q[perm[i] * bn]), bn));
217:     } else {
218:       tetar[i] = PetscRealPart(wr[perm[bn - nb + i]]);
219:       tetai[i] = PetscImaginaryPart(wr[perm[bn - nb + i]]);
220:       PetscCall(PetscArraycpy(&SR[i * bn], &(Q[perm[bn - nb + i] * bn]), bn)); /* permute columns of Q */
221:     }
222:   }
223: #else
224:   /* count the number of extracted (Harmonic) Ritz pairs (with complex conjugates) */
225:   if (small) {
226:     while (nb < NbrRitz) {
227:       if (!wi[perm[nb]]) nb += 1;
228:       else {
229:         if (nb < NbrRitz - 1) nb += 2;
230:         else break;
231:       }
232:     }
233:     PetscCall(PetscMalloc1(nb * bn, &SR));
234:     for (i = 0; i < nb; i++) {
235:       tetar[i] = wr[perm[i]];
236:       tetai[i] = wi[perm[i]];
237:       PetscCall(PetscArraycpy(&SR[i * bn], &(Q[perm[i] * bn]), bn));
238:     }
239:   } else {
240:     while (nb < NbrRitz) {
241:       if (wi[perm[bn - nb - 1]] == 0) nb += 1;
242:       else {
243:         if (nb < NbrRitz - 1) nb += 2;
244:         else break;
245:       }
246:     }
247:     PetscCall(PetscMalloc1(nb * bn, &SR)); /* bn rows, nb columns */
248:     for (i = 0; i < nb; i++) {
249:       tetar[i] = wr[perm[bn - nb + i]];
250:       tetai[i] = wi[perm[bn - nb + i]];
251:       PetscCall(PetscArraycpy(&SR[i * bn], &(Q[perm[bn - nb + i] * bn]), bn)); /* permute columns of Q */
252:     }
253:   }
254: #endif
255:   PetscCall(PetscFree2(modul, perm));
256:   PetscCall(PetscFree4(H, Q, wr, wi));

258:   /* Form the (Harmonic) Ritz vectors S = SR*V, columns of VV correspond to the basis of the Krylov subspace */
259:   for (j = 0; j < nb; j++) {
260:     PetscCall(VecZeroEntries(S[j]));
261:     PetscCall(VecMAXPY(S[j], bn, &SR[j * bn], gmres->fullcycle ? gmres->vecb : &VEC_VV(0)));
262:   }

264:   PetscCall(PetscFree(SR));
265:   *nrit = nb;
266:   PetscFunctionReturn(PETSC_SUCCESS);
267: }