Actual source code: dsghiep.c

slepc-main 2024-11-22
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: */

 11: #include <slepc/private/dsimpl.h>
 12: #include <slepcblaslapack.h>

 14: static PetscErrorCode DSAllocate_GHIEP(DS ds,PetscInt ld)
 15: {
 16:   PetscFunctionBegin;
 17:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
 18:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
 19:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q));
 20:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
 21:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
 22:   PetscCall(PetscFree(ds->perm));
 23:   PetscCall(PetscMalloc1(ld,&ds->perm));
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

 27: PetscErrorCode DSSwitchFormat_GHIEP(DS ds,PetscBool tocompact)
 28: {
 29:   PetscReal      *T,*S;
 30:   PetscScalar    *A,*B;
 31:   PetscInt       i,n,ld;

 33:   PetscFunctionBegin;
 34:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
 35:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
 36:   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
 37:   PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
 38:   n = ds->n;
 39:   ld = ds->ld;
 40:   if (tocompact) { /* switch from dense (arrow) to compact storage */
 41:     PetscCall(PetscArrayzero(T,n));
 42:     PetscCall(PetscArrayzero(T+ld,n));
 43:     PetscCall(PetscArrayzero(T+2*ld,n));
 44:     PetscCall(PetscArrayzero(S,n));
 45:     for (i=0;i<n-1;i++) {
 46:       T[i]    = PetscRealPart(A[i+i*ld]);
 47:       T[ld+i] = PetscRealPart(A[i+1+i*ld]);
 48:       S[i]    = PetscRealPart(B[i+i*ld]);
 49:     }
 50:     T[n-1] = PetscRealPart(A[n-1+(n-1)*ld]);
 51:     S[n-1] = PetscRealPart(B[n-1+(n-1)*ld]);
 52:     for (i=ds->l;i<ds->k;i++) T[2*ld+i] = PetscRealPart(A[ds->k+i*ld]);
 53:   } else { /* switch from compact (arrow) to dense storage */
 54:     for (i=0;i<n;i++) {
 55:       PetscCall(PetscArrayzero(A+i*ld,n));
 56:       PetscCall(PetscArrayzero(B+i*ld,n));
 57:     }
 58:     for (i=0;i<n-1;i++) {
 59:       A[i+i*ld]     = T[i];
 60:       A[i+1+i*ld]   = T[ld+i];
 61:       A[i+(i+1)*ld] = T[ld+i];
 62:       B[i+i*ld]     = S[i];
 63:     }
 64:     A[n-1+(n-1)*ld] = T[n-1];
 65:     B[n-1+(n-1)*ld] = S[n-1];
 66:     for (i=ds->l;i<ds->k;i++) {
 67:       A[ds->k+i*ld] = T[2*ld+i];
 68:       A[i+ds->k*ld] = T[2*ld+i];
 69:     }
 70:   }
 71:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
 72:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
 73:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
 74:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }

 78: static PetscErrorCode DSView_GHIEP(DS ds,PetscViewer viewer)
 79: {
 80:   PetscViewerFormat format;
 81:   PetscInt          i,j;
 82:   PetscReal         *T,*S,value;
 83:   const char        *methodname[] = {
 84:                      "QR + Inverse Iteration",
 85:                      "HZ method",
 86:                      "QR"
 87:   };
 88:   const int         nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);

 90:   PetscFunctionBegin;
 91:   PetscCall(PetscViewerGetFormat(viewer,&format));
 92:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
 93:     if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
 94:     PetscFunctionReturn(PETSC_SUCCESS);
 95:   }
 96:   if (ds->compact) {
 97:     PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
 98:     PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
 99:     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
100:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
101:       PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",ds->n,ds->n));
102:       PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",3*ds->n));
103:       PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
104:       for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "  %18.16e\n",i+1,i+1,(double)T[i]));
105:       for (i=0;i<ds->n-1;i++) {
106:         if (T[i+ds->ld] !=0 && i!=ds->k-1) {
107:           PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "  %18.16e\n",i+2,i+1,(double)T[i+ds->ld]));
108:           PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "  %18.16e\n",i+1,i+2,(double)T[i+ds->ld]));
109:         }
110:       }
111:       for (i = ds->l;i<ds->k;i++) {
112:         if (T[i+2*ds->ld]) {
113:           PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "  %18.16e\n",ds->k+1,i+1,(double)T[i+2*ds->ld]));
114:           PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "  %18.16e\n",i+1,ds->k+1,(double)T[i+2*ds->ld]));
115:         }
116:       }
117:       PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_A]));

119:       PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",ds->n,ds->n));
120:       PetscCall(PetscViewerASCIIPrintf(viewer,"omega = zeros(%" PetscInt_FMT ",3);\n",3*ds->n));
121:       PetscCall(PetscViewerASCIIPrintf(viewer,"omega = [\n"));
122:       for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "  %18.16e\n",i+1,i+1,(double)S[i]));
123:       PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(omega);\n",DSMatName[DS_MAT_B]));

125:     } else {
126:       PetscCall(PetscViewerASCIIPrintf(viewer,"T\n"));
127:       for (i=0;i<ds->n;i++) {
128:         for (j=0;j<ds->n;j++) {
129:           if (i==j) value = T[i];
130:           else if (i==j+1 || j==i+1) value = T[PetscMin(i,j)+ds->ld];
131:           else if ((i<ds->k && j==ds->k) || (i==ds->k && j<ds->k)) value = T[PetscMin(i,j)+2*ds->ld];
132:           else value = 0.0;
133:           PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
134:         }
135:         PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
136:       }
137:       PetscCall(PetscViewerASCIIPrintf(viewer,"omega\n"));
138:       for (i=0;i<ds->n;i++) {
139:         for (j=0;j<ds->n;j++) {
140:           if (i==j) value = S[i];
141:           else value = 0.0;
142:           PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
143:         }
144:         PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
145:       }
146:     }
147:     PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
148:     PetscCall(PetscViewerFlush(viewer));
149:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
150:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
151:   } else {
152:     PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
153:     PetscCall(DSViewMat(ds,viewer,DS_MAT_B));
154:   }
155:   if (ds->state>DS_STATE_INTERMEDIATE) PetscCall(DSViewMat(ds,viewer,DS_MAT_Q));
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

159: static PetscErrorCode DSVectors_GHIEP_Eigen_Some(DS ds,PetscInt *idx,PetscReal *rnorm)
160: {
161:   PetscReal         b[4],M[4],*T,*S,d1,d2,s1,s2,e,scal1,scal2,wr1,wr2,wi,ep,norm;
162:   PetscScalar       *X,Y[4],alpha,szero=0.0;
163:   const PetscScalar *A,*B,*Q;
164:   PetscInt          k;
165:   PetscBLASInt      two=2,n_,ld,one=1;
166: #if !defined(PETSC_USE_COMPLEX)
167:   PetscBLASInt      four=4;
168: #endif

170:   PetscFunctionBegin;
171:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
172:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
173:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
174:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
175:   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
176:   PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
177:   k = *idx;
178:   PetscCall(PetscBLASIntCast(ds->n,&n_));
179:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
180:   if (k < ds->n-1) e = (ds->compact)?T[k+ld]:PetscRealPart(A[(k+1)+ld*k]);
181:   else e = 0.0;
182:   if (e == 0.0) { /* Real */
183:     if (ds->state>=DS_STATE_CONDENSED) PetscCall(PetscArraycpy(X+k*ld,Q+k*ld,ld));
184:     else {
185:       PetscCall(PetscArrayzero(X+k*ds->ld,ds->ld));
186:       X[k+k*ds->ld] = 1.0;
187:     }
188:     if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
189:   } else { /* 2x2 block */
190:     if (ds->compact) {
191:       s1 = S[k];
192:       d1 = T[k];
193:       s2 = S[k+1];
194:       d2 = T[k+1];
195:     } else {
196:       s1 = PetscRealPart(B[k*ld+k]);
197:       d1 = PetscRealPart(A[k+k*ld]);
198:       s2 = PetscRealPart(B[(k+1)*ld+k+1]);
199:       d2 = PetscRealPart(A[k+1+(k+1)*ld]);
200:     }
201:     M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
202:     b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
203:     ep = LAPACKlamch_("S");
204:     /* Compute eigenvalues of the block */
205:     PetscCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
206:     PetscCheck(wi!=0.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Real block in DSVectors_GHIEP");
207:     /* Complex eigenvalues */
208:     PetscCheck(scal1>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
209:     wr1 /= scal1;
210:     wi  /= scal1;
211: #if !defined(PETSC_USE_COMPLEX)
212:     if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
213:       Y[0] = wr1-s2*d2; Y[1] = s2*e; Y[2] = wi; Y[3] = 0.0;
214:     } else {
215:       Y[0] = s1*e; Y[1] = wr1-s1*d1; Y[2] = 0.0; Y[3] = wi;
216:     }
217:     norm = BLASnrm2_(&four,Y,&one);
218:     norm = 1.0/norm;
219:     if (ds->state >= DS_STATE_CONDENSED) {
220:       alpha = norm;
221:       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&two,&two,&alpha,Q+k*ld,&ld,Y,&two,&szero,X+k*ld,&ld));
222:       if (rnorm) *rnorm = SlepcAbsEigenvalue(X[ds->n-1+k*ld],X[ds->n-1+(k+1)*ld]);
223:     } else {
224:       PetscCall(PetscArrayzero(X+k*ld,2*ld));
225:       X[k*ld+k]       = Y[0]*norm;
226:       X[k*ld+k+1]     = Y[1]*norm;
227:       X[(k+1)*ld+k]   = Y[2]*norm;
228:       X[(k+1)*ld+k+1] = Y[3]*norm;
229:     }
230: #else
231:     if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
232:       Y[0] = PetscCMPLX(wr1-s2*d2,wi);
233:       Y[1] = s2*e;
234:     } else {
235:       Y[0] = s1*e;
236:       Y[1] = PetscCMPLX(wr1-s1*d1,wi);
237:     }
238:     norm = BLASnrm2_(&two,Y,&one);
239:     norm = 1.0/norm;
240:     if (ds->state >= DS_STATE_CONDENSED) {
241:       alpha = norm;
242:       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&two,&alpha,Q+k*ld,&ld,Y,&one,&szero,X+k*ld,&one));
243:       if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
244:     } else {
245:       PetscCall(PetscArrayzero(X+k*ld,2*ld));
246:       X[k*ld+k]   = Y[0]*norm;
247:       X[k*ld+k+1] = Y[1]*norm;
248:     }
249:     X[(k+1)*ld+k]   = PetscConj(X[k*ld+k]);
250:     X[(k+1)*ld+k+1] = PetscConj(X[k*ld+k+1]);
251: #endif
252:     (*idx)++;
253:   }
254:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
255:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
256:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
257:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
258:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
259:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: static PetscErrorCode DSVectors_GHIEP(DS ds,DSMatType mat,PetscInt *k,PetscReal *rnorm)
264: {
265:   PetscScalar       *Z;
266:   const PetscScalar *A,*Q;
267:   PetscInt          i;
268:   PetscReal         e,*T;

270:   PetscFunctionBegin;
271:   switch (mat) {
272:     case DS_MAT_X:
273:     case DS_MAT_Y:
274:       if (k) PetscCall(DSVectors_GHIEP_Eigen_Some(ds,k,rnorm));
275:       else {
276:         PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
277:         PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
278:         PetscCall(MatDenseGetArray(ds->omat[mat],&Z));
279:         PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
280:         for (i=0; i<ds->n; i++) {
281:           e = (ds->compact)?T[i+ds->ld]:PetscRealPart(A[(i+1)+ds->ld*i]);
282:           if (e == 0.0) { /* real */
283:             if (ds->state >= DS_STATE_CONDENSED) PetscCall(PetscArraycpy(Z+i*ds->ld,Q+i*ds->ld,ds->ld));
284:             else {
285:               PetscCall(PetscArrayzero(Z+i*ds->ld,ds->ld));
286:               Z[i+i*ds->ld] = 1.0;
287:             }
288:           } else PetscCall(DSVectors_GHIEP_Eigen_Some(ds,&i,rnorm));
289:         }
290:         PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
291:         PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
292:         PetscCall(MatDenseRestoreArray(ds->omat[mat],&Z));
293:         PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
294:       }
295:       break;
296:     case DS_MAT_U:
297:     case DS_MAT_V:
298:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
299:     default:
300:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
301:   }
302:   PetscFunctionReturn(PETSC_SUCCESS);
303: }

305: /*
306:   Extract the eigenvalues contained in the block-diagonal of the indefinite problem.
307:   Only the index range n0..n1 is processed.
308: */
309: PetscErrorCode DSGHIEPComplexEigs(DS ds,PetscInt n0,PetscInt n1,PetscScalar *wr,PetscScalar *wi)
310: {
311:   PetscInt          k,ld;
312:   PetscBLASInt      two=2;
313:   const PetscScalar *A,*B;
314:   PetscReal         *D,*T,b[4],M[4],d1,d2,s1,s2,e,scal1,scal2,ep,wr1,wr2,wi1;

316:   PetscFunctionBegin;
317:   ld = ds->ld;
318:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
319:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
320:   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
321:   PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
322:   for (k=n0;k<n1;k++) {
323:     if (k < n1-1) e = (ds->compact)?T[ld+k]:PetscRealPart(A[(k+1)+ld*k]);
324:     else e = 0.0;
325:     if (e==0.0) { /* real eigenvalue */
326:       wr[k] = (ds->compact)?T[k]/D[k]:A[k+k*ld]/B[k+k*ld];
327: #if !defined(PETSC_USE_COMPLEX)
328:       wi[k] = 0.0 ;
329: #endif
330:     } else { /* diagonal block */
331:       if (ds->compact) {
332:         s1 = D[k];
333:         d1 = T[k];
334:         s2 = D[k+1];
335:         d2 = T[k+1];
336:       } else {
337:         s1 = PetscRealPart(B[k*ld+k]);
338:         d1 = PetscRealPart(A[k+k*ld]);
339:         s2 = PetscRealPart(B[(k+1)*ld+k+1]);
340:         d2 = PetscRealPart(A[k+1+(k+1)*ld]);
341:       }
342:       M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
343:       b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
344:       ep = LAPACKlamch_("S");
345:       /* Compute eigenvalues of the block */
346:       PetscCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi1));
347:       PetscCheck(scal1>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
348:       if (wi1==0.0) { /* Real eigenvalues */
349:         PetscCheck(scal2>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
350:         wr[k] = wr1/scal1; wr[k+1] = wr2/scal2;
351: #if !defined(PETSC_USE_COMPLEX)
352:         wi[k] = wi[k+1] = 0.0;
353: #endif
354:       } else { /* Complex eigenvalues */
355: #if !defined(PETSC_USE_COMPLEX)
356:         wr[k]   = wr1/scal1;
357:         wr[k+1] = wr[k];
358:         wi[k]   = wi1/scal1;
359:         wi[k+1] = -wi[k];
360: #else
361:         wr[k]   = PetscCMPLX(wr1,wi1)/scal1;
362:         wr[k+1] = PetscConj(wr[k]);
363: #endif
364:       }
365:       k++;
366:     }
367:   }
368: #if defined(PETSC_USE_COMPLEX)
369:   if (wi) {
370:     for (k=n0;k<n1;k++) wi[k] = 0.0;
371:   }
372: #endif
373:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
374:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
375:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
376:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
377:   PetscFunctionReturn(PETSC_SUCCESS);
378: }

380: static PetscErrorCode DSSort_GHIEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
381: {
382:   PetscInt       n,i,*perm;
383:   PetscReal      *d,*e,*s;

385:   PetscFunctionBegin;
386: #if !defined(PETSC_USE_COMPLEX)
387:   PetscAssertPointer(wi,3);
388: #endif
389:   n = ds->n;
390:   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
391:   e = d + ds->ld;
392:   PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
393:   PetscCall(DSAllocateWork_Private(ds,ds->ld,ds->ld,0));
394:   perm = ds->perm;
395:   if (!rr) {
396:     rr = wr;
397:     ri = wi;
398:   }
399:   PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_TRUE));
400:   if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_TRUE));
401:   PetscCall(PetscArraycpy(ds->work,wr,n));
402:   for (i=ds->l;i<n;i++) wr[i] = *(ds->work+perm[i]);
403: #if !defined(PETSC_USE_COMPLEX)
404:   PetscCall(PetscArraycpy(ds->work,wi,n));
405:   for (i=ds->l;i<n;i++) wi[i] = *(ds->work+perm[i]);
406: #endif
407:   PetscCall(PetscArraycpy(ds->rwork,s,n));
408:   for (i=ds->l;i<n;i++) s[i] = *(ds->rwork+perm[i]);
409:   PetscCall(PetscArraycpy(ds->rwork,d,n));
410:   for (i=ds->l;i<n;i++) d[i] = *(ds->rwork+perm[i]);
411:   PetscCall(PetscArraycpy(ds->rwork,e,n-1));
412:   PetscCall(PetscArrayzero(e+ds->l,n-1-ds->l));
413:   for (i=ds->l;i<n-1;i++) {
414:     if (perm[i]<n-1) e[i] = *(ds->rwork+perm[i]);
415:   }
416:   if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_FALSE));
417:   PetscCall(DSPermuteColumns_Private(ds,ds->l,n,n,DS_MAT_Q,perm));
418:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
419:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
420:   PetscFunctionReturn(PETSC_SUCCESS);
421: }

423: static PetscErrorCode DSUpdateExtraRow_GHIEP(DS ds)
424: {
425:   PetscInt          i;
426:   PetscBLASInt      n,ld,incx=1;
427:   PetscScalar       *A,*x,*y,one=1.0,zero=0.0;
428:   const PetscScalar *Q;
429:   PetscReal         *T,*b,*r,beta;

431:   PetscFunctionBegin;
432:   PetscCall(PetscBLASIntCast(ds->n,&n));
433:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
434:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
435:   if (ds->compact) {
436:     PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
437:     b = T+ld;
438:     r = T+2*ld;
439:     beta = b[n-1];   /* in compact, we assume all entries are zero except the last one */
440:     for (i=0;i<n;i++) r[i] = PetscRealPart(beta*Q[n-1+i*ld]);
441:     ds->k = n;
442:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
443:   } else {
444:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
445:     PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
446:     x = ds->work;
447:     y = ds->work+ld;
448:     for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
449:     PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
450:     for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
451:     ds->k = n;
452:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
453:   }
454:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
455:   PetscFunctionReturn(PETSC_SUCCESS);
456: }

458: /*
459:   Get eigenvectors with inverse iteration.
460:   The system matrix is in Hessenberg form.
461: */
462: PetscErrorCode DSGHIEPInverseIteration(DS ds,PetscScalar *wr,PetscScalar *wi)
463: {
464:   PetscInt          i,off;
465:   PetscBLASInt      *select,*infoC,ld,n1,mout,info;
466:   const PetscScalar *A,*B;
467:   PetscScalar       *H,*X;
468:   PetscReal         *s,*d,*e;
469: #if defined(PETSC_USE_COMPLEX)
470:   PetscInt          j;
471: #endif

473:   PetscFunctionBegin;
474:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
475:   PetscCall(PetscBLASIntCast(ds->n-ds->l,&n1));
476:   PetscCall(DSAllocateWork_Private(ds,ld*ld+2*ld,ld,2*ld));
477:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
478:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
479:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
480:   PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&H));
481:   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
482:   PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
483:   e = d + ld;
484:   select = ds->iwork;
485:   infoC = ds->iwork + ld;
486:   off = ds->l+ds->l*ld;
487:   if (ds->compact) {
488:     H[off] = d[ds->l]*s[ds->l];
489:     H[off+ld] = e[ds->l]*s[ds->l];
490:     for (i=ds->l+1;i<ds->n-1;i++) {
491:       H[i+(i-1)*ld] = e[i-1]*s[i];
492:       H[i+i*ld] = d[i]*s[i];
493:       H[i+(i+1)*ld] = e[i]*s[i];
494:     }
495:     H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
496:     H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
497:   } else {
498:     s[ds->l]  = PetscRealPart(B[off]);
499:     H[off]    = A[off]*s[ds->l];
500:     H[off+ld] = A[off+ld]*s[ds->l];
501:     for (i=ds->l+1;i<ds->n-1;i++) {
502:       s[i] = PetscRealPart(B[i+i*ld]);
503:       H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
504:       H[i+i*ld]     = A[i+i*ld]*s[i];
505:       H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
506:     }
507:     s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
508:     H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
509:     H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
510:   }
511:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
512:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
513:   for (i=0;i<n1;i++) select[i] = 1;
514: #if !defined(PETSC_USE_COMPLEX)
515:   PetscCallBLAS("LAPACKhsein",LAPACKhsein_("R","N","N",select,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,X+off,&ld,&n1,&mout,ds->work,NULL,infoC,&info));
516: #else
517:   PetscCallBLAS("LAPACKhsein",LAPACKhsein_("R","N","N",select,&n1,H+off,&ld,wr+ds->l,NULL,&ld,X+off,&ld,&n1,&mout,ds->work,ds->rwork,NULL,infoC,&info));

519:   /* Separate real and imaginary part of complex eigenvectors */
520:   for (j=ds->l;j<ds->n;j++) {
521:     if (PetscAbsReal(PetscImaginaryPart(wr[j])) > PetscAbsScalar(wr[j])*PETSC_SQRT_MACHINE_EPSILON) {
522:       for (i=ds->l;i<ds->n;i++) {
523:         X[i+(j+1)*ds->ld] = PetscImaginaryPart(X[i+j*ds->ld]);
524:         X[i+j*ds->ld] = PetscRealPart(X[i+j*ds->ld]);
525:       }
526:       j++;
527:     }
528:   }
529: #endif
530:   SlepcCheckLapackInfo("hsein",info);
531:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
532:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
533:   PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&H));
534:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
535:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
536:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
537:   PetscCall(DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_TRUE));
538:   PetscFunctionReturn(PETSC_SUCCESS);
539: }

541: /*
542:    Undo 2x2 blocks that have real eigenvalues.
543: */
544: PetscErrorCode DSGHIEPRealBlocks(DS ds)
545: {
546:   PetscInt       i;
547:   PetscReal      e,d1,d2,s1,s2,ss1,ss2,t,dd,ss;
548:   PetscReal      maxy,ep,scal1,scal2,snorm;
549:   PetscReal      *T,*D,b[4],M[4],wr1,wr2,wi;
550:   PetscScalar    *A,*B,*Q,Y[4],sone=1.0,szero=0.0;
551:   PetscBLASInt   m,two=2,ld;
552:   PetscBool      isreal;

554:   PetscFunctionBegin;
555:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
556:   PetscCall(PetscBLASIntCast(ds->n-ds->l,&m));
557:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
558:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
559:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
560:   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
561:   PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
562:   PetscCall(DSAllocateWork_Private(ds,2*m,0,0));
563:   for (i=ds->l;i<ds->n-1;i++) {
564:     e = (ds->compact)?T[ld+i]:PetscRealPart(A[(i+1)+ld*i]);
565:     if (e != 0.0) { /* 2x2 block */
566:       if (ds->compact) {
567:         s1 = D[i];
568:         d1 = T[i];
569:         s2 = D[i+1];
570:         d2 = T[i+1];
571:       } else {
572:         s1 = PetscRealPart(B[i*ld+i]);
573:         d1 = PetscRealPart(A[i*ld+i]);
574:         s2 = PetscRealPart(B[(i+1)*ld+i+1]);
575:         d2 = PetscRealPart(A[(i+1)*ld+i+1]);
576:       }
577:       isreal = PETSC_FALSE;
578:       if (s1==s2) { /* apply a Jacobi rotation to compute the eigendecomposition */
579:         dd = d1-d2;
580:         if (2*PetscAbsReal(e) <= dd) {
581:           t = 2*e/dd;
582:           t = t/(1 + PetscSqrtReal(1+t*t));
583:         } else {
584:           t = dd/(2*e);
585:           ss = (t>=0)?1.0:-1.0;
586:           t = ss/(PetscAbsReal(t)+PetscSqrtReal(1+t*t));
587:         }
588:         Y[0] = 1/PetscSqrtReal(1 + t*t); Y[3] = Y[0]; /* c */
589:         Y[1] = Y[0]*t; Y[2] = -Y[1]; /* s */
590:         wr1 = d1+t*e; wr2 = d2-t*e;
591:         ss1 = s1; ss2 = s2;
592:         isreal = PETSC_TRUE;
593:       } else {
594:         ss1 = 1.0; ss2 = 1.0,
595:         M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
596:         b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
597:         ep = LAPACKlamch_("S");

599:         /* Compute eigenvalues of the block */
600:         PetscCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
601:         if (wi==0.0) { /* Real eigenvalues */
602:           isreal = PETSC_TRUE;
603:           PetscCheck(scal1>=ep && scal2>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
604:           wr1 /= scal1;
605:           wr2 /= scal2;
606:           if (PetscAbsReal(s1*d1-wr1)<PetscAbsReal(s2*d2-wr1)) {
607:             Y[0] = wr1-s2*d2;
608:             Y[1] = s2*e;
609:           } else {
610:             Y[0] = s1*e;
611:             Y[1] = wr1-s1*d1;
612:           }
613:           /* normalize with a signature*/
614:           maxy = PetscMax(PetscAbsScalar(Y[0]),PetscAbsScalar(Y[1]));
615:           scal1 = PetscRealPart(Y[0])/maxy;
616:           scal2 = PetscRealPart(Y[1])/maxy;
617:           snorm = scal1*scal1*s1 + scal2*scal2*s2;
618:           if (snorm<0) { ss1 = -1.0; snorm = -snorm; }
619:           snorm = maxy*PetscSqrtReal(snorm);
620:           Y[0] = Y[0]/snorm;
621:           Y[1] = Y[1]/snorm;
622:           if (PetscAbsReal(s1*d1-wr2)<PetscAbsReal(s2*d2-wr2)) {
623:             Y[2] = wr2-s2*d2;
624:             Y[3] = s2*e;
625:           } else {
626:             Y[2] = s1*e;
627:             Y[3] = wr2-s1*d1;
628:           }
629:           maxy = PetscMax(PetscAbsScalar(Y[2]),PetscAbsScalar(Y[3]));
630:           scal1 = PetscRealPart(Y[2])/maxy;
631:           scal2 = PetscRealPart(Y[3])/maxy;
632:           snorm = scal1*scal1*s1 + scal2*scal2*s2;
633:           if (snorm<0) { ss2 = -1.0; snorm = -snorm; }
634:           snorm = maxy*PetscSqrtReal(snorm); Y[2] = Y[2]/snorm; Y[3] = Y[3]/snorm;
635:         }
636:         wr1 *= ss1; wr2 *= ss2;
637:       }
638:       if (isreal) {
639:         if (ds->compact) {
640:           D[i]    = ss1;
641:           T[i]    = wr1;
642:           D[i+1]  = ss2;
643:           T[i+1]  = wr2;
644:           T[ld+i] = 0.0;
645:         } else {
646:           B[i*ld+i]       = ss1;
647:           A[i*ld+i]       = wr1;
648:           B[(i+1)*ld+i+1] = ss2;
649:           A[(i+1)*ld+i+1] = wr2;
650:           A[(i+1)+ld*i]   = 0.0;
651:           A[i+ld*(i+1)]   = 0.0;
652:         }
653:         PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&two,&two,&sone,Q+ds->l+i*ld,&ld,Y,&two,&szero,ds->work,&m));
654:         PetscCall(PetscArraycpy(Q+ds->l+i*ld,ds->work,m));
655:         PetscCall(PetscArraycpy(Q+ds->l+(i+1)*ld,ds->work+m,m));
656:       }
657:       i++;
658:     }
659:   }
660:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
661:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
662:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
663:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
664:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
665:   PetscFunctionReturn(PETSC_SUCCESS);
666: }

668: static PetscErrorCode DSSolve_GHIEP_QR_II(DS ds,PetscScalar *wr,PetscScalar *wi)
669: {
670:   PetscInt          i,off;
671:   PetscBLASInt      n1,ld,one=1,info,lwork;
672:   const PetscScalar *A,*B;
673:   PetscScalar       *H,*Q;
674:   PetscReal         *d,*e,*s;
675: #if defined(PETSC_USE_COMPLEX)
676:   PetscInt          j;
677: #endif

679:   PetscFunctionBegin;
680: #if !defined(PETSC_USE_COMPLEX)
681:   PetscAssertPointer(wi,3);
682: #endif
683:   PetscCall(PetscBLASIntCast(ds->n-ds->l,&n1));
684:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
685:   off = ds->l + ds->l*ld;
686:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
687:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
688:   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
689:   PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
690:   e = d + ld;
691: #if defined(PETSC_USE_DEBUG)
692:   /* Check signature */
693:   for (i=0;i<ds->n;i++) {
694:     PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
695:     PetscCheck(de==1.0 || de==-1.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diagonal elements of the signature matrix must be 1 or -1");
696:   }
697: #endif

699:   /* Quick return if possible */
700:   if (n1 == 1) {
701:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
702:     for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
703:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
704:     PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
705:     if (!ds->compact) {
706:       d[ds->l] = PetscRealPart(A[off]);
707:       s[ds->l] = PetscRealPart(B[off]);
708:     }
709:     PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
710:     PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
711:     wr[ds->l] = d[ds->l]/s[ds->l];
712:     if (wi) wi[ds->l] = 0.0;
713:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
714:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
715:     PetscFunctionReturn(PETSC_SUCCESS);
716:   }

718:   PetscCall(DSAllocateWork_Private(ds,ld*ld,2*ld,ld*2));
719:   lwork = ld*ld;

721:   /* Reduce to pseudotriadiagonal form */
722:   PetscCall(DSIntermediate_GHIEP(ds));

724:   /* Compute Eigenvalues (QR) */
725:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
726:   PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&H));
727:   if (ds->compact) {
728:     H[off]    = d[ds->l]*s[ds->l];
729:     H[off+ld] = e[ds->l]*s[ds->l];
730:     for (i=ds->l+1;i<ds->n-1;i++) {
731:       H[i+(i-1)*ld] = e[i-1]*s[i];
732:       H[i+i*ld]     = d[i]*s[i];
733:       H[i+(i+1)*ld] = e[i]*s[i];
734:     }
735:     H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
736:     H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
737:   } else {
738:     s[ds->l]  = PetscRealPart(B[off]);
739:     H[off]    = A[off]*s[ds->l];
740:     H[off+ld] = A[off+ld]*s[ds->l];
741:     for (i=ds->l+1;i<ds->n-1;i++) {
742:       s[i] = PetscRealPart(B[i+i*ld]);
743:       H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
744:       H[i+i*ld]     = A[i+i*ld]*s[i];
745:       H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
746:     }
747:     s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
748:     H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
749:     H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
750:   }

752: #if !defined(PETSC_USE_COMPLEX)
753:   PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,ds->work,&lwork,&info));
754: #else
755:   PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,NULL,&ld,ds->work,&lwork,&info));
756:   for (i=ds->l;i<ds->n;i++) if (PetscAbsReal(PetscImaginaryPart(wr[i]))<10*PETSC_MACHINE_EPSILON) wr[i] = PetscRealPart(wr[i]);
757:   /* Sort to have consecutive conjugate pairs */
758:   for (i=ds->l;i<ds->n;i++) {
759:     j=i+1;
760:     while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
761:     if (j==ds->n) {
762:       PetscCheck(PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
763:       wr[i] = PetscRealPart(wr[i]);
764:     } else { /* complex eigenvalue */
765:       wr[j] = wr[i+1];
766:       if (PetscImaginaryPart(wr[i])<0) wr[i] = PetscConj(wr[i]);
767:       wr[i+1] = PetscConj(wr[i]);
768:       i++;
769:     }
770:   }
771: #endif
772:   SlepcCheckLapackInfo("hseqr",info);
773:   PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&H));
774:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
775:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
776:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
777:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));

779:   /* Compute Eigenvectors with Inverse Iteration */
780:   PetscCall(DSGHIEPInverseIteration(ds,wr,wi));

782:   /* Recover eigenvalues from diagonal */
783:   PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
784: #if defined(PETSC_USE_COMPLEX)
785:   if (wi) {
786:     for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
787:   }
788: #endif
789:   PetscFunctionReturn(PETSC_SUCCESS);
790: }

792: static PetscErrorCode DSSolve_GHIEP_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
793: {
794:   PetscInt          i,j,off,nwu=0,n,lw,lwr,nwru=0;
795:   PetscBLASInt      n_,ld,info,lwork,ilo,ihi;
796:   const PetscScalar *A,*B;
797:   PetscScalar       *H,*Q,*X;
798:   PetscReal         *d,*s,*scale,nrm,*rcde,*rcdv;
799: #if defined(PETSC_USE_COMPLEX)
800:   PetscInt          k;
801: #endif

803:   PetscFunctionBegin;
804: #if !defined(PETSC_USE_COMPLEX)
805:   PetscAssertPointer(wi,3);
806: #endif
807:   n = ds->n-ds->l;
808:   PetscCall(PetscBLASIntCast(n,&n_));
809:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
810:   off = ds->l + ds->l*ld;
811:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
812:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
813:   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
814:   PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
815: #if defined(PETSC_USE_DEBUG)
816:   /* Check signature */
817:   for (i=0;i<ds->n;i++) {
818:     PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
819:     PetscCheck(de==1.0 || de==-1.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diagonal elements of the signature matrix must be 1 or -1");
820:   }
821: #endif

823:   /* Quick return if possible */
824:   if (n_ == 1) {
825:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
826:     for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
827:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
828:     PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
829:     if (!ds->compact) {
830:       d[ds->l] = PetscRealPart(A[off]);
831:       s[ds->l] = PetscRealPart(B[off]);
832:     }
833:     PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
834:     PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
835:     wr[ds->l] = d[ds->l]/s[ds->l];
836:     if (wi) wi[ds->l] = 0.0;
837:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
838:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
839:     PetscFunctionReturn(PETSC_SUCCESS);
840:   }

842:   lw = 14*ld+ld*ld;
843:   lwr = 7*ld;
844:   PetscCall(DSAllocateWork_Private(ds,lw,lwr,0));
845:   scale = ds->rwork+nwru;
846:   nwru += ld;
847:   rcde = ds->rwork+nwru;
848:   nwru += ld;
849:   rcdv = ds->rwork+nwru;

851:   /* Form pseudo-symmetric matrix */
852:   H =  ds->work+nwu;
853:   nwu += n*n;
854:   PetscCall(PetscArrayzero(H,n*n));
855:   if (ds->compact) {
856:     for (i=0;i<n-1;i++) {
857:       H[i+i*n]     = s[ds->l+i]*d[ds->l+i];
858:       H[i+1+i*n]   = s[ds->l+i+1]*d[ld+ds->l+i];
859:       H[i+(i+1)*n] = s[ds->l+i]*d[ld+ds->l+i];
860:     }
861:     H[n-1+(n-1)*n] = s[ds->l+n-1]*d[ds->l+n-1];
862:     for (i=0;i<ds->k-ds->l;i++) {
863:       H[ds->k-ds->l+i*n] = s[ds->k]*d[2*ld+ds->l+i];
864:       H[i+(ds->k-ds->l)*n] = s[i+ds->l]*d[2*ld+ds->l+i];
865:     }
866:   } else {
867:     for (j=0;j<n;j++) {
868:       for (i=0;i<n;i++) H[i+j*n] = B[off+i+i*ld]*A[off+i+j*ld];
869:     }
870:   }

872:   /* Compute eigenpairs */
873:   PetscCall(PetscBLASIntCast(lw-nwu,&lwork));
874:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
875:   PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_X],&X));
876: #if !defined(PETSC_USE_COMPLEX)
877:   PetscCallBLAS("LAPACKgeevx",LAPACKgeevx_("B","N","V","N",&n_,H,&n_,wr+ds->l,wi+ds->l,NULL,&ld,X+off,&ld,&ilo,&ihi,scale,&nrm,rcde,rcdv,ds->work+nwu,&lwork,NULL,&info));
878: #else
879:   PetscCallBLAS("LAPACKgeevx",LAPACKgeevx_("B","N","V","N",&n_,H,&n_,wr+ds->l,NULL,&ld,X+off,&ld,&ilo,&ihi,scale,&nrm,rcde,rcdv,ds->work+nwu,&lwork,ds->rwork+nwru,&info));

881:   /* Sort to have consecutive conjugate pairs
882:      Separate real and imaginary part of complex eigenvectors*/
883:   for (i=ds->l;i<ds->n;i++) {
884:     j=i+1;
885:     while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
886:     if (j==ds->n) {
887:       PetscCheck(PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
888:       wr[i]=PetscRealPart(wr[i]); /* real eigenvalue */
889:       for (k=ds->l;k<ds->n;k++) {
890:         X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
891:       }
892:     } else { /* complex eigenvalue */
893:       if (j!=i+1) {
894:         wr[j] = wr[i+1];
895:         PetscCall(PetscArraycpy(X+j*ds->ld,X+(i+1)*ds->ld,ds->ld));
896:       }
897:       if (PetscImaginaryPart(wr[i])<0) {
898:         wr[i] = PetscConj(wr[i]);
899:         for (k=ds->l;k<ds->n;k++) {
900:           X[k+(i+1)*ds->ld] = -PetscImaginaryPart(X[k+i*ds->ld]);
901:           X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
902:         }
903:       } else {
904:         for (k=ds->l;k<ds->n;k++) {
905:           X[k+(i+1)*ds->ld] = PetscImaginaryPart(X[k+i*ds->ld]);
906:           X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
907:         }
908:       }
909:       wr[i+1] = PetscConj(wr[i]);
910:       i++;
911:     }
912:   }
913: #endif
914:   SlepcCheckLapackInfo("geevx",info);
915:   PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_X],&X));
916:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
917:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
918:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
919:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));

921:   /* Compute real s-orthonormal basis */
922:   PetscCall(DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_FALSE));

924:   /* Recover eigenvalues from diagonal */
925:   PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
926: #if defined(PETSC_USE_COMPLEX)
927:   if (wi) {
928:     for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
929:   }
930: #endif
931:   PetscFunctionReturn(PETSC_SUCCESS);
932: }

934: static PetscErrorCode DSGetTruncateSize_GHIEP(DS ds,PetscInt l,PetscInt n,PetscInt *k)
935: {
936:   PetscReal *T;

938:   PetscFunctionBegin;
939:   PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
940:   if (T[l+(*k)-1+ds->ld] !=0.0) {
941:     if (l+(*k)<n-1) (*k)++;
942:     else (*k)--;
943:   }
944:   PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
945:   PetscFunctionReturn(PETSC_SUCCESS);
946: }

948: static PetscErrorCode DSTruncate_GHIEP(DS ds,PetscInt n,PetscBool trim)
949: {
950:   PetscInt    i,ld=ds->ld,l=ds->l;
951:   PetscScalar *A;
952:   PetscReal   *T,*b,*r,*omega;

954:   PetscFunctionBegin;
955:   if (ds->compact) {
956:     PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
957:     PetscCall(DSGetArrayReal(ds,DS_MAT_D,&omega));
958: #if defined(PETSC_USE_DEBUG)
959:     /* make sure diagonal 2x2 block is not broken */
960:     PetscCheck(ds->state<DS_STATE_CONDENSED || n==0 || n==ds->n || T[n-1+ld]==0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"The given size would break a 2x2 block, call DSGetTruncateSize() first");
961: #endif
962:   }
963:   if (trim) {
964:     if (!ds->compact && ds->extrarow) {   /* clean extra row */
965:       PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
966:       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
967:       PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
968:     }
969:     ds->l = 0;
970:     ds->k = 0;
971:     ds->n = n;
972:     ds->t = ds->n;   /* truncated length equal to the new dimension */
973:   } else {
974:     if (!ds->compact && ds->extrarow && ds->k==ds->n) {
975:       /* copy entries of extra row to the new position, then clean last row */
976:       PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
977:       for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
978:       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
979:       PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
980:     }
981:     if (ds->compact) {
982:       b = T+ld;
983:       r = T+2*ld;
984:       b[n-1] = r[n-1];
985:       b[n] = b[ds->n];
986:       omega[n] = omega[ds->n];
987:     }
988:     ds->k = (ds->extrarow)? n: 0;
989:     ds->t = ds->n;   /* truncated length equal to previous dimension */
990:     ds->n = n;
991:   }
992:   if (ds->compact) {
993:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
994:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&omega));
995:   }
996:   PetscFunctionReturn(PETSC_SUCCESS);
997: }

999: #if !defined(PETSC_HAVE_MPIUNI)
1000: static PetscErrorCode DSSynchronize_GHIEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
1001: {
1002:   PetscScalar    *A,*B,*Q;
1003:   PetscReal      *T,*D;
1004:   PetscInt       ld=ds->ld,l=ds->l,k=0,kr=0;
1005:   PetscMPIInt    n,rank,off=0,size,ldn,ld3,ld_;

1007:   PetscFunctionBegin;
1008:   if (ds->compact) kr = 4*ld;
1009:   else k = 2*(ds->n-l)*ld;
1010:   if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
1011:   if (eigr) k += (ds->n-l);
1012:   if (eigi) k += (ds->n-l);
1013:   PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
1014:   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
1015:   PetscCall(PetscMPIIntCast(ds->n-l,&n));
1016:   PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
1017:   PetscCall(PetscMPIIntCast(ld*3,&ld3));
1018:   PetscCall(PetscMPIIntCast(ld,&ld_));
1019:   if (ds->compact) {
1020:     PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
1021:     PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
1022:   } else {
1023:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
1024:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
1025:   }
1026:   if (ds->state>DS_STATE_RAW) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
1027:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
1028:   if (!rank) {
1029:     if (ds->compact) {
1030:       PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1031:       PetscCallMPI(MPI_Pack(D,ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1032:     } else {
1033:       PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1034:       PetscCallMPI(MPI_Pack(B+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1035:     }
1036:     if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Pack(Q+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1037:     if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1038: #if !defined(PETSC_USE_COMPLEX)
1039:     if (eigi) PetscCallMPI(MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1040: #endif
1041:   }
1042:   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
1043:   if (rank) {
1044:     if (ds->compact) {
1045:       PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
1046:       PetscCallMPI(MPI_Unpack(ds->work,size,&off,D,ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
1047:     } else {
1048:       PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1049:       PetscCallMPI(MPI_Unpack(ds->work,size,&off,B+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1050:     }
1051:     if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Unpack(ds->work,size,&off,Q+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1052:     if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1053: #if !defined(PETSC_USE_COMPLEX)
1054:     if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1055: #endif
1056:   }
1057:   if (ds->compact) {
1058:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
1059:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
1060:   } else {
1061:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
1062:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
1063:   }
1064:   if (ds->state>DS_STATE_RAW) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
1065:   PetscFunctionReturn(PETSC_SUCCESS);
1066: }
1067: #endif

1069: static PetscErrorCode DSHermitian_GHIEP(DS ds,DSMatType m,PetscBool *flg)
1070: {
1071:   PetscFunctionBegin;
1072:   if ((m==DS_MAT_A && !ds->extrarow) || m==DS_MAT_B) *flg = PETSC_TRUE;
1073:   else *flg = PETSC_FALSE;
1074:   PetscFunctionReturn(PETSC_SUCCESS);
1075: }

1077: /*MC
1078:    DSGHIEP - Dense Generalized Hermitian Indefinite Eigenvalue Problem.

1080:    Level: beginner

1082:    Notes:
1083:    The problem is expressed as A*X = B*X*Lambda, where both A and B are
1084:    real symmetric (or complex Hermitian) and possibly indefinite. Lambda
1085:    is a diagonal matrix whose diagonal elements are the arguments of DSSolve().
1086:    After solve, A is overwritten with Lambda. Note that in the case of real
1087:    scalars, A is overwritten with a real representation of Lambda, i.e.,
1088:    complex conjugate eigenvalue pairs are stored as a 2x2 block in the
1089:    quasi-diagonal matrix.

1091:    In the intermediate state A is reduced to tridiagonal form and B is
1092:    transformed into a signature matrix. In compact storage format, these
1093:    matrices are stored in T and D, respectively.

1095:    Used DS matrices:
1096: +  DS_MAT_A - first problem matrix
1097: .  DS_MAT_B - second problem matrix
1098: .  DS_MAT_T - symmetric tridiagonal matrix of the reduced pencil
1099: .  DS_MAT_D - diagonal matrix (signature) of the reduced pencil
1100: -  DS_MAT_Q - pseudo-orthogonal transformation that reduces (A,B) to
1101:    tridiagonal-diagonal form (intermediate step) or a real basis of eigenvectors

1103:    Implemented methods:
1104: +  0 - QR iteration plus inverse iteration for the eigenvectors
1105: .  1 - HZ iteration
1106: -  2 - QR iteration plus pseudo-orthogonalization for the eigenvectors

1108:    References:
1109: .  1. - C. Campos and J. E. Roman, "Restarted Q-Arnoldi-type methods exploiting
1110:    symmetry in quadratic eigenvalue problems", BIT Numer. Math. 56(4):1213-1236, 2016.

1112: .seealso: DSCreate(), DSSetType(), DSType
1113: M*/
1114: SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS ds)
1115: {
1116:   PetscFunctionBegin;
1117:   ds->ops->allocate        = DSAllocate_GHIEP;
1118:   ds->ops->view            = DSView_GHIEP;
1119:   ds->ops->vectors         = DSVectors_GHIEP;
1120:   ds->ops->solve[0]        = DSSolve_GHIEP_QR_II;
1121:   ds->ops->solve[1]        = DSSolve_GHIEP_HZ;
1122:   ds->ops->solve[2]        = DSSolve_GHIEP_QR;
1123:   ds->ops->sort            = DSSort_GHIEP;
1124: #if !defined(PETSC_HAVE_MPIUNI)
1125:   ds->ops->synchronize     = DSSynchronize_GHIEP;
1126: #endif
1127:   ds->ops->gettruncatesize = DSGetTruncateSize_GHIEP;
1128:   ds->ops->truncate        = DSTruncate_GHIEP;
1129:   ds->ops->update          = DSUpdateExtraRow_GHIEP;
1130:   ds->ops->hermitian       = DSHermitian_GHIEP;
1131:   PetscFunctionReturn(PETSC_SUCCESS);
1132: }