Actual source code: dsutil.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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,DSMatType mA,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
22: {
23: PetscScalar *work,*tau,*A,*Q;
24: PetscInt i,j;
25: PetscBLASInt ilo,lwork,n,k,ld;
27: PetscFunctionBegin;
28: PetscCall(MatDenseGetArray(ds->omat[mA],&A));
29: PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
30: PetscCall(PetscBLASIntCast(ds->n,&n));
31: PetscCall(PetscBLASIntCast(ds->ld,&ld));
32: PetscCall(PetscBLASIntCast(ds->l+1,&ilo));
33: PetscCall(PetscBLASIntCast(ds->k,&k));
34: PetscCall(DSAllocateWork_Private(ds,ld+6*ld,0,0));
35: tau = ds->work;
36: work = ds->work+ld;
37: lwork = 6*ld;
39: /* initialize orthogonal matrix */
40: PetscCall(PetscArrayzero(Q,ld*ld));
41: for (i=0;i<n;i++) Q[i+i*ld] = 1.0;
42: if (n==1) { /* quick return */
43: wr[0] = A[0];
44: if (wi) wi[0] = 0.0;
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
48: /* reduce to upper Hessenberg form */
49: if (ds->state<DS_STATE_INTERMEDIATE) {
50: PetscCallLAPACKInfo("LAPACKgehrd",LAPACKgehrd_(&n,&ilo,&n,A,&ld,tau,work,&lwork,&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: PetscCallLAPACKInfo("LAPACKorghr",LAPACKorghr_(&n,&ilo,&n,Q,&ld,tau,work,&lwork,&info));
58: }
60: /* compute the (real) Schur form */
61: #if !PetscDefined(USE_COMPLEX)
62: PetscCallLAPACKInfo("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,wi,Q,&ld,work,&lwork,&info));
63: for (j=0;j<ds->l;j++) {
64: if (j==n-1 || A[j+1+j*ld] == 0.0) {
65: /* real eigenvalue */
66: wr[j] = A[j+j*ld];
67: wi[j] = 0.0;
68: } else {
69: /* complex eigenvalue */
70: wr[j] = A[j+j*ld];
71: wr[j+1] = A[j+j*ld];
72: wi[j] = PetscSqrtReal(PetscAbsReal(A[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(A[j+(j+1)*ld]));
73: wi[j+1] = -wi[j];
74: j++;
75: }
76: }
77: #else
78: PetscCallLAPACKInfo("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,Q,&ld,work,&lwork,&info));
79: if (wi) for (i=ds->l;i<n;i++) wi[i] = 0.0;
80: #endif
81: PetscCall(MatDenseRestoreArray(ds->omat[mA],&A));
82: PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
83: PetscFunctionReturn(PETSC_SUCCESS);
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,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
91: {
92: PetscScalar re,*T,*Q;
93: PetscInt i,j,pos,result;
94: PetscBLASInt ifst,ilst,n,ld;
95: #if !PetscDefined(USE_COMPLEX)
96: PetscScalar *work,im;
97: #endif
99: PetscFunctionBegin;
100: PetscCall(MatDenseGetArray(ds->omat[mT],&T));
101: PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
102: PetscCall(PetscBLASIntCast(ds->n,&n));
103: PetscCall(PetscBLASIntCast(ds->ld,&ld));
104: #if !PetscDefined(USE_COMPLEX)
105: PetscCall(DSAllocateWork_Private(ds,ld,0,0));
106: work = ds->work;
107: #endif
108: /* selection sort */
109: for (i=ds->l;i<n-1;i++) {
110: re = wr[i];
111: #if !PetscDefined(USE_COMPLEX)
112: im = wi[i];
113: #endif
114: pos = 0;
115: j=i+1; /* j points to the next eigenvalue */
116: #if !PetscDefined(USE_COMPLEX)
117: if (im != 0) j=i+2;
118: #endif
119: /* find minimum eigenvalue */
120: for (;j<n;j++) {
121: #if !PetscDefined(USE_COMPLEX)
122: PetscCall(SlepcSCCompare(ds->sc,re,im,wr[j],wi[j],&result));
123: #else
124: PetscCall(SlepcSCCompare(ds->sc,re,0.0,wr[j],0.0,&result));
125: #endif
126: if (result > 0) {
127: re = wr[j];
128: #if !PetscDefined(USE_COMPLEX)
129: im = wi[j];
130: #endif
131: pos = j;
132: }
133: #if !PetscDefined(USE_COMPLEX)
134: if (wi[j] != 0) j++;
135: #endif
136: }
137: if (pos) {
138: /* interchange blocks */
139: PetscCall(PetscBLASIntCast(pos+1,&ifst));
140: PetscCall(PetscBLASIntCast(i+1,&ilst));
141: #if !PetscDefined(USE_COMPLEX)
142: PetscCallLAPACKInfo("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
143: #else
144: PetscCallLAPACKInfo("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
145: #endif
146: /* recover original eigenvalues from T matrix */
147: for (j=i;j<n;j++) {
148: wr[j] = T[j+j*ld];
149: #if !PetscDefined(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 !PetscDefined(USE_COMPLEX)
161: if (wi[i] != 0) i++;
162: #endif
163: }
164: PetscCall(MatDenseRestoreArray(ds->omat[mT],&T));
165: PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: /*
170: Reorder a Schur form represented by T,Q according to a permutation perm,
171: and return the sorted eigenvalues in wr,wi
172: */
173: PetscErrorCode DSSortWithPermutation_NHEP_Private(DS ds,PetscInt *perm,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
174: {
175: PetscInt i,j,pos,inc=1;
176: PetscBLASInt ifst,ilst,n,ld;
177: PetscScalar *T,*Q;
178: #if !PetscDefined(USE_COMPLEX)
179: PetscScalar *work;
180: #endif
182: PetscFunctionBegin;
183: PetscCall(MatDenseGetArray(ds->omat[mT],&T));
184: PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
185: PetscCall(PetscBLASIntCast(ds->n,&n));
186: PetscCall(PetscBLASIntCast(ds->ld,&ld));
187: #if !PetscDefined(USE_COMPLEX)
188: PetscCall(DSAllocateWork_Private(ds,ld,0,0));
189: work = ds->work;
190: #endif
191: for (i=ds->l;i<n-1;i++) {
192: pos = perm[i];
193: #if !PetscDefined(USE_COMPLEX)
194: inc = (pos<n-1 && T[pos+1+pos*ld] != 0.0)? 2: 1;
195: #endif
196: if (pos!=i) {
197: #if !PetscDefined(USE_COMPLEX)
198: PetscCheck((T[pos+(pos-1)*ld]==0.0 || perm[i+1]==pos-1) && (pos==n-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);
199: #endif
200: /* interchange blocks */
201: PetscCall(PetscBLASIntCast(pos+1,&ifst));
202: PetscCall(PetscBLASIntCast(i+1,&ilst));
203: #if !PetscDefined(USE_COMPLEX)
204: PetscCallLAPACKInfo("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
205: #else
206: PetscCallLAPACKInfo("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
207: #endif
208: for (j=i+1;j<n;j++) {
209: if (perm[j]>=i && perm[j]<pos) perm[j]+=inc;
210: }
211: perm[i] = i;
212: if (inc==2) perm[i+1] = i+1;
213: }
214: if (inc==2) i++;
215: }
216: /* recover original eigenvalues from T matrix */
217: for (j=ds->l;j<n;j++) {
218: wr[j] = T[j+j*ld];
219: #if !PetscDefined(USE_COMPLEX)
220: if (j<n-1 && T[j+1+j*ld] != 0.0) {
221: /* complex conjugate eigenvalue */
222: wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
223: wr[j+1] = wr[j];
224: wi[j+1] = -wi[j];
225: j++;
226: } else wi[j] = 0.0;
227: #endif
228: }
229: PetscCall(MatDenseRestoreArray(ds->omat[mT],&T));
230: PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }