Actual source code: dsghiep.c

slepc-3.18.2 2023-01-26
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: PetscErrorCode DSAllocate_GHIEP(DS ds,PetscInt ld)
 15: {
 16:   DSAllocateMat_Private(ds,DS_MAT_A);
 17:   DSAllocateMat_Private(ds,DS_MAT_B);
 18:   DSAllocateMat_Private(ds,DS_MAT_Q);
 19:   DSAllocateMat_Private(ds,DS_MAT_T);
 20:   DSAllocateMat_Private(ds,DS_MAT_D);
 21:   PetscFree(ds->perm);
 22:   PetscMalloc1(ld,&ds->perm);
 23:   return 0;
 24: }

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

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

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

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

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

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

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

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

259: PetscErrorCode DSVectors_GHIEP(DS ds,DSMatType mat,PetscInt *k,PetscReal *rnorm)
260: {
261:   PetscScalar       *Z;
262:   const PetscScalar *A,*Q;
263:   PetscInt          i;
264:   PetscReal         e,*T;

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

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

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

374: PetscErrorCode DSSort_GHIEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
375: {
376:   PetscInt       n,i,*perm;
377:   PetscReal      *d,*e,*s;

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

416: PetscErrorCode DSUpdateExtraRow_GHIEP(DS ds)
417: {
418:   PetscInt          i;
419:   PetscBLASInt      n,ld,incx=1;
420:   PetscScalar       *A,*x,*y,one=1.0,zero=0.0;
421:   const PetscScalar *Q;
422:   PetscReal         *T,*b,*r,beta;

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

450: /*
451:   Get eigenvectors with inverse iteration.
452:   The system matrix is in Hessenberg form.
453: */
454: PetscErrorCode DSGHIEPInverseIteration(DS ds,PetscScalar *wr,PetscScalar *wi)
455: {
456:   PetscInt          i,off;
457:   PetscBLASInt      *select,*infoC,ld,n1,mout,info;
458:   const PetscScalar *A,*B;
459:   PetscScalar       *H,*X;
460:   PetscReal         *s,*d,*e;
461: #if defined(PETSC_USE_COMPLEX)
462:   PetscInt          j;
463: #endif

465:   PetscBLASIntCast(ds->ld,&ld);
466:   PetscBLASIntCast(ds->n-ds->l,&n1);
467:   DSAllocateWork_Private(ds,ld*ld+2*ld,ld,2*ld);
468:   DSAllocateMat_Private(ds,DS_MAT_W);
469:   MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A);
470:   MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B);
471:   MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&H);
472:   DSGetArrayReal(ds,DS_MAT_T,&d);
473:   DSGetArrayReal(ds,DS_MAT_D,&s);
474:   e = d + ld;
475:   select = ds->iwork;
476:   infoC = ds->iwork + ld;
477:   off = ds->l+ds->l*ld;
478:   if (ds->compact) {
479:     H[off] = d[ds->l]*s[ds->l];
480:     H[off+ld] = e[ds->l]*s[ds->l];
481:     for (i=ds->l+1;i<ds->n-1;i++) {
482:       H[i+(i-1)*ld] = e[i-1]*s[i];
483:       H[i+i*ld] = d[i]*s[i];
484:       H[i+(i+1)*ld] = e[i]*s[i];
485:     }
486:     H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
487:     H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
488:   } else {
489:     s[ds->l]  = PetscRealPart(B[off]);
490:     H[off]    = A[off]*s[ds->l];
491:     H[off+ld] = A[off+ld]*s[ds->l];
492:     for (i=ds->l+1;i<ds->n-1;i++) {
493:       s[i] = PetscRealPart(B[i+i*ld]);
494:       H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
495:       H[i+i*ld]     = A[i+i*ld]*s[i];
496:       H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
497:     }
498:     s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
499:     H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
500:     H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
501:   }
502:   DSAllocateMat_Private(ds,DS_MAT_X);
503:   MatDenseGetArray(ds->omat[DS_MAT_X],&X);
504:   for (i=0;i<n1;i++) select[i] = 1;
505: #if !defined(PETSC_USE_COMPLEX)
506:   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));
507: #else
508:   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));

510:   /* Separate real and imaginary part of complex eigenvectors */
511:   for (j=ds->l;j<ds->n;j++) {
512:     if (PetscAbsReal(PetscImaginaryPart(wr[j])) > PetscAbsScalar(wr[j])*PETSC_SQRT_MACHINE_EPSILON) {
513:       for (i=ds->l;i<ds->n;i++) {
514:         X[i+(j+1)*ds->ld] = PetscImaginaryPart(X[i+j*ds->ld]);
515:         X[i+j*ds->ld] = PetscRealPart(X[i+j*ds->ld]);
516:       }
517:       j++;
518:     }
519:   }
520: #endif
521:   SlepcCheckLapackInfo("hsein",info);
522:   MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A);
523:   MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B);
524:   MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&H);
525:   MatDenseRestoreArray(ds->omat[DS_MAT_X],&X);
526:   DSRestoreArrayReal(ds,DS_MAT_T,&d);
527:   DSRestoreArrayReal(ds,DS_MAT_D,&s);
528:   DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_TRUE);
529:   return 0;
530: }

532: /*
533:    Undo 2x2 blocks that have real eigenvalues.
534: */
535: PetscErrorCode DSGHIEPRealBlocks(DS ds)
536: {
537:   PetscInt       i;
538:   PetscReal      e,d1,d2,s1,s2,ss1,ss2,t,dd,ss;
539:   PetscReal      maxy,ep,scal1,scal2,snorm;
540:   PetscReal      *T,*D,b[4],M[4],wr1,wr2,wi;
541:   PetscScalar    *A,*B,*Q,Y[4],sone=1.0,szero=0.0;
542:   PetscBLASInt   m,two=2,ld;
543:   PetscBool      isreal;

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

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

658: PetscErrorCode DSSolve_GHIEP_QR_II(DS ds,PetscScalar *wr,PetscScalar *wi)
659: {
660:   PetscInt          i,off;
661:   PetscBLASInt      n1,ld,one=1,info,lwork;
662:   const PetscScalar *A,*B;
663:   PetscScalar       *H,*Q;
664:   PetscReal         *d,*e,*s;
665: #if defined(PETSC_USE_COMPLEX)
666:   PetscInt          j;
667: #endif

669: #if !defined(PETSC_USE_COMPLEX)
671: #endif
672:   PetscBLASIntCast(ds->n-ds->l,&n1);
673:   PetscBLASIntCast(ds->ld,&ld);
674:   off = ds->l + ds->l*ld;
675:   MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A);
676:   MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B);
677:   DSGetArrayReal(ds,DS_MAT_T,&d);
678:   DSGetArrayReal(ds,DS_MAT_D,&s);
679:   e = d + ld;
680: #if defined(PETSC_USE_DEBUG)
681:   /* Check signature */
682:   for (i=0;i<ds->n;i++) {
683:     PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
685:   }
686: #endif

688:   /* Quick return if possible */
689:   if (n1 == 1) {
690:     MatDenseGetArray(ds->omat[DS_MAT_Q],&Q);
691:     for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
692:     MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q);
693:     DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
694:     if (!ds->compact) {
695:       d[ds->l] = PetscRealPart(A[off]);
696:       s[ds->l] = PetscRealPart(B[off]);
697:     }
698:     MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A);
699:     MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B);
700:     wr[ds->l] = d[ds->l]/s[ds->l];
701:     if (wi) wi[ds->l] = 0.0;
702:     DSRestoreArrayReal(ds,DS_MAT_T,&d);
703:     DSRestoreArrayReal(ds,DS_MAT_D,&s);
704:     return 0;
705:   }

707:   DSAllocateWork_Private(ds,ld*ld,2*ld,ld*2);
708:   lwork = ld*ld;

710:   /* Reduce to pseudotriadiagonal form */
711:   DSIntermediate_GHIEP(ds);

713:   /* Compute Eigenvalues (QR) */
714:   DSAllocateMat_Private(ds,DS_MAT_W);
715:   MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&H);
716:   if (ds->compact) {
717:     H[off]    = d[ds->l]*s[ds->l];
718:     H[off+ld] = e[ds->l]*s[ds->l];
719:     for (i=ds->l+1;i<ds->n-1;i++) {
720:       H[i+(i-1)*ld] = e[i-1]*s[i];
721:       H[i+i*ld]     = d[i]*s[i];
722:       H[i+(i+1)*ld] = e[i]*s[i];
723:     }
724:     H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
725:     H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
726:   } else {
727:     s[ds->l]  = PetscRealPart(B[off]);
728:     H[off]    = A[off]*s[ds->l];
729:     H[off+ld] = A[off+ld]*s[ds->l];
730:     for (i=ds->l+1;i<ds->n-1;i++) {
731:       s[i] = PetscRealPart(B[i+i*ld]);
732:       H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
733:       H[i+i*ld]     = A[i+i*ld]*s[i];
734:       H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
735:     }
736:     s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
737:     H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
738:     H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
739:   }

741: #if !defined(PETSC_USE_COMPLEX)
742:   PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,ds->work,&lwork,&info));
743: #else
744:   PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,NULL,&ld,ds->work,&lwork,&info));
745:   for (i=ds->l;i<ds->n;i++) if (PetscAbsReal(PetscImaginaryPart(wr[i]))<10*PETSC_MACHINE_EPSILON) wr[i] = PetscRealPart(wr[i]);
746:   /* Sort to have consecutive conjugate pairs */
747:   for (i=ds->l;i<ds->n;i++) {
748:     j=i+1;
749:     while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
750:     if (j==ds->n) {
752:       wr[i] = PetscRealPart(wr[i]);
753:     } else { /* complex eigenvalue */
754:       wr[j] = wr[i+1];
755:       if (PetscImaginaryPart(wr[i])<0) wr[i] = PetscConj(wr[i]);
756:       wr[i+1] = PetscConj(wr[i]);
757:       i++;
758:     }
759:   }
760: #endif
761:   SlepcCheckLapackInfo("hseqr",info);
762:   MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&H);
763:   MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A);
764:   MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B);
765:   DSRestoreArrayReal(ds,DS_MAT_T,&d);
766:   DSRestoreArrayReal(ds,DS_MAT_D,&s);

768:   /* Compute Eigenvectors with Inverse Iteration */
769:   DSGHIEPInverseIteration(ds,wr,wi);

771:   /* Recover eigenvalues from diagonal */
772:   DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
773: #if defined(PETSC_USE_COMPLEX)
774:   if (wi) {
775:     for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
776:   }
777: #endif
778:   return 0;
779: }

781: PetscErrorCode DSSolve_GHIEP_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
782: {
783:   PetscInt          i,j,off,nwu=0,n,lw,lwr,nwru=0;
784:   PetscBLASInt      n_,ld,info,lwork,ilo,ihi;
785:   const PetscScalar *A,*B;
786:   PetscScalar       *H,*Q,*X;
787:   PetscReal         *d,*s,*scale,nrm,*rcde,*rcdv;
788: #if defined(PETSC_USE_COMPLEX)
789:   PetscInt          k;
790: #endif

792: #if !defined(PETSC_USE_COMPLEX)
794: #endif
795:   n = ds->n-ds->l;
796:   PetscBLASIntCast(n,&n_);
797:   PetscBLASIntCast(ds->ld,&ld);
798:   off = ds->l + ds->l*ld;
799:   MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A);
800:   MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B);
801:   DSGetArrayReal(ds,DS_MAT_T,&d);
802:   DSGetArrayReal(ds,DS_MAT_D,&s);
803: #if defined(PETSC_USE_DEBUG)
804:   /* Check signature */
805:   for (i=0;i<ds->n;i++) {
806:     PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
808:   }
809: #endif

811:   /* Quick return if possible */
812:   if (n_ == 1) {
813:     MatDenseGetArray(ds->omat[DS_MAT_Q],&Q);
814:     for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
815:     MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q);
816:     DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
817:     if (!ds->compact) {
818:       d[ds->l] = PetscRealPart(A[off]);
819:       s[ds->l] = PetscRealPart(B[off]);
820:     }
821:     MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A);
822:     MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B);
823:     wr[ds->l] = d[ds->l]/s[ds->l];
824:     if (wi) wi[ds->l] = 0.0;
825:     DSRestoreArrayReal(ds,DS_MAT_T,&d);
826:     DSRestoreArrayReal(ds,DS_MAT_D,&s);
827:     return 0;
828:   }

830:   lw = 14*ld+ld*ld;
831:   lwr = 7*ld;
832:   DSAllocateWork_Private(ds,lw,lwr,0);
833:   scale = ds->rwork+nwru;
834:   nwru += ld;
835:   rcde = ds->rwork+nwru;
836:   nwru += ld;
837:   rcdv = ds->rwork+nwru;

839:   /* Form pseudo-symmetric matrix */
840:   H =  ds->work+nwu;
841:   nwu += n*n;
842:   PetscArrayzero(H,n*n);
843:   if (ds->compact) {
844:     for (i=0;i<n-1;i++) {
845:       H[i+i*n]     = s[ds->l+i]*d[ds->l+i];
846:       H[i+1+i*n]   = s[ds->l+i+1]*d[ld+ds->l+i];
847:       H[i+(i+1)*n] = s[ds->l+i]*d[ld+ds->l+i];
848:     }
849:     H[n-1+(n-1)*n] = s[ds->l+n-1]*d[ds->l+n-1];
850:     for (i=0;i<ds->k-ds->l;i++) {
851:       H[ds->k-ds->l+i*n] = s[ds->k]*d[2*ld+ds->l+i];
852:       H[i+(ds->k-ds->l)*n] = s[i+ds->l]*d[2*ld+ds->l+i];
853:     }
854:   } else {
855:     for (j=0;j<n;j++) {
856:       for (i=0;i<n;i++) H[i+j*n] = B[off+i+i*ld]*A[off+i+j*ld];
857:     }
858:   }

860:   /* Compute eigenpairs */
861:   PetscBLASIntCast(lw-nwu,&lwork);
862:   DSAllocateMat_Private(ds,DS_MAT_X);
863:   MatDenseGetArrayWrite(ds->omat[DS_MAT_X],&X);
864: #if !defined(PETSC_USE_COMPLEX)
865:   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));
866: #else
867:   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));

869:   /* Sort to have consecutive conjugate pairs
870:      Separate real and imaginary part of complex eigenvectors*/
871:   for (i=ds->l;i<ds->n;i++) {
872:     j=i+1;
873:     while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
874:     if (j==ds->n) {
876:       wr[i]=PetscRealPart(wr[i]); /* real eigenvalue */
877:       for (k=ds->l;k<ds->n;k++) {
878:         X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
879:       }
880:     } else { /* complex eigenvalue */
881:       if (j!=i+1) {
882:         wr[j] = wr[i+1];
883:         PetscArraycpy(X+j*ds->ld,X+(i+1)*ds->ld,ds->ld);
884:       }
885:       if (PetscImaginaryPart(wr[i])<0) {
886:         wr[i] = PetscConj(wr[i]);
887:         for (k=ds->l;k<ds->n;k++) {
888:           X[k+(i+1)*ds->ld] = -PetscImaginaryPart(X[k+i*ds->ld]);
889:           X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
890:         }
891:       } else {
892:         for (k=ds->l;k<ds->n;k++) {
893:           X[k+(i+1)*ds->ld] = PetscImaginaryPart(X[k+i*ds->ld]);
894:           X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
895:         }
896:       }
897:       wr[i+1] = PetscConj(wr[i]);
898:       i++;
899:     }
900:   }
901: #endif
902:   SlepcCheckLapackInfo("geevx",info);
903:   MatDenseRestoreArrayWrite(ds->omat[DS_MAT_X],&X);
904:   MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A);
905:   MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B);
906:   DSRestoreArrayReal(ds,DS_MAT_T,&d);
907:   DSRestoreArrayReal(ds,DS_MAT_D,&s);

909:   /* Compute real s-orthonormal basis */
910:   DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_FALSE);

912:   /* Recover eigenvalues from diagonal */
913:   DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
914: #if defined(PETSC_USE_COMPLEX)
915:   if (wi) {
916:     for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
917:   }
918: #endif
919:   return 0;
920: }

922: PetscErrorCode DSGetTruncateSize_GHIEP(DS ds,PetscInt l,PetscInt n,PetscInt *k)
923: {
924:   PetscReal *T;

926:   DSGetArrayReal(ds,DS_MAT_T,&T);
927:   if (T[l+(*k)-1+ds->ld] !=0.0) {
928:     if (l+(*k)<n-1) (*k)++;
929:     else (*k)--;
930:   }
931:   DSRestoreArrayReal(ds,DS_MAT_T,&T);
932:   return 0;
933: }

935: PetscErrorCode DSTruncate_GHIEP(DS ds,PetscInt n,PetscBool trim)
936: {
937:   PetscInt    i,ld=ds->ld,l=ds->l;
938:   PetscScalar *A;
939:   PetscReal   *T,*b,*r,*omega;

941:   if (ds->compact) {
942:     DSGetArrayReal(ds,DS_MAT_T,&T);
943:     DSGetArrayReal(ds,DS_MAT_D,&omega);
944: #if defined(PETSC_USE_DEBUG)
945:     /* make sure diagonal 2x2 block is not broken */
947: #endif
948:   }
949:   if (trim) {
950:     if (!ds->compact && ds->extrarow) {   /* clean extra row */
951:       MatDenseGetArray(ds->omat[DS_MAT_A],&A);
952:       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
953:       MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
954:     }
955:     ds->l = 0;
956:     ds->k = 0;
957:     ds->n = n;
958:     ds->t = ds->n;   /* truncated length equal to the new dimension */
959:   } else {
960:     if (!ds->compact && ds->extrarow && ds->k==ds->n) {
961:       /* copy entries of extra row to the new position, then clean last row */
962:       MatDenseGetArray(ds->omat[DS_MAT_A],&A);
963:       for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
964:       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
965:       MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
966:     }
967:     if (ds->compact) {
968:       b = T+ld;
969:       r = T+2*ld;
970:       b[n-1] = r[n-1];
971:       b[n] = b[ds->n];
972:       omega[n] = omega[ds->n];
973:     }
974:     ds->k = (ds->extrarow)? n: 0;
975:     ds->t = ds->n;   /* truncated length equal to previous dimension */
976:     ds->n = n;
977:   }
978:   if (ds->compact) {
979:     DSRestoreArrayReal(ds,DS_MAT_T,&T);
980:     DSRestoreArrayReal(ds,DS_MAT_D,&omega);
981:   }
982:   return 0;
983: }

985: #if !defined(PETSC_HAVE_MPIUNI)
986: PetscErrorCode DSSynchronize_GHIEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
987: {
988:   PetscScalar    *A,*B,*Q;
989:   PetscReal      *T,*D;
990:   PetscInt       ld=ds->ld,l=ds->l,k=0,kr=0;
991:   PetscMPIInt    n,rank,off=0,size,ldn,ld3,ld_;

993:   if (ds->compact) kr = 4*ld;
994:   else k = 2*(ds->n-l)*ld;
995:   if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
996:   if (eigr) k += (ds->n-l);
997:   if (eigi) k += (ds->n-l);
998:   DSAllocateWork_Private(ds,k+kr,0,0);
999:   PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size);
1000:   PetscMPIIntCast(ds->n-l,&n);
1001:   PetscMPIIntCast(ld*(ds->n-l),&ldn);
1002:   PetscMPIIntCast(ld*3,&ld3);
1003:   PetscMPIIntCast(ld,&ld_);
1004:   if (ds->compact) {
1005:     DSGetArrayReal(ds,DS_MAT_T,&T);
1006:     DSGetArrayReal(ds,DS_MAT_D,&D);
1007:   } else {
1008:     MatDenseGetArray(ds->omat[DS_MAT_A],&A);
1009:     MatDenseGetArray(ds->omat[DS_MAT_B],&B);
1010:   }
1011:   if (ds->state>DS_STATE_RAW) MatDenseGetArray(ds->omat[DS_MAT_Q],&Q);
1012:   MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
1013:   if (!rank) {
1014:     if (ds->compact) {
1015:       MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
1016:       MPI_Pack(D,ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
1017:     } else {
1018:       MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
1019:       MPI_Pack(B+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
1020:     }
1021:     if (ds->state>DS_STATE_RAW) MPI_Pack(Q+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
1022:     if (eigr) MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
1023: #if !defined(PETSC_USE_COMPLEX)
1024:     if (eigi) MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
1025: #endif
1026:   }
1027:   MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
1028:   if (rank) {
1029:     if (ds->compact) {
1030:       MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds));
1031:       MPI_Unpack(ds->work,size,&off,D,ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds));
1032:     } else {
1033:       MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
1034:       MPI_Unpack(ds->work,size,&off,B+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
1035:     }
1036:     if (ds->state>DS_STATE_RAW) MPI_Unpack(ds->work,size,&off,Q+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
1037:     if (eigr) MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
1038: #if !defined(PETSC_USE_COMPLEX)
1039:     if (eigi) MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
1040: #endif
1041:   }
1042:   if (ds->compact) {
1043:     DSRestoreArrayReal(ds,DS_MAT_T,&T);
1044:     DSRestoreArrayReal(ds,DS_MAT_D,&D);
1045:   } else {
1046:     MatDenseRestoreArray(ds->omat[DS_MAT_A],&A);
1047:     MatDenseRestoreArray(ds->omat[DS_MAT_B],&B);
1048:   }
1049:   if (ds->state>DS_STATE_RAW) MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q);
1050:   return 0;
1051: }
1052: #endif

1054: PetscErrorCode DSHermitian_GHIEP(DS ds,DSMatType m,PetscBool *flg)
1055: {
1056:   if ((m==DS_MAT_A && !ds->extrarow) || m==DS_MAT_B) *flg = PETSC_TRUE;
1057:   else *flg = PETSC_FALSE;
1058:   return 0;
1059: }

1061: /*MC
1062:    DSGHIEP - Dense Generalized Hermitian Indefinite Eigenvalue Problem.

1064:    Level: beginner

1066:    Notes:
1067:    The problem is expressed as A*X = B*X*Lambda, where both A and B are
1068:    real symmetric (or complex Hermitian) and possibly indefinite. Lambda
1069:    is a diagonal matrix whose diagonal elements are the arguments of DSSolve().
1070:    After solve, A is overwritten with Lambda. Note that in the case of real
1071:    scalars, A is overwritten with a real representation of Lambda, i.e.,
1072:    complex conjugate eigenvalue pairs are stored as a 2x2 block in the
1073:    quasi-diagonal matrix.

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

1079:    Used DS matrices:
1080: +  DS_MAT_A - first problem matrix
1081: .  DS_MAT_B - second problem matrix
1082: .  DS_MAT_T - symmetric tridiagonal matrix of the reduced pencil
1083: .  DS_MAT_D - diagonal matrix (signature) of the reduced pencil
1084: -  DS_MAT_Q - pseudo-orthogonal transformation that reduces (A,B) to
1085:    tridiagonal-diagonal form (intermediate step) or a real basis of eigenvectors

1087:    Implemented methods:
1088: +  0 - QR iteration plus inverse iteration for the eigenvectors
1089: .  1 - HZ iteration
1090: -  2 - QR iteration plus pseudo-orthogonalization for the eigenvectors

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

1096: .seealso: DSCreate(), DSSetType(), DSType
1097: M*/
1098: SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS ds)
1099: {
1100:   ds->ops->allocate        = DSAllocate_GHIEP;
1101:   ds->ops->view            = DSView_GHIEP;
1102:   ds->ops->vectors         = DSVectors_GHIEP;
1103:   ds->ops->solve[0]        = DSSolve_GHIEP_QR_II;
1104:   ds->ops->solve[1]        = DSSolve_GHIEP_HZ;
1105:   ds->ops->solve[2]        = DSSolve_GHIEP_QR;
1106:   ds->ops->sort            = DSSort_GHIEP;
1107: #if !defined(PETSC_HAVE_MPIUNI)
1108:   ds->ops->synchronize     = DSSynchronize_GHIEP;
1109: #endif
1110:   ds->ops->gettruncatesize = DSGetTruncateSize_GHIEP;
1111:   ds->ops->truncate        = DSTruncate_GHIEP;
1112:   ds->ops->update          = DSUpdateExtraRow_GHIEP;
1113:   ds->ops->hermitian       = DSHermitian_GHIEP;
1114:   return 0;
1115: }