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, lierr;
      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 !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());

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    260:   PetscCall(PetscFree(SR));
    261:   *nrit = nb;
    262:   PetscFunctionReturn(PETSC_SUCCESS);
    263: }