Actual source code: dsghiep.c

slepc-3.17.2 2022-08-09
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:   DSAllocateMatReal_Private(ds,DS_MAT_T);
 20:   DSAllocateMatReal_Private(ds,DS_MAT_D);
 21:   PetscFree(ds->perm);
 22:   PetscMalloc1(ld,&ds->perm);
 23:   PetscLogObjectMemory((PetscObject)ds,ld*sizeof(PetscInt));
 24:   PetscFunctionReturn(0);
 25: }

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

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

 73: PetscErrorCode DSView_GHIEP(DS ds,PetscViewer viewer)
 74: {
 75:   PetscViewerFormat format;
 76:   PetscInt          i,j;
 77:   PetscReal         value;
 78:   const char        *methodname[] = {
 79:                      "QR + Inverse Iteration",
 80:                      "HZ method",
 81:                      "QR"
 82:   };
 83:   const int         nmeth=sizeof(methodname)/sizeof(methodname[0]);

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

111:       PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",ds->n,ds->n);
112:       PetscViewerASCIIPrintf(viewer,"omega = zeros(%" PetscInt_FMT ",3);\n",3*ds->n);
113:       PetscViewerASCIIPrintf(viewer,"omega = [\n");
114:       for (i=0;i<ds->n;i++) PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "  %18.16e\n",i+1,i+1,(double)*(ds->rmat[DS_MAT_D]+i));
115:       PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(omega);\n",DSMatName[DS_MAT_B]);

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

149: static PetscErrorCode DSVectors_GHIEP_Eigen_Some(DS ds,PetscInt *idx,PetscReal *rnorm)
150: {
151:   PetscReal      b[4],M[4],d1,d2,s1,s2,e;
152:   PetscReal      scal1,scal2,wr1,wr2,wi,ep,norm;
153:   PetscScalar    *Q,*X,Y[4],alpha,zeroS = 0.0;
154:   PetscInt       k;
155:   PetscBLASInt   two = 2,n_,ld,one=1;
156: #if !defined(PETSC_USE_COMPLEX)
157:   PetscBLASInt   four=4;
158: #endif

160:   X = ds->mat[DS_MAT_X];
161:   Q = ds->mat[DS_MAT_Q];
162:   k = *idx;
163:   PetscBLASIntCast(ds->n,&n_);
164:   PetscBLASIntCast(ds->ld,&ld);
165:   if (k < ds->n-1) e = (ds->compact)?*(ds->rmat[DS_MAT_T]+ld+k):PetscRealPart(*(ds->mat[DS_MAT_A]+(k+1)+ld*k));
166:   else e = 0.0;
167:   if (e == 0.0) { /* Real */
168:     if (ds->state>=DS_STATE_CONDENSED) PetscArraycpy(X+k*ld,Q+k*ld,ld);
169:     else {
170:       PetscArrayzero(X+k*ds->ld,ds->ld);
171:       X[k+k*ds->ld] = 1.0;
172:     }
173:     if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
174:   } else { /* 2x2 block */
175:     if (ds->compact) {
176:       s1 = *(ds->rmat[DS_MAT_D]+k);
177:       d1 = *(ds->rmat[DS_MAT_T]+k);
178:       s2 = *(ds->rmat[DS_MAT_D]+k+1);
179:       d2 = *(ds->rmat[DS_MAT_T]+k+1);
180:     } else {
181:       s1 = PetscRealPart(*(ds->mat[DS_MAT_B]+k*ld+k));
182:       d1 = PetscRealPart(*(ds->mat[DS_MAT_A]+k+k*ld));
183:       s2 = PetscRealPart(*(ds->mat[DS_MAT_B]+(k+1)*ld+k+1));
184:       d2 = PetscRealPart(*(ds->mat[DS_MAT_A]+k+1+(k+1)*ld));
185:     }
186:     M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
187:     b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
188:     ep = LAPACKlamch_("S");
189:     /* Compute eigenvalues of the block */
190:     PetscStackCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
192:     /* Complex eigenvalues */
194:     wr1 /= scal1;
195:     wi  /= scal1;
196: #if !defined(PETSC_USE_COMPLEX)
197:     if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
198:       Y[0] = wr1-s2*d2; Y[1] = s2*e; Y[2] = wi; Y[3] = 0.0;
199:     } else {
200:       Y[0] = s1*e; Y[1] = wr1-s1*d1; Y[2] = 0.0; Y[3] = wi;
201:     }
202:     norm = BLASnrm2_(&four,Y,&one);
203:     norm = 1.0/norm;
204:     if (ds->state >= DS_STATE_CONDENSED) {
205:       alpha = norm;
206:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&two,&two,&alpha,ds->mat[DS_MAT_Q]+k*ld,&ld,Y,&two,&zeroS,X+k*ld,&ld));
207:       if (rnorm) *rnorm = SlepcAbsEigenvalue(X[ds->n-1+k*ld],X[ds->n-1+(k+1)*ld]);
208:     } else {
209:       PetscArrayzero(X+k*ld,2*ld);
210:       X[k*ld+k]       = Y[0]*norm;
211:       X[k*ld+k+1]     = Y[1]*norm;
212:       X[(k+1)*ld+k]   = Y[2]*norm;
213:       X[(k+1)*ld+k+1] = Y[3]*norm;
214:     }
215: #else
216:     if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
217:       Y[0] = PetscCMPLX(wr1-s2*d2,wi);
218:       Y[1] = s2*e;
219:     } else {
220:       Y[0] = s1*e;
221:       Y[1] = PetscCMPLX(wr1-s1*d1,wi);
222:     }
223:     norm = BLASnrm2_(&two,Y,&one);
224:     norm = 1.0/norm;
225:     if (ds->state >= DS_STATE_CONDENSED) {
226:       alpha = norm;
227:       PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n_,&two,&alpha,ds->mat[DS_MAT_Q]+k*ld,&ld,Y,&one,&zeroS,X+k*ld,&one));
228:       if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
229:     } else {
230:       PetscArrayzero(X+k*ld,2*ld);
231:       X[k*ld+k]   = Y[0]*norm;
232:       X[k*ld+k+1] = Y[1]*norm;
233:     }
234:     X[(k+1)*ld+k]   = PetscConj(X[k*ld+k]);
235:     X[(k+1)*ld+k+1] = PetscConj(X[k*ld+k+1]);
236: #endif
237:     (*idx)++;
238:   }
239:   PetscFunctionReturn(0);
240: }

242: PetscErrorCode DSVectors_GHIEP(DS ds,DSMatType mat,PetscInt *k,PetscReal *rnorm)
243: {
244:   PetscInt       i;
245:   PetscReal      e;

247:   switch (mat) {
248:     case DS_MAT_X:
249:     case DS_MAT_Y:
250:       if (k) DSVectors_GHIEP_Eigen_Some(ds,k,rnorm);
251:       else {
252:         for (i=0; i<ds->n; i++) {
253:           e = (ds->compact)?*(ds->rmat[DS_MAT_T]+ds->ld+i):PetscRealPart(*(ds->mat[DS_MAT_A]+(i+1)+ds->ld*i));
254:           if (e == 0.0) { /* real */
255:             if (ds->state >= DS_STATE_CONDENSED) PetscArraycpy(ds->mat[mat]+i*ds->ld,ds->mat[DS_MAT_Q]+i*ds->ld,ds->ld);
256:             else {
257:               PetscArrayzero(ds->mat[mat]+i*ds->ld,ds->ld);
258:               *(ds->mat[mat]+i+i*ds->ld) = 1.0;
259:             }
260:           } else DSVectors_GHIEP_Eigen_Some(ds,&i,rnorm);
261:         }
262:       }
263:       break;
264:     case DS_MAT_U:
265:     case DS_MAT_V:
266:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
267:     default:
268:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
269:   }
270:   PetscFunctionReturn(0);
271: }

273: /*
274:   Extract the eigenvalues contained in the block-diagonal of the indefinite problem.
275:   Only the index range n0..n1 is processed.
276: */
277: PetscErrorCode DSGHIEPComplexEigs(DS ds,PetscInt n0,PetscInt n1,PetscScalar *wr,PetscScalar *wi)
278: {
279:   PetscInt     k,ld;
280:   PetscBLASInt two=2;
281:   PetscScalar  *A,*B;
282:   PetscReal    *D,*T;
283:   PetscReal    b[4],M[4],d1,d2,s1,s2,e;
284:   PetscReal    scal1,scal2,ep,wr1,wr2,wi1;

286:   ld = ds->ld;
287:   A = ds->mat[DS_MAT_A];
288:   B = ds->mat[DS_MAT_B];
289:   D = ds->rmat[DS_MAT_D];
290:   T = ds->rmat[DS_MAT_T];
291:   for (k=n0;k<n1;k++) {
292:     if (k < n1-1) e = (ds->compact)?T[ld+k]:PetscRealPart(A[(k+1)+ld*k]);
293:     else e = 0.0;
294:     if (e==0.0) { /* real eigenvalue */
295:       wr[k] = (ds->compact)?T[k]/D[k]:A[k+k*ld]/B[k+k*ld];
296: #if !defined(PETSC_USE_COMPLEX)
297:       wi[k] = 0.0 ;
298: #endif
299:     } else { /* diagonal block */
300:       if (ds->compact) {
301:         s1 = D[k];
302:         d1 = T[k];
303:         s2 = D[k+1];
304:         d2 = T[k+1];
305:       } else {
306:         s1 = PetscRealPart(B[k*ld+k]);
307:         d1 = PetscRealPart(A[k+k*ld]);
308:         s2 = PetscRealPart(B[(k+1)*ld+k+1]);
309:         d2 = PetscRealPart(A[k+1+(k+1)*ld]);
310:       }
311:       M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
312:       b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
313:       ep = LAPACKlamch_("S");
314:       /* Compute eigenvalues of the block */
315:       PetscStackCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi1));
317:       if (wi1==0.0) { /* Real eigenvalues */
319:         wr[k] = wr1/scal1; wr[k+1] = wr2/scal2;
320: #if !defined(PETSC_USE_COMPLEX)
321:         wi[k] = wi[k+1] = 0.0;
322: #endif
323:       } else { /* Complex eigenvalues */
324: #if !defined(PETSC_USE_COMPLEX)
325:         wr[k]   = wr1/scal1;
326:         wr[k+1] = wr[k];
327:         wi[k]   = wi1/scal1;
328:         wi[k+1] = -wi[k];
329: #else
330:         wr[k]   = PetscCMPLX(wr1,wi1)/scal1;
331:         wr[k+1] = PetscConj(wr[k]);
332: #endif
333:       }
334:       k++;
335:     }
336:   }
337: #if defined(PETSC_USE_COMPLEX)
338:   if (wi) {
339:     for (k=n0;k<n1;k++) wi[k] = 0.0;
340:   }
341: #endif
342:   PetscFunctionReturn(0);
343: }

345: PetscErrorCode DSSort_GHIEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
346: {
347:   PetscInt       n,i,*perm;
348:   PetscReal      *d,*e,*s;

350: #if !defined(PETSC_USE_COMPLEX)
352: #endif
353:   n = ds->n;
354:   d = ds->rmat[DS_MAT_T];
355:   e = d + ds->ld;
356:   s = ds->rmat[DS_MAT_D];
357:   DSAllocateWork_Private(ds,ds->ld,ds->ld,0);
358:   perm = ds->perm;
359:   if (!rr) {
360:     rr = wr;
361:     ri = wi;
362:   }
363:   DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_TRUE);
364:   if (!ds->compact) DSSwitchFormat_GHIEP(ds,PETSC_TRUE);
365:   PetscArraycpy(ds->work,wr,n);
366:   for (i=ds->l;i<n;i++) wr[i] = *(ds->work+perm[i]);
367: #if !defined(PETSC_USE_COMPLEX)
368:   PetscArraycpy(ds->work,wi,n);
369:   for (i=ds->l;i<n;i++) wi[i] = *(ds->work+perm[i]);
370: #endif
371:   PetscArraycpy(ds->rwork,s,n);
372:   for (i=ds->l;i<n;i++) s[i] = *(ds->rwork+perm[i]);
373:   PetscArraycpy(ds->rwork,d,n);
374:   for (i=ds->l;i<n;i++) d[i] = *(ds->rwork+perm[i]);
375:   PetscArraycpy(ds->rwork,e,n-1);
376:   PetscArrayzero(e+ds->l,n-1-ds->l);
377:   for (i=ds->l;i<n-1;i++) {
378:     if (perm[i]<n-1) e[i] = *(ds->rwork+perm[i]);
379:   }
380:   if (!ds->compact) DSSwitchFormat_GHIEP(ds,PETSC_FALSE);
381:   DSPermuteColumns_Private(ds,ds->l,n,n,DS_MAT_Q,perm);
382:   PetscFunctionReturn(0);
383: }

385: PetscErrorCode DSUpdateExtraRow_GHIEP(DS ds)
386: {
387:   PetscInt       i;
388:   PetscBLASInt   n,ld,incx=1;
389:   PetscScalar    *A,*Q,*x,*y,one=1.0,zero=0.0;
390:   PetscReal      *b,*r,beta;

392:   PetscBLASIntCast(ds->n,&n);
393:   PetscBLASIntCast(ds->ld,&ld);
394:   A  = ds->mat[DS_MAT_A];
395:   Q  = ds->mat[DS_MAT_Q];
396:   b  = ds->rmat[DS_MAT_T]+ld;
397:   r  = ds->rmat[DS_MAT_T]+2*ld;

399:   if (ds->compact) {
400:     beta = b[n-1];   /* in compact, we assume all entries are zero except the last one */
401:     for (i=0;i<n;i++) r[i] = PetscRealPart(beta*Q[n-1+i*ld]);
402:     ds->k = n;
403:   } else {
404:     DSAllocateWork_Private(ds,2*ld,0,0);
405:     x = ds->work;
406:     y = ds->work+ld;
407:     for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
408:     PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
409:     for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
410:     ds->k = n;
411:   }
412:   PetscFunctionReturn(0);
413: }

415: /*
416:   Get eigenvectors with inverse iteration.
417:   The system matrix is in Hessenberg form.
418: */
419: PetscErrorCode DSGHIEPInverseIteration(DS ds,PetscScalar *wr,PetscScalar *wi)
420: {
421:   PetscInt       i,off;
422:   PetscBLASInt   *select,*infoC,ld,n1,mout,info;
423:   PetscScalar    *A,*B,*H,*X;
424:   PetscReal      *s,*d,*e;
425: #if defined(PETSC_USE_COMPLEX)
426:   PetscInt       j;
427: #endif

429:   PetscBLASIntCast(ds->ld,&ld);
430:   PetscBLASIntCast(ds->n-ds->l,&n1);
431:   DSAllocateWork_Private(ds,ld*ld+2*ld,ld,2*ld);
432:   DSAllocateMat_Private(ds,DS_MAT_W);
433:   A = ds->mat[DS_MAT_A];
434:   B = ds->mat[DS_MAT_B];
435:   H = ds->mat[DS_MAT_W];
436:   s = ds->rmat[DS_MAT_D];
437:   d = ds->rmat[DS_MAT_T];
438:   e = d + ld;
439:   select = ds->iwork;
440:   infoC = ds->iwork + ld;
441:   off = ds->l+ds->l*ld;
442:   if (ds->compact) {
443:     H[off] = d[ds->l]*s[ds->l];
444:     H[off+ld] = e[ds->l]*s[ds->l];
445:     for (i=ds->l+1;i<ds->n-1;i++) {
446:       H[i+(i-1)*ld] = e[i-1]*s[i];
447:       H[i+i*ld] = d[i]*s[i];
448:       H[i+(i+1)*ld] = e[i]*s[i];
449:     }
450:     H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
451:     H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
452:   } else {
453:     s[ds->l]  = PetscRealPart(B[off]);
454:     H[off]    = A[off]*s[ds->l];
455:     H[off+ld] = A[off+ld]*s[ds->l];
456:     for (i=ds->l+1;i<ds->n-1;i++) {
457:       s[i] = PetscRealPart(B[i+i*ld]);
458:       H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
459:       H[i+i*ld]     = A[i+i*ld]*s[i];
460:       H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
461:     }
462:     s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
463:     H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
464:     H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
465:   }
466:   DSAllocateMat_Private(ds,DS_MAT_X);
467:   X = ds->mat[DS_MAT_X];
468:   for (i=0;i<n1;i++) select[i] = 1;
469: #if !defined(PETSC_USE_COMPLEX)
470:   PetscStackCallBLAS("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));
471: #else
472:   PetscStackCallBLAS("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));

474:   /* Separate real and imaginary part of complex eigenvectors */
475:   for (j=ds->l;j<ds->n;j++) {
476:     if (PetscAbsReal(PetscImaginaryPart(wr[j])) > PetscAbsScalar(wr[j])*PETSC_SQRT_MACHINE_EPSILON) {
477:       for (i=ds->l;i<ds->n;i++) {
478:         X[i+(j+1)*ds->ld] = PetscImaginaryPart(X[i+j*ds->ld]);
479:         X[i+j*ds->ld] = PetscRealPart(X[i+j*ds->ld]);
480:       }
481:       j++;
482:     }
483:   }
484: #endif
485:   SlepcCheckLapackInfo("hsein",info);
486:   DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_TRUE);
487:   PetscFunctionReturn(0);
488: }

490: /*
491:    Undo 2x2 blocks that have real eigenvalues.
492: */
493: PetscErrorCode DSGHIEPRealBlocks(DS ds)
494: {
495:   PetscInt       i;
496:   PetscReal      e,d1,d2,s1,s2,ss1,ss2,t,dd,ss;
497:   PetscReal      maxy,ep,scal1,scal2,snorm;
498:   PetscReal      *T,*D,b[4],M[4],wr1,wr2,wi;
499:   PetscScalar    *A,*B,Y[4],oneS = 1.0,zeroS = 0.0;
500:   PetscBLASInt   m,two=2,ld;
501:   PetscBool      isreal;

503:   PetscBLASIntCast(ds->ld,&ld);
504:   PetscBLASIntCast(ds->n-ds->l,&m);
505:   A = ds->mat[DS_MAT_A];
506:   B = ds->mat[DS_MAT_B];
507:   T = ds->rmat[DS_MAT_T];
508:   D = ds->rmat[DS_MAT_D];
509:   DSAllocateWork_Private(ds,2*m,0,0);
510:   for (i=ds->l;i<ds->n-1;i++) {
511:     e = (ds->compact)?T[ld+i]:PetscRealPart(A[(i+1)+ld*i]);
512:     if (e != 0.0) { /* 2x2 block */
513:       if (ds->compact) {
514:         s1 = D[i];
515:         d1 = T[i];
516:         s2 = D[i+1];
517:         d2 = T[i+1];
518:       } else {
519:         s1 = PetscRealPart(B[i*ld+i]);
520:         d1 = PetscRealPart(A[i*ld+i]);
521:         s2 = PetscRealPart(B[(i+1)*ld+i+1]);
522:         d2 = PetscRealPart(A[(i+1)*ld+i+1]);
523:       }
524:       isreal = PETSC_FALSE;
525:       if (s1==s2) { /* apply a Jacobi rotation to compute the eigendecomposition */
526:         dd = d1-d2;
527:         if (2*PetscAbsReal(e) <= dd) {
528:           t = 2*e/dd;
529:           t = t/(1 + PetscSqrtReal(1+t*t));
530:         } else {
531:           t = dd/(2*e);
532:           ss = (t>=0)?1.0:-1.0;
533:           t = ss/(PetscAbsReal(t)+PetscSqrtReal(1+t*t));
534:         }
535:         Y[0] = 1/PetscSqrtReal(1 + t*t); Y[3] = Y[0]; /* c */
536:         Y[1] = Y[0]*t; Y[2] = -Y[1]; /* s */
537:         wr1 = d1+t*e; wr2 = d2-t*e;
538:         ss1 = s1; ss2 = s2;
539:         isreal = PETSC_TRUE;
540:       } else {
541:         ss1 = 1.0; ss2 = 1.0,
542:         M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
543:         b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
544:         ep = LAPACKlamch_("S");

546:         /* Compute eigenvalues of the block */
547:         PetscStackCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
548:         if (wi==0.0) { /* Real eigenvalues */
549:           isreal = PETSC_TRUE;
551:           wr1 /= scal1;
552:           wr2 /= scal2;
553:           if (PetscAbsReal(s1*d1-wr1)<PetscAbsReal(s2*d2-wr1)) {
554:             Y[0] = wr1-s2*d2;
555:             Y[1] = s2*e;
556:           } else {
557:             Y[0] = s1*e;
558:             Y[1] = wr1-s1*d1;
559:           }
560:           /* normalize with a signature*/
561:           maxy = PetscMax(PetscAbsScalar(Y[0]),PetscAbsScalar(Y[1]));
562:           scal1 = PetscRealPart(Y[0])/maxy;
563:           scal2 = PetscRealPart(Y[1])/maxy;
564:           snorm = scal1*scal1*s1 + scal2*scal2*s2;
565:           if (snorm<0) { ss1 = -1.0; snorm = -snorm; }
566:           snorm = maxy*PetscSqrtReal(snorm);
567:           Y[0] = Y[0]/snorm;
568:           Y[1] = Y[1]/snorm;
569:           if (PetscAbsReal(s1*d1-wr2)<PetscAbsReal(s2*d2-wr2)) {
570:             Y[2] = wr2-s2*d2;
571:             Y[3] = s2*e;
572:           } else {
573:             Y[2] = s1*e;
574:             Y[3] = wr2-s1*d1;
575:           }
576:           maxy = PetscMax(PetscAbsScalar(Y[2]),PetscAbsScalar(Y[3]));
577:           scal1 = PetscRealPart(Y[2])/maxy;
578:           scal2 = PetscRealPart(Y[3])/maxy;
579:           snorm = scal1*scal1*s1 + scal2*scal2*s2;
580:           if (snorm<0) { ss2 = -1.0; snorm = -snorm; }
581:           snorm = maxy*PetscSqrtReal(snorm); Y[2] = Y[2]/snorm; Y[3] = Y[3]/snorm;
582:         }
583:         wr1 *= ss1; wr2 *= ss2;
584:       }
585:       if (isreal) {
586:         if (ds->compact) {
587:           D[i]    = ss1;
588:           T[i]    = wr1;
589:           D[i+1]  = ss2;
590:           T[i+1]  = wr2;
591:           T[ld+i] = 0.0;
592:         } else {
593:           B[i*ld+i]       = ss1;
594:           A[i*ld+i]       = wr1;
595:           B[(i+1)*ld+i+1] = ss2;
596:           A[(i+1)*ld+i+1] = wr2;
597:           A[(i+1)+ld*i]   = 0.0;
598:           A[i+ld*(i+1)]   = 0.0;
599:         }
600:         PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&two,&two,&oneS,ds->mat[DS_MAT_Q]+ds->l+i*ld,&ld,Y,&two,&zeroS,ds->work,&m));
601:         PetscArraycpy(ds->mat[DS_MAT_Q]+ds->l+i*ld,ds->work,m);
602:         PetscArraycpy(ds->mat[DS_MAT_Q]+ds->l+(i+1)*ld,ds->work+m,m);
603:       }
604:       i++;
605:     }
606:   }
607:   PetscFunctionReturn(0);
608: }

610: PetscErrorCode DSSolve_GHIEP_QR_II(DS ds,PetscScalar *wr,PetscScalar *wi)
611: {
612:   PetscInt       i,off;
613:   PetscBLASInt   n1,ld,one,info,lwork;
614:   PetscScalar    *H,*A,*B,*Q;
615:   PetscReal      *d,*e,*s;
616: #if defined(PETSC_USE_COMPLEX)
617:   PetscInt       j;
618: #endif

620: #if !defined(PETSC_USE_COMPLEX)
622: #endif
623:   one = 1;
624:   PetscBLASIntCast(ds->n-ds->l,&n1);
625:   PetscBLASIntCast(ds->ld,&ld);
626:   off = ds->l + ds->l*ld;
627:   A = ds->mat[DS_MAT_A];
628:   B = ds->mat[DS_MAT_B];
629:   Q = ds->mat[DS_MAT_Q];
630:   d = ds->rmat[DS_MAT_T];
631:   e = ds->rmat[DS_MAT_T] + ld;
632:   s = ds->rmat[DS_MAT_D];
633: #if defined(PETSC_USE_DEBUG)
634:   /* Check signature */
635:   for (i=0;i<ds->n;i++) {
636:     PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
638:   }
639: #endif
640:   DSAllocateWork_Private(ds,ld*ld,2*ld,ld*2);
641:   lwork = ld*ld;

643:   /* Quick return if possible */
644:   if (n1 == 1) {
645:     for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
646:     DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
647:     if (!ds->compact) {
648:       d[ds->l] = PetscRealPart(A[off]);
649:       s[ds->l] = PetscRealPart(B[off]);
650:     }
651:     wr[ds->l] = d[ds->l]/s[ds->l];
652:     if (wi) wi[ds->l] = 0.0;
653:     PetscFunctionReturn(0);
654:   }
655:   /* Reduce to pseudotriadiagonal form */
656:   DSIntermediate_GHIEP(ds);

658:   /* Compute Eigenvalues (QR) */
659:   DSAllocateMat_Private(ds,DS_MAT_W);
660:   H = ds->mat[DS_MAT_W];
661:   if (ds->compact) {
662:     H[off]    = d[ds->l]*s[ds->l];
663:     H[off+ld] = e[ds->l]*s[ds->l];
664:     for (i=ds->l+1;i<ds->n-1;i++) {
665:       H[i+(i-1)*ld] = e[i-1]*s[i];
666:       H[i+i*ld]     = d[i]*s[i];
667:       H[i+(i+1)*ld] = e[i]*s[i];
668:     }
669:     H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
670:     H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
671:   } else {
672:     s[ds->l]  = PetscRealPart(B[off]);
673:     H[off]    = A[off]*s[ds->l];
674:     H[off+ld] = A[off+ld]*s[ds->l];
675:     for (i=ds->l+1;i<ds->n-1;i++) {
676:       s[i] = PetscRealPart(B[i+i*ld]);
677:       H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
678:       H[i+i*ld]     = A[i+i*ld]*s[i];
679:       H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
680:     }
681:     s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
682:     H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
683:     H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
684:   }

686: #if !defined(PETSC_USE_COMPLEX)
687:   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,ds->work,&lwork,&info));
688: #else
689:   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,NULL,&ld,ds->work,&lwork,&info));
690:   for (i=ds->l;i<ds->n;i++) if (PetscAbsReal(PetscImaginaryPart(wr[i]))<10*PETSC_MACHINE_EPSILON) wr[i] = PetscRealPart(wr[i]);
691:   /* Sort to have consecutive conjugate pairs */
692:   for (i=ds->l;i<ds->n;i++) {
693:       j=i+1;
694:       while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
695:       if (j==ds->n) {
697:         wr[i]=PetscRealPart(wr[i]);
698:       } else { /* complex eigenvalue */
699:         wr[j] = wr[i+1];
700:         if (PetscImaginaryPart(wr[i])<0) wr[i] = PetscConj(wr[i]);
701:         wr[i+1] = PetscConj(wr[i]);
702:         i++;
703:       }
704:   }
705: #endif
706:   SlepcCheckLapackInfo("hseqr",info);
707:   /* Compute Eigenvectors with Inverse Iteration */
708:   DSGHIEPInverseIteration(ds,wr,wi);

710:   /* Recover eigenvalues from diagonal */
711:   DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
712: #if defined(PETSC_USE_COMPLEX)
713:   if (wi) {
714:     for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
715:   }
716: #endif
717:   PetscFunctionReturn(0);
718: }

720: PetscErrorCode DSSolve_GHIEP_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
721: {
722:   PetscInt       i,j,off,nwu=0,n,lw,lwr,nwru=0;
723:   PetscBLASInt   n_,ld,info,lwork,ilo,ihi;
724:   PetscScalar    *H,*A,*B,*Q,*X;
725:   PetscReal      *d,*s,*scale,nrm,*rcde,*rcdv;
726: #if defined(PETSC_USE_COMPLEX)
727:   PetscInt       k;
728: #endif

730: #if !defined(PETSC_USE_COMPLEX)
732: #endif
733:   n = ds->n-ds->l;
734:   PetscBLASIntCast(n,&n_);
735:   PetscBLASIntCast(ds->ld,&ld);
736:   off = ds->l + ds->l*ld;
737:   A = ds->mat[DS_MAT_A];
738:   B = ds->mat[DS_MAT_B];
739:   Q = ds->mat[DS_MAT_Q];
740:   d = ds->rmat[DS_MAT_T];
741:   s = ds->rmat[DS_MAT_D];
742: #if defined(PETSC_USE_DEBUG)
743:   /* Check signature */
744:   for (i=0;i<ds->n;i++) {
745:     PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
747:   }
748: #endif
749:   lw = 14*ld+ld*ld;
750:   lwr = 7*ld;
751:   DSAllocateWork_Private(ds,lw,lwr,0);
752:   scale = ds->rwork+nwru;
753:   nwru += ld;
754:   rcde = ds->rwork+nwru;
755:   nwru += ld;
756:   rcdv = ds->rwork+nwru;
757:   /* Quick return if possible */
758:   if (n_ == 1) {
759:     for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
760:     DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
761:     if (!ds->compact) {
762:       d[ds->l] = PetscRealPart(A[off]);
763:       s[ds->l] = PetscRealPart(B[off]);
764:     }
765:     wr[ds->l] = d[ds->l]/s[ds->l];
766:     if (wi) wi[ds->l] = 0.0;
767:     PetscFunctionReturn(0);
768:   }

770:   /* Form pseudo-symmetric matrix */
771:   H =  ds->work+nwu;
772:   nwu += n*n;
773:   PetscArrayzero(H,n*n);
774:   if (ds->compact) {
775:     for (i=0;i<n-1;i++) {
776:       H[i+i*n]     = s[ds->l+i]*d[ds->l+i];
777:       H[i+1+i*n]   = s[ds->l+i+1]*d[ld+ds->l+i];
778:       H[i+(i+1)*n] = s[ds->l+i]*d[ld+ds->l+i];
779:     }
780:     H[n-1+(n-1)*n] = s[ds->l+n-1]*d[ds->l+n-1];
781:     for (i=0;i<ds->k-ds->l;i++) {
782:       H[ds->k-ds->l+i*n] = s[ds->k]*d[2*ld+ds->l+i];
783:       H[i+(ds->k-ds->l)*n] = s[i+ds->l]*d[2*ld+ds->l+i];
784:     }
785:   } else {
786:     for (j=0;j<n;j++) {
787:       for (i=0;i<n;i++) H[i+j*n] = B[off+i+i*ld]*A[off+i+j*ld];
788:     }
789:   }

791:   /* Compute eigenpairs */
792:   PetscBLASIntCast(lw-nwu,&lwork);
793:   DSAllocateMat_Private(ds,DS_MAT_X);
794:   X = ds->mat[DS_MAT_X];
795: #if !defined(PETSC_USE_COMPLEX)
796:   PetscStackCallBLAS("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));
797: #else
798:   PetscStackCallBLAS("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));

800:   /* Sort to have consecutive conjugate pairs
801:      Separate real and imaginary part of complex eigenvectors*/
802:   for (i=ds->l;i<ds->n;i++) {
803:     j=i+1;
804:     while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
805:     if (j==ds->n) {
807:       wr[i]=PetscRealPart(wr[i]); /* real eigenvalue */
808:       for (k=ds->l;k<ds->n;k++) {
809:         X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
810:       }
811:     } else { /* complex eigenvalue */
812:       if (j!=i+1) {
813:         wr[j] = wr[i+1];
814:         PetscArraycpy(X+j*ds->ld,X+(i+1)*ds->ld,ds->ld);
815:       }
816:       if (PetscImaginaryPart(wr[i])<0) {
817:         wr[i] = PetscConj(wr[i]);
818:         for (k=ds->l;k<ds->n;k++) {
819:           X[k+(i+1)*ds->ld] = -PetscImaginaryPart(X[k+i*ds->ld]);
820:           X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
821:         }
822:       } else {
823:         for (k=ds->l;k<ds->n;k++) {
824:           X[k+(i+1)*ds->ld] = PetscImaginaryPart(X[k+i*ds->ld]);
825:           X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
826:         }
827:       }
828:       wr[i+1] = PetscConj(wr[i]);
829:       i++;
830:     }
831:   }
832: #endif
833:   SlepcCheckLapackInfo("geevx",info);

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

838:   /* Recover eigenvalues from diagonal */
839:   DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
840: #if defined(PETSC_USE_COMPLEX)
841:   if (wi) {
842:     for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
843:   }
844: #endif
845:   PetscFunctionReturn(0);
846: }

848: PetscErrorCode DSGetTruncateSize_GHIEP(DS ds,PetscInt l,PetscInt n,PetscInt *k)
849: {
850:   PetscReal *T = ds->rmat[DS_MAT_T];

852:   if (T[l+(*k)-1+ds->ld] !=0.0) {
853:     if (l+(*k)<n-1) (*k)++;
854:     else (*k)--;
855:   }
856:   PetscFunctionReturn(0);
857: }

859: PetscErrorCode DSTruncate_GHIEP(DS ds,PetscInt n,PetscBool trim)
860: {
861:   PetscInt    i,ld=ds->ld,l=ds->l;
862:   PetscScalar *A = ds->mat[DS_MAT_A];
863:   PetscReal   *T = ds->rmat[DS_MAT_T],*b,*r,*omega;

865: #if defined(PETSC_USE_DEBUG)
866:   /* make sure diagonal 2x2 block is not broken */
868: #endif
869:   if (trim) {
870:     if (!ds->compact && ds->extrarow) {   /* clean extra row */
871:       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
872:     }
873:     ds->l = 0;
874:     ds->k = 0;
875:     ds->n = n;
876:     ds->t = ds->n;   /* truncated length equal to the new dimension */
877:   } else {
878:     if (!ds->compact && ds->extrarow && ds->k==ds->n) {
879:       /* copy entries of extra row to the new position, then clean last row */
880:       for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
881:       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
882:     }
883:     if (ds->compact) {
884:       b = T+ld;
885:       r = T+2*ld;
886:       omega = ds->rmat[DS_MAT_D];
887:       b[n-1] = r[n-1];
888:       b[n] = b[ds->n];
889:       omega[n] = omega[ds->n];
890:     }
891:     ds->k = (ds->extrarow)? n: 0;
892:     ds->t = ds->n;   /* truncated length equal to previous dimension */
893:     ds->n = n;
894:   }
895:   PetscFunctionReturn(0);
896: }

898: PetscErrorCode DSSynchronize_GHIEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
899: {
900:   PetscInt       ld=ds->ld,l=ds->l,k=0,kr=0;
901:   PetscMPIInt    n,rank,off=0,size,ldn,ld3,ld_;

903:   if (ds->compact) kr = 4*ld;
904:   else k = 2*(ds->n-l)*ld;
905:   if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
906:   if (eigr) k += (ds->n-l);
907:   if (eigi) k += (ds->n-l);
908:   DSAllocateWork_Private(ds,k+kr,0,0);
909:   PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size);
910:   PetscMPIIntCast(ds->n-l,&n);
911:   PetscMPIIntCast(ld*(ds->n-l),&ldn);
912:   PetscMPIIntCast(ld*3,&ld3);
913:   PetscMPIIntCast(ld,&ld_);
914:   MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
915:   if (!rank) {
916:     if (ds->compact) {
917:       MPI_Pack(ds->rmat[DS_MAT_T],ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
918:       MPI_Pack(ds->rmat[DS_MAT_D],ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
919:     } else {
920:       MPI_Pack(ds->mat[DS_MAT_A]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
921:       MPI_Pack(ds->mat[DS_MAT_B]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
922:     }
923:     if (ds->state>DS_STATE_RAW) MPI_Pack(ds->mat[DS_MAT_Q]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
924:     if (eigr) MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
925: #if !defined(PETSC_USE_COMPLEX)
926:     if (eigi) MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
927: #endif
928:   }
929:   MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
930:   if (rank) {
931:     if (ds->compact) {
932:       MPI_Unpack(ds->work,size,&off,ds->rmat[DS_MAT_T],ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds));
933:       MPI_Unpack(ds->work,size,&off,ds->rmat[DS_MAT_D],ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds));
934:     } else {
935:       MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_A]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
936:       MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_B]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
937:     }
938:     if (ds->state>DS_STATE_RAW) MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_Q]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
939:     if (eigr) MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
940: #if !defined(PETSC_USE_COMPLEX)
941:     if (eigi) MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
942: #endif
943:   }
944:   PetscFunctionReturn(0);
945: }

947: PetscErrorCode DSHermitian_GHIEP(DS ds,DSMatType m,PetscBool *flg)
948: {
949:   if (m==DS_MAT_A || m==DS_MAT_B) *flg = PETSC_TRUE;
950:   else *flg = PETSC_FALSE;
951:   PetscFunctionReturn(0);
952: }

954: /*MC
955:    DSGHIEP - Dense Generalized Hermitian Indefinite Eigenvalue Problem.

957:    Level: beginner

959:    Notes:
960:    The problem is expressed as A*X = B*X*Lambda, where both A and B are
961:    real symmetric (or complex Hermitian) and possibly indefinite. Lambda
962:    is a diagonal matrix whose diagonal elements are the arguments of DSSolve().
963:    After solve, A is overwritten with Lambda. Note that in the case of real
964:    scalars, A is overwritten with a real representation of Lambda, i.e.,
965:    complex conjugate eigenvalue pairs are stored as a 2x2 block in the
966:    quasi-diagonal matrix.

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

972:    Used DS matrices:
973: +  DS_MAT_A - first problem matrix
974: .  DS_MAT_B - second problem matrix
975: .  DS_MAT_T - symmetric tridiagonal matrix of the reduced pencil
976: .  DS_MAT_D - diagonal matrix (signature) of the reduced pencil
977: -  DS_MAT_Q - pseudo-orthogonal transformation that reduces (A,B) to
978:    tridiagonal-diagonal form (intermediate step) or a real basis of eigenvectors

980:    Implemented methods:
981: +  0 - QR iteration plus inverse iteration for the eigenvectors
982: .  1 - HZ iteration
983: -  2 - QR iteration plus pseudo-orthogonalization for the eigenvectors

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

989: .seealso: DSCreate(), DSSetType(), DSType
990: M*/
991: SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS ds)
992: {
993:   ds->ops->allocate        = DSAllocate_GHIEP;
994:   ds->ops->view            = DSView_GHIEP;
995:   ds->ops->vectors         = DSVectors_GHIEP;
996:   ds->ops->solve[0]        = DSSolve_GHIEP_QR_II;
997:   ds->ops->solve[1]        = DSSolve_GHIEP_HZ;
998:   ds->ops->solve[2]        = DSSolve_GHIEP_QR;
999:   ds->ops->sort            = DSSort_GHIEP;
1000:   ds->ops->synchronize     = DSSynchronize_GHIEP;
1001:   ds->ops->gettruncatesize = DSGetTruncateSize_GHIEP;
1002:   ds->ops->truncate        = DSTruncate_GHIEP;
1003:   ds->ops->update          = DSUpdateExtraRow_GHIEP;
1004:   ds->ops->hermitian       = DSHermitian_GHIEP;
1005:   PetscFunctionReturn(0);
1006: }