Actual source code: lmekrylov.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */
     10: /*
     11:    SLEPc matrix equation solver: "krylov"

     13:    Method: Arnoldi with Eiermann-Ernst restart

     15:    Algorithm:

     17:        Project the equation onto the Arnoldi basis and solve the compressed
     18:        equation the Hessenberg matrix H, restart by discarding the Krylov
     19:        basis but keeping H.

     21:    References:

     23:        [1] Y. Saad, "Numerical solution of large Lyapunov equations", in
     24:            Signal processing, scattering and operator theory, and numerical
     25:            methods, vol. 5 of Progr. Systems Control Theory, pages 503-511,
     26:            1990.

     28:        [2] D. Kressner, "Memory-efficient Krylov subspace techniques for
     29:            solving large-scale Lyapunov equations", in 2008 IEEE Int. Conf.
     30:            Computer-Aided Control Systems, pages 613-618, 2008.
     31: */

     33: #include <slepc/private/lmeimpl.h>
     34: #include <slepcblaslapack.h>

     36: static PetscErrorCode LMESetUp_Krylov(LME lme)
     37: {
     38:   PetscInt       N;

     40:   PetscFunctionBegin;
     41:   PetscCall(MatGetSize(lme->A,&N,NULL));
     42:   if (lme->ncv==PETSC_DETERMINE) lme->ncv = PetscMin(30,N);
     43:   if (lme->max_it==PETSC_DETERMINE) lme->max_it = 100;
     44:   PetscCall(LMEAllocateSolution(lme,1));
     45:   PetscFunctionReturn(PETSC_SUCCESS);
     46: }

     48: static PetscErrorCode LMESolve_Krylov_Lyapunov_Vec(LME lme,Vec b,PetscBool fixed,PetscInt rrank,BV C1,BV *X1,PetscInt *col,PetscBool *fail,PetscInt *totalits)
     49: {
     50:   PetscInt       n=0,m,ldh,ldg=0,i,j,rank=0,lrank,pass,nouter=0,its;
     51:   PetscReal      bnorm,beta,errest;
     52:   PetscBool      breakdown;
     53:   PetscScalar    *Harray,*G=NULL,*Gnew=NULL,*L=NULL,*U=NULL,*r,*Qarray,sone=1.0,zero=0.0;
     54:   PetscBLASInt   n_,m_,rk_;
     55:   Mat            Q,H;

     57:   PetscFunctionBegin;
     58:   *fail = PETSC_FALSE;
     59:   its = 0;
     60:   m  = lme->ncv;
     61:   ldh = m+1;
     62:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ldh,m,NULL,&H));
     63:   PetscCall(MatDenseGetArray(H,&Harray));

     65:   PetscCall(VecNorm(b,NORM_2,&bnorm));
     66:   PetscCheck(bnorm,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Cannot process a zero vector in the right-hand side");

     68:   for (pass=0;pass<2;pass++) {

     70:     /* set initial vector to b/||b|| */
     71:     PetscCall(BVInsertVec(lme->V,0,b));
     72:     PetscCall(BVScaleColumn(lme->V,0,1.0/bnorm));

     74:     /* Restart loop */
     75:     while ((pass==0 && !*fail) || (pass==1 && its+1<nouter)) {
     76:       its++;

     78:       /* compute Arnoldi factorization */
     79:       PetscCall(BVMatArnoldi(lme->V,lme->A,H,0,&m,&beta,&breakdown));
     80:       PetscCall(BVSetActiveColumns(lme->V,0,m));

     82:       if (pass==0) {
     83:         /* glue together the previous H and the new H obtained with Arnoldi */
     84:         ldg = n+m+1;
     85:         PetscCall(PetscCalloc1(ldg*(n+m),&Gnew));
     86:         for (j=0;j<m;j++) PetscCall(PetscArraycpy(Gnew+n+(j+n)*ldg,Harray+j*ldh,m));
     87:         Gnew[n+m+(n+m-1)*ldg] = beta;
     88:         if (G) {
     89:           for (j=0;j<n;j++) PetscCall(PetscArraycpy(Gnew+j*ldg,G+j*(n+1),n+1));
     90:           PetscCall(PetscFree(G));
     91:         }
     92:         G = Gnew;
     93:         n += m;
     94:       } else {
     95:         /* update Z = Z + V(:,1:m)*Q    with   Q=U(blk,:)*P(1:nrk,:)'  */
     96:         PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,m,*col+rank,NULL,&Q));
     97:         PetscCall(MatDenseGetArray(Q,&Qarray));
     98:         PetscCall(PetscBLASIntCast(m,&m_));
     99:         PetscCall(PetscBLASIntCast(n,&n_));
    100:         PetscCall(PetscBLASIntCast(rank,&rk_));
    101:         PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&m_,&rk_,&rk_,&sone,U+its*m,&n_,L,&n_,&zero,Qarray+(*col)*m,&m_));
    102:         PetscCall(MatDenseRestoreArray(Q,&Qarray));
    103:         PetscCall(BVSetActiveColumns(*X1,*col,*col+rank));
    104:         PetscCall(BVMult(*X1,1.0,1.0,lme->V,Q));
    105:         PetscCall(MatDestroy(&Q));
    106:       }

    108:       if (pass==0) {
    109:         /* solve compressed Lyapunov equation */
    110:         PetscCall(PetscCalloc1(n,&r));
    111:         PetscCall(PetscCalloc1(n*n,&L));
    112:         r[0] = bnorm;
    113:         errest = PetscAbsScalar(G[n+(n-1)*ldg]);
    114:         PetscCall(LMEDenseHessLyapunovChol(lme,n,G,ldg,1,r,n,L,n,&errest));
    115:         PetscCall(LMEMonitor(lme,*totalits+its,errest));
    116:         PetscCall(PetscFree(r));

    118:         /* check convergence */
    119:         if (errest<lme->tol) {
    120:           lme->errest += errest;
    121:           PetscCall(PetscMalloc1(n*n,&U));
    122:           /* transpose L */
    123:           for (j=0;j<n;j++) {
    124:             for (i=j+1;i<n;i++) {
    125:               L[i+j*n] = PetscConj(L[j+i*n]);
    126:               L[j+i*n] = 0.0;
    127:             }
    128:           }
    129:           PetscCall(LMEDenseRankSVD(lme,n,L,n,U,n,&lrank));
    130:           PetscCall(PetscInfo(lme,"Rank of the Cholesky factor = %" PetscInt_FMT "\n",lrank));
    131:           nouter = its;
    132:           its = -1;
    133:           if (!fixed) {  /* X1 was not set by user, allocate it with rank columns */
    134:             rank = lrank;
    135:             if (*col) PetscCall(BVResize(*X1,*col+rank,PETSC_TRUE));
    136:             else PetscCall(BVDuplicateResize(C1,rank,X1));
    137:           } else rank = PetscMin(lrank,rrank);
    138:           PetscCall(PetscFree(G));
    139:           break;
    140:         } else {
    141:           PetscCall(PetscFree(L));
    142:           if (*totalits+its>=lme->max_it) *fail = PETSC_TRUE;
    143:         }
    144:       }

    146:       /* restart with vector v_{m+1} */
    147:       if (!*fail) PetscCall(BVCopyColumn(lme->V,m,0));
    148:     }
    149:   }

    151:   *col += rank;
    152:   *totalits += its+1;
    153:   PetscCall(MatDenseRestoreArray(H,&Harray));
    154:   PetscCall(MatDestroy(&H));
    155:   PetscCall(PetscFree(L));
    156:   PetscCall(PetscFree(U));
    157:   PetscFunctionReturn(PETSC_SUCCESS);
    158: }

    160: static PetscErrorCode LMESolve_Krylov_Lyapunov(LME lme)
    161: {
    162:   PetscBool      fail,fixed = lme->X? PETSC_TRUE: PETSC_FALSE;
    163:   PetscInt       i,k,rank=0,col=0;
    164:   Vec            b;
    165:   BV             X1=NULL,C1;
    166:   Mat            X1m,X1t,C1m;

    168:   PetscFunctionBegin;
    169:   PetscCall(MatLRCGetMats(lme->C,NULL,&C1m,NULL,NULL));
    170:   PetscCall(BVCreateFromMat(C1m,&C1));
    171:   PetscCall(BVSetFromOptions(C1));
    172:   PetscCall(BVGetActiveColumns(C1,NULL,&k));
    173:   if (fixed) {
    174:     PetscCall(MatLRCGetMats(lme->X,NULL,&X1m,NULL,NULL));
    175:     PetscCall(BVCreateFromMat(X1m,&X1));
    176:     PetscCall(BVSetFromOptions(X1));
    177:     PetscCall(BVGetActiveColumns(X1,NULL,&rank));
    178:     rank = rank/k;
    179:   }
    180:   for (i=0;i<k;i++) {
    181:     PetscCall(BVGetColumn(C1,i,&b));
    182:     PetscCall(LMESolve_Krylov_Lyapunov_Vec(lme,b,fixed,rank,C1,&X1,&col,&fail,&lme->its));
    183:     PetscCall(BVRestoreColumn(C1,i,&b));
    184:     if (fail) {
    185:       lme->reason = LME_DIVERGED_ITS;
    186:       break;
    187:     }
    188:   }
    189:   if (lme->reason==LME_CONVERGED_ITERATING) lme->reason = LME_CONVERGED_TOL;
    190:   PetscCall(BVCreateMat(X1,&X1t));
    191:   if (fixed) PetscCall(MatCopy(X1t,X1m,SAME_NONZERO_PATTERN));
    192:   else PetscCall(MatCreateLRC(NULL,X1t,NULL,NULL,&lme->X));
    193:   PetscCall(MatDestroy(&X1t));
    194:   PetscCall(BVDestroy(&C1));
    195:   PetscCall(BVDestroy(&X1));
    196:   PetscFunctionReturn(PETSC_SUCCESS);
    197: }

    199: /*MC
    200:    LMEKRYLOV - LMEKRYLOV = "krylov" - A restarted Krylov method for the solution of
    201:    linear matrix equations.

    203:    Note:
    204:    A description of the method is provided in chapter [](ch:lme)
    205:    of the Users Manual.

    207:    Level: beginner

    209: .seealso: [](ch:lme), `LME`, `LMEType`, `LMESetType()`
    210: M*/
    211: SLEPC_EXTERN PetscErrorCode LMECreate_Krylov(LME lme)
    212: {
    213:   PetscFunctionBegin;
    214:   lme->ops->solve[LME_LYAPUNOV]      = LMESolve_Krylov_Lyapunov;
    215:   lme->ops->setup                    = LMESetUp_Krylov;
    216:   PetscFunctionReturn(PETSC_SUCCESS);
    217: }