Actual source code: agmresdeflation.c

      1: /*
      2:   This file computes data for the deflated restarting in the Newton-basis GMRES.
      3:   At each restart (or at each detected stagnation in the adaptive strategy), a basis of an
      4:   (approximated)invariant subspace corresponding to the smallest eigenvalues is extracted from the Krylov subspace.
      5:   It is then used to augment the Newton basis.
      6: */

      8: #include <../src/ksp/ksp/impls/gmres/agmres/agmresimpl.h>

     10: /* Quicksort algorithm to sort  the eigenvalues in increasing orders
     11:    val_r - real part of eigenvalues, unchanged on exit.
     12:    val_i - Imaginary part of eigenvalues unchanged on exit.
     13:    size - Number of eigenvalues (with complex conjugates)
     14:    perm - contains on exit the permutation vector to reorder the vectors val_r and val_i.
     15: */
     16: #define DEPTH 500
     17: static PetscErrorCode KSPAGMRESQuickSort(PetscScalar *val_r, PetscScalar *val_i, PetscInt size, PetscInt *perm)
     18: {
     19:   PetscInt    deb[DEPTH], fin[DEPTH];
     20:   PetscInt    ipivot;
     21:   PetscScalar pivot_r, pivot_i;
     22:   PetscInt    i, L, R, j;
     23:   PetscScalar abs_pivot;
     24:   PetscScalar abs_val;

     26:   PetscFunctionBegin;
     27:   /* initialize perm vector */
     28:   for (j = 0; j < size; j++) perm[j] = j;

     30:   deb[0] = 0;
     31:   fin[0] = size;
     32:   i      = 0;
     33:   while (i >= 0) {
     34:     L = deb[i];
     35:     R = fin[i] - 1;
     36:     if (L < R) {
     37:       pivot_r   = val_r[L];
     38:       pivot_i   = val_i[L];
     39:       abs_pivot = PetscSqrtReal(pivot_r * pivot_r + pivot_i * pivot_i);
     40:       ipivot    = perm[L];
     41:       PetscCheck(i != DEPTH - 1, PETSC_COMM_SELF, PETSC_ERR_MEM, "Could cause stack overflow: Try to increase the value of DEPTH ");
     42:       while (L < R) {
     43:         abs_val = PetscSqrtReal(val_r[R] * val_r[R] + val_i[R] * val_i[R]);
     44:         while (abs_val >= abs_pivot && L < R) {
     45:           R--;
     46:           abs_val = PetscSqrtReal(val_r[R] * val_r[R] + val_i[R] * val_i[R]);
     47:         }
     48:         if (L < R) {
     49:           val_r[L] = val_r[R];
     50:           val_i[L] = val_i[R];
     51:           perm[L]  = perm[R];
     52:           L += 1;
     53:         }
     54:         abs_val = PetscSqrtReal(val_r[L] * val_r[L] + val_i[L] * val_i[L]);
     55:         while (abs_val <= abs_pivot && L < R) {
     56:           L++;
     57:           abs_val = PetscSqrtReal(val_r[L] * val_r[L] + val_i[L] * val_i[L]);
     58:         }
     59:         if (L < R) {
     60:           val_r[R] = val_r[L];
     61:           val_i[R] = val_i[L];
     62:           perm[R]  = perm[L];
     63:           R -= 1;
     64:         }
     65:       }
     66:       val_r[L]   = pivot_r;
     67:       val_i[L]   = pivot_i;
     68:       perm[L]    = ipivot;
     69:       deb[i + 1] = L + 1;
     70:       fin[i + 1] = fin[i];
     71:       fin[i]     = L;
     72:       i += 1;
     73:       PetscCheck(i != DEPTH - 1, PETSC_COMM_SELF, PETSC_ERR_MEM, "Could cause stack overflow: Try to increase the value of DEPTH ");
     74:     } else i--;
     75:   }
     76:   PetscFunctionReturn(PETSC_SUCCESS);
     77: }

     79: /*
     80:  Compute the Schur vectors from the generalized eigenvalue problem A.u =\lambda.B.u
     81:  KspSize -  rank of the matrices A and B, size of the current Krylov basis
     82:  A - Left matrix
     83:  B - Right matrix
     84:  ldA - first dimension of A as declared  in the calling program
     85:  ldB - first dimension of B as declared  in the calling program
     86:  IsReduced - specifies if the matrices are already in the reduced form,
     87:  i.e A is a Hessenberg matrix and B is upper triangular.
     88:  Sr - on exit, the extracted Schur vectors corresponding
     89:  the smallest eigenvalues (with complex conjugates)
     90:  CurNeig - Number of extracted eigenvalues
     91: */
     92: static PetscErrorCode KSPAGMRESSchurForm(KSP ksp, PetscBLASInt KspSize, PetscScalar *A, PetscBLASInt ldA, PetscScalar *B, PetscBLASInt ldB, PetscBool IsReduced, PetscScalar *Sr, PetscInt *CurNeig)
     93: {
     94:   KSP_AGMRES   *agmres = (KSP_AGMRES *)ksp->data;
     95:   PetscInt      max_k  = agmres->max_k;
     96:   PetscBLASInt  r;
     97:   PetscInt      neig   = agmres->neig;
     98:   PetscScalar  *wr     = agmres->wr;
     99:   PetscScalar  *wi     = agmres->wi;
    100:   PetscScalar  *beta   = agmres->beta;
    101:   PetscScalar  *Q      = agmres->Q;
    102:   PetscScalar  *Z      = agmres->Z;
    103:   PetscScalar  *work   = agmres->work;
    104:   PetscBLASInt *select = agmres->select;
    105:   PetscInt     *perm   = agmres->perm;
    106:   PetscBLASInt  sdim   = 0;
    107:   PetscInt      i, j;
    108:   PetscBLASInt  info;
    109:   PetscBLASInt *iwork = agmres->iwork;
    110:   PetscBLASInt  N     = (PetscBLASInt)MAXKSPSIZE;
    111:   PetscBLASInt  lwork, liwork;
    112:   PetscBLASInt  ilo;
    113:   PetscBLASInt  ijob, wantQ, wantZ;
    114:   PetscScalar   Dif[2];

    116:   PetscFunctionBegin;
    117:   ijob  = 2;
    118:   wantQ = 1;
    119:   wantZ = 1;
    120:   PetscCall(PetscBLASIntCast(PetscMax(8 * N + 16, 4 * neig * (N - neig)), &lwork));
    121:   PetscCall(PetscBLASIntCast(2 * N * neig, &liwork));
    122:   ilo = 1;

    124:   /* Compute the Schur form */
    125:   if (IsReduced) { /* The eigenvalue problem is already in reduced form, meaning that A is upper Hessenberg and B is triangular */
    126:     PetscCallBLAS("LAPACKhgeqz", LAPACKhgeqz_("S", "I", "I", &KspSize, &ilo, &KspSize, A, &ldA, B, &ldB, wr, wi, beta, Q, &N, Z, &N, work, &lwork, &info));
    127:     PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Error while calling LAPACK routine xhgeqz_");
    128:   } else {
    129:     PetscCallBLAS("LAPACKgges", LAPACKgges_("V", "V", "N", NULL, &KspSize, A, &ldA, B, &ldB, &sdim, wr, wi, beta, Q, &N, Z, &N, work, &lwork, NULL, &info));
    130:     PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Error while calling LAPACK routine xgges_");
    131:   }

    133:   /* We should avoid computing these ratio...  */
    134:   for (i = 0; i < KspSize; i++) {
    135:     if (beta[i] != 0.0) {
    136:       wr[i] /= beta[i];
    137:       wi[i] /= beta[i];
    138:     }
    139:   }

    141:   /* Sort the eigenvalues to extract the smallest ones */
    142:   PetscCall(KSPAGMRESQuickSort(wr, wi, KspSize, perm));

    144:   /* Count the number of extracted eigenvalues (with complex conjugates) */
    145:   r = 0;
    146:   while (r < neig) {
    147:     if (wi[r] != 0) r += 2;
    148:     else r += 1;
    149:   }
    150:   /* Reorder the Schur decomposition so that the cluster of smallest/largest eigenvalues appears in the leading diagonal blocks of A (and B)*/
    151:   PetscCall(PetscArrayzero(select, N));
    152:   if (!agmres->GreatestEig) {
    153:     for (j = 0; j < r; j++) select[perm[j]] = 1;
    154:   } else {
    155:     for (j = 0; j < r; j++) select[perm[KspSize - j - 1]] = 1;
    156:   }
    157:   PetscCallBLAS("LAPACKtgsen", LAPACKtgsen_(&ijob, &wantQ, &wantZ, select, &KspSize, A, &ldA, B, &ldB, wr, wi, beta, Q, &N, Z, &N, &r, NULL, NULL, &Dif[0], work, &lwork, iwork, &liwork, &info));
    158:   PetscCheck(info != 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "UNABLE TO REORDER THE EIGENVALUES WITH THE LAPACK ROUTINE : ILL-CONDITIONED PROBLEM");
    159:   /* Extract the Schur vectors associated to the r smallest eigenvalues */
    160:   PetscCall(PetscArrayzero(Sr, (N + 1) * r));
    161:   for (j = 0; j < r; j++) {
    162:     for (i = 0; i < KspSize; i++) Sr[j * (N + 1) + i] = Z[j * N + i];
    163:   }

    165:   /* Broadcast Sr to all other processes to have consistent data;
    166:    * FIXME should investigate how to get unique Schur vectors (unique QR factorization, probably the sign of rotations) */
    167:   PetscCallMPI(MPI_Bcast(Sr, (N + 1) * r, MPIU_SCALAR, agmres->First, PetscObjectComm((PetscObject)ksp)));
    168:   /* Update the Shift values for the Newton basis. This is surely necessary when applying the DeflationPrecond */
    169:   if (agmres->DeflPrecond) PetscCall(KSPAGMRESLejaOrdering(wr, wi, agmres->Rshift, agmres->Ishift, max_k));
    170:   *CurNeig = r; /* Number of extracted eigenvalues */
    171:   PetscFunctionReturn(PETSC_SUCCESS);
    172: }

    174: /*
    175:   Forms the matrices for the generalized eigenvalue problem,
    176:   it then compute the Schur vectors needed to augment the Newton basis.
    177: */
    178: PetscErrorCode KSPAGMRESComputeDeflationData(KSP ksp)
    179: {
    180:   KSP_AGMRES  *agmres  = (KSP_AGMRES *)ksp->data;
    181:   Vec         *U       = agmres->U;
    182:   Vec         *TmpU    = agmres->TmpU;
    183:   PetscScalar *MatEigL = agmres->MatEigL;
    184:   PetscScalar *MatEigR = agmres->MatEigR;
    185:   PetscScalar *Sr      = agmres->Sr;
    186:   PetscScalar  alpha, beta;
    187:   PetscInt     i, j;
    188:   PetscInt     max_k = agmres->max_k; /* size of the non - augmented subspace */
    189:   PetscInt     CurNeig;               /* Current number of extracted eigenvalues */
    190:   PetscInt     N = MAXKSPSIZE;
    191:   PetscBLASInt bN, iKspSize;
    192:   PetscInt     lC      = N + 1;
    193:   PetscInt     KspSize = KSPSIZE;
    194:   PetscBLASInt blC, bKspSize;
    195:   PetscInt     PrevNeig = agmres->r;

    197:   PetscFunctionBegin;
    198:   PetscCall(PetscLogEventBegin(KSP_AGMRESComputeDeflationData, ksp, 0, 0, 0));
    199:   if (agmres->neig <= 1) PetscFunctionReturn(PETSC_SUCCESS);
    200:   /* Explicitly form MatEigL = H^T*H, It can also be formed as H^T+h_{N+1,N}H^-1e^T */
    201:   alpha = 1.0;
    202:   beta  = 0.0;
    203:   PetscCall(PetscBLASIntCast(KspSize, &bKspSize));
    204:   PetscCall(PetscBLASIntCast(lC, &blC));
    205:   PetscCall(PetscBLASIntCast(N, &bN));
    206:   PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &bKspSize, &bKspSize, &blC, &alpha, agmres->hes_origin, &blC, agmres->hes_origin, &blC, &beta, MatEigL, &bN));
    207:   if (!agmres->ritz) {
    208:     /* Form TmpU = V*H where V is the Newton basis orthogonalized  with roddec*/
    209:     for (j = 0; j < KspSize; j++) {
    210:       /* Apply the elementary reflectors (stored in Qloc) on H */
    211:       PetscCall(KSPAGMRESRodvec(ksp, KspSize + 1, &agmres->hes_origin[j * lC], TmpU[j]));
    212:     }
    213:     /* Now form MatEigR = TmpU^T*W where W is [VEC_V(1:max_k); U] */
    214:     for (j = 0; j < max_k; j++) PetscCall(VecMDot(VEC_V(j), KspSize, TmpU, &MatEigR[j * N]));
    215:     for (j = max_k; j < KspSize; j++) PetscCall(VecMDot(U[j - max_k], KspSize, TmpU, &MatEigR[j * N]));
    216:   } else { /* Form H^T */
    217:     for (j = 0; j < N; j++) {
    218:       for (i = 0; i < N; i++) MatEigR[j * N + i] = agmres->hes_origin[i * lC + j];
    219:     }
    220:   }
    221:   /* Obtain the Schur form of  the generalized eigenvalue problem MatEigL*y = \lambda*MatEigR*y */
    222:   PetscCall(PetscBLASIntCast(KspSize, &iKspSize));
    223:   if (agmres->DeflPrecond) {
    224:     PetscCall(KSPAGMRESSchurForm(ksp, iKspSize, agmres->hes_origin, blC, agmres->Rloc, blC, PETSC_TRUE, Sr, &CurNeig));
    225:   } else {
    226:     PetscCall(KSPAGMRESSchurForm(ksp, iKspSize, MatEigL, bN, MatEigR, bN, PETSC_FALSE, Sr, &CurNeig));
    227:   }

    229:   if (agmres->DeflPrecond) { /* Switch to DGMRES to improve the basis of the invariant subspace associated to the deflation */
    230:     agmres->HasSchur = PETSC_TRUE;
    231:     PetscCall(KSPDGMRESComputeDeflationData(ksp, &CurNeig));
    232:     PetscFunctionReturn(PETSC_SUCCESS);
    233:   }
    234:   /* Form the Schur vectors in the entire subspace: U = W * Sr where W = [VEC_V(1:max_k); U]*/
    235:   for (j = 0; j < PrevNeig; j++) { /* First, copy U to a temporary place */
    236:     PetscCall(VecCopy(U[j], TmpU[j]));
    237:   }

    239:   for (j = 0; j < CurNeig; j++) {
    240:     PetscCall(VecMAXPBY(U[j], max_k, &Sr[j * (N + 1)], 0, &VEC_V(0)));
    241:     PetscCall(VecMAXPY(U[j], PrevNeig, &Sr[j * (N + 1) + max_k], TmpU));
    242:   }
    243:   agmres->r = CurNeig;
    244:   PetscCall(PetscLogEventEnd(KSP_AGMRESComputeDeflationData, ksp, 0, 0, 0));
    245:   PetscFunctionReturn(PETSC_SUCCESS);
    246: }