Actual source code: lmekrylov.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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: PetscErrorCode LMESetUp_Krylov(LME lme)
 37: {
 39:   PetscInt       N;

 42:   MatGetSize(lme->A,&N,NULL);
 43:   if (!lme->ncv) lme->ncv = PetscMin(30,N);
 44:   if (!lme->max_it) lme->max_it = 100;
 45:   LMEAllocateSolution(lme,1);
 46:   return(0);
 47: }

 49: PetscErrorCode LMEBasicArnoldi(LME lme,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,PetscReal *beta,PetscBool *breakdown)
 50: {
 52:   PetscScalar    *a;
 53:   PetscInt       j,nc,n,m = *M;
 54:   Vec            vj,vj1,buf;

 57:   BVSetActiveColumns(lme->V,0,m);
 58:   for (j=k;j<m;j++) {
 59:     BVGetColumn(lme->V,j,&vj);
 60:     BVGetColumn(lme->V,j+1,&vj1);
 61:     MatMult(lme->A,vj,vj1);
 62:     BVRestoreColumn(lme->V,j,&vj);
 63:     BVRestoreColumn(lme->V,j+1,&vj1);
 64:     BVOrthonormalizeColumn(lme->V,j+1,PETSC_FALSE,beta,breakdown);
 65:     if (*breakdown) {
 66:       *M = j+1;
 67:       break;
 68:     }
 69:   }
 70:   /* extract Hessenberg matrix from the BV object */
 71:   BVGetNumConstraints(lme->V,&nc);
 72:   BVGetSizes(lme->V,NULL,NULL,&n);
 73:   BVGetBufferVec(lme->V,&buf);
 74:   VecGetArray(buf,&a);
 75:   for (j=k;j<*M;j++) {
 76:     PetscMemcpy(H+j*ldh,a+nc+(j+1)*(nc+n),(j+2)*sizeof(PetscScalar));
 77:   }
 78:   VecRestoreArray(buf,&a);
 79:   return(0);
 80: }

 82: PetscErrorCode LMESolve_Krylov_Lyapunov_Vec(LME lme,Vec b,PetscBool fixed,PetscInt rrank,BV C1,BV *X1,PetscInt *col,PetscBool *fail,PetscInt *totalits)
 83: {
 85:   PetscInt       n=0,m,ldh,ldg,j,rank=0,lrank,pass,nouter=0,its;
 86:   PetscReal      bnorm,beta,errest;
 87:   PetscBool      breakdown;
 88:   PetscScalar    *H,*G=NULL,*Gnew=NULL,*Gcopy,*L,*U,*r,*Qarray,sone=1.0,zero=0.0;
 89:   PetscBLASInt   n_,m_,rk_;
 90:   Mat            Q;

 93:   *fail = PETSC_FALSE;
 94:   its = 0;
 95:   m  = lme->ncv;
 96:   ldh = m+1;
 97:   PetscCalloc1(ldh*m,&H);

 99:   VecNorm(b,NORM_2,&bnorm);
100:   if (!bnorm) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Cannot process a zero vector in the right-hand side");

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

104:     /* set initial vector to b/||b|| */
105:     BVInsertVec(lme->V,0,b);
106:     BVScaleColumn(lme->V,0,1.0/bnorm);

108:     /* Restart loop */
109:     while ((pass==0 && !*fail) || (pass==1 && its+1<nouter)) {
110:       its++;

112:       /* compute Arnoldi factorization */
113:       LMEBasicArnoldi(lme,H,ldh,0,&m,&beta,&breakdown);

115:       if (pass==0) {
116:         /* glue together the previous H and the new H obtained with Arnoldi */
117:         ldg = n+m+1;
118:         PetscCalloc1(ldg*(n+m),&Gnew);
119:         for (j=0;j<m;j++) {
120:           PetscMemcpy(Gnew+n+(j+n)*ldg,H+j*ldh,m*sizeof(PetscScalar));
121:         }
122:         Gnew[n+m+(n+m-1)*ldg] = beta;
123:         if (G) {
124:           for (j=0;j<n;j++) {
125:             PetscMemcpy(Gnew+j*ldg,G+j*(n+1),(n+1)*sizeof(PetscScalar));
126:           }
127:           PetscFree(G);
128:         }
129:         G = Gnew;
130:         n += m;
131:       } else {
132:         /* update Z = Z + V(:,1:m)*Q    with   Q=U(blk,:)*P(1:nrk,:)'  */
133:         MatCreateDense(PETSC_COMM_SELF,m,*col+rank,m,*col+rank,NULL,&Q);
134:         MatDenseGetArray(Q,&Qarray);
135:         PetscBLASIntCast(m,&m_);
136:         PetscBLASIntCast(n,&n_);
137:         PetscBLASIntCast(rank,&rk_);
138:         PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m_,&rk_,&n_,&sone,U+its*m,&n_,L,&n_,&zero,Qarray+(*col)*m,&m_));
139:         MatDenseRestoreArray(Q,&Qarray);
140:         BVSetActiveColumns(*X1,*col,*col+rank);
141:         BVMult(*X1,1.0,1.0,lme->V,Q);
142:         MatDestroy(&Q);
143:       }

