Actual source code: dsutil.c
slepc-3.21.1 2024-04-26
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,info,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: PetscCallBLAS("LAPACKgehrd",LAPACKgehrd_(&n,&ilo,&n,A,&ld,tau,work,&lwork,&info));
51: SlepcCheckLapackInfo("gehrd",info);
52: for (j=0;j<n-1;j++) {
53: for (i=j+2;i<n;i++) {
54: Q[i+j*ld] = A[i+j*ld];
55: A[i+j*ld] = 0.0;
56: }
57: }
58: PetscCallBLAS("LAPACKorghr",LAPACKorghr_(&n,&ilo,&n,Q,&ld,tau,work,&lwork,&info));
59: SlepcCheckLapackInfo("orghr",info);
60: }
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: 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: SlepcCheckLapackInfo("hseqr",info);
84: PetscCall(MatDenseRestoreArray(ds->omat[mA],&A));
85: PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
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: PetscErrorCode DSSort_NHEP_Total(DS ds,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
94: {
95: PetscScalar re,*T,*Q;
96: PetscInt i,j,pos,result;
97: PetscBLASInt ifst,ilst,info,n,ld;
98: #if !defined(PETSC_USE_COMPLEX)
99: PetscScalar *work,im;
100: #endif
102: PetscFunctionBegin;
103: PetscCall(MatDenseGetArray(ds->omat[mT],&T));
104: PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
105: PetscCall(PetscBLASIntCast(ds->n,&n));
106: 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: for (i=ds->l;i<n-1;i++) {
113: re = wr[i];
114: #if !defined(PETSC_USE_COMPLEX)
115: im = wi[i];
116: #endif
117: pos = 0;
118: 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: 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: PetscCall(SlepcSCCompare(ds->sc,re,0.0,wr[j],0.0,&result));
128: #endif
129: if (result > 0) {
130: re = wr[j];
131: #if !defined(PETSC_USE_COMPLEX)
132: im = wi[j];
133: #endif
134: pos = j;
135: }
136: #if !defined(PETSC_USE_COMPLEX)
137: if (wi[j] != 0) j++;
138: #endif
139: }
140: if (pos) {
141: /* interchange blocks */
142: PetscCall(PetscBLASIntCast(pos+1,&ifst));
143: 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: PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
148: #endif
149: SlepcCheckLapackInfo("trexc",info);
150: /* recover original eigenvalues from T matrix */
151: for (j=i;j<n;j++) {
152: 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: PetscCall(MatDenseRestoreArray(ds->omat[mT],&T));
169: PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
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: PetscErrorCode DSSortWithPermutation_NHEP_Private(DS ds,PetscInt *perm,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
178: {
179: PetscInt i,j,pos,inc=1;
180: PetscBLASInt ifst,ilst,info,n,ld;
181: PetscScalar *T,*Q;
182: #if !defined(PETSC_USE_COMPLEX)
183: PetscScalar *work;
184: #endif
186: PetscFunctionBegin;
187: PetscCall(MatDenseGetArray(ds->omat[mT],&T));
188: PetscCall(MatDenseGetArray(ds->omat[mQ],&Q));
189: PetscCall(PetscBLASIntCast(ds->n,&n));
190: 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: for (i=ds->l;i<n-1;i++) {
196: 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: 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: PetscCall(PetscBLASIntCast(pos+1,&ifst));
206: 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: PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
211: #endif
212: SlepcCheckLapackInfo("trexc",info);
213: for (j=i+1;j<n;j++) {
214: if (perm[j]>=i && perm[j]<pos) perm[j]+=inc;
215: }
216: perm[i] = i;
217: if (inc==2) perm[i+1] = i+1;
218: }
219: if (inc==2) i++;
220: }
221: /* recover original eigenvalues from T matrix */
222: for (j=ds->l;j<n;j++) {
223: 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: PetscCall(MatDenseRestoreArray(ds->omat[mT],&T));
235: PetscCall(MatDenseRestoreArray(ds->omat[mQ],&Q));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }