Actual source code: dsutil.c

slepc-3.18.0 2022-10-01
Report Typos and Errors
  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:   MatDenseGetArray(ds->omat[mA],&A);
 28:   MatDenseGetArray(ds->omat[mQ],&Q);
 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:     PetscCallBLAS("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:     PetscCallBLAS("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:   PetscCallBLAS("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:   PetscCallBLAS("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:   MatDenseRestoreArray(ds->omat[mA],&A);
 84:   MatDenseRestoreArray(ds->omat[mQ],&Q);
 85:   return 0;
 86: }

 88: /*
 89:    Sort a Schur form represented by the (quasi-)triangular matrix T and
 90:    the unitary matrix Q, and return the sorted eigenvalues in wr,wi
 91: */
 92: PetscErrorCode DSSort_NHEP_Total(DS ds,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
 93: {
 94:   PetscScalar    re,*T,*Q;
 95:   PetscInt       i,j,pos,result;
 96:   PetscBLASInt   ifst,ilst,info,n,ld;
 97: #if !defined(PETSC_USE_COMPLEX)
 98:   PetscScalar    *work,im;
 99: #endif

101:   MatDenseGetArray(ds->omat[mT],&T);
102:   MatDenseGetArray(ds->omat[mQ],&Q);
103:   PetscBLASIntCast(ds->n,&n);
104:   PetscBLASIntCast(ds->ld,&ld);
105: #if !defined(PETSC_USE_COMPLEX)
106:   DSAllocateWork_Private(ds,ld,0,0);
107:   work = ds->work;
108: #endif
109:   /* selection sort */
110:   for (i=ds->l;i<n-1;i++) {
111:     re = wr[i];
112: #if !defined(PETSC_USE_COMPLEX)
113:     im = wi[i];
114: #endif
115:     pos = 0;
116:     j=i+1; /* j points to the next eigenvalue */
117: #if !defined(PETSC_USE_COMPLEX)
118:     if (im != 0) j=i+2;
119: #endif
120:     /* find minimum eigenvalue */
121:     for (;j<n;j++) {
122: #if !defined(PETSC_USE_COMPLEX)
123:       SlepcSCCompare(ds->sc,re,im,wr[j],wi[j],&result);
124: #else
125:       SlepcSCCompare(ds->sc,re,0.0,wr[j],0.0,&result);
126: #endif
127:       if (result > 0) {
128:         re = wr[j];
129: #if !defined(PETSC_USE_COMPLEX)
130:         im = wi[j];
131: #endif
132:         pos = j;
133:       }
134: #if !defined(PETSC_USE_COMPLEX)
135:       if (wi[j] != 0) j++;
136: #endif
137:     }
138:     if (pos) {
139:       /* interchange blocks */
140:       PetscBLASIntCast(pos+1,&ifst);
141:       PetscBLASIntCast(i+1,&ilst);
142: #if !defined(PETSC_USE_COMPLEX)
143:       PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
144: #else
145:       PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
146: #endif
147:       SlepcCheckLapackInfo("trexc",info);
148:       /* recover original eigenvalues from T matrix */
149:       for (j=i;j<n;j++) {
150:         wr[j] = T[j+j*ld];
151: #if !defined(PETSC_USE_COMPLEX)
152:         if (j<n-1 && T[j+1+j*ld] != 0.0) {
153:           /* complex conjugate eigenvalue */
154:           wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
155:           wr[j+1] = wr[j];
156:           wi[j+1] = -wi[j];
157:           j++;
158:         } else wi[j] = 0.0;
159: #endif
160:       }
161:     }
162: #if !defined(PETSC_USE_COMPLEX)
163:     if (wi[i] != 0) i++;
164: #endif
165:   }
166:   MatDenseRestoreArray(ds->omat[mT],&T);
167:   MatDenseRestoreArray(ds->omat[mQ],&Q);
168:   return 0;
169: }

171: /*
172:    Reorder a Schur form represented by T,Q according to a permutation perm,
173:    and return the sorted eigenvalues in wr,wi
174: */
175: PetscErrorCode DSSortWithPermutation_NHEP_Private(DS ds,PetscInt *perm,DSMatType mT,DSMatType mQ,PetscScalar *wr,PetscScalar *wi)
176: {
177:   PetscInt       i,j,pos,inc=1;
178:   PetscBLASInt   ifst,ilst,info,n,ld;
179:   PetscScalar    *T,*Q;
180: #if !defined(PETSC_USE_COMPLEX)
181:   PetscScalar    *work;
182: #endif

184:   MatDenseGetArray(ds->omat[mT],&T);
185:   MatDenseGetArray(ds->omat[mQ],&Q);
186:   PetscBLASIntCast(ds->n,&n);
187:   PetscBLASIntCast(ds->ld,&ld);
188: #if !defined(PETSC_USE_COMPLEX)
189:   DSAllocateWork_Private(ds,ld,0,0);
190:   work = ds->work;
191: #endif
192:   for (i=ds->l;i<n-1;i++) {
193:     pos = perm[i];
194: #if !defined(PETSC_USE_COMPLEX)
195:     inc = (pos<n-1 && T[pos+1+pos*ld] != 0.0)? 2: 1;
196: #endif
197:     if (pos!=i) {
198: #if !defined(PETSC_USE_COMPLEX)
200: #endif
201:       /* interchange blocks */
202:       PetscBLASIntCast(pos+1,&ifst);
203:       PetscBLASIntCast(i+1,&ilst);
204: #if !defined(PETSC_USE_COMPLEX)
205:       PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
206: #else
207:       PetscCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
208: #endif
209:       SlepcCheckLapackInfo("trexc",info);
210:       for (j=i+1;j<n;j++) {
211:         if (perm[j]>=i && perm[j]<pos) perm[j]+=inc;
212:       }
213:       perm[i] = i;
214:       if (inc==2) perm[i+1] = i+1;
215:     }
216:     if (inc==2) i++;
217:   }
218:   /* recover original eigenvalues from T matrix */
219:   for (j=ds->l;j<n;j++) {
220:     wr[j] = T[j+j*ld];
221: #if !defined(PETSC_USE_COMPLEX)
222:     if (j<n-1 && T[j+1+j*ld] != 0.0) {
223:       /* complex conjugate eigenvalue */
224:       wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld]))*PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
225:       wr[j+1] = wr[j];
226:       wi[j+1] = -wi[j];
227:       j++;
228:     } else wi[j] = 0.0;
229: #endif
230:   }
231:   MatDenseRestoreArray(ds->omat[mT],&T);
232:   MatDenseRestoreArray(ds->omat[mQ],&Q);
233:   return 0;
234: }