Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : #include <slepc/private/dsimpl.h>
      12             : #include <slepcblaslapack.h>
      13             : 
      14             : typedef struct {
      15             :   PetscScalar *wr,*wi;     /* eigenvalues of B */
      16             : } DS_NHEPTS;
      17             : 
      18          13 : static PetscErrorCode DSAllocate_NHEPTS(DS ds,PetscInt ld)
      19             : {
      20          13 :   DS_NHEPTS      *ctx = (DS_NHEPTS*)ds->data;
      21             : 
      22          13 :   PetscFunctionBegin;
      23          13 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
      24          13 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
      25          13 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q));
      26          13 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_Z));
      27          13 :   PetscCall(PetscFree(ds->perm));
      28          13 :   PetscCall(PetscMalloc1(ld,&ds->perm));
      29          13 :   PetscCall(PetscMalloc1(ld,&ctx->wr));
      30             : #if !defined(PETSC_USE_COMPLEX)
      31             :   PetscCall(PetscMalloc1(ld,&ctx->wi));
      32             : #endif
      33          13 :   PetscFunctionReturn(PETSC_SUCCESS);
      34             : }
      35             : 
      36           0 : static PetscErrorCode DSView_NHEPTS(DS ds,PetscViewer viewer)
      37             : {
      38           0 :   PetscViewerFormat format;
      39             : 
      40           0 :   PetscFunctionBegin;
      41           0 :   PetscCall(PetscViewerGetFormat(viewer,&format));
      42           0 :   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
      43           0 :   PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
      44           0 :   PetscCall(DSViewMat(ds,viewer,DS_MAT_B));
      45           0 :   if (ds->state>DS_STATE_INTERMEDIATE) {
      46           0 :     PetscCall(DSViewMat(ds,viewer,DS_MAT_Q));
      47           0 :     PetscCall(DSViewMat(ds,viewer,DS_MAT_Z));
      48             :   }
      49           0 :   if (ds->omat[DS_MAT_X]) PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
      50           0 :   if (ds->omat[DS_MAT_Y]) PetscCall(DSViewMat(ds,viewer,DS_MAT_Y));
      51           0 :   PetscFunctionReturn(PETSC_SUCCESS);
      52             : }
      53             : 
      54         322 : static PetscErrorCode DSVectors_NHEPTS_Eigen_Some(DS ds,PetscInt *k,PetscReal *rnorm,PetscBool left)
      55             : {
      56         322 :   PetscInt          i;
      57         322 :   PetscBLASInt      mm=1,mout,info,ld,n,*select,inc=1,cols=1,zero=0;
      58         322 :   PetscScalar       sone=1.0,szero=0.0;
      59         322 :   PetscReal         norm,done=1.0;
      60         322 :   PetscBool         iscomplex = PETSC_FALSE;
      61         322 :   PetscScalar       *X,*Y;
      62         322 :   const PetscScalar *A,*Q;
      63             : 
      64         322 :   PetscFunctionBegin;
      65         322 :   PetscCall(PetscBLASIntCast(ds->n,&n));
      66         322 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
      67         322 :   PetscCall(DSAllocateWork_Private(ds,0,0,ld));
      68         322 :   select = ds->iwork;
      69        6492 :   for (i=0;i<n;i++) select[i] = (PetscBLASInt)PETSC_FALSE;
      70             : 
      71             :   /* compute k-th eigenvector Y of A */
      72         322 :   PetscCall(MatDenseGetArrayRead(ds->omat[left?DS_MAT_B:DS_MAT_A],&A));
      73         483 :   PetscCall(MatDenseGetArray(ds->omat[left?DS_MAT_Y:DS_MAT_X],&X));
      74         322 :   Y = X+(*k)*ld;
      75         322 :   select[*k] = (PetscBLASInt)PETSC_TRUE;
      76             : #if !defined(PETSC_USE_COMPLEX)
      77             :   if ((*k)<n-1 && A[(*k)+1+(*k)*ld]!=0.0) iscomplex = PETSC_TRUE;
      78             :   mm = iscomplex? 2: 1;
      79             :   if (iscomplex) select[(*k)+1] = (PetscBLASInt)PETSC_TRUE;
      80             :   PetscCall(DSAllocateWork_Private(ds,3*ld,0,0));
      81             :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","S",select,&n,(PetscScalar*)A,&ld,Y,&ld,Y,&ld,&mm,&mout,ds->work,&info));
      82             : #else
      83         322 :   PetscCall(DSAllocateWork_Private(ds,2*ld,ld,0));
      84         322 :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R","S",select,&n,(PetscScalar*)A,&ld,Y,&ld,Y,&ld,&mm,&mout,ds->work,ds->rwork,&info));
      85             : #endif
      86         322 :   SlepcCheckLapackInfo("trevc",info);
      87         322 :   PetscCheck(mout==mm,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Inconsistent arguments");
      88         322 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[left?DS_MAT_B:DS_MAT_A],&A));
      89             : 
      90             :   /* accumulate and normalize eigenvectors */
      91         322 :   if (ds->state>=DS_STATE_CONDENSED) {
      92         483 :     PetscCall(MatDenseGetArrayRead(ds->omat[left?DS_MAT_Z:DS_MAT_Q],&Q));
      93         322 :     PetscCall(PetscArraycpy(ds->work,Y,mout*ld));
      94         322 :     PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,Q,&ld,ds->work,&inc,&szero,Y,&inc));
      95             : #if !defined(PETSC_USE_COMPLEX)
      96             :     if (iscomplex) PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,Q,&ld,ds->work+ld,&inc,&szero,Y+ld,&inc));
      97             : #endif
      98         322 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[left?DS_MAT_Z:DS_MAT_Q],&Q));
      99         322 :     cols = 1;
     100         322 :     norm = BLASnrm2_(&n,Y,&inc);
     101             : #if !defined(PETSC_USE_COMPLEX)
     102             :     if (iscomplex) {
     103             :       norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,Y+ld,&inc));
     104             :       cols = 2;
     105             :     }
     106             : #endif
     107         322 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Y,&ld,&info));
     108         322 :     SlepcCheckLapackInfo("lascl",info);
     109             :   }
     110             : 
     111             :   /* set output arguments */
     112         322 :   if (iscomplex) (*k)++;
     113         322 :   if (rnorm) {
     114         322 :     if (iscomplex) *rnorm = SlepcAbsEigenvalue(Y[n-1],Y[n-1+ld]);
     115         322 :     else *rnorm = PetscAbsScalar(Y[n-1]);
     116             :   }
     117         322 :   PetscCall(MatDenseRestoreArray(ds->omat[left?DS_MAT_Y:DS_MAT_X],&X));
     118         322 :   PetscFunctionReturn(PETSC_SUCCESS);
     119             : }
     120             : 
     121          26 : static PetscErrorCode DSVectors_NHEPTS_Eigen_All(DS ds,PetscBool left)
     122             : {
     123          26 :   PetscInt          i;
     124          26 :   PetscBLASInt      n,ld,mout,info,inc=1,cols,zero=0;
     125          26 :   PetscBool         iscomplex;
     126          26 :   PetscScalar       *X;
     127          26 :   const PetscScalar *A;
     128          26 :   PetscReal         norm,done=1.0;
     129          26 :   const char        *back;
     130             : 
     131          26 :   PetscFunctionBegin;
     132          26 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     133          26 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     134          26 :   PetscCall(MatDenseGetArrayRead(ds->omat[left?DS_MAT_B:DS_MAT_A],&A));
     135          39 :   PetscCall(MatDenseGetArrayWrite(ds->omat[left?DS_MAT_Y:DS_MAT_X],&X));
     136          26 :   if (ds->state>=DS_STATE_CONDENSED) {
     137             :     /* DSSolve() has been called, backtransform with matrix Q */
     138           0 :     back = "B";
     139           0 :     PetscCall(MatCopy(ds->omat[left?DS_MAT_Z:DS_MAT_Q],ds->omat[left?DS_MAT_Y:DS_MAT_X],SAME_NONZERO_PATTERN));
     140             :   } else back = "A";
     141             : #if !defined(PETSC_USE_COMPLEX)
     142             :   PetscCall(DSAllocateWork_Private(ds,3*ld,0,0));
     143             :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R",back,NULL,&n,(PetscScalar*)A,&ld,X,&ld,X,&ld,&n,&mout,ds->work,&info));
     144             : #else
     145          26 :   PetscCall(DSAllocateWork_Private(ds,2*ld,ld,0));
     146          26 :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_("R",back,NULL,&n,(PetscScalar*)A,&ld,X,&ld,X,&ld,&n,&mout,ds->work,ds->rwork,&info));
     147             : #endif
     148          26 :   SlepcCheckLapackInfo("trevc",info);
     149             : 
     150             :   /* normalize eigenvectors */
     151         152 :   for (i=0;i<n;i++) {
     152         126 :     iscomplex = (i<n-1 && A[i+1+i*ld]!=0.0)? PETSC_TRUE: PETSC_FALSE;
     153         126 :     cols = 1;
     154         126 :     norm = BLASnrm2_(&n,X+i*ld,&inc);
     155             : #if !defined(PETSC_USE_COMPLEX)
     156             :     if (iscomplex) {
     157             :       norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,X+(i+1)*ld,&inc));
     158             :       cols = 2;
     159             :     }
     160             : #endif
     161         126 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,X+i*ld,&ld,&info));
     162         126 :     SlepcCheckLapackInfo("lascl",info);
     163         126 :     if (iscomplex) i++;
     164             :   }
     165          26 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[left?DS_MAT_B:DS_MAT_A],&A));
     166          26 :   PetscCall(MatDenseRestoreArrayWrite(ds->omat[left?DS_MAT_Y:DS_MAT_X],&X));
     167          26 :   PetscFunctionReturn(PETSC_SUCCESS);
     168             : }
     169             : 
     170         348 : static PetscErrorCode DSVectors_NHEPTS(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
     171             : {
     172         348 :   PetscFunctionBegin;
     173         348 :   switch (mat) {
     174         174 :     case DS_MAT_X:
     175         174 :       PetscCheck(!ds->refined,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
     176         174 :       if (j) PetscCall(DSVectors_NHEPTS_Eigen_Some(ds,j,rnorm,PETSC_FALSE));
     177          13 :       else PetscCall(DSVectors_NHEPTS_Eigen_All(ds,PETSC_FALSE));
     178             :       break;
     179         174 :     case DS_MAT_Y:
     180         174 :       PetscCheck(!ds->refined,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
     181         174 :       if (j) PetscCall(DSVectors_NHEPTS_Eigen_Some(ds,j,rnorm,PETSC_TRUE));
     182          13 :       else PetscCall(DSVectors_NHEPTS_Eigen_All(ds,PETSC_TRUE));
     183             :       break;
     184           0 :     case DS_MAT_U:
     185             :     case DS_MAT_V:
     186           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
     187           0 :     default:
     188           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
     189             :   }
     190         348 :   PetscFunctionReturn(PETSC_SUCCESS);
     191             : }
     192             : 
     193          82 : static PetscErrorCode DSSort_NHEPTS(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
     194             : {
     195          82 :   DS_NHEPTS      *ctx = (DS_NHEPTS*)ds->data;
     196          82 :   PetscInt       i,j,cont,id=0,*p,*idx,*idx2;
     197          82 :   PetscReal      s,t;
     198             : #if defined(PETSC_USE_COMPLEX)
     199          82 :   Mat            A,U;
     200             : #endif
     201             : 
     202          82 :   PetscFunctionBegin;
     203          82 :   PetscCheck(!rr || wr==rr,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
     204          82 :   PetscCall(PetscMalloc3(ds->ld,&idx,ds->ld,&idx2,ds->ld,&p));
     205          82 :   PetscCall(DSSort_NHEP_Total(ds,DS_MAT_A,DS_MAT_Q,wr,wi));
     206             : #if defined(PETSC_USE_COMPLEX)
     207          82 :   PetscCall(DSGetMat(ds,DS_MAT_B,&A));
     208          82 :   PetscCall(MatConjugate(A));
     209          82 :   PetscCall(DSRestoreMat(ds,DS_MAT_B,&A));
     210          82 :   PetscCall(DSGetMat(ds,DS_MAT_Z,&U));
     211          82 :   PetscCall(MatConjugate(U));
     212          82 :   PetscCall(DSRestoreMat(ds,DS_MAT_Z,&U));
     213        1650 :   for (i=0;i<ds->n;i++) ctx->wr[i] = PetscConj(ctx->wr[i]);
     214             : #endif
     215          82 :   PetscCall(DSSort_NHEP_Total(ds,DS_MAT_B,DS_MAT_Z,ctx->wr,ctx->wi));
     216             :   /* check correct eigenvalue correspondence */
     217             :   cont = 0;
     218        1650 :   for (i=0;i<ds->n;i++) {
     219        1568 :     if (SlepcAbsEigenvalue(ctx->wr[i]-wr[i],ctx->wi[i]-wi[i])>PETSC_SQRT_MACHINE_EPSILON) {idx2[cont] = i; idx[cont++] = i;}
     220        1568 :     p[i] = -1;
     221             :   }
     222          82 :   if (cont) {
     223          33 :     for (i=0;i<cont;i++) {
     224             :       t = PETSC_MAX_REAL;
     225         102 :       for (j=0;j<cont;j++) if (idx2[j]!=-1 && (s=SlepcAbsEigenvalue(ctx->wr[idx[j]]-wr[idx[i]],ctx->wi[idx[j]]-wi[idx[i]]))<t) { id = j; t = s; }
     226          24 :       p[idx[i]] = idx[id];
     227          24 :       idx2[id] = -1;
     228             :     }
     229         188 :     for (i=0;i<ds->n;i++) if (p[i]==-1) p[i] = i;
     230           9 :     PetscCall(DSSortWithPermutation_NHEP_Private(ds,p,DS_MAT_B,DS_MAT_Z,ctx->wr,ctx->wi));
     231             :   }
     232             : #if defined(PETSC_USE_COMPLEX)
     233          82 :   PetscCall(DSGetMat(ds,DS_MAT_B,&A));
     234          82 :   PetscCall(MatConjugate(A));
     235          82 :   PetscCall(DSRestoreMat(ds,DS_MAT_B,&A));
     236          82 :   PetscCall(DSGetMat(ds,DS_MAT_Z,&U));
     237          82 :   PetscCall(MatConjugate(U));
     238          82 :   PetscCall(DSRestoreMat(ds,DS_MAT_Z,&U));
     239             : #endif
     240          82 :   PetscCall(PetscFree3(idx,idx2,p));
     241          82 :   PetscFunctionReturn(PETSC_SUCCESS);
     242             : }
     243             : 
     244          82 : static PetscErrorCode DSUpdateExtraRow_NHEPTS(DS ds)
     245             : {
     246          82 :   PetscInt          i;
     247          82 :   PetscBLASInt      n,ld,incx=1;
     248          82 :   PetscScalar       *A,*x,*y,one=1.0,zero=0.0;
     249          82 :   const PetscScalar *Q;
     250             : 
     251          82 :   PetscFunctionBegin;
     252          82 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     253          82 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     254          82 :   PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
     255          82 :   x = ds->work;
     256          82 :   y = ds->work+ld;
     257          82 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     258          82 :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
     259        1650 :   for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
     260          82 :   PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
     261        1650 :   for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
     262          82 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     263          82 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
     264          82 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&A));
     265          82 :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Z],&Q));
     266        1650 :   for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
     267          82 :   PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
     268        1650 :   for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
     269          82 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&A));
     270          82 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Z],&Q));
     271          82 :   ds->k = n;
     272          82 :   PetscFunctionReturn(PETSC_SUCCESS);
     273             : }
     274             : 
     275          82 : static PetscErrorCode DSSolve_NHEPTS(DS ds,PetscScalar *wr,PetscScalar *wi)
     276             : {
     277          82 :   DS_NHEPTS      *ctx = (DS_NHEPTS*)ds->data;
     278             : 
     279          82 :   PetscFunctionBegin;
     280             : #if !defined(PETSC_USE_COMPLEX)
     281             :   PetscAssertPointer(wi,3);
     282             : #endif
     283          82 :   PetscCall(DSSolve_NHEP_Private(ds,DS_MAT_A,DS_MAT_Q,wr,wi));
     284          82 :   PetscCall(DSSolve_NHEP_Private(ds,DS_MAT_B,DS_MAT_Z,ctx->wr,ctx->wi));
     285          82 :   PetscFunctionReturn(PETSC_SUCCESS);
     286             : }
     287             : 
     288             : #if !defined(PETSC_HAVE_MPIUNI)
     289          28 : static PetscErrorCode DSSynchronize_NHEPTS(DS ds,PetscScalar eigr[],PetscScalar eigi[])
     290             : {
     291          28 :   PetscInt       ld=ds->ld,l=ds->l,k;
     292          28 :   PetscMPIInt    n,rank,off=0,size,ldn;
     293          28 :   DS_NHEPTS      *ctx = (DS_NHEPTS*)ds->data;
     294          28 :   PetscScalar    *A,*B,*Q,*Z;
     295             : 
     296          28 :   PetscFunctionBegin;
     297          28 :   k = 2*(ds->n-l)*ld;
     298          28 :   if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
     299          28 :   if (eigr) k += ds->n-l;
     300          28 :   if (eigi) k += ds->n-l;
     301          28 :   if (ctx->wr) k += ds->n-l;
     302          28 :   if (ctx->wi) k += ds->n-l;
     303          28 :   PetscCall(DSAllocateWork_Private(ds,k,0,0));
     304          28 :   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
     305          28 :   PetscCall(PetscMPIIntCast(ds->n-l,&n));
     306          28 :   PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
     307          28 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     308          28 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
     309          28 :   if (ds->state>DS_STATE_RAW) {
     310          28 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
     311          28 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Z],&Z));
     312             :   }
     313          28 :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
     314          28 :   if (!rank) {
     315          14 :     PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     316          14 :     PetscCallMPI(MPI_Pack(B+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     317          14 :     if (ds->state>DS_STATE_RAW) {
     318          14 :       PetscCallMPI(MPI_Pack(Q+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     319          14 :       PetscCallMPI(MPI_Pack(Z+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     320             :     }
     321          14 :     if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     322             : #if !defined(PETSC_USE_COMPLEX)
     323             :     if (eigi) PetscCallMPI(MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     324             : #endif
     325          14 :     if (ctx->wr) PetscCallMPI(MPI_Pack(ctx->wr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     326          14 :     if (ctx->wi) PetscCallMPI(MPI_Pack(ctx->wi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     327             :   }
     328          56 :   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
     329          28 :   if (rank) {
     330          14 :     PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     331          14 :     PetscCallMPI(MPI_Unpack(ds->work,size,&off,B+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     332          14 :     if (ds->state>DS_STATE_RAW) {
     333          14 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,Q+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     334          14 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,Z+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     335             :     }
     336          14 :     if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     337             : #if !defined(PETSC_USE_COMPLEX)
     338             :     if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     339             : #endif
     340          14 :     if (ctx->wr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,ctx->wr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     341          14 :     if (ctx->wi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,ctx->wi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     342             :   }
     343          28 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     344          28 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
     345          28 :   if (ds->state>DS_STATE_RAW) {
     346          28 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
     347          28 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Z],&Z));
     348             :   }
     349          28 :   PetscFunctionReturn(PETSC_SUCCESS);
     350             : }
     351             : #endif
     352             : 
     353          69 : static PetscErrorCode DSGetTruncateSize_NHEPTS(DS ds,PetscInt l,PetscInt n,PetscInt *k)
     354             : {
     355             : #if !defined(PETSC_USE_COMPLEX)
     356             :   const PetscScalar *A,*B;
     357             : #endif
     358             : 
     359          69 :   PetscFunctionBegin;
     360             : #if !defined(PETSC_USE_COMPLEX)
     361             :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
     362             :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
     363             :   if (A[l+(*k)+(l+(*k)-1)*ds->ld] != 0.0 || B[l+(*k)+(l+(*k)-1)*ds->ld] != 0.0) {
     364             :     if (l+(*k)<n-1) (*k)++;
     365             :     else (*k)--;
     366             :   }
     367             :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
     368             :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
     369             : #endif
     370          69 :   PetscFunctionReturn(PETSC_SUCCESS);
     371             : }
     372             : 
     373          82 : static PetscErrorCode DSTruncate_NHEPTS(DS ds,PetscInt n,PetscBool trim)
     374             : {
     375          82 :   PetscInt    i,ld=ds->ld,l=ds->l;
     376          82 :   PetscScalar *A,*B;
     377             : 
     378          82 :   PetscFunctionBegin;
     379          82 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     380          82 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
     381             : #if defined(PETSC_USE_DEBUG)
     382             :   /* make sure diagonal 2x2 block is not broken */
     383          82 :   PetscCheck(ds->state<DS_STATE_CONDENSED || n==0 || n==ds->n || A[n+(n-1)*ld]==0.0 || B[n+(n-1)*ld]==0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"The given size would break a 2x2 block, call DSGetTruncateSize() first");
     384             : #endif
     385          82 :   if (trim) {
     386          13 :     if (ds->extrarow) {   /* clean extra row */
     387         233 :       for (i=l;i<ds->n;i++) { A[ds->n+i*ld] = 0.0; B[ds->n+i*ld] = 0.0; }
     388             :     }
     389          13 :     ds->l = 0;
     390          13 :     ds->k = 0;
     391          13 :     ds->n = n;
     392          13 :     ds->t = ds->n;   /* truncated length equal to the new dimension */
     393             :   } else {
     394          69 :     if (ds->extrarow && ds->k==ds->n) {
     395             :       /* copy entries of extra row to the new position, then clean last row */
     396         686 :       for (i=l;i<n;i++) { A[n+i*ld] = A[ds->n+i*ld]; B[n+i*ld] = B[ds->n+i*ld]; }
     397        1343 :       for (i=l;i<ds->n;i++) { A[ds->n+i*ld] = 0.0; B[ds->n+i*ld] = 0.0; }
     398             :     }
     399          69 :     ds->k = (ds->extrarow)? n: 0;
     400          69 :     ds->t = ds->n;   /* truncated length equal to previous dimension */
     401          69 :     ds->n = n;
     402             :   }
     403          82 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     404          82 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
     405          82 :   PetscFunctionReturn(PETSC_SUCCESS);
     406             : }
     407             : 
     408          13 : static PetscErrorCode DSDestroy_NHEPTS(DS ds)
     409             : {
     410          13 :   DS_NHEPTS      *ctx = (DS_NHEPTS*)ds->data;
     411             : 
     412          13 :   PetscFunctionBegin;
     413          13 :   if (ctx->wr) PetscCall(PetscFree(ctx->wr));
     414          13 :   if (ctx->wi) PetscCall(PetscFree(ctx->wi));
     415          13 :   PetscCall(PetscFree(ds->data));
     416          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     417             : }
     418             : 
     419         682 : static PetscErrorCode DSMatGetSize_NHEPTS(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
     420             : {
     421         682 :   PetscFunctionBegin;
     422         682 :   *rows = ((t==DS_MAT_A || t==DS_MAT_B) && ds->extrarow)? ds->n+1: ds->n;
     423         682 :   *cols = ds->n;
     424         682 :   PetscFunctionReturn(PETSC_SUCCESS);
     425             : }
     426             : 
     427             : /*MC
     428             :    DSNHEPTS - Dense Non-Hermitian Eigenvalue Problem (special variant intended
     429             :    for two-sided Krylov solvers).
     430             : 
     431             :    Level: beginner
     432             : 
     433             :    Notes:
     434             :    Two related problems are solved, A*X = X*Lambda and B*Y = Y*Lambda', where A and
     435             :    B are supposed to come from the Arnoldi factorizations of a certain matrix and its
     436             :    (conjugate) transpose, respectively. Hence, in exact arithmetic the columns of Y
     437             :    are equal to the left eigenvectors of A. Lambda is a diagonal matrix whose diagonal
     438             :    elements are the arguments of DSSolve(). After solve, A is overwritten with the
     439             :    upper quasi-triangular matrix T of the (real) Schur form, A*Q = Q*T, and similarly
     440             :    another (real) Schur relation is computed, B*Z = Z*S, overwriting B.
     441             : 
     442             :    In the intermediate state A and B are reduced to upper Hessenberg form.
     443             : 
     444             :    When left eigenvectors DS_MAT_Y are requested, right eigenvectors of B are returned,
     445             :    while DS_MAT_X contains right eigenvectors of A.
     446             : 
     447             :    Used DS matrices:
     448             : +  DS_MAT_A - first problem matrix obtained from Arnoldi
     449             : .  DS_MAT_B - second problem matrix obtained from Arnoldi on the transpose
     450             : .  DS_MAT_Q - orthogonal/unitary transformation that reduces A to Hessenberg form
     451             :    (intermediate step) or matrix of orthogonal Schur vectors of A
     452             : -  DS_MAT_Z - orthogonal/unitary transformation that reduces B to Hessenberg form
     453             :    (intermediate step) or matrix of orthogonal Schur vectors of B
     454             : 
     455             :    Implemented methods:
     456             : .  0 - Implicit QR (_hseqr)
     457             : 
     458             : .seealso: DSCreate(), DSSetType(), DSType
     459             : M*/
     460          13 : SLEPC_EXTERN PetscErrorCode DSCreate_NHEPTS(DS ds)
     461             : {
     462          13 :   DS_NHEPTS      *ctx;
     463             : 
     464          13 :   PetscFunctionBegin;
     465          13 :   PetscCall(PetscNew(&ctx));
     466          13 :   ds->data = (void*)ctx;
     467             : 
     468          13 :   ds->ops->allocate        = DSAllocate_NHEPTS;
     469          13 :   ds->ops->view            = DSView_NHEPTS;
     470          13 :   ds->ops->vectors         = DSVectors_NHEPTS;
     471          13 :   ds->ops->solve[0]        = DSSolve_NHEPTS;
     472          13 :   ds->ops->sort            = DSSort_NHEPTS;
     473             : #if !defined(PETSC_HAVE_MPIUNI)
     474          13 :   ds->ops->synchronize     = DSSynchronize_NHEPTS;
     475             : #endif
     476          13 :   ds->ops->gettruncatesize = DSGetTruncateSize_NHEPTS;
     477          13 :   ds->ops->truncate        = DSTruncate_NHEPTS;
     478          13 :   ds->ops->update          = DSUpdateExtraRow_NHEPTS;
     479          13 :   ds->ops->destroy         = DSDestroy_NHEPTS;
     480          13 :   ds->ops->matgetsize      = DSMatGetSize_NHEPTS;
     481          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     482             : }