Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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"
12 :
13 : Method: Arnoldi with Eiermann-Ernst restart
14 :
15 : Algorithm:
16 :
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.
20 :
21 : References:
22 :
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.
27 :
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 : */
32 :
33 : #include <slepc/private/lmeimpl.h>
34 : #include <slepcblaslapack.h>
35 :
36 7 : static PetscErrorCode LMESetUp_Krylov(LME lme)
37 : {
38 7 : PetscInt N;
39 :
40 7 : PetscFunctionBegin;
41 7 : PetscCall(MatGetSize(lme->A,&N,NULL));
42 7 : if (lme->ncv==PETSC_DETERMINE) lme->ncv = PetscMin(30,N);
43 7 : if (lme->max_it==PETSC_DETERMINE) lme->max_it = 100;
44 7 : PetscCall(LMEAllocateSolution(lme,1));
45 7 : PetscFunctionReturn(PETSC_SUCCESS);
46 : }
47 :
48 55 : 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 55 : PetscInt n=0,m,ldh,ldg=0,i,j,rank=0,lrank,pass,nouter=0,its;
51 55 : PetscReal bnorm,beta,errest;
52 55 : PetscBool breakdown;
53 55 : PetscScalar *Harray,*G=NULL,*Gnew=NULL,*L=NULL,*U=NULL,*r,*Qarray,sone=1.0,zero=0.0;
54 55 : PetscBLASInt n_,m_,rk_;
55 55 : Mat Q,H;
56 :
57 55 : PetscFunctionBegin;
58 55 : *fail = PETSC_FALSE;
59 55 : its = 0;
60 55 : m = lme->ncv;
61 55 : ldh = m+1;
62 55 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ldh,m,NULL,&H));
63 55 : PetscCall(MatDenseGetArray(H,&Harray));
64 :
65 55 : PetscCall(VecNorm(b,NORM_2,&bnorm));
66 55 : PetscCheck(bnorm,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Cannot process a zero vector in the right-hand side");
67 :
68 165 : for (pass=0;pass<2;pass++) {
69 :
70 : /* set initial vector to b/||b|| */
71 110 : PetscCall(BVInsertVec(lme->V,0,b));
72 110 : PetscCall(BVScaleColumn(lme->V,0,1.0/bnorm));
73 :
74 : /* Restart loop */
75 173 : while ((pass==0 && !*fail) || (pass==1 && its+1<nouter)) {
76 118 : its++;
77 :
78 : /* compute Arnoldi factorization */
79 118 : PetscCall(BVMatArnoldi(lme->V,lme->A,H,0,&m,&beta,&breakdown));
80 118 : PetscCall(BVSetActiveColumns(lme->V,0,m));
81 :
82 118 : if (pass==0) {
83 : /* glue together the previous H and the new H obtained with Arnoldi */
84 59 : ldg = n+m+1;
85 59 : PetscCall(PetscCalloc1(ldg*(n+m),&Gnew));
86 1757 : for (j=0;j<m;j++) PetscCall(PetscArraycpy(Gnew+n+(j+n)*ldg,Harray+j*ldh,m));
87 59 : Gnew[n+m+(n+m-1)*ldg] = beta;
88 59 : if (G) {
89 88 : for (j=0;j<n;j++) PetscCall(PetscArraycpy(Gnew+j*ldg,G+j*(n+1),n+1));
90 4 : PetscCall(PetscFree(G));
91 : }
92 59 : G = Gnew;
93 59 : n += m;
94 : } else {
95 : /* update Z = Z + V(:,1:m)*Q with Q=U(blk,:)*P(1:nrk,:)' */
96 59 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,m,*col+rank,NULL,&Q));
97 59 : PetscCall(MatDenseGetArray(Q,&Qarray));
98 59 : PetscCall(PetscBLASIntCast(m,&m_));
99 59 : PetscCall(PetscBLASIntCast(n,&n_));
100 59 : PetscCall(PetscBLASIntCast(rank,&rk_));
101 59 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&m_,&rk_,&rk_,&sone,U+its*m,&n_,L,&n_,&zero,Qarray+(*col)*m,&m_));
102 59 : PetscCall(MatDenseRestoreArray(Q,&Qarray));
103 59 : PetscCall(BVSetActiveColumns(*X1,*col,*col+rank));
104 59 : PetscCall(BVMult(*X1,1.0,1.0,lme->V,Q));
105 59 : PetscCall(MatDestroy(&Q));
106 : }
107 :
108 118 : if (pass==0) {
109 : /* solve compressed Lyapunov equation */
110 59 : PetscCall(PetscCalloc1(n,&r));
111 59 : PetscCall(PetscCalloc1(n*n,&L));
112 59 : r[0] = bnorm;
113 59 : errest = PetscAbsScalar(G[n+(n-1)*ldg]);
114 59 : PetscCall(LMEDenseHessLyapunovChol(lme,n,G,ldg,1,r,n,L,n,&errest));
115 59 : PetscCall(LMEMonitor(lme,*totalits+its,errest));
116 59 : PetscCall(PetscFree(r));
117 :
118 : /* check convergence */
119 59 : if (errest<lme->tol) {
120 55 : lme->errest += errest;
121 55 : PetscCall(PetscMalloc1(n*n,&U));
122 : /* transpose L */
123 1753 : for (j=0;j<n;j++) {
124 27327 : for (i=j+1;i<n;i++) {
125 25629 : L[i+j*n] = PetscConj(L[j+i*n]);
126 25629 : L[j+i*n] = 0.0;
127 : }
128 : }
129 55 : PetscCall(LMEDenseRankSVD(lme,n,L,n,U,n,&lrank));
130 55 : PetscCall(PetscInfo(lme,"Rank of the Cholesky factor = %" PetscInt_FMT "\n",lrank));
131 55 : nouter = its;
132 55 : its = -1;
133 55 : if (!fixed) { /* X1 was not set by user, allocate it with rank columns */
134 6 : rank = lrank;
135 6 : if (*col) PetscCall(BVResize(*X1,*col+rank,PETSC_TRUE));
136 3 : else PetscCall(BVDuplicateResize(C1,rank,X1));
137 49 : } else rank = PetscMin(lrank,rrank);
138 55 : PetscCall(PetscFree(G));
139 : break;
140 : } else {
141 4 : PetscCall(PetscFree(L));
142 4 : if (*totalits+its>=lme->max_it) *fail = PETSC_TRUE;
143 : }
144 : }
145 :
146 : /* restart with vector v_{m+1} */
147 236 : if (!*fail) PetscCall(BVCopyColumn(lme->V,m,0));
148 : }
149 : }
150 :
151 55 : *col += rank;
152 55 : *totalits += its+1;
153 55 : PetscCall(MatDenseRestoreArray(H,&Harray));
154 55 : PetscCall(MatDestroy(&H));
155 55 : if (L) PetscCall(PetscFree(L));
156 55 : if (U) PetscCall(PetscFree(U));
157 55 : PetscFunctionReturn(PETSC_SUCCESS);
158 : }
159 :
160 32 : static PetscErrorCode LMESolve_Krylov_Lyapunov(LME lme)
161 : {
162 32 : PetscBool fail,fixed = lme->X? PETSC_TRUE: PETSC_FALSE;
163 32 : PetscInt i,k,rank=0,col=0;
164 32 : Vec b;
165 32 : BV X1=NULL,C1;
166 32 : Mat X1m,X1t,C1m;
167 :
168 32 : PetscFunctionBegin;
169 32 : PetscCall(MatLRCGetMats(lme->C,NULL,&C1m,NULL,NULL));
170 32 : PetscCall(BVCreateFromMat(C1m,&C1));
171 32 : PetscCall(BVSetFromOptions(C1));
172 32 : PetscCall(BVGetActiveColumns(C1,NULL,&k));
173 32 : if (fixed) {
174 29 : PetscCall(MatLRCGetMats(lme->X,NULL,&X1m,NULL,NULL));
175 29 : PetscCall(BVCreateFromMat(X1m,&X1));
176 29 : PetscCall(BVSetFromOptions(X1));
177 29 : PetscCall(BVGetActiveColumns(X1,NULL,&rank));
178 29 : rank = rank/k;
179 : }
180 87 : for (i=0;i<k;i++) {
181 55 : PetscCall(BVGetColumn(C1,i,&b));
182 55 : PetscCall(LMESolve_Krylov_Lyapunov_Vec(lme,b,fixed,rank,C1,&X1,&col,&fail,&lme->its));
183 55 : PetscCall(BVRestoreColumn(C1,i,&b));
184 55 : if (fail) {
185 0 : lme->reason = LME_DIVERGED_ITS;
186 0 : break;
187 : }
188 : }
189 32 : if (lme->reason==LME_CONVERGED_ITERATING) lme->reason = LME_CONVERGED_TOL;
190 32 : PetscCall(BVCreateMat(X1,&X1t));
191 32 : if (fixed) PetscCall(MatCopy(X1t,X1m,SAME_NONZERO_PATTERN));
192 3 : else PetscCall(MatCreateLRC(NULL,X1t,NULL,NULL,&lme->X));
193 32 : PetscCall(MatDestroy(&X1t));
194 32 : PetscCall(BVDestroy(&C1));
195 32 : PetscCall(BVDestroy(&X1));
196 32 : PetscFunctionReturn(PETSC_SUCCESS);
197 : }
198 :
199 7 : SLEPC_EXTERN PetscErrorCode LMECreate_Krylov(LME lme)
200 : {
201 7 : PetscFunctionBegin;
202 7 : lme->ops->solve[LME_LYAPUNOV] = LMESolve_Krylov_Lyapunov;
203 7 : lme->ops->setup = LMESetUp_Krylov;
204 7 : PetscFunctionReturn(PETSC_SUCCESS);
205 : }
|
