Actual source code: agmresorthog.c

      1: #define PETSCKSP_DLL

      3: #include <../src/ksp/ksp/impls/gmres/agmres/agmresimpl.h>
      4: /*
      5:   This file implements the RODDEC algorithm : its purpose is to orthogonalize a set of vectors distributed across several processes.
      6:   These processes are organized in a virtual ring.

      8:   References : [1] Sidje, Roger B. Alternatives for parallel Krylov subspace basis computation. Numer. Linear Algebra Appl. 4 (1997), no. 4, 305-331

     10:   initial author R. B. SIDJE,
     11:   modified : G.-A Atenekeng-Kahou, 2008
     12:   modified : D. NUENTSA WAKAM, 2011

     14: */

     16: /*
     17: * Take the processes that own the vectors and organize them on a virtual ring.
     18: */
     19: static PetscErrorCode KSPAGMRESRoddecGivens(PetscReal *, PetscReal *, PetscReal *, PetscInt);

     21: PetscErrorCode KSPAGMRESRoddecInitNeighboor(KSP ksp)
     22: {
     23:   MPI_Comm    comm;
     24:   KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
     25:   PetscMPIInt First, Last, rank, size;

     27:   PetscFunctionBegin;
     28:   PetscCall(PetscObjectGetComm((PetscObject)agmres->vecs[0], &comm));
     29:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
     30:   PetscCallMPI(MPI_Comm_size(comm, &size));
     31:   PetscCallMPI(MPIU_Allreduce(&rank, &First, 1, MPI_INT, MPI_MIN, comm));
     32:   PetscCallMPI(MPIU_Allreduce(&rank, &Last, 1, MPI_INT, MPI_MAX, comm));

     34:   if ((rank != Last) && (rank == 0)) {
     35:     agmres->Ileft  = rank - 1;
     36:     agmres->Iright = rank + 1;
     37:   } else {
     38:     if (rank == Last) {
     39:       agmres->Ileft  = rank - 1;
     40:       agmres->Iright = First;
     41:     } else {
     42:       {
     43:         agmres->Ileft  = Last;
     44:         agmres->Iright = rank + 1;
     45:       }
     46:     }
     47:   }
     48:   agmres->rank  = rank;
     49:   agmres->size  = size;
     50:   agmres->First = First;
     51:   agmres->Last  = Last;
     52:   PetscFunctionReturn(PETSC_SUCCESS);
     53: }

     55: static PetscErrorCode KSPAGMRESRoddecGivens(PetscReal *c, PetscReal *s, PetscReal *r, PetscInt make_r)
     56: {
     57:   PetscReal a, b, t;

     59:   PetscFunctionBegin;
     60:   if (make_r == 1) {
     61:     a = *c;
     62:     b = *s;
     63:     if (b == 0.e0) {
     64:       *c = 1.e0;
     65:       *s = 0.e0;
     66:     } else {
     67:       if (PetscAbsReal(b) > PetscAbsReal(a)) {
     68:         t  = -a / b;
     69:         *s = 1.e0 / PetscSqrtReal(1.e0 + t * t);
     70:         *c = (*s) * t;
     71:       } else {
     72:         t  = -b / a;
     73:         *c = 1.e0 / PetscSqrtReal(1.e0 + t * t);
     74:         *s = (*c) * t;
     75:       }
     76:     }
     77:     if (*c == 0.e0) {
     78:       *r = 1.e0;
     79:     } else {
     80:       if (PetscAbsReal(*s) < PetscAbsReal(*c)) {
     81:         *r = PetscSign(*c) * (*s) / 2.e0;
     82:       } else {
     83:         *r = PetscSign(*s) * 2.e0 / (*c);
     84:       }
     85:     }
     86:   }

     88:   if (*r == 1.e0) {
     89:     *c = 0.e0;
     90:     *s = 1.e0;
     91:   } else {
     92:     if (PetscAbsReal(*r) < 1.e0) {
     93:       *s = 2.e0 * (*r);
     94:       *c = PetscSqrtReal(1.e0 - (*s) * (*s));
     95:     } else {
     96:       *c = 2.e0 / (*r);
     97:       *s = PetscSqrtReal(1.e0 - (*c) * (*c));
     98:     }
     99:   }
    100:   PetscFunctionReturn(PETSC_SUCCESS);
    101: }

    103: /*
    104:   Compute the QR factorization of the Krylov basis vectors
    105:   Input :
    106:    - the vectors are available in agmres->vecs (alias VEC_V)
    107:    - nvec :  the number of vectors
    108:   Output :
    109:    - agmres->Qloc : product of elementary reflectors for the QR of the (local part) of the vectors.
    110:    - agmres->sgn :  Sign of the rotations
    111:    - agmres->tloc : scalar factors of the elementary reflectors.

    113: */
    114: PetscErrorCode KSPAGMRESRoddec(KSP ksp, PetscInt nvec)
    115: {
    116:   KSP_AGMRES  *agmres = (KSP_AGMRES *)ksp->data;
    117:   MPI_Comm     comm;
    118:   PetscScalar *Qloc    = agmres->Qloc;
    119:   PetscScalar *sgn     = agmres->sgn;
    120:   PetscScalar *tloc    = agmres->tloc;
    121:   PetscReal   *wbufptr = agmres->wbufptr;
    122:   PetscMPIInt  rank    = agmres->rank;
    123:   PetscMPIInt  First   = agmres->First;
    124:   PetscMPIInt  Last    = agmres->Last;
    125:   PetscBLASInt pas, len, bnloc, bpos;
    126:   PetscInt     nloc, d, i, j, k;
    127:   PetscInt     pos;
    128:   PetscReal    c, s, rho, Ajj, val, tt, old;
    129:   PetscScalar *col;
    130:   MPI_Status   status;
    131:   PetscBLASInt N = (PetscBLASInt)(MAXKSPSIZE + 1);

    133:   PetscFunctionBegin;
    134:   PetscCall(PetscObjectGetComm((PetscObject)ksp, &comm));
    135:   PetscCall(PetscLogEventBegin(KSP_AGMRESRoddec, ksp, 0, 0, 0));
    136:   PetscCall(PetscArrayzero(agmres->Rloc, N * N));
    137:   /* check input arguments */
    138:   PetscCheck(nvec >= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "The number of input vectors should be positive");
    139:   PetscCall(VecGetLocalSize(VEC_V(0), &nloc));
    140:   PetscCall(PetscBLASIntCast(nloc, &bnloc));
    141:   PetscCheck(nvec <= nloc, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_WRONG, "In QR factorization, the number of local rows should be greater or equal to the number of columns");
    142:   pas = 1;
    143:   /* Copy the vectors of the basis */
    144:   for (j = 0; j < nvec; j++) {
    145:     PetscCall(VecGetArray(VEC_V(j), &col));
    146:     PetscCallBLAS("BLAScopy", BLAScopy_(&bnloc, col, &pas, &Qloc[j * nloc], &pas));
    147:     PetscCall(VecRestoreArray(VEC_V(j), &col));
    148:   }
    149:   /* Each process performs a local QR on its own block */
    150:   for (j = 0; j < nvec; j++) {
    151:     PetscCall(PetscBLASIntCast(nloc - j, &len));
    152:     Ajj = Qloc[j * nloc + j];
    153:     PetscCallBLAS("BLASnrm2", rho = -PetscSign(Ajj) * BLASnrm2_(&len, &Qloc[j * nloc + j], &pas));
    154:     if (rho == 0.0) tloc[j] = 0.0;
    155:     else {
    156:       tloc[j] = (Ajj - rho) / rho;
    157:       len     = len - 1;
    158:       val     = 1.0 / (Ajj - rho);
    159:       PetscCallBLAS("BLASscal", BLASscal_(&len, &val, &Qloc[j * nloc + j + 1], &pas));
    160:       Qloc[j * nloc + j] = 1.0;
    161:       len                = len + 1;
    162:       for (k = j + 1; k < nvec; k++) {
    163:         PetscCallBLAS("BLASdot", tt = tloc[j] * BLASdot_(&len, &Qloc[j * nloc + j], &pas, &Qloc[k * nloc + j], &pas));
    164:         PetscCallBLAS("BLASaxpy", BLASaxpy_(&len, &tt, &Qloc[j * nloc + j], &pas, &Qloc[k * nloc + j], &pas));
    165:       }
    166:       Qloc[j * nloc + j] = rho;
    167:     }
    168:   }
    169:   /* annihilate undesirable Rloc, diagonal by diagonal*/
    170:   for (d = 0; d < nvec; d++) {
    171:     PetscCall(PetscBLASIntCast(nloc - j, &len));
    172:     if (rank == First) {
    173:       PetscCallBLAS("BLAScopy", BLAScopy_(&len, &Qloc[d * nloc + d], &bnloc, &wbufptr[d], &pas));
    174:       PetscCallMPI(MPI_Send(&wbufptr[d], len, MPIU_SCALAR, rank + 1, agmres->tag, comm));
    175:     } else {
    176:       PetscCallMPI(MPI_Recv(&wbufptr[d], len, MPIU_SCALAR, rank - 1, agmres->tag, comm, &status));
    177:       /* Elimination of Rloc(1,d)*/
    178:       c = wbufptr[d];
    179:       s = Qloc[d * nloc];
    180:       PetscCall(KSPAGMRESRoddecGivens(&c, &s, &rho, 1));
    181:       /* Apply Givens Rotation*/
    182:       for (k = d; k < nvec; k++) {
    183:         old            = wbufptr[k];
    184:         wbufptr[k]     = c * old - s * Qloc[k * nloc];
    185:         Qloc[k * nloc] = s * old + c * Qloc[k * nloc];
    186:       }
    187:       Qloc[d * nloc] = rho;
    188:       if (rank != Last) PetscCallMPI(MPI_Send(&wbufptr[d], len, MPIU_SCALAR, rank + 1, agmres->tag, comm));
    189:       /* zero-out the d-th diagonal of Rloc ...*/
    190:       for (j = d + 1; j < nvec; j++) {
    191:         /* elimination of Rloc[i][j]*/
    192:         i = j - d;
    193:         c = Qloc[j * nloc + i - 1];
    194:         s = Qloc[j * nloc + i];
    195:         PetscCall(KSPAGMRESRoddecGivens(&c, &s, &rho, 1));
    196:         for (k = j; k < nvec; k++) {
    197:           old                    = Qloc[k * nloc + i - 1];
    198:           Qloc[k * nloc + i - 1] = c * old - s * Qloc[k * nloc + i];
    199:           Qloc[k * nloc + i]     = s * old + c * Qloc[k * nloc + i];
    200:         }
    201:         Qloc[j * nloc + i] = rho;
    202:       }
    203:       if (rank == Last) {
    204:         PetscCallBLAS("BLAScopy", BLAScopy_(&len, &wbufptr[d], &pas, RLOC(d, d), &N));
    205:         for (k = d + 1; k < nvec; k++) *RLOC(k, d) = 0.0;
    206:       }
    207:     }
    208:   }

    210:   if (rank == Last) {
    211:     for (d = 0; d < nvec; d++) {
    212:       pos = nvec - d;
    213:       PetscCall(PetscBLASIntCast(pos, &bpos));
    214:       sgn[d] = PetscSign(*RLOC(d, d));
    215:       PetscCallBLAS("BLASscal", BLASscal_(&bpos, &sgn[d], RLOC(d, d), &N));
    216:     }
    217:   }
    218:   /* BroadCast Rloc to all other processes
    219:    * NWD : should not be needed
    220:    */
    221:   PetscCallMPI(MPI_Bcast(agmres->Rloc, N * N, MPIU_SCALAR, Last, comm));
    222:   PetscCall(PetscLogEventEnd(KSP_AGMRESRoddec, ksp, 0, 0, 0));
    223:   PetscFunctionReturn(PETSC_SUCCESS);
    224: }

    226: /*
    227:    Computes Out <-- Q * In where Q is the orthogonal matrix from AGMRESRoddec
    228:    Input
    229:     - Qloc, sgn, tloc, nvec (see AGMRESRoddec above)
    230:     - In : input vector (size nvec)
    231:    Output :
    232:     - Out : Petsc vector (distributed as the basis vectors)
    233: */
    234: PetscErrorCode KSPAGMRESRodvec(KSP ksp, PetscInt nvec, PetscScalar *In, Vec Out)
    235: {
    236:   KSP_AGMRES  *agmres = (KSP_AGMRES *)ksp->data;
    237:   MPI_Comm     comm;
    238:   PetscScalar *Qloc  = agmres->Qloc;
    239:   PetscScalar *sgn   = agmres->sgn;
    240:   PetscScalar *tloc  = agmres->tloc;
    241:   PetscMPIInt  rank  = agmres->rank;
    242:   PetscMPIInt  First = agmres->First, Last = agmres->Last;
    243:   PetscMPIInt  Iright = agmres->Iright, Ileft = agmres->Ileft;
    244:   PetscScalar *y, *zloc;
    245:   PetscInt     nloc, d, len, i, j;
    246:   PetscBLASInt bnvec, pas, blen;
    247:   PetscInt     dpt;
    248:   PetscReal    c, s, rho, zp, zq, yd = 0.0, tt;
    249:   MPI_Status   status;

    251:   PetscFunctionBegin;
    252:   PetscCall(PetscBLASIntCast(nvec, &bnvec));
    253:   PetscCall(PetscObjectGetComm((PetscObject)ksp, &comm));
    254:   pas = 1;
    255:   PetscCall(VecGetLocalSize(VEC_V(0), &nloc));
    256:   PetscCall(PetscMalloc1(nvec, &y));
    257:   PetscCall(PetscArraycpy(y, In, nvec));
    258:   PetscCall(VecGetArray(Out, &zloc));

    260:   if (rank == Last) {
    261:     for (i = 0; i < nvec; i++) y[i] = sgn[i] * y[i];
    262:   }
    263:   for (i = 0; i < nloc; i++) zloc[i] = 0.0;
    264:   if (agmres->size == 1) PetscCallBLAS("BLAScopy", BLAScopy_(&bnvec, y, &pas, &zloc[0], &pas));
    265:   else {
    266:     for (d = nvec - 1; d >= 0; d--) {
    267:       if (rank == First) {
    268:         PetscCallMPI(MPI_Recv(&zloc[d], 1, MPIU_SCALAR, Iright, agmres->tag, comm, &status));
    269:       } else {
    270:         for (j = nvec - 1; j >= d + 1; j--) {
    271:           i = j - d;
    272:           PetscCall(KSPAGMRESRoddecGivens(&c, &s, &Qloc[j * nloc + i], 0));
    273:           zp          = zloc[i - 1];
    274:           zq          = zloc[i];
    275:           zloc[i - 1] = c * zp + s * zq;
    276:           zloc[i]     = -s * zp + c * zq;
    277:         }
    278:         PetscCall(KSPAGMRESRoddecGivens(&c, &s, &Qloc[d * nloc], 0));
    279:         if (rank == Last) {
    280:           zp      = y[d];
    281:           zq      = zloc[0];
    282:           y[d]    = c * zp + s * zq;
    283:           zloc[0] = -s * zp + c * zq;
    284:           PetscCallMPI(MPI_Send(&y[d], 1, MPIU_SCALAR, Ileft, agmres->tag, comm));
    285:         } else {
    286:           PetscCallMPI(MPI_Recv(&yd, 1, MPIU_SCALAR, Iright, agmres->tag, comm, &status));
    287:           zp      = yd;
    288:           zq      = zloc[0];
    289:           yd      = c * zp + s * zq;
    290:           zloc[0] = -s * zp + c * zq;
    291:           PetscCallMPI(MPI_Send(&yd, 1, MPIU_SCALAR, Ileft, agmres->tag, comm));
    292:         }
    293:       }
    294:     }
    295:   }
    296:   for (j = nvec - 1; j >= 0; j--) {
    297:     dpt = j * nloc + j;
    298:     if (tloc[j] != 0.0) {
    299:       len = nloc - j;
    300:       PetscCall(PetscBLASIntCast(len, &blen));
    301:       rho       = Qloc[dpt];
    302:       Qloc[dpt] = 1.0;
    303:       tt        = tloc[j] * (BLASdot_(&blen, &Qloc[dpt], &pas, &zloc[j], &pas));
    304:       PetscCallBLAS("BLASaxpy", BLASaxpy_(&blen, &tt, &Qloc[dpt], &pas, &zloc[j], &pas));
    305:       Qloc[dpt] = rho;
    306:     }
    307:   }
    308:   PetscCall(VecRestoreArray(Out, &zloc));
    309:   PetscCall(PetscFree(y));
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }