Actual source code: dsghiep.c

slepc-3.15.2 2021-09-20
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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: {

 19:   DSAllocateMat_Private(ds,DS_MAT_A);
 20:   DSAllocateMat_Private(ds,DS_MAT_B);
 21:   DSAllocateMat_Private(ds,DS_MAT_Q);
 22:   DSAllocateMatReal_Private(ds,DS_MAT_T);
 23:   DSAllocateMatReal_Private(ds,DS_MAT_D);
 24:   PetscFree(ds->perm);
 25:   PetscMalloc1(ld,&ds->perm);
 26:   PetscLogObjectMemory((PetscObject)ds,ld*sizeof(PetscInt));
 27:   return(0);
 28: }

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

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

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

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

122:       PetscViewerASCIIPrintf(viewer,"%% Size = %D %D\n",ds->n,ds->n);
123:       PetscViewerASCIIPrintf(viewer,"omega = zeros(%D,3);\n",3*ds->n);
124:       PetscViewerASCIIPrintf(viewer,"omega = [\n");
125:       for (i=0;i<ds->n;i++) {
126:         PetscViewerASCIIPrintf(viewer,"%D %D  %18.16e\n",i+1,i+1,(double)*(ds->rmat[DS_MAT_D]+i));
127:       }
128:       PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(omega);\n",DSMatName[DS_MAT_B]);

130:     } else {
131:       PetscViewerASCIIPrintf(viewer,"T\n");
132:       for (i=0;i<ds->n;i++) {
133:         for (j=0;j<ds->n;j++) {
134:           if (i==j) value = *(ds->rmat[DS_MAT_T]+i);
135:           else if (i==j+1 || j==i+1) value = *(ds->rmat[DS_MAT_T]+ds->ld+PetscMin(i,j));
136:           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));
137:           else value = 0.0;
138:           PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value);
139:         }
140:         PetscViewerASCIIPrintf(viewer,"\n");
141:       }
142:       PetscViewerASCIIPrintf(viewer,"omega\n");
143:       for (i=0;i<ds->n;i++) {
144:         for (j=0;j<ds->n;j++) {
145:           if (i==j) value = *(ds->rmat[DS_MAT_D]+i);
146:           else value = 0.0;
147:           PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value);
148:         }
149:         PetscViewerASCIIPrintf(viewer,"\n");
150:       }
151:     }
152:     PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
153:     PetscViewerFlush(viewer);
154:   } else {
155:     DSViewMat(ds,viewer,DS_MAT_A);
156:     DSViewMat(ds,viewer,DS_MAT_B);
157:   }
158:   if (ds->state>DS_STATE_INTERMEDIATE) { DSViewMat(ds,viewer,DS_MAT_Q); }
159:   return(0);
160: }

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

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

259: PetscErrorCode DSVectors_GHIEP(DS ds,DSMatType mat,PetscInt *k,PetscReal *rnorm)
260: {
261:   PetscInt       i;
262:   PetscReal      e;

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

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

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

369: PetscErrorCode DSSort_GHIEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
370: {
372:   PetscInt       n,i,*perm;
373:   PetscReal      *d,*e,*s;

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

411: PetscErrorCode DSUpdateExtraRow_GHIEP(DS ds)
412: {
414:   PetscInt       i;
415:   PetscBLASInt   n,ld,incx=1;
416:   PetscScalar    *A,*Q,*x,*y,one=1.0,zero=0.0;
417:   PetscReal      *b,*r,beta;

420:   PetscBLASIntCast(ds->n,&n);
421:   PetscBLASIntCast(ds->ld,&ld);
422:   A  = ds->mat[DS_MAT_A];
423:   Q  = ds->mat[DS_MAT_Q];
424:   b  = ds->rmat[DS_MAT_T]+ld;
425:   r  = ds->rmat[DS_MAT_T]+2*ld;

427:   if (ds->compact) {
428:     beta = b[n-1];   /* in compact, we assume all entries are zero except the last one */
429:     for (i=0;i<n;i++) r[i] = PetscRealPart(beta*Q[n-1+i*ld]);
430:     ds->k = n;
431:   } else {
432:     DSAllocateWork_Private(ds,2*ld,0,0);
433:     x = ds->work;
434:     y = ds->work+ld;
435:     for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
436:     PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
437:     for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
438:     ds->k = n;
439:   }
440:   return(0);
441: }

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

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

504:   /* Separate real and imaginary part of complex eigenvectors */
505:   for (j=ds->l;j<ds->n;j++) {
506:     if (PetscAbsReal(PetscImaginaryPart(wr[j])) > PetscAbsScalar(wr[j])*PETSC_SQRT_MACHINE_EPSILON) {
507:       for (i=ds->l;i<ds->n;i++) {
508:         X[i+(j+1)*ds->ld] = PetscImaginaryPart(X[i+j*ds->ld]);
509:         X[i+j*ds->ld] = PetscRealPart(X[i+j*ds->ld]);
510:       }
511:       j++;
512:     }
513:   }
514: #endif
515:   SlepcCheckLapackInfo("hsein",info);
516:   DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_TRUE);
517:   return(0);
518: }

520: /*
521:    Undo 2x2 blocks that have real eigenvalues.
522: */
523: PetscErrorCode DSGHIEPRealBlocks(DS ds)
524: {
526:   PetscInt       i;
527:   PetscReal      e,d1,d2,s1,s2,ss1,ss2,t,dd,ss;
528:   PetscReal      maxy,ep,scal1,scal2,snorm;
529:   PetscReal      *T,*D,b[4],M[4],wr1,wr2,wi;
530:   PetscScalar    *A,*B,Y[4],oneS = 1.0,zeroS = 0.0;
531:   PetscBLASInt   m,two=2,ld;
532:   PetscBool      isreal;

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

578:         /* Compute eigenvalues of the block */
579:         PetscStackCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
580:         if (wi==0.0) { /* Real eigenvalues */
581:           isreal = PETSC_TRUE;
582:           if (scal1<ep||scal2<ep) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
583:           wr1 /= scal1;
584:           wr2 /= scal2;
585:           if (PetscAbsReal(s1*d1-wr1)<PetscAbsReal(s2*d2-wr1)) {
586:             Y[0] = wr1-s2*d2;
587:             Y[1] = s2*e;
588:           } else {
589:             Y[0] = s1*e;
590:             Y[1] = wr1-s1*d1;
591:           }
592:           /* normalize with a signature*/
593:           maxy = PetscMax(PetscAbsScalar(Y[0]),PetscAbsScalar(Y[1]));
594:           scal1 = PetscRealPart(Y[0])/maxy;
595:           scal2 = PetscRealPart(Y[1])/maxy;
596:           snorm = scal1*scal1*s1 + scal2*scal2*s2;
597:           if (snorm<0) { ss1 = -1.0; snorm = -snorm; }
598:           snorm = maxy*PetscSqrtReal(snorm);
599:           Y[0] = Y[0]/snorm;
600:           Y[1] = Y[1]/snorm;
601:           if (PetscAbsReal(s1*d1-wr2)<PetscAbsReal(s2*d2-wr2)) {
602:             Y[2] = wr2-s2*d2;
603:             Y[3] = s2*e;
604:           } else {
605:             Y[2] = s1*e;
606:             Y[3] = wr2-s1*d1;
607:           }
608:           maxy = PetscMax(PetscAbsScalar(Y[2]),PetscAbsScalar(Y[3]));
609:           scal1 = PetscRealPart(Y[2])/maxy;
610:           scal2 = PetscRealPart(Y[3])/maxy;
611:           snorm = scal1*scal1*s1 + scal2*scal2*s2;
612:           if (snorm<0) { ss2 = -1.0; snorm = -snorm; }
613:           snorm = maxy*PetscSqrtReal(snorm); Y[2] = Y[2]/snorm; Y[3] = Y[3]/snorm;
614:         }
615:         wr1 *= ss1; wr2 *= ss2;
616:       }
617:       if (isreal) {
618:         if (ds->compact) {
619:           D[i]    = ss1;
620:           T[i]    = wr1;
621:           D[i+1]  = ss2;
622:           T[i+1]  = wr2;
623:           T[ld+i] = 0.0;
624:         } else {
625:           B[i*ld+i]       = ss1;
626:           A[i*ld+i]       = wr1;
627:           B[(i+1)*ld+i+1] = ss2;
628:           A[(i+1)*ld+i+1] = wr2;
629:           A[(i+1)+ld*i]   = 0.0;
630:           A[i+ld*(i+1)]   = 0.0;
631:         }
632:         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));
633:         PetscArraycpy(ds->mat[DS_MAT_Q]+ds->l+i*ld,ds->work,m);
634:         PetscArraycpy(ds->mat[DS_MAT_Q]+ds->l+(i+1)*ld,ds->work+m,m);
635:       }
636:       i++;
637:     }
638:   }
639:   return(0);
640: }

642: PetscErrorCode DSSolve_GHIEP_QR_II(DS ds,PetscScalar *wr,PetscScalar *wi)
643: {
645:   PetscInt       i,off;
646:   PetscBLASInt   n1,ld,one,info,lwork;
647:   PetscScalar    *H,*A,*B,*Q;
648:   PetscReal      *d,*e,*s;
649: #if defined(PETSC_USE_COMPLEX)
650:   PetscInt       j;
651: #endif

654: #if !defined(PETSC_USE_COMPLEX)
656: #endif
657:   one = 1;
658:   PetscBLASIntCast(ds->n-ds->l,&n1);
659:   PetscBLASIntCast(ds->ld,&ld);
660:   off = ds->l + ds->l*ld;
661:   A = ds->mat[DS_MAT_A];
662:   B = ds->mat[DS_MAT_B];
663:   Q = ds->mat[DS_MAT_Q];
664:   d = ds->rmat[DS_MAT_T];
665:   e = ds->rmat[DS_MAT_T] + ld;
666:   s = ds->rmat[DS_MAT_D];
667: #if defined(PETSC_USE_DEBUG)
668:   /* Check signature */
669:   for (i=0;i<ds->n;i++) {
670:     PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
671:     if (de != 1.0 && de != -1.0) SETERRQ(PETSC_COMM_SELF,1,"Diagonal elements of the signature matrix must be 1 or -1");
672:   }
673: #endif
674:   DSAllocateWork_Private(ds,ld*ld,2*ld,ld*2);
675:   lwork = ld*ld;

677:   /* Quick return if possible */
678:   if (n1 == 1) {
679:     for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
680:     DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
681:     if (!ds->compact) {
682:       d[ds->l] = PetscRealPart(A[off]);
683:       s[ds->l] = PetscRealPart(B[off]);
684:     }
685:     wr[ds->l] = d[ds->l]/s[ds->l];
686:     if (wi) wi[ds->l] = 0.0;
687:     return(0);
688:   }
689:   /* Reduce to pseudotriadiagonal form */
690:   DSIntermediate_GHIEP(ds);

692:   /* Compute Eigenvalues (QR)*/
693:   DSAllocateMat_Private(ds,DS_MAT_W);
694:   H = ds->mat[DS_MAT_W];
695:   if (ds->compact) {
696:     H[off]    = d[ds->l]*s[ds->l];
697:     H[off+ld] = e[ds->l]*s[ds->l];
698:     for (i=ds->l+1;i<ds->n-1;i++) {
699:       H[i+(i-1)*ld] = e[i-1]*s[i];
700:       H[i+i*ld]     = d[i]*s[i];
701:       H[i+(i+1)*ld] = e[i]*s[i];
702:     }
703:     H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
704:     H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
705:   } else {
706:     s[ds->l]  = PetscRealPart(B[off]);
707:     H[off]    = A[off]*s[ds->l];
708:     H[off+ld] = A[off+ld]*s[ds->l];
709:     for (i=ds->l+1;i<ds->n-1;i++) {
710:       s[i] = PetscRealPart(B[i+i*ld]);
711:       H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
712:       H[i+i*ld]     = A[i+i*ld]*s[i];
713:       H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
714:     }
715:     s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
716:     H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
717:     H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
718:   }

720: #if !defined(PETSC_USE_COMPLEX)
721:   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,ds->work,&lwork,&info));
722: #else
723:   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,NULL,&ld,ds->work,&lwork,&info));
724:   for (i=ds->l;i<ds->n;i++) if (PetscAbsReal(PetscImaginaryPart(wr[i]))<10*PETSC_MACHINE_EPSILON) wr[i] = PetscRealPart(wr[i]);
725:   /* Sort to have consecutive conjugate pairs */
726:   for (i=ds->l;i<ds->n;i++) {
727:       j=i+1;
728:       while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
729:       if (j==ds->n) {
730:         if (PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON) wr[i]=PetscRealPart(wr[i]);
731:         else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
732:       } else { /* complex eigenvalue */
733:         wr[j] = wr[i+1];
734:         if (PetscImaginaryPart(wr[i])<0) wr[i] = PetscConj(wr[i]);
735:         wr[i+1] = PetscConj(wr[i]);
736:         i++;
737:       }
738:   }
739: #endif
740:   SlepcCheckLapackInfo("hseqr",info);
741:   /* Compute Eigenvectors with Inverse Iteration */
742:   DSGHIEPInverseIteration(ds,wr,wi);

744:   /* Recover eigenvalues from diagonal */
745:   DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
746: #if defined(PETSC_USE_COMPLEX)
747:   if (wi) {
748:     for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
749:   }
750: #endif
751:   return(0);
752: }

754: PetscErrorCode DSSolve_GHIEP_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
755: {
757:   PetscInt       i,j,off,nwu=0,n,lw,lwr,nwru=0;
758:   PetscBLASInt   n_,ld,info,lwork,ilo,ihi;
759:   PetscScalar    *H,*A,*B,*Q,*X;
760:   PetscReal      *d,*s,*scale,nrm,*rcde,*rcdv;
761: #if defined(PETSC_USE_COMPLEX)
762:   PetscInt       k;
763: #endif

766: #if !defined(PETSC_USE_COMPLEX)
768: #endif
769:   n = ds->n-ds->l;
770:   PetscBLASIntCast(n,&n_);
771:   PetscBLASIntCast(ds->ld,&ld);
772:   off = ds->l + ds->l*ld;
773:   A = ds->mat[DS_MAT_A];
774:   B = ds->mat[DS_MAT_B];
775:   Q = ds->mat[DS_MAT_Q];
776:   d = ds->rmat[DS_MAT_T];
777:   s = ds->rmat[DS_MAT_D];
778: #if defined(PETSC_USE_DEBUG)
779:   /* Check signature */
780:   for (i=0;i<ds->n;i++) {
781:     PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
782:     if (de != 1.0 && de != -1.0) SETERRQ(PETSC_COMM_SELF,1,"Diagonal elements of the signature matrix must be 1 or -1");
783:   }
784: #endif
785:   lw = 14*ld+ld*ld;
786:   lwr = 7*ld;
787:   DSAllocateWork_Private(ds,lw,lwr,0);
788:   scale = ds->rwork+nwru;
789:   nwru += ld;
790:   rcde = ds->rwork+nwru;
791:   nwru += ld;
792:   rcdv = ds->rwork+nwru;
793:   nwru += ld;
794:   /* Quick return if possible */
795:   if (n_ == 1) {
796:     for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
797:     DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
798:     if (!ds->compact) {
799:       d[ds->l] = PetscRealPart(A[off]);
800:       s[ds->l] = PetscRealPart(B[off]);
801:     }
802:     wr[ds->l] = d[ds->l]/s[ds->l];
803:     if (wi) wi[ds->l] = 0.0;
804:     return(0);
805:   }

807:   /* Form pseudo-symmetric matrix */
808:   H =  ds->work+nwu;
809:   nwu += n*n;
810:   PetscArrayzero(H,n*n);
811:   if (ds->compact) {
812:     for (i=0;i<n-1;i++) {
813:       H[i+i*n]     = s[ds->l+i]*d[ds->l+i];
814:       H[i+1+i*n]   = s[ds->l+i+1]*d[ld+ds->l+i];
815:       H[i+(i+1)*n] = s[ds->l+i]*d[ld+ds->l+i];
816:     }
817:     H[n-1+(n-1)*n] = s[ds->l+n-1]*d[ds->l+n-1];
818:     for (i=0;i<ds->k-ds->l;i++) {
819:       H[ds->k-ds->l+i*n] = s[ds->k]*d[2*ld+ds->l+i];
820:       H[i+(ds->k-ds->l)*n] = s[i+ds->l]*d[2*ld+ds->l+i];
821:     }
822:   } else {
823:     for (j=0;j<n;j++) {
824:       for (i=0;i<n;i++) H[i+j*n] = B[off+i+i*ld]*A[off+i+j*ld];
825:     }
826:   }

828:   /* Compute eigenpairs */
829:   PetscBLASIntCast(lw-nwu,&lwork);
830:   DSAllocateMat_Private(ds,DS_MAT_X);
831:   X = ds->mat[DS_MAT_X];
832: #if !defined(PETSC_USE_COMPLEX)
833:   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));
834: #else
835:   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));

837:   /* Sort to have consecutive conjugate pairs
838:      Separate real and imaginary part of complex eigenvectors*/
839:   for (i=ds->l;i<ds->n;i++) {
840:     j=i+1;
841:     while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
842:     if (j==ds->n) {
843:       if (PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON) {
844:         wr[i]=PetscRealPart(wr[i]); /* real eigenvalue */
845:         for (k=ds->l;k<ds->n;k++) {
846:           X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
847:         }
848:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
849:     } else { /* complex eigenvalue */
850:       if (j!=i+1) {
851:         wr[j] = wr[i+1];
852:         PetscArraycpy(X+j*ds->ld,X+(i+1)*ds->ld,ds->ld);
853:       }
854:       if (PetscImaginaryPart(wr[i])<0) {
855:         wr[i] = PetscConj(wr[i]);
856:         for (k=ds->l;k<ds->n;k++) {
857:           X[k+(i+1)*ds->ld] = -PetscImaginaryPart(X[k+i*ds->ld]);
858:           X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
859:         }
860:       } else {
861:         for (k=ds->l;k<ds->n;k++) {
862:           X[k+(i+1)*ds->ld] = PetscImaginaryPart(X[k+i*ds->ld]);
863:           X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
864:         }
865:       }
866:       wr[i+1] = PetscConj(wr[i]);
867:       i++;
868:     }
869:   }
870: #endif
871:   SlepcCheckLapackInfo("geevx",info);

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

876:   /* Recover eigenvalues from diagonal */
877:   DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
878: #if defined(PETSC_USE_COMPLEX)
879:   if (wi) {
880:     for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
881:   }
882: #endif
883:   return(0);
884: }

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

891:   if (T[l+(*k)-1+ds->ld] !=0.0) {
892:     if (l+(*k)<n-1) (*k)++;
893:     else (*k)--;
894:   }
895:   return(0);
896: }

898: PetscErrorCode DSTruncate_GHIEP(DS ds,PetscInt n,PetscBool trim)
899: {
900:   PetscInt    i,ld=ds->ld,l=ds->l;
901:   PetscScalar *A = ds->mat[DS_MAT_A];
902:   PetscReal   *T = ds->rmat[DS_MAT_T],*b,*r,*omega;

905: #if defined(PETSC_USE_DEBUG)
906:   /* make sure diagonal 2x2 block is not broken */
907:   if (ds->state>=DS_STATE_CONDENSED && n>0 && n<ds->n && T[n-1+ld]!=0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"The given size would break a 2x2 block, call DSGetTruncateSize() first");
908: #endif
909:   if (trim) {
910:     if (!ds->compact && ds->extrarow) {   /* clean extra row */
911:       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
912:     }
913:     ds->l = 0;
914:     ds->k = 0;
915:     ds->n = n;
916:     ds->t = ds->n;   /* truncated length equal to the new dimension */
917:   } else {
918:     if (!ds->compact && ds->extrarow && ds->k==ds->n) {
919:       /* copy entries of extra row to the new position, then clean last row */
920:       for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
921:       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
922:     }
923:     if (ds->compact) {
924:       b = T+ld;
925:       r = T+2*ld;
926:       omega = ds->rmat[DS_MAT_D];
927:       b[n-1] = r[n-1];
928:       b[n] = b[ds->n];
929:       omega[n] = omega[ds->n];
930:     }
931:     ds->k = (ds->extrarow)? n: 0;
932:     ds->t = ds->n;   /* truncated length equal to previous dimension */
933:     ds->n = n;
934:   }
935:   return(0);
936: }

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

945:   if (ds->compact) kr = 4*ld;
946:   else k = 2*(ds->n-l)*ld;
947:   if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
948:   if (eigr) k += (ds->n-l);
949:   if (eigi) k += (ds->n-l);
950:   DSAllocateWork_Private(ds,k+kr,0,0);
951:   PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size);
952:   PetscMPIIntCast(ds->n-l,&n);
953:   PetscMPIIntCast(ld*(ds->n-l),&ldn);
954:   PetscMPIIntCast(ld*3,&ld3);
955:   PetscMPIIntCast(ld,&ld_);
956:   MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);CHKERRMPI(ierr);
957:   if (!rank) {
958:     if (ds->compact) {
959:       MPI_Pack(ds->rmat[DS_MAT_T],ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
960:       MPI_Pack(ds->rmat[DS_MAT_D],ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
961:     } else {
962:       MPI_Pack(ds->mat[DS_MAT_A]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
963:       MPI_Pack(ds->mat[DS_MAT_B]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
964:     }
965:     if (ds->state>DS_STATE_RAW) {
966:       MPI_Pack(ds->mat[DS_MAT_Q]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
967:     }
968:     if (eigr) {
969:       MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
970:     }
971:     if (eigi) {
972:       MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
973:     }
974:   }
975:   MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
976:   if (rank) {
977:     if (ds->compact) {
978:       MPI_Unpack(ds->work,size,&off,ds->rmat[DS_MAT_T],ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
979:       MPI_Unpack(ds->work,size,&off,ds->rmat[DS_MAT_D],ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
980:     } else {
981:       MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_A]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
982:       MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_B]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
983:     }
984:     if (ds->state>DS_STATE_RAW) {
985:       MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_Q]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
986:     }
987:     if (eigr) {
988:       MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
989:     }
990:     if (eigi) {
991:       MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));CHKERRMPI(ierr);
992:     }
993:   }
994:   return(0);
995: }

997: PetscErrorCode DSHermitian_GHIEP(DS ds,DSMatType m,PetscBool *flg)
998: {
1000:   if (m==DS_MAT_A || m==DS_MAT_B) *flg = PETSC_TRUE;
1001:   else *flg = PETSC_FALSE;
1002:   return(0);
1003: }

1005: SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS ds)
1006: {
1008:   ds->ops->allocate        = DSAllocate_GHIEP;
1009:   ds->ops->view            = DSView_GHIEP;
1010:   ds->ops->vectors         = DSVectors_GHIEP;
1011:   ds->ops->solve[0]        = DSSolve_GHIEP_QR_II;
1012:   ds->ops->solve[1]        = DSSolve_GHIEP_HZ;
1013:   ds->ops->solve[2]        = DSSolve_GHIEP_QR;
1014:   ds->ops->sort            = DSSort_GHIEP;
1015:   ds->ops->synchronize     = DSSynchronize_GHIEP;
1016:   ds->ops->gettruncatesize = DSGetTruncateSize_GHIEP;
1017:   ds->ops->truncate        = DSTruncate_GHIEP;
1018:   ds->ops->update          = DSUpdateExtraRow_GHIEP;
1019:   ds->ops->hermitian       = DSHermitian_GHIEP;
1020:   return(0);
1021: }