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 : Utility subroutines common to several impls
12 : */
13 :
14 : #include <slepc/private/dsimpl.h> /*I "slepcds.h" I*/
15 : #include <slepcblaslapack.h>
16 :
17 : /*
18 : Compute the (real) Schur form of A. At the end, A is (quasi-)triangular and Q
19 : contains the unitary matrix of Schur vectors. Eigenvalues are returned in wr,wi
20 : */
21 3551 : PetscErrorCode DSSolve_NHEP_Private(DS ds,DSMatType mA,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
22 : {
23 3551 : PetscScalar *work,*tau,*A,*Q;
24 3551 : PetscInt i,j;
25 3551 : PetscBLASInt ilo,lwork,info,n,k,ld;
26 :
27 3551 : PetscFunctionBegin;
28 3551 : PetscCall(MatDenseGetArray(ds->omat[mA],&A));
29 3551 : PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
30 3551 : PetscCall(PetscBLASIntCast(ds->n,&n));
31 3551 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
32 3551 : PetscCall(PetscBLASIntCast(ds->l+1,&ilo));
33 3551 : PetscCall(PetscBLASIntCast(ds->k,&k));
34 3551 : PetscCall(DSAllocateWork_Private(ds,ld+6*ld,0,0));
35 3551 : tau = ds->work;
36 3551 : work = ds->work+ld;
37 3551 : lwork = 6*ld;
38 :
39 : /* initialize orthogonal matrix */
40 3551 : PetscCall(PetscArrayzero(Q,ld*ld));
41 67782 : for (i=0;i<n;i++) Q[i+i*ld] = 1.0;
42 3551 : if (n==1) { /* quick return */
43 33 : wr[0] = A[0];
44 33 : if (wi) wi[0] = 0.0;
45 33 : PetscFunctionReturn(PETSC_SUCCESS);
46 : }
47 :
48 : /* reduce to upper Hessenberg form */
49 3518 : if (ds->state<DS_STATE_INTERMEDIATE) {
50 2346 : PetscCallBLAS("LAPACKgehrd",LAPACKgehrd_(&n,&ilo,&n,A,&ld,tau,work,&lwork,&info));
51 2346 : SlepcCheckLapackInfo("gehrd",info);
52 43923 : for (j=0;j<n-1;j++) {
53 533465 : for (i=j+2;i<n;i++) {
54 491888 : Q[i+j*ld] = A[i+j*ld];
55 491888 : A[i+j*ld] = 0.0;
56 : }
57 : }
58 2346 : PetscCallBLAS("LAPACKorghr",LAPACKorghr_(&n,&ilo,&n,Q,&ld,tau,work,&lwork,&info));
59 2346 : SlepcCheckLapackInfo("orghr",info);
60 : }
61 :
62 : /* compute the (real) Schur form */
63 : #if !defined(PETSC_USE_COMPLEX)
64 3518 : PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,wi,Q,&ld,work,&lwork,&info));
65 7188 : for (j=0;j<ds->l;j++) {
66 3670 : if (j==n-1 || A[j+1+j*ld] == 0.0) {
67 : /* real eigenvalue */
68 3380 : wr[j] = A[j+j*ld];
69 3380 : wi[j] = 0.0;
70 : } else {
71 : /* complex eigenvalue */
72 290 : wr[j] = A[j+j*ld];
73 290 : wr[j+1] = A[j+j*ld];
74 290 : wi[j] = PetscSqrtReal(PetscAbsReal(A[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(A[j+(j+1)*ld]));
75 290 : wi[j+1] = -wi[j];
76 290 : j++;
77 : }
78 : }
79 : #else
80 : PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,Q,&ld,work,&lwork,&info));
81 : if (wi) for (i=ds->l;i<n;i++) wi[i] = 0.0;
82 : #endif
83 3518 : SlepcCheckLapackInfo("hseqr",info);
84 3518 : PetscCall(MatDenseRestoreArray(ds->omat[mA],&A));
85 3518 : PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
86 3518 : PetscFunctionReturn(PETSC_SUCCESS);
87 : }
88 :
89 : /*
90 : Sort a Schur form represented by the (quasi-)triangular matrix T and
91 : the unitary matrix Q, and return the sorted eigenvalues in wr,wi
92 : */
93 3933 : PetscErrorCode DSSort_NHEP_Total(DS ds,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
94 : {
95 3933 : PetscScalar re,*T,*Q;
96 3933 : PetscInt i,j,pos,result;
97 3933 : PetscBLASInt ifst,ilst,info,n,ld;
98 : #if !defined(PETSC_USE_COMPLEX)
99 3933 : PetscScalar *work,im;
100 : #endif
101 :
102 3933 : PetscFunctionBegin;
103 3933 : PetscCall(MatDenseGetArray(ds->omat[mT],&T));
104 3933 : PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
105 3933 : PetscCall(PetscBLASIntCast(ds->n,&n));
106 3933 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
107 : #if !defined(PETSC_USE_COMPLEX)
108 3933 : PetscCall(DSAllocateWork_Private(ds,ld,0,0));
109 3933 : work = ds->work;
110 : #endif
111 : /* selection sort */
112 57110 : for (i=ds->l;i<n-1;i++) {
113 53177 : re = wr[i];
114 : #if !defined(PETSC_USE_COMPLEX)
115 53177 : im = wi[i];
116 : #endif
117 53177 : pos = 0;
118 53177 : j=i+1; /* j points to the next eigenvalue */
119 : #if !defined(PETSC_USE_COMPLEX)
120 53177 : if (im != 0) j=i+2;
121 : #endif
122 : /* find minimum eigenvalue */
123 597892 : for (;j<n;j++) {
124 : #if !defined(PETSC_USE_COMPLEX)
125 544715 : PetscCall(SlepcSCCompare(ds->sc,re,im,wr[j],wi[j],&result));
126 : #else
127 : PetscCall(SlepcSCCompare(ds->sc,re,0.0,wr[j],0.0,&result));
128 : #endif
129 544715 : if (result > 0) {
130 213519 : re = wr[j];
131 : #if !defined(PETSC_USE_COMPLEX)
132 213519 : im = wi[j];
133 : #endif
134 213519 : pos = j;
135 : }
136 : #if !defined(PETSC_USE_COMPLEX)
137 544715 : if (wi[j] != 0) j++;
138 : #endif
139 : }
140 53177 : if (pos) {
141 : /* interchange blocks */
142 35702 : PetscCall(PetscBLASIntCast(pos+1,&ifst));
143 35702 : PetscCall(PetscBLASIntCast(i+1,&ilst));
144 : #if !defined(PETSC_USE_COMPLEX)
145 35702 : PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
146 : #else
147 : PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
148 : #endif
149 35702 : SlepcCheckLapackInfo("trexc",info);
150 : /* recover original eigenvalues from T matrix */
151 465195 : for (j=i;j<n;j++) {
152 429493 : wr[j] = T[j+j*ld];
153 : #if !defined(PETSC_USE_COMPLEX)
154 429493 : if (j<n-1 && T[j+1+j*ld] != 0.0) {
155 : /* complex conjugate eigenvalue */
156 43419 : wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
157 43419 : wr[j+1] = wr[j];
158 43419 : wi[j+1] = -wi[j];
159 43419 : j++;
160 386074 : } else wi[j] = 0.0;
161 : #endif
162 : }
163 : }
164 : #if !defined(PETSC_USE_COMPLEX)
165 53177 : if (wi[i] != 0) i++;
166 : #endif
167 : }
168 3933 : PetscCall(MatDenseRestoreArray(ds->omat[mT],&T));
169 3933 : PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
170 3933 : PetscFunctionReturn(PETSC_SUCCESS);
171 : }
172 :
173 : /*
174 : Reorder a Schur form represented by T,Q according to a permutation perm,
175 : and return the sorted eigenvalues in wr,wi
176 : */
177 1 : PetscErrorCode DSSortWithPermutation_NHEP_Private(DS ds,PetscInt *perm,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
178 : {
179 1 : PetscInt i,j,pos,inc=1;
180 1 : PetscBLASInt ifst,ilst,info,n,ld;
181 1 : PetscScalar *T,*Q;
182 : #if !defined(PETSC_USE_COMPLEX)
183 1 : PetscScalar *work;
184 : #endif
185 :
186 1 : PetscFunctionBegin;
187 1 : PetscCall(MatDenseGetArray(ds->omat[mT],&T));
188 1 : PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
189 1 : PetscCall(PetscBLASIntCast(ds->n,&n));
190 1 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
191 : #if !defined(PETSC_USE_COMPLEX)
192 1 : PetscCall(DSAllocateWork_Private(ds,ld,0,0));
193 1 : work = ds->work;
194 : #endif
195 7 : for (i=ds->l;i<n-1;i++) {
196 6 : pos = perm[i];
197 : #if !defined(PETSC_USE_COMPLEX)
198 6 : inc = (pos<n-1 && T[pos+1+pos*ld] != 0.0)? 2: 1;
199 : #endif
200 6 : if (pos!=i) {
201 : #if !defined(PETSC_USE_COMPLEX)
202 3 : PetscCheck((T[pos+(pos-1)*ld]==0.0 || perm[i+1]==pos-1) && (T[pos+1+pos*ld]==0.0 || perm[i+1]==pos+1),PETSC_COMM_SELF,PETSC_ERR_FP,"Invalid permutation due to a 2x2 block at position %" PetscInt_FMT,pos);
203 : #endif
204 : /* interchange blocks */
205 3 : PetscCall(PetscBLASIntCast(pos+1,&ifst));
206 3 : PetscCall(PetscBLASIntCast(i+1,&ilst));
207 : #if !defined(PETSC_USE_COMPLEX)
208 3 : PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
209 : #else
210 : PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
211 : #endif
212 3 : SlepcCheckLapackInfo("trexc",info);
213 30 : for (j=i+1;j<n;j++) {
214 27 : if (perm[j]>=i && perm[j]<pos) perm[j]+=inc;
215 : }
216 3 : perm[i] = i;
217 3 : if (inc==2) perm[i+1] = i+1;
218 : }
219 6 : if (inc==2) i++;
220 : }
221 : /* recover original eigenvalues from T matrix */
222 7 : for (j=ds->l;j<n;j++) {
223 6 : wr[j] = T[j+j*ld];
224 : #if !defined(PETSC_USE_COMPLEX)
225 6 : if (j<n-1 && T[j+1+j*ld] != 0.0) {
226 : /* complex conjugate eigenvalue */
227 6 : wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
228 6 : wr[j+1] = wr[j];
229 6 : wi[j+1] = -wi[j];
230 6 : j++;
231 0 : } else wi[j] = 0.0;
232 : #endif
233 : }
234 1 : PetscCall(MatDenseRestoreArray(ds->omat[mT],&T));
235 1 : PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
236 1 : PetscFunctionReturn(PETSC_SUCCESS);
237 : }
|
