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 : #include <slepc/private/dsimpl.h>
12 : #include <slepcblaslapack.h>
13 :
14 : /*
15 : 1) Patterns of A and B
16 : DS_STATE_RAW: DS_STATE_INTERM/CONDENSED
17 : 0 n-1 0 n-1
18 : ------------- -------------
19 : 0 |* * * * * *| 0 |* * * * * *|
20 : |* * * * * *| | * * * * *|
21 : |* * * * * *| | * * * *|
22 : |* * * * * *| | * * * *|
23 : |* * * * * *| | * *|
24 : n-1 |* * * * * *| n-1 | *|
25 : ------------- -------------
26 :
27 : 2) Moreover, P and Q are assumed to be the identity in DS_STATE_INTERMEDIATE.
28 : */
29 :
30 : static PetscErrorCode CleanDenseSchur(PetscInt n,PetscInt k,PetscScalar *S,PetscInt ldS,PetscScalar *T,PetscInt ldT,PetscScalar *X,PetscInt ldX,PetscScalar *Y,PetscInt ldY);
31 :
32 29 : static PetscErrorCode DSAllocate_GNHEP(DS ds,PetscInt ld)
33 : {
34 29 : PetscFunctionBegin;
35 29 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
36 29 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
37 29 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_Z));
38 29 : PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q));
39 29 : PetscCall(PetscFree(ds->perm));
40 29 : PetscCall(PetscMalloc1(ld,&ds->perm));
41 29 : PetscFunctionReturn(PETSC_SUCCESS);
42 : }
43 :
44 4 : static PetscErrorCode DSView_GNHEP(DS ds,PetscViewer viewer)
45 : {
46 4 : PetscViewerFormat format;
47 :
48 4 : PetscFunctionBegin;
49 4 : PetscCall(PetscViewerGetFormat(viewer,&format));
50 4 : if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
51 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
52 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_B));
53 0 : if (ds->state>DS_STATE_INTERMEDIATE) {
54 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_Z));
55 0 : PetscCall(DSViewMat(ds,viewer,DS_MAT_Q));
56 : }
57 0 : if (ds->omat[DS_MAT_X]) PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
58 0 : if (ds->omat[DS_MAT_Y]) PetscCall(DSViewMat(ds,viewer,DS_MAT_Y));
59 0 : PetscFunctionReturn(PETSC_SUCCESS);
60 : }
61 :
62 11736 : static PetscErrorCode DSVectors_GNHEP_Eigen_Some(DS ds,PetscInt *k,PetscReal *rnorm,PetscBool left)
63 : {
64 11736 : PetscInt i;
65 11736 : PetscBLASInt n,ld,mout,info,*select,mm,inc=1,cols=1,zero=0;
66 11736 : PetscScalar *X,*Y,*XY,*Z,*Q,*A,*B,fone=1.0,fzero=0.0;
67 11736 : PetscReal norm,done=1.0;
68 11736 : PetscBool iscomplex = PETSC_FALSE;
69 11736 : const char *side;
70 :
71 11736 : PetscFunctionBegin;
72 11736 : PetscCall(PetscBLASIntCast(ds->n,&n));
73 11736 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
74 11736 : if (left) {
75 615 : X = NULL;
76 615 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
77 : side = "L";
78 : } else {
79 11121 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
80 11121 : Y = NULL;
81 11121 : side = "R";
82 : }
83 11736 : XY = left? Y: X;
84 11736 : PetscCall(DSAllocateWork_Private(ds,0,0,ld));
85 11736 : select = ds->iwork;
86 156940 : for (i=0;i<n;i++) select[i] = (PetscBLASInt)PETSC_FALSE;
87 11736 : if (ds->state <= DS_STATE_INTERMEDIATE) {
88 1 : PetscCall(DSSetIdentity(ds,DS_MAT_Q));
89 1 : PetscCall(DSSetIdentity(ds,DS_MAT_Z));
90 : }
91 11736 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
92 11736 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
93 11736 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
94 11736 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Z],&Z));
95 11736 : PetscCall(CleanDenseSchur(n,0,A,ld,B,ld,Q,ld,Z,ld));
96 11736 : if (ds->state < DS_STATE_CONDENSED) PetscCall(DSSetState(ds,DS_STATE_CONDENSED));
97 :
98 : /* compute k-th eigenvector */
99 11736 : select[*k] = (PetscBLASInt)PETSC_TRUE;
100 : #if defined(PETSC_USE_COMPLEX)
101 : mm = 1;
102 : PetscCall(DSAllocateWork_Private(ds,2*ld,2*ld,0));
103 : PetscCallBLAS("LAPACKtgevc",LAPACKtgevc_(side,"S",select,&n,A,&ld,B,&ld,PetscSafePointerPlusOffset(Y,(*k)*ld),&ld,PetscSafePointerPlusOffset(X,(*k)*ld),&ld,&mm,&mout,ds->work,ds->rwork,&info));
104 : #else
105 11736 : if ((*k)<n-1 && (A[ld*(*k)+(*k)+1] != 0.0 || B[ld*(*k)+(*k)+1] != 0.0)) iscomplex = PETSC_TRUE;
106 11736 : mm = iscomplex? 2: 1;
107 11736 : if (iscomplex) select[(*k)+1] = (PetscBLASInt)PETSC_TRUE;
108 11736 : PetscCall(DSAllocateWork_Private(ds,6*ld,0,0));
109 11736 : PetscCallBLAS("LAPACKtgevc",LAPACKtgevc_(side,"S",select,&n,A,&ld,B,&ld,PetscSafePointerPlusOffset(Y,(*k)*ld),&ld,PetscSafePointerPlusOffset(X,(*k)*ld),&ld,&mm,&mout,ds->work,&info));
110 : #endif
111 11736 : SlepcCheckLapackInfo("tgevc",info);
112 11736 : PetscCheck(select[*k] && mout==mm,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Wrong arguments in call to Lapack xTGEVC");
113 11736 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
114 11736 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
115 :
116 : /* accumulate and normalize eigenvectors */
117 11736 : PetscCall(PetscArraycpy(ds->work,XY+(*k)*ld,mm*ld));
118 11736 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&mm,&n,&fone,left?Z:Q,&ld,ds->work,&ld,&fzero,XY+(*k)*ld,&ld));
119 11736 : norm = BLASnrm2_(&n,XY+(*k)*ld,&inc);
120 : #if !defined(PETSC_USE_COMPLEX)
121 11736 : if (iscomplex) {
122 100 : norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,XY+(*k+1)*ld,&inc));
123 100 : cols = 2;
124 : }
125 : #endif
126 11736 : PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,XY+(*k)*ld,&ld,&info));
127 11736 : SlepcCheckLapackInfo("lascl",info);
128 11736 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
129 11736 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Z],&Z));
130 :
131 : /* set output arguments */
132 11736 : if (rnorm) {
133 18 : if (iscomplex) *rnorm = SlepcAbsEigenvalue(XY[n-1+(*k)*ld],XY[n-1+(*k+1)*ld]);
134 16 : else *rnorm = PetscAbsScalar(XY[n-1+(*k)*ld]);
135 : }
136 11736 : if (iscomplex) (*k)++;
137 22857 : PetscCall(MatDenseRestoreArray(ds->omat[left?DS_MAT_Y:DS_MAT_X],&XY));
138 11736 : PetscFunctionReturn(PETSC_SUCCESS);
139 : }
140 :
141 44 : static PetscErrorCode DSVectors_GNHEP_Eigen_All(DS ds,PetscBool left)
142 : {
143 44 : PetscInt i;
144 44 : PetscBLASInt n,ld,mout,info,inc = 1;
145 44 : PetscBool iscomplex;
146 44 : PetscScalar *X,*Y,*XY,*Q,*Z,*A,*B,tmp;
147 44 : PetscReal norm;
148 44 : const char *side,*back;
149 :
150 44 : PetscFunctionBegin;
151 44 : PetscCall(PetscBLASIntCast(ds->n,&n));
152 44 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
153 44 : if (left) {
154 0 : X = NULL;
155 0 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
156 : side = "L";
157 : } else {
158 44 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
159 44 : Y = NULL;
160 44 : side = "R";
161 : }
162 44 : XY = left? Y: X;
163 44 : if (ds->state <= DS_STATE_INTERMEDIATE) {
164 3 : PetscCall(DSSetIdentity(ds,DS_MAT_Q));
165 3 : PetscCall(DSSetIdentity(ds,DS_MAT_Z));
166 : }
167 44 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
168 44 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
169 44 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
170 44 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Z],&Z));
171 44 : PetscCall(CleanDenseSchur(n,0,A,ld,B,ld,Q,ld,Z,ld));
172 44 : if (ds->state>=DS_STATE_CONDENSED) {
173 : /* DSSolve() has been called, backtransform with matrix Q */
174 41 : back = "B";
175 41 : PetscCall(PetscArraycpy(left?Y:X,left?Z:Q,ld*ld));
176 : } else {
177 3 : back = "A";
178 3 : PetscCall(DSSetState(ds,DS_STATE_CONDENSED));
179 : }
180 : #if defined(PETSC_USE_COMPLEX)
181 : PetscCall(DSAllocateWork_Private(ds,2*ld,2*ld,0));
182 : PetscCallBLAS("LAPACKtgevc",LAPACKtgevc_(side,back,NULL,&n,A,&ld,B,&ld,Y,&ld,X,&ld,&n,&mout,ds->work,ds->rwork,&info));
183 : #else
184 44 : PetscCall(DSAllocateWork_Private(ds,6*ld,0,0));
185 44 : PetscCallBLAS("LAPACKtgevc",LAPACKtgevc_(side,back,NULL,&n,A,&ld,B,&ld,Y,&ld,X,&ld,&n,&mout,ds->work,&info));
186 : #endif
187 44 : SlepcCheckLapackInfo("tgevc",info);
188 44 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
189 44 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Z],&Z));
190 :
191 : /* normalize eigenvectors */
192 213 : for (i=0;i<n;i++) {
193 169 : iscomplex = (i<n-1 && (A[i+1+i*ld]!=0.0 || B[i+1+i*ld]!=0.0))? PETSC_TRUE: PETSC_FALSE;
194 169 : norm = BLASnrm2_(&n,XY+i*ld,&inc);
195 : #if !defined(PETSC_USE_COMPLEX)
196 169 : if (iscomplex) {
197 9 : tmp = BLASnrm2_(&n,XY+(i+1)*ld,&inc);
198 9 : norm = SlepcAbsEigenvalue(norm,tmp);
199 : }
200 : #endif
201 169 : tmp = 1.0 / norm;
202 169 : PetscCallBLAS("BLASscal",BLASscal_(&n,&tmp,XY+i*ld,&inc));
203 : #if !defined(PETSC_USE_COMPLEX)
204 169 : if (iscomplex) PetscCallBLAS("BLASscal",BLASscal_(&n,&tmp,XY+(i+1)*ld,&inc));
205 : #endif
206 169 : if (iscomplex) i++;
207 : }
208 44 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
209 44 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
210 88 : PetscCall(MatDenseRestoreArray(ds->omat[left?DS_MAT_Y:DS_MAT_X],&XY));
211 44 : PetscFunctionReturn(PETSC_SUCCESS);
212 : }
213 :
214 11780 : static PetscErrorCode DSVectors_GNHEP(DS ds,DSMatType mat,PetscInt *k,PetscReal *rnorm)
215 : {
216 11780 : PetscFunctionBegin;
217 11780 : switch (mat) {
218 11780 : case DS_MAT_X:
219 : case DS_MAT_Y:
220 11780 : if (k) PetscCall(DSVectors_GNHEP_Eigen_Some(ds,k,rnorm,mat == DS_MAT_Y?PETSC_TRUE:PETSC_FALSE));
221 44 : else PetscCall(DSVectors_GNHEP_Eigen_All(ds,mat == DS_MAT_Y?PETSC_TRUE:PETSC_FALSE));
222 11780 : break;
223 0 : default:
224 0 : SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
225 : }
226 11780 : PetscFunctionReturn(PETSC_SUCCESS);
227 : }
228 :
229 2230 : static PetscErrorCode DSSort_GNHEP_Arbitrary(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
230 : {
231 2230 : PetscInt i;
232 2230 : PetscBLASInt info,n,ld,mout,lwork,liwork,*iwork,*selection,zero_=0,true_=1;
233 2230 : PetscScalar *S,*T,*Q,*Z,*work,*beta;
234 :
235 2230 : PetscFunctionBegin;
236 2230 : if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
237 2230 : PetscCall(PetscBLASIntCast(ds->n,&n));
238 2230 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
239 : #if !defined(PETSC_USE_COMPLEX)
240 2230 : lwork = 4*n+16;
241 : #else
242 : lwork = 1;
243 : #endif
244 2230 : liwork = 1;
245 2230 : PetscCall(DSAllocateWork_Private(ds,lwork+2*n,0,liwork+n));
246 2230 : beta = ds->work;
247 2230 : work = ds->work + n;
248 2230 : PetscCall(PetscBLASIntCast(ds->lwork-n,&lwork));
249 2230 : selection = ds->iwork;
250 2230 : iwork = ds->iwork + n;
251 2230 : PetscCall(PetscBLASIntCast(ds->liwork-n,&liwork));
252 : /* Compute the selected eigenvalue to be in the leading position */
253 2230 : PetscCall(DSSortEigenvalues_Private(ds,rr,ri,ds->perm,PETSC_FALSE));
254 2230 : PetscCall(PetscArrayzero(selection,n));
255 9673 : for (i=0; i<*k; i++) selection[ds->perm[i]] = 1;
256 2230 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&S));
257 2230 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&T));
258 2230 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
259 2230 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Z],&Z));
260 : #if !defined(PETSC_USE_COMPLEX)
261 2230 : PetscCallBLAS("LAPACKtgsen",LAPACKtgsen_(&zero_,&true_,&true_,selection,&n,S,&ld,T,&ld,wr,wi,beta,Z,&ld,Q,&ld,&mout,NULL,NULL,NULL,work,&lwork,iwork,&liwork,&info));
262 : #else
263 : PetscCallBLAS("LAPACKtgsen",LAPACKtgsen_(&zero_,&true_,&true_,selection,&n,S,&ld,T,&ld,wr,beta,Z,&ld,Q,&ld,&mout,NULL,NULL,NULL,work,&lwork,iwork,&liwork,&info));
264 : #endif
265 2230 : SlepcCheckLapackInfo("tgsen",info);
266 2230 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&S));
267 2230 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&T));
268 2230 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
269 2230 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Z],&Z));
270 2230 : *k = mout;
271 26834 : for (i=0;i<n;i++) {
272 24604 : if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
273 24604 : else wr[i] /= beta[i];
274 : #if !defined(PETSC_USE_COMPLEX)
275 24604 : if (beta[i]==0.0) wi[i] = (wi[i]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
276 24604 : else wi[i] /= beta[i];
277 : #endif
278 : }
279 2230 : PetscFunctionReturn(PETSC_SUCCESS);
280 : }
281 :
282 6 : static PetscErrorCode DSSort_GNHEP_Total(DS ds,PetscScalar *wr,PetscScalar *wi)
283 : {
284 6 : PetscScalar re;
285 6 : PetscInt i,j,pos,result;
286 6 : PetscBLASInt ifst,ilst,info,n,ld,one=1;
287 6 : PetscScalar *S,*T,*Z,*Q;
288 : #if !defined(PETSC_USE_COMPLEX)
289 6 : PetscBLASInt lwork;
290 6 : PetscScalar *work,a,safmin,scale1,scale2,im;
291 : #endif
292 :
293 6 : PetscFunctionBegin;
294 6 : if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
295 6 : PetscCall(PetscBLASIntCast(ds->n,&n));
296 6 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
297 6 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&S));
298 6 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&T));
299 6 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
300 6 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Z],&Z));
301 : #if !defined(PETSC_USE_COMPLEX)
302 6 : lwork = -1;
303 6 : PetscCallBLAS("LAPACKtgexc",LAPACKtgexc_(&one,&one,&ld,NULL,&ld,NULL,&ld,NULL,&ld,NULL,&ld,&one,&one,&a,&lwork,&info));
304 6 : SlepcCheckLapackInfo("tgexc",info);
305 6 : safmin = LAPACKlamch_("S");
306 6 : PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
307 6 : PetscCall(DSAllocateWork_Private(ds,lwork,0,0));
308 6 : work = ds->work;
309 : #endif
310 : /* selection sort */
311 84 : for (i=ds->l;i<n-1;i++) {
312 78 : re = wr[i];
313 : #if !defined(PETSC_USE_COMPLEX)
314 78 : im = wi[i];
315 : #endif
316 78 : pos = 0;
317 78 : j = i+1; /* j points to the next eigenvalue */
318 : #if !defined(PETSC_USE_COMPLEX)
319 78 : if (im != 0) j=i+2;
320 : #endif
321 : /* find minimum eigenvalue */
322 618 : for (;j<n;j++) {
323 : #if !defined(PETSC_USE_COMPLEX)
324 540 : PetscCall(SlepcSCCompare(ds->sc,re,im,wr[j],wi[j],&result));
325 : #else
326 : PetscCall(SlepcSCCompare(ds->sc,re,0.0,wr[j],0.0,&result));
327 : #endif
328 540 : if (result > 0) {
329 443 : re = wr[j];
330 : #if !defined(PETSC_USE_COMPLEX)
331 443 : im = wi[j];
332 : #endif
333 443 : pos = j;
334 : }
335 : #if !defined(PETSC_USE_COMPLEX)
336 540 : if (wi[j] != 0) j++;
337 : #endif
338 : }
339 78 : if (pos) {
340 : /* interchange blocks */
341 72 : PetscCall(PetscBLASIntCast(pos+1,&ifst));
342 72 : PetscCall(PetscBLASIntCast(i+1,&ilst));
343 : #if !defined(PETSC_USE_COMPLEX)
344 72 : PetscCallBLAS("LAPACKtgexc",LAPACKtgexc_(&one,&one,&n,S,&ld,T,&ld,Z,&ld,Q,&ld,&ifst,&ilst,work,&lwork,&info));
345 : #else
346 : PetscCallBLAS("LAPACKtgexc",LAPACKtgexc_(&one,&one,&n,S,&ld,T,&ld,Z,&ld,Q,&ld,&ifst,&ilst,&info));
347 : #endif
348 72 : SlepcCheckLapackInfo("tgexc",info);
349 : /* recover original eigenvalues from T and S matrices */
350 672 : for (j=i;j<n;j++) {
351 : #if !defined(PETSC_USE_COMPLEX)
352 600 : if (j<n-1 && S[j*ld+j+1] != 0.0) {
353 : /* complex conjugate eigenvalue */
354 151 : PetscCallBLAS("LAPACKlag2",LAPACKlag2_(S+j*ld+j,&ld,T+j*ld+j,&ld,&safmin,&scale1,&scale2,&re,&a,&im));
355 151 : wr[j] = re / scale1;
356 151 : wi[j] = im / scale1;
357 151 : wr[j+1] = a / scale2;
358 151 : wi[j+1] = -wi[j];
359 151 : j++;
360 : } else
361 : #endif
362 : {
363 449 : if (T[j*ld+j] == 0.0) wr[j] = (PetscRealPart(S[j*ld+j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
364 449 : else wr[j] = S[j*ld+j] / T[j*ld+j];
365 : #if !defined(PETSC_USE_COMPLEX)
366 449 : wi[j] = 0.0;
367 : #endif
368 : }
369 : }
370 : }
371 : #if !defined(PETSC_USE_COMPLEX)
372 78 : if (wi[i] != 0.0) i++;
373 : #endif
374 : }
375 6 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&S));
376 6 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&T));
377 6 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
378 6 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Z],&Z));
379 6 : PetscFunctionReturn(PETSC_SUCCESS);
380 : }
381 :
382 2236 : static PetscErrorCode DSSort_GNHEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
383 : {
384 2236 : PetscFunctionBegin;
385 2236 : if (!rr || wr == rr) PetscCall(DSSort_GNHEP_Total(ds,wr,wi));
386 2230 : else PetscCall(DSSort_GNHEP_Arbitrary(ds,wr,wi,rr,ri,k));
387 2236 : PetscFunctionReturn(PETSC_SUCCESS);
388 : }
389 :
390 5 : static PetscErrorCode DSUpdateExtraRow_GNHEP(DS ds)
391 : {
392 5 : PetscInt i;
393 5 : PetscBLASInt n,ld,incx=1;
394 5 : PetscScalar *A,*B,*x,*y,one=1.0,zero=0.0;
395 5 : const PetscScalar *Q;
396 :
397 5 : PetscFunctionBegin;
398 5 : PetscCall(PetscBLASIntCast(ds->n,&n));
399 5 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
400 5 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
401 5 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
402 5 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
403 5 : PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
404 5 : x = ds->work;
405 5 : y = ds->work+ld;
406 95 : for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
407 5 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
408 95 : for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
409 95 : for (i=0;i<n;i++) x[i] = PetscConj(B[n+i*ld]);
410 5 : PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
411 95 : for (i=0;i<n;i++) B[n+i*ld] = PetscConj(y[i]);
412 5 : ds->k = n;
413 5 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
414 5 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
415 5 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
416 5 : PetscFunctionReturn(PETSC_SUCCESS);
417 : }
418 :
419 : /*
420 : Write zeros from the column k to n in the lower triangular part of the
421 : matrices S and T, and inside 2-by-2 diagonal blocks of T in order to
422 : make (S,T) a valid Schur decompositon.
423 : */
424 11780 : static PetscErrorCode CleanDenseSchur(PetscInt n,PetscInt k,PetscScalar *S,PetscInt ldS,PetscScalar *T,PetscInt ldT,PetscScalar *X,PetscInt ldX,PetscScalar *Y,PetscInt ldY)
425 : {
426 11780 : PetscInt i;
427 : #if defined(PETSC_USE_COMPLEX)
428 : PetscInt j;
429 : PetscScalar s;
430 : #else
431 11780 : PetscBLASInt ldS_,ldT_,n_i,n_i_2,one=1,n_,i_2,i_;
432 11780 : PetscScalar b11,b22,sr,cr,sl,cl;
433 : #endif
434 :
435 11780 : PetscFunctionBegin;
436 : #if defined(PETSC_USE_COMPLEX)
437 : for (i=k; i<n; i++) {
438 : /* Some functions need the diagonal elements in T be real */
439 : if (T && PetscImaginaryPart(T[ldT*i+i]) != 0.0) {
440 : s = PetscConj(T[ldT*i+i])/PetscAbsScalar(T[ldT*i+i]);
441 : for (j=0;j<=i;j++) {
442 : T[ldT*i+j] *= s;
443 : S[ldS*i+j] *= s;
444 : }
445 : T[ldT*i+i] = PetscRealPart(T[ldT*i+i]);
446 : if (X) for (j=0;j<n;j++) X[ldX*i+j] *= s;
447 : }
448 : j = i+1;
449 : if (j<n) {
450 : S[ldS*i+j] = 0.0;
451 : if (T) T[ldT*i+j] = 0.0;
452 : }
453 : }
454 : #else
455 11780 : PetscCall(PetscBLASIntCast(ldS,&ldS_));
456 11780 : PetscCall(PetscBLASIntCast(ldT,&ldT_));
457 11780 : PetscCall(PetscBLASIntCast(n,&n_));
458 144776 : for (i=k;i<n-1;i++) {
459 132996 : if (S[ldS*i+i+1] != 0.0) {
460 : /* Check if T(i+1,i) and T(i,i+1) are zero */
461 700 : if (T[ldT*(i+1)+i] != 0.0 || T[ldT*i+i+1] != 0.0) {
462 : /* Check if T(i+1,i) and T(i,i+1) are negligible */
463 1 : if (PetscAbs(T[ldT*(i+1)+i])+PetscAbs(T[ldT*i+i+1]) < (PetscAbs(T[ldT*i+i])+PetscAbs(T[ldT*(i+1)+i+1]))*PETSC_MACHINE_EPSILON) {
464 0 : T[ldT*i+i+1] = 0.0;
465 0 : T[ldT*(i+1)+i] = 0.0;
466 : } else {
467 : /* If one of T(i+1,i) or T(i,i+1) is negligible, we make zero the other element */
468 1 : if (PetscAbs(T[ldT*i+i+1]) < (PetscAbs(T[ldT*i+i])+PetscAbs(T[ldT*(i+1)+i+1])+PetscAbs(T[ldT*(i+1)+i]))*PETSC_MACHINE_EPSILON) {
469 1 : PetscCallBLAS("LAPACKlasv2",LAPACKlasv2_(&T[ldT*i+i],&T[ldT*(i+1)+i],&T[ldT*(i+1)+i+1],&b22,&b11,&sl,&cl,&sr,&cr));
470 0 : } else if (PetscAbs(T[ldT*(i+1)+i]) < (PetscAbs(T[ldT*i+i])+PetscAbs(T[ldT*(i+1)+i+1])+PetscAbs(T[ldT*i+i+1]))*PETSC_MACHINE_EPSILON) {
471 0 : PetscCallBLAS("LAPACKlasv2",LAPACKlasv2_(&T[ldT*i+i],&T[ldT*i+i+1],&T[ldT*(i+1)+i+1],&b22,&b11,&sr,&cr,&sl,&cl));
472 0 : } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported format. Call DSSolve before this function");
473 1 : PetscCall(PetscBLASIntCast(n-i,&n_i));
474 1 : n_i_2 = n_i - 2;
475 1 : PetscCall(PetscBLASIntCast(i+2,&i_2));
476 1 : PetscCall(PetscBLASIntCast(i,&i_));
477 1 : if (b11 < 0.0) {
478 0 : cr = -cr; sr = -sr;
479 0 : b11 = -b11; b22 = -b22;
480 : }
481 1 : PetscCallBLAS("BLASrot",BLASrot_(&n_i,&S[ldS*i+i],&ldS_,&S[ldS*i+i+1],&ldS_,&cl,&sl));
482 1 : PetscCallBLAS("BLASrot",BLASrot_(&i_2,&S[ldS*i],&one,&S[ldS*(i+1)],&one,&cr,&sr));
483 1 : PetscCallBLAS("BLASrot",BLASrot_(&n_i_2,&T[ldT*(i+2)+i],&ldT_,&T[ldT*(i+2)+i+1],&ldT_,&cl,&sl));
484 1 : PetscCallBLAS("BLASrot",BLASrot_(&i_,&T[ldT*i],&one,&T[ldT*(i+1)],&one,&cr,&sr));
485 1 : if (X) PetscCallBLAS("BLASrot",BLASrot_(&n_,&X[ldX*i],&one,&X[ldX*(i+1)],&one,&cr,&sr));
486 1 : if (Y) PetscCallBLAS("BLASrot",BLASrot_(&n_,&Y[ldY*i],&one,&Y[ldY*(i+1)],&one,&cl,&sl));
487 1 : T[ldT*i+i] = b11; T[ldT*i+i+1] = 0.0;
488 1 : T[ldT*(i+1)+i] = 0.0; T[ldT*(i+1)+i+1] = b22;
489 : }
490 : }
491 : i++;
492 : }
493 : }
494 : #endif
495 11780 : PetscFunctionReturn(PETSC_SUCCESS);
496 : }
497 :
498 1170 : static PetscErrorCode DSSolve_GNHEP(DS ds,PetscScalar *wr,PetscScalar *wi)
499 : {
500 1170 : PetscScalar *work,*beta,a;
501 1170 : PetscInt i;
502 1170 : PetscBLASInt lwork,info,n,ld,iaux;
503 1170 : PetscScalar *A,*B,*Z,*Q;
504 1170 : PetscBool usegges3=(ds->method==1)?PETSC_TRUE:PETSC_FALSE;
505 :
506 1170 : PetscFunctionBegin;
507 : #if !defined(PETSC_USE_COMPLEX)
508 1170 : PetscAssertPointer(wi,3);
509 : #endif
510 1170 : PetscCall(PetscBLASIntCast(ds->n,&n));
511 1170 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
512 1170 : lwork = -1;
513 1170 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
514 1170 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
515 1170 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
516 1170 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Z],&Z));
517 : #if !defined(PETSC_USE_COMPLEX)
518 1170 : if (usegges3) PetscCallBLAS("LAPACKgges3",LAPACKgges3_("V","V","N",NULL,&n,A,&ld,B,&ld,&iaux,wr,wi,NULL,Z,&ld,Q,&ld,&a,&lwork,NULL,&info));
519 1170 : else PetscCallBLAS("LAPACKgges",LAPACKgges_("V","V","N",NULL,&n,A,&ld,B,&ld,&iaux,wr,wi,NULL,Z,&ld,Q,&ld,&a,&lwork,NULL,&info));
520 1170 : PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
521 1170 : PetscCall(DSAllocateWork_Private(ds,lwork+ld,0,0));
522 1170 : beta = ds->work;
523 1170 : work = beta+ds->n;
524 1170 : PetscCall(PetscBLASIntCast(ds->lwork-ds->n,&lwork));
525 1170 : if (usegges3) PetscCallBLAS("LAPACKgges3",LAPACKgges3_("V","V","N",NULL,&n,A,&ld,B,&ld,&iaux,wr,wi,beta,Z,&ld,Q,&ld,work,&lwork,NULL,&info));
526 1170 : else PetscCallBLAS("LAPACKgges",LAPACKgges_("V","V","N",NULL,&n,A,&ld,B,&ld,&iaux,wr,wi,beta,Z,&ld,Q,&ld,work,&lwork,NULL,&info));
527 : #else
528 : if (usegges3) PetscCallBLAS("LAPACKgges3",LAPACKgges3_("V","V","N",NULL,&n,A,&ld,B,&ld,&iaux,wr,NULL,Z,&ld,Q,&ld,&a,&lwork,NULL,NULL,&info));
529 : else PetscCallBLAS("LAPACKgges",LAPACKgges_("V","V","N",NULL,&n,A,&ld,B,&ld,&iaux,wr,NULL,Z,&ld,Q,&ld,&a,&lwork,NULL,NULL,&info));
530 : PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
531 : PetscCall(DSAllocateWork_Private(ds,lwork+ld,8*ld,0));
532 : beta = ds->work;
533 : work = beta+ds->n;
534 : PetscCall(PetscBLASIntCast(ds->lwork-ds->n,&lwork));
535 : if (usegges3) PetscCallBLAS("LAPACKgges3",LAPACKgges3_("V","V","N",NULL,&n,A,&ld,B,&ld,&iaux,wr,beta,Z,&ld,Q,&ld,work,&lwork,ds->rwork,NULL,&info));
536 : else PetscCallBLAS("LAPACKgges",LAPACKgges_("V","V","N",NULL,&n,A,&ld,B,&ld,&iaux,wr,beta,Z,&ld,Q,&ld,work,&lwork,ds->rwork,NULL,&info));
537 : #endif
538 1170 : SlepcCheckLapackInfo(usegges3?"gges3":"gges",info);
539 13845 : for (i=0;i<n;i++) {
540 12675 : if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
541 12675 : else wr[i] /= beta[i];
542 : #if !defined(PETSC_USE_COMPLEX)
543 12675 : if (beta[i]==0.0) wi[i] = (wi[i]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
544 12675 : else wi[i] /= beta[i];
545 : #else
546 : if (wi) wi[i] = 0.0;
547 : #endif
548 : }
549 1170 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
550 1170 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
551 1170 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
552 1170 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Z],&Z));
553 1170 : PetscFunctionReturn(PETSC_SUCCESS);
554 : }
555 :
556 : #if !defined(PETSC_HAVE_MPIUNI)
557 56 : static PetscErrorCode DSSynchronize_GNHEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
558 : {
559 56 : PetscInt ld=ds->ld,l=ds->l,k;
560 56 : PetscMPIInt n,rank,off=0,size,ldn;
561 56 : PetscScalar *A,*B,*Q,*Z;
562 :
563 56 : PetscFunctionBegin;
564 56 : k = 2*(ds->n-l)*ld;
565 56 : if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
566 56 : if (eigr) k += (ds->n-l);
567 56 : if (eigi) k += (ds->n-l);
568 56 : PetscCall(DSAllocateWork_Private(ds,k,0,0));
569 56 : PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
570 56 : PetscCall(PetscMPIIntCast(ds->n-l,&n));
571 56 : PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
572 56 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
573 56 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
574 56 : if (ds->state>DS_STATE_RAW) {
575 56 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
576 56 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Z],&Z));
577 : }
578 56 : PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
579 56 : if (!rank) {
580 28 : PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
581 28 : PetscCallMPI(MPI_Pack(B+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
582 28 : if (ds->state>DS_STATE_RAW) {
583 28 : PetscCallMPI(MPI_Pack(Q+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
584 28 : PetscCallMPI(MPI_Pack(Z+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
585 : }
586 28 : if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
587 : #if !defined(PETSC_USE_COMPLEX)
588 28 : if (eigi) PetscCallMPI(MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
589 : #endif
590 : }
591 112 : PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
592 56 : if (rank) {
593 28 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
594 28 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,B+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
595 28 : if (ds->state>DS_STATE_RAW) {
596 28 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,Q+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
597 28 : PetscCallMPI(MPI_Unpack(ds->work,size,&off,Z+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
598 : }
599 28 : if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
600 : #if !defined(PETSC_USE_COMPLEX)
601 28 : if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
602 : #endif
603 : }
604 56 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
605 56 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
606 56 : if (ds->state>DS_STATE_RAW) {
607 56 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
608 56 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Z],&Z));
609 : }
610 56 : PetscFunctionReturn(PETSC_SUCCESS);
611 : }
612 : #endif
613 :
614 5 : static PetscErrorCode DSTruncate_GNHEP(DS ds,PetscInt n,PetscBool trim)
615 : {
616 5 : PetscInt i,ld=ds->ld,l=ds->l;
617 5 : PetscScalar *A,*B;
618 :
619 5 : PetscFunctionBegin;
620 5 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
621 5 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
622 : #if defined(PETSC_USE_DEBUG)
623 : /* make sure diagonal 2x2 block is not broken */
624 5 : PetscCheck(ds->state<DS_STATE_CONDENSED || n==0 || n==ds->n || (A[n+(n-1)*ld]==0.0 && B[n+(n-1)*ld]==0.0),PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"The given size would break a 2x2 block, call DSGetTruncateSize() first");
625 : #endif
626 5 : if (trim) {
627 3 : if (ds->extrarow) { /* clean extra row */
628 53 : for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
629 53 : for (i=l;i<ds->n;i++) B[ds->n+i*ld] = 0.0;
630 : }
631 3 : ds->l = 0;
632 3 : ds->k = 0;
633 3 : ds->n = n;
634 3 : ds->t = ds->n; /* truncated length equal to the new dimension */
635 : } else {
636 2 : if (ds->extrarow && ds->k==ds->n) {
637 : /* copy entries of extra row to the new position, then clean last row */
638 18 : for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
639 34 : for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
640 18 : for (i=l;i<n;i++) B[n+i*ld] = B[ds->n+i*ld];
641 34 : for (i=l;i<ds->n;i++) B[ds->n+i*ld] = 0.0;
642 : }
643 2 : ds->k = (ds->extrarow)? n: 0;
644 2 : ds->t = ds->n; /* truncated length equal to previous dimension */
645 2 : ds->n = n;
646 : }
647 5 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
648 5 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
649 5 : PetscFunctionReturn(PETSC_SUCCESS);
650 : }
651 :
652 : /*MC
653 : DSGNHEP - Dense Generalized Non-Hermitian Eigenvalue Problem.
654 :
655 : Level: beginner
656 :
657 : Notes:
658 : The problem is expressed as A*X = B*X*Lambda, where (A,B) is the input
659 : matrix pencil. Lambda is a diagonal matrix whose diagonal elements are the
660 : arguments of DSSolve(). After solve, (A,B) is overwritten with the
661 : generalized (real) Schur form (S,T) = (Z'*A*Q,Z'*B*Q), with the first
662 : matrix being upper quasi-triangular and the second one triangular.
663 :
664 : Used DS matrices:
665 : + DS_MAT_A - first problem matrix
666 : . DS_MAT_B - second problem matrix
667 : . DS_MAT_Q - first orthogonal/unitary transformation that reduces to
668 : generalized (real) Schur form
669 : - DS_MAT_Z - second orthogonal/unitary transformation that reduces to
670 : generalized (real) Schur form
671 :
672 : Implemented methods:
673 : + 0 - QZ iteration (_gges)
674 : - 1 - blocked QZ iteration (_gges3, if available)
675 :
676 : .seealso: DSCreate(), DSSetType(), DSType
677 : M*/
678 29 : SLEPC_EXTERN PetscErrorCode DSCreate_GNHEP(DS ds)
679 : {
680 29 : PetscFunctionBegin;
681 29 : ds->ops->allocate = DSAllocate_GNHEP;
682 29 : ds->ops->view = DSView_GNHEP;
683 29 : ds->ops->vectors = DSVectors_GNHEP;
684 29 : ds->ops->solve[0] = DSSolve_GNHEP;
685 : #if !defined(SLEPC_MISSING_LAPACK_GGES3)
686 29 : ds->ops->solve[1] = DSSolve_GNHEP;
687 : #endif
688 29 : ds->ops->sort = DSSort_GNHEP;
689 : #if !defined(PETSC_HAVE_MPIUNI)
690 29 : ds->ops->synchronize = DSSynchronize_GNHEP;
691 : #endif
692 29 : ds->ops->gettruncatesize = DSGetTruncateSize_Default;
693 29 : ds->ops->truncate = DSTruncate_GNHEP;
694 29 : ds->ops->update = DSUpdateExtraRow_GNHEP;
695 29 : PetscFunctionReturn(PETSC_SUCCESS);
696 : }
|
