Actual source code: lmekrylov.c
slepc-3.16.2 2022-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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==PETSC_DEFAULT) lme->ncv = PetscMin(30,N);
44: if (lme->max_it==PETSC_DEFAULT) lme->max_it = 100;
45: LMEAllocateSolution(lme,1);
46: return(0);
47: }
49: PetscErrorCode LMESolve_Krylov_Lyapunov_Vec(LME lme,Vec b,PetscBool fixed,PetscInt rrank,BV C1,BV *X1,PetscInt *col,PetscBool *fail,PetscInt *totalits)
50: {
52: PetscInt n=0,m,ldh,ldg=0,i,j,rank=0,lrank,pass,nouter=0,its;
53: PetscReal bnorm,beta,errest;
54: PetscBool breakdown;
55: PetscScalar *Harray,*G=NULL,*Gnew=NULL,*L,*U,*r,*Qarray,sone=1.0,zero=0.0;
56: PetscBLASInt n_,m_,rk_;
57: Mat Q,H;
60: *fail = PETSC_FALSE;
61: its = 0;
62: m = lme->ncv;
63: ldh = m+1;
64: MatCreateSeqDense(PETSC_COMM_SELF,ldh,m,NULL,&H);
65: MatDenseGetArray(H,&Harray);
67: VecNorm(b,NORM_2,&bnorm);
68: if (!bnorm) SETERRQ(PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Cannot process a zero vector in the right-hand side");
70: for (pass=0;pass<2;pass++) {
72: /* set initial vector to b/||b|| */
73: BVInsertVec(lme->V,0,b);
74: BVScaleColumn(lme->V,0,1.0/bnorm);
76: /* Restart loop */
77: while ((pass==0 && !*fail) || (pass==1 && its+1<nouter)) {
78: its++;
80: /* compute Arnoldi factorization */
81: BVMatArnoldi(lme->V,lme->A,H,0,&m,&beta,&breakdown);
82: BVSetActiveColumns(lme->V,0,m);
84: if (pass==0) {
85: /* glue together the previous H and the new H obtained with Arnoldi */
86: ldg = n+m+1;
87: PetscCalloc1(ldg*(n+m),&Gnew);
88: for (j=0;j<m;j++) {
89: PetscArraycpy(Gnew+n+(j+n)*ldg,Harray+j*ldh,m);
90: }
91: Gnew[n+m+(n+m-1)*ldg] = beta;
92: if (G) {
93: for (j=0;j<n;j++) {
94: PetscArraycpy(Gnew+j*ldg,G+j*(n+1),n+1);
95: }
96: PetscFree(G);
97: }
98: G = Gnew;
99: n += m;
100: } else {
101: /* update Z = Z + V(:,1:m)*Q with Q=U(blk,:)*P(1:nrk,:)' */
102: MatCreateSeqDense(PETSC_COMM_SELF,m,*col+rank,NULL,&Q);
103: MatDenseGetArray(Q,&Qarray);
104: PetscBLASIntCast(m,&m_);
105: PetscBLASIntCast(n,&n_);
106: PetscBLASIntCast(rank,&rk_);
107: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m_,&rk_,&rk_,&sone,U+its*m,&n_,L,&n_,&zero,Qarray+(*col)*m,&m_));
108: MatDenseRestoreArray(Q,&Qarray);
109: BVSetActiveColumns(*X1,*col,*col+rank);
110: BVMult(*X1,1.0,1.0,lme->V,Q);
111: MatDestroy(&Q);
112: }
114: if (pass==0) {
115: /* solve compressed Lyapunov equation */
116: PetscCalloc1(n,&r);
117: PetscCalloc1(n*n,&L);
118: r[0] = bnorm;
119: errest = PetscAbsScalar(G[n+(n-1)*ldg]);
120: LMEDenseHessLyapunovChol(lme,n,G,ldg,1,r,n,L,n,&errest);
121: LMEMonitor(lme,*totalits+its,errest);
122: PetscFree(r);
124: /* check convergence */
125: if (errest<lme->tol) {
126: lme->errest += errest;
127: PetscMalloc1(n*n,&U);
128: /* transpose L */
129: for (j=0;j<n;j++) {
130: for (i=j+1;i<n;i++) {
131: L[i+j*n] = PetscConj(L[j+i*n]);
132: L[j+i*n] = 0.0;
133: }
134: }
135: LMEDenseRankSVD(lme,n,L,n,U,n,&lrank);
136: PetscInfo1(lme,"Rank of the Cholesky factor = %D\n",lrank);
137: nouter = its;
138: its = -1;
139: if (!fixed) { /* X1 was not set by user, allocate it with rank columns */
140: rank = lrank;
141: if (*col) {
142: BVResize(*X1,*col+rank,PETSC_TRUE);
143: } else {
144: BVDuplicateResize(C1,rank,X1);
145: }
146: } else rank = PetscMin(lrank,rrank);
147: PetscFree(G);
148: break;
149: } else {
150: PetscFree(L);
151: if (*totalits+its>=lme->max_it) *fail = PETSC_TRUE;
152: }
153: }
155: /* restart with vector v_{m+1} */
156: if (!*fail) {
157: BVCopyColumn(lme->V,m,0);
158: }
159: }
160: }
162: *col += rank;
163: *totalits += its+1;
164: MatDenseRestoreArray(H,&Harray);
165: MatDestroy(&H);
166: if (L) { PetscFree(L); }
167: if (U) { PetscFree(U); }
168: return(0);
169: }
171: PetscErrorCode LMESolve_Krylov_Lyapunov(LME lme)
172: {
174: PetscBool fail,fixed = lme->X? PETSC_TRUE: PETSC_FALSE;
175: PetscInt i,k,rank=0,col=0;
176: Vec b;
177: BV X1=NULL,C1;
178: Mat X1m,X1t,C1m;
181: MatLRCGetMats(lme->C,NULL,&C1m,NULL,NULL);
182: BVCreateFromMat(C1m,&C1);
183: BVSetFromOptions(C1);
184: BVGetActiveColumns(C1,NULL,&k);
185: if (fixed) {
186: MatLRCGetMats(lme->X,NULL,&X1m,NULL,NULL);
187: BVCreateFromMat(X1m,&X1);
188: BVSetFromOptions(X1);
189: BVGetActiveColumns(X1,NULL,&rank);
190: rank = rank/k;
191: }
192: for (i=0;i<k;i++) {
193: BVGetColumn(C1,i,&b);
194: LMESolve_Krylov_Lyapunov_Vec(lme,b,fixed,rank,C1,&X1,&col,&fail,&lme->its);
195: BVRestoreColumn(C1,i,&b);
196: if (fail) {
197: lme->reason = LME_DIVERGED_ITS;
198: break;
199: }
200: }
201: if (lme->reason==LME_CONVERGED_ITERATING) lme->reason = LME_CONVERGED_TOL;
202: BVCreateMat(X1,&X1t);
203: if (fixed) {
204: MatCopy(X1t,X1m,SAME_NONZERO_PATTERN);
205: } else {
206: MatCreateLRC(NULL,X1t,NULL,NULL,&lme->X);
207: }
208: MatDestroy(&X1t);
209: BVDestroy(&C1);
210: BVDestroy(&X1);
211: return(0);
212: }
214: SLEPC_EXTERN PetscErrorCode LMECreate_Krylov(LME lme)
215: {
217: lme->ops->solve[LME_LYAPUNOV] = LMESolve_Krylov_Lyapunov;
218: lme->ops->setup = LMESetUp_Krylov;
219: return(0);
220: }