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 3760 : PetscErrorCode DSSolve_NHEP_Private(DS ds,DSMatType mA,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
22 : {
23 3760 : PetscScalar *work,*tau,*A,*Q;
24 3760 : PetscInt i,j;
25 3760 : PetscBLASInt ilo,lwork,info,n,k,ld;
26 :
27 3760 : PetscFunctionBegin;
28 3760 : PetscCall(MatDenseGetArray(ds->omat[mA],&A));
29 3760 : PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
30 3760 : PetscCall(PetscBLASIntCast(ds->n,&n));
31 3760 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
32 3760 : PetscCall(PetscBLASIntCast(ds->l+1,&ilo));
33 3760 : PetscCall(PetscBLASIntCast(ds->k,&k));
34 3760 : PetscCall(DSAllocateWork_Private(ds,ld+6*ld,0,0));
35 3760 : tau = ds->work;
36 3760 : work = ds->work+ld;
37 3760 : lwork = 6*ld;
38 :
39 : /* initialize orthogonal matrix */
40 3760 : PetscCall(PetscArrayzero(Q,ld*ld));
41 72233 : for (i=0;i<n;i++) Q[i+i*ld] = 1.0;
42 3760 : 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 3727 : if (ds->state<DS_STATE_INTERMEDIATE) {
50 2544 : PetscCallBLAS("LAPACKgehrd",LAPACKgehrd_(&n,&ilo,&n,A,&ld,tau,work,&lwork,&info));
51 2544 : SlepcCheckLapackInfo("gehrd",info);
52 47885 : for (j=0;j<n-1;j++) {
53 583127 : for (i=j+2;i<n;i++) {
54 537786 : Q[i+j*ld] = A[i+j*ld];
55 537786 : A[i+j*ld] = 0.0;
56 : }
57 : }
58 2544 : PetscCallBLAS("LAPACKorghr",LAPACKorghr_(&n,&ilo,&n,Q,&ld,tau,work,&lwork,&info));
59 2544 : SlepcCheckLapackInfo("orghr",info);
60 : }
61 :
62 : /* compute the (real) Schur form */
63 : #if !defined(PETSC_USE_COMPLEX)
64 : PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,wi,Q,&ld,work,&lwork,&info));
65 : for (j=0;j<ds->l;j++) {
66 : if (j==n-1 || A[j+1+j*ld] == 0.0) {
67 : /* real eigenvalue */
68 : wr[j] = A[j+j*ld];
69 : wi[j] = 0.0;
70 : } else {
71 : /* complex eigenvalue */
72 : wr[j] = A[j+j*ld];
73 : wr[j+1] = A[j+j*ld];
74 : wi[j] = PetscSqrtReal(PetscAbsReal(A[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(A[j+(j+1)*ld]));
75 : wi[j+1] = -wi[j];
76 : j++;
77 : }
78 : }
79 : #else
80 3727 : PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,Q,&ld,work,&lwork,&info));
81 66112 : if (wi) for (i=ds->l;i<n;i++) wi[i] = 0.0;
82 : #endif
83 3727 : SlepcCheckLapackInfo("hseqr",info);
84 3727 : PetscCall(MatDenseRestoreArray(ds->omat[mA],&A));
85 3727 : PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
86 3727 : 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 4170 : PetscErrorCode DSSort_NHEP_Total(DS ds,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
94 : {
95 4170 : PetscScalar re,*T,*Q;
96 4170 : PetscInt i,j,pos,result;
97 4170 : PetscBLASInt ifst,ilst,info,n,ld;
98 : #if !defined(PETSC_USE_COMPLEX)
99 : PetscScalar *work,im;
100 : #endif
101 :
102 4170 : PetscFunctionBegin;
103 4170 : PetscCall(MatDenseGetArray(ds->omat[mT],&T));
104 4170 : PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
105 4170 : PetscCall(PetscBLASIntCast(ds->n,&n));
106 4170 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
107 : #if !defined(PETSC_USE_COMPLEX)
108 : PetscCall(DSAllocateWork_Private(ds,ld,0,0));
109 : work = ds->work;
110 : #endif
111 : /* selection sort */
112 69515 : for (i=ds->l;i<n-1;i++) {
113 65345 : re = wr[i];
114 : #if !defined(PETSC_USE_COMPLEX)
115 : im = wi[i];
116 : #endif
117 65345 : pos = 0;
118 65345 : j=i+1; /* j points to the next eigenvalue */
119 : #if !defined(PETSC_USE_COMPLEX)
120 : if (im != 0) j=i+2;
121 : #endif
122 : /* find minimum eigenvalue */
123 792529 : for (;j<n;j++) {
124 : #if !defined(PETSC_USE_COMPLEX)
125 : PetscCall(SlepcSCCompare(ds->sc,re,im,wr[j],wi[j],&result));
126 : #else
127 727184 : PetscCall(SlepcSCCompare(ds->sc,re,0.0,wr[j],0.0,&result));
128 : #endif
129 727184 : if (result > 0) {
130 266936 : re = wr[j];
131 : #if !defined(PETSC_USE_COMPLEX)
132 : im = wi[j];
133 : #endif
134 266936 : pos = j;
135 : }
136 : #if !defined(PETSC_USE_COMPLEX)
137 : if (wi[j] != 0) j++;
138 : #endif
139 : }
140 65345 : if (pos) {
141 : /* interchange blocks */
142 44409 : PetscCall(PetscBLASIntCast(pos+1,&ifst));
143 44409 : PetscCall(PetscBLASIntCast(i+1,&ilst));
144 : #if !defined(PETSC_USE_COMPLEX)
145 : PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
146 : #else
147 44409 : PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
148 : #endif
149 44409 : SlepcCheckLapackInfo("trexc",info);
150 : /* recover original eigenvalues from T matrix */
151 607967 : for (j=i;j<n;j++) {
152 563558 : wr[j] = T[j+j*ld];
153 : #if !defined(PETSC_USE_COMPLEX)
154 : if (j<n-1 && T[j+1+j*ld] != 0.0) {
155 : /* complex conjugate eigenvalue */
156 : wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
157 : wr[j+1] = wr[j];
158 : wi[j+1] = -wi[j];
159 : j++;
160 : } else wi[j] = 0.0;
161 : #endif
162 : }
163 : }
164 : #if !defined(PETSC_USE_COMPLEX)
165 : if (wi[i] != 0) i++;
166 : #endif
167 : }
168 4170 : PetscCall(MatDenseRestoreArray(ds->omat[mT],&T));
169 4170 : PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
170 4170 : 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 10 : PetscErrorCode DSSortWithPermutation_NHEP_Private(DS ds,PetscInt *perm,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
178 : {
179 10 : PetscInt i,j,pos,inc=1;
180 10 : PetscBLASInt ifst,ilst,info,n,ld;
181 10 : PetscScalar *T,*Q;
182 : #if !defined(PETSC_USE_COMPLEX)
183 : PetscScalar *work;
184 : #endif
185 :
186 10 : PetscFunctionBegin;
187 10 : PetscCall(MatDenseGetArray(ds->omat[mT],&T));
188 10 : PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
189 10 : PetscCall(PetscBLASIntCast(ds->n,&n));
190 10 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
191 : #if !defined(PETSC_USE_COMPLEX)
192 : PetscCall(DSAllocateWork_Private(ds,ld,0,0));
193 : work = ds->work;
194 : #endif
195 181 : for (i=ds->l;i<n-1;i++) {
196 171 : pos = perm[i];
197 : #if !defined(PETSC_USE_COMPLEX)
198 : inc = (pos<n-1 && T[pos+1+pos*ld] != 0.0)? 2: 1;
199 : #endif
200 171 : if (pos!=i) {
201 : #if !defined(PETSC_USE_COMPLEX)
202 : 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 15 : PetscCall(PetscBLASIntCast(pos+1,&ifst));
206 15 : PetscCall(PetscBLASIntCast(i+1,&ilst));
207 : #if !defined(PETSC_USE_COMPLEX)
208 : PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
209 : #else
210 15 : PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
211 : #endif
212 15 : SlepcCheckLapackInfo("trexc",info);
213 206 : for (j=i+1;j<n;j++) {
214 191 : if (perm[j]>=i && perm[j]<pos) perm[j]+=inc;
215 : }
216 15 : perm[i] = i;
217 15 : if (inc==2) perm[i+1] = i+1;
218 : }
219 171 : if (inc==2) i++;
220 : }
221 : /* recover original eigenvalues from T matrix */
222 191 : for (j=ds->l;j<n;j++) {
223 181 : wr[j] = T[j+j*ld];
224 : #if !defined(PETSC_USE_COMPLEX)
225 : if (j<n-1 && T[j+1+j*ld] != 0.0) {
226 : /* complex conjugate eigenvalue */
227 : wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
228 : wr[j+1] = wr[j];
229 : wi[j+1] = -wi[j];
230 : j++;
231 : } else wi[j] = 0.0;
232 : #endif
233 : }
234 10 : PetscCall(MatDenseRestoreArray(ds->omat[mT],&T));
235 10 : PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
236 10 : PetscFunctionReturn(PETSC_SUCCESS);
237 : }
|