145:       if (pass==0) {
146:         /* solve compressed Lyapunov equation */
147:         PetscCalloc2(n,&r,ldg*n,&Gcopy);
148:         PetscCalloc1(n*n,&L);
149:         r[0] = bnorm;
150:         PetscMemcpy(Gcopy,G,ldg*n*sizeof(PetscScalar));
151:         LMEDenseLyapunovChol(lme,Gcopy,n,ldg,r,L,n,&errest);
152:         LMEMonitor(lme,*totalits+its,errest);
153:         PetscFree2(r,Gcopy);

155:         /* check convergence */
156:         if (errest<lme->tol) {
157:           lme->errest += errest;
158:           PetscMalloc1(n*n,&U);
159:           LMERankSVD(lme,n,L,U,&lrank);
160:           nouter = its;
161:           its = -1;
162:           if (!fixed) {  /* X1 was not set by user, allocate it with rank columns */
163:             rank = lrank;
164:             if (*col) {
165:               BVResize(*X1,*col+rank,PETSC_TRUE);
166:             } else {
167:               BVDuplicateResize(C1,rank,X1);
168:             }
169:           } else rank = PetscMin(lrank,rrank);
170:           PetscFree(G);
171:           break;
172:         } else {
173:           PetscFree(L);
174:           if (*totalits+its>=lme->max_it) *fail = PETSC_TRUE;
175:         }
176:       }

178:       /* restart with vector v_{m+1} */
179:       if (!*fail) {
180:         BVCopyColumn(lme->V,m,0);
181:       }
182:     }
183:   }

185:   *col += rank;
186:   *totalits += its+1;
187:   PetscFree(H);
188:   if (L) { PetscFree(L); }
189:   if (U) { PetscFree(U); }
190:   return(0);
191: }

193: PetscErrorCode LMESolve_Krylov_Lyapunov(LME lme)
194: {
196:   PetscBool      fail,fixed = lme->X? PETSC_TRUE: PETSC_FALSE;
197:   PetscInt       i,k,rank,col=0;
198:   Vec            b;
199:   BV             X1=NULL,C1;
200:   Mat            X1m,X1t,C1m;

203:   MatLRCGetMats(lme->C,NULL,&C1m,NULL,NULL);
204:   BVCreateFromMat(C1m,&C1);
205:   BVSetFromOptions(C1);
206:   BVGetActiveColumns(C1,NULL,&k);
207:   if (fixed) {
208:     MatLRCGetMats(lme->X,NULL,&X1m,NULL,NULL);
209:     BVCreateFromMat(X1m,&X1);
210:     BVSetFromOptions(X1);
211:     BVGetActiveColumns(X1,NULL,&rank);
212:     rank = rank/k;
213:   }
214:   for (i=0;i<k;i++) {
215:     BVGetColumn(C1,i,&b);
216:     LMESolve_Krylov_Lyapunov_Vec(lme,b,fixed,rank,C1,&X1,&col,&fail,&lme->its);
217:     BVRestoreColumn(C1,i,&b);
218:     if (fail) {
219:       lme->reason = LME_DIVERGED_ITS;
220:       break;
221:     }
222:   }
223:   if (lme->reason==LME_CONVERGED_ITERATING) lme->reason = LME_CONVERGED_TOL;
224:   BVCreateMat(X1,&X1t);
225:   if (fixed) {
226:     MatCopy(X1t,X1m,SAME_NONZERO_PATTERN);
227:   } else {
228:     MatCreateLRC(NULL,X1t,NULL,NULL,&lme->X);
229:   }
230:   MatDestroy(&X1t);
231:   BVDestroy(&C1);
232:   BVDestroy(&X1);
233:   return(0);
234: }

236: SLEPC_EXTERN PetscErrorCode LMECreate_Krylov(LME lme)
237: {
239:   lme->ops->solve[LME_LYAPUNOV]      = LMESolve_Krylov_Lyapunov;
240:   lme->ops->setup                    = LMESetUp_Krylov;
241:   return(0);
242: }