Actual source code: dsutil.c
slepc-3.15.0 2021-03-31
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: Utility subroutines common to several impls
12: */
14: #include <slepc/private/dsimpl.h>
15: #include <slepcblaslapack.h>
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: PetscErrorCode DSSolve_NHEP_Private(DS ds,PetscScalar *A,PetscScalar *Q,PetscScalar *wr,PetscScalar *wi)
22: {
24: PetscScalar *work,*tau;
25: PetscInt i,j;
26: PetscBLASInt ilo,lwork,info,n,k,ld;
29: PetscBLASIntCast(ds->n,&n);
30: PetscBLASIntCast(ds->ld,&ld);
31: PetscBLASIntCast(ds->l+1,&ilo);
32: PetscBLASIntCast(ds->k,&k);
33: DSAllocateWork_Private(ds,ld+6*ld,0,0);
34: tau = ds->work;
35: work = ds->work+ld;
36: lwork = 6*ld;
38: /* initialize orthogonal matrix */
39: PetscArrayzero(Q,ld*ld);
40: for (i=0;i<n;i++) Q[i+i*ld] = 1.0;
41: if (n==1) { /* quick return */
42: wr[0] = A[0];
43: if (wi) wi[0] = 0.0;
44: return(0);
45: }
47: /* reduce to upper Hessenberg form */
48: if (ds->state<DS_STATE_INTERMEDIATE) {
49: PetscStackCallBLAS("LAPACKgehrd",LAPACKgehrd_(&n,&ilo,&n,A,&ld,tau,work,&lwork,&info));
50: SlepcCheckLapackInfo("gehrd",info);
51: for (j=0;j<n-1;j++) {
52: for (i=j+2;i<n;i++) {
53: Q[i+j*ld] = A[i+j*ld];
54: A[i+j*ld] = 0.0;
55: }
56: }
57: PetscStackCallBLAS("LAPACKorghr",LAPACKorghr_(&n,&ilo,&n,Q,&ld,tau,work,&lwork,&info));
58: SlepcCheckLapackInfo("orghr",info);
59: }
61: /* compute the (real) Schur form */
62: #if !defined(PETSC_USE_COMPLEX)
63: PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,wi,Q,&ld,work,&lwork,&info));
64: for (j=0;j<ds->l;j++) {
65: if (j==n-1 || A[j+1+j*ld] == 0.0) {
66: /* real eigenvalue */
67: wr[j] = A[j+j*ld];
68: wi[j] = 0.0;
69: } else {
70: /* complex eigenvalue */
71: wr[j] = A[j+j*ld];
72: wr[j+1] = A[j+j*ld];
73: wi[j] = PetscSqrtReal(PetscAbsReal(A[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(A[j+(j+1)*ld]));
74: wi[j+1] = -wi[j];
75: j++;
76: }
77: }
78: #else
79: PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,Q,&ld,work,&lwork,&info));
80: if (wi) for (i=ds->l;i<n;i++) wi[i] = 0.0;
81: #endif
82: SlepcCheckLapackInfo("hseqr",info);
83: return(0);
84: }
86: /*
87: Sort a Schur form represented by the (quasi-)triangular matrix T and
88: the unitary matrix Q, and return the sorted eigenvalues in wr,wi
89: */
90: PetscErrorCode DSSort_NHEP_Total(DS ds,PetscScalar *T,PetscScalar *Q,PetscScalar *wr,PetscScalar *wi)
91: {
93: PetscScalar re;
94: PetscInt i,j,pos,result;
95: PetscBLASInt ifst,ilst,info,n,ld;
96: #if !defined(PETSC_USE_COMPLEX)
97: PetscScalar *work,im;
98: #endif
101: PetscBLASIntCast(ds->n,&n);
102: PetscBLASIntCast(ds->ld,&ld);
103: #if !defined(PETSC_USE_COMPLEX)
104: DSAllocateWork_Private(ds,ld,0,0);
105: work = ds->work;
106: #endif
107: /* selection sort */
108: for (i=ds->l;i<n-1;i++) {
109: re = wr[i];
110: #if !defined(PETSC_USE_COMPLEX)
111: im = wi[i];
112: #endif
113: pos = 0;
114: j=i+1; /* j points to the next eigenvalue */
115: #if !defined(PETSC_USE_COMPLEX)
116: if (im != 0) j=i+2;
117: #endif
118: /* find minimum eigenvalue */
119: for (;j<n;j++) {
120: #if !defined(PETSC_USE_COMPLEX)
121: SlepcSCCompare(ds->sc,re,im,wr[j],wi[j],&result);
122: #else
123: SlepcSCCompare(ds->sc,re,0.0,wr[j],0.0,&result);
124: #endif
125: if (result > 0) {
126: re = wr[j];
127: #if !defined(PETSC_USE_COMPLEX)
128: im = wi[j];
129: #endif
130: pos = j;
131: }
132: #if !defined(PETSC_USE_COMPLEX)
133: if (wi[j] != 0) j++;
134: #endif
135: }
136: if (pos) {
137: /* interchange blocks */
138: PetscBLASIntCast(pos+1,&ifst);
139: PetscBLASIntCast(i+1,&ilst);
140: #if !defined(PETSC_USE_COMPLEX)
141: PetscStackCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
142: #else
143: PetscStackCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
144: #endif
145: SlepcCheckLapackInfo("trexc",info);
146: /* recover original eigenvalues from T matrix */
147: for (j=i;j<n;j++) {
148: wr[j] = T[j+j*ld];
149: #if !defined(PETSC_USE_COMPLEX)
150: if (j<n-1 && T[j+1+j*ld] != 0.0) {
151: /* complex conjugate eigenvalue */
152: wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
153: wr[j+1] = wr[j];
154: wi[j+1] = -wi[j];
155: j++;
156: } else wi[j] = 0.0;
157: #endif
158: }
159: }
160: #if !defined(PETSC_USE_COMPLEX)
161: if (wi[i] != 0) i++;
162: #endif
163: }
164: return(0);
165: }
167: /*
168: Reorder a Schur form represented by T,Q according to a permutation perm,
169: and return the sorted eigenvalues in wr,wi
170: */
171: PetscErrorCode DSSortWithPermutation_NHEP_Private(DS ds,PetscInt *perm,PetscScalar *T,PetscScalar *Q,PetscScalar *wr,PetscScalar *wi)
172: {
174: PetscInt i,j,pos,inc=1;
175: PetscBLASInt ifst,ilst,info,n,ld;
176: #if !defined(PETSC_USE_COMPLEX)
177: PetscScalar *work;
178: #endif
181: PetscBLASIntCast(ds->n,&n);
182: PetscBLASIntCast(ds->ld,&ld);
183: #if !defined(PETSC_USE_COMPLEX)
184: DSAllocateWork_Private(ds,ld,0,0);
185: work = ds->work;
186: #endif
187: for (i=ds->l;i<n-1;i++) {
188: pos = perm[i];
189: #if !defined(PETSC_USE_COMPLEX)
190: inc = (pos<n-1 && T[pos+1+pos*ld] != 0.0)? 2: 1;
191: #endif
192: if (pos!=i) {
193: #if !defined(PETSC_USE_COMPLEX)
194: if ((T[pos+(pos-1)*ld] != 0.0 && perm[i+1]!=pos-1) || (T[pos+1+pos*ld] != 0.0 && perm[i+1]!=pos+1))
195: SETERRQ1(PETSC_COMM_SELF,1,"Invalid permutation due to a 2x2 block at position %D",pos);
196: #endif
197: /* interchange blocks */
198: PetscBLASIntCast(pos+1,&ifst);
199: PetscBLASIntCast(i+1,&ilst);
200: #if !defined(PETSC_USE_COMPLEX)
201: PetscStackCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
202: #else
203: PetscStackCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
204: #endif
205: SlepcCheckLapackInfo("trexc",info);
206: for (j=i+1;j<n;j++) {
207: if (perm[j]>=i && perm[j]<pos) perm[j]+=inc;
208: }
209: perm[i] = i;
210: if (inc==2) perm[i+1] = i+1;
211: }
212: if (inc==2) i++;
213: }
214: /* recover original eigenvalues from T matrix */
215: for (j=ds->l;j<n;j++) {
216: wr[j] = T[j+j*ld];
217: #if !defined(PETSC_USE_COMPLEX)
218: if (j<n-1 && T[j+1+j*ld] != 0.0) {
219: /* complex conjugate eigenvalue */
220: wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
221: wr[j+1] = wr[j];
222: wi[j+1] = -wi[j];
223: j++;
224: } else wi[j] = 0.0;
225: #endif
226: }
227: return(0);
228: }