Actual source code: lmekrylov.c
slepc-3.11.2 2019-07-30
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: }