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         399 : static PetscErrorCode DSAllocate_NHEP(DS ds,PetscInt ld)
      15             : {
      16         399 :   PetscFunctionBegin;
      17         399 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
      18         399 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q));
      19         399 :   PetscCall(PetscFree(ds->perm));
      20         399 :   PetscCall(PetscMalloc1(ld,&ds->perm));
      21         399 :   PetscFunctionReturn(PETSC_SUCCESS);
      22             : }
      23             : 
      24          10 : static PetscErrorCode DSView_NHEP(DS ds,PetscViewer viewer)
      25             : {
      26          10 :   PetscViewerFormat format;
      27             : 
      28          10 :   PetscFunctionBegin;
      29          10 :   PetscCall(PetscViewerGetFormat(viewer,&format));
      30          10 :   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
      31           0 :   PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
      32           0 :   if (ds->state>DS_STATE_INTERMEDIATE) PetscCall(DSViewMat(ds,viewer,DS_MAT_Q));
      33           0 :   if (ds->omat[DS_MAT_X]) PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
      34           0 :   if (ds->omat[DS_MAT_Y]) PetscCall(DSViewMat(ds,viewer,DS_MAT_Y));
      35           0 :   PetscFunctionReturn(PETSC_SUCCESS);
      36             : }
      37             : 
      38          51 : static PetscErrorCode DSVectors_NHEP_Refined_Some(DS ds,PetscInt *k,PetscReal *rnorm,PetscBool left)
      39             : {
      40          51 :   PetscInt          i,j;
      41          51 :   PetscBLASInt      info,ld,n,n1,lwork,inc=1;
      42          51 :   PetscScalar       sdummy,done=1.0,zero=0.0;
      43          51 :   PetscReal         *sigma;
      44          51 :   PetscBool         iscomplex = PETSC_FALSE;
      45          51 :   PetscScalar       *X,*W;
      46          51 :   const PetscScalar *A,*Q;
      47             : 
      48          51 :   PetscFunctionBegin;
      49          51 :   PetscCheck(!left,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for left vectors");
      50          51 :   PetscCall(PetscBLASIntCast(ds->n,&n));
      51          51 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
      52          51 :   n1 = n+1;
      53          51 :   PetscCall(DSAllocateWork_Private(ds,5*ld,6*ld,0));
      54          51 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
      55          51 :   lwork = 5*ld;
      56          51 :   sigma = ds->rwork+5*ld;
      57             : 
      58             :   /* build A-w*I in W */
      59          51 :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
      60          51 :   PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&W));
      61          51 :   if ((*k)<n-1 && A[(*k)+1+(*k)*ld]!=0.0) iscomplex = PETSC_TRUE;
      62          51 :   PetscCheck(!iscomplex,PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for complex eigenvalues yet");
      63        2152 :   for (j=0;j<n;j++)
      64       92403 :     for (i=0;i<=n;i++)
      65       90302 :       W[i+j*ld] = A[i+j*ld];
      66        2152 :   for (i=0;i<n;i++)
      67        2101 :     W[i+i*ld] -= A[(*k)+(*k)*ld];
      68          51 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
      69             : 
      70             :   /* compute SVD of W */
      71             : #if !defined(PETSC_USE_COMPLEX)
      72             :   PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&n1,&n,W,&ld,sigma,&sdummy,&ld,&sdummy,&ld,ds->work,&lwork,&info));
      73             : #else
      74          51 :   PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&n1,&n,W,&ld,sigma,&sdummy,&ld,&sdummy,&ld,ds->work,&lwork,ds->rwork,&info));
      75             : #endif
      76          51 :   SlepcCheckLapackInfo("gesvd",info);
      77             : 
      78             :   /* the smallest singular value is the new error estimate */
      79          51 :   if (rnorm) *rnorm = sigma[n-1];
      80             : 
      81             :   /* update vector with right singular vector associated to smallest singular value,
      82             :      accumulating the transformation matrix Q */
      83          51 :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
      84         102 :   PetscCall(MatDenseGetArray(ds->omat[left?DS_MAT_Y:DS_MAT_X],&X));
      85          51 :   PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&done,Q,&ld,W+n-1,&ld,&zero,X+(*k)*ld,&inc));
      86          51 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
      87          51 :   PetscCall(MatDenseRestoreArray(ds->omat[left?DS_MAT_Y:DS_MAT_X],&X));
      88          51 :   PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&W));
      89          51 :   PetscFunctionReturn(PETSC_SUCCESS);
      90             : }
      91             : 
      92           1 : static PetscErrorCode DSVectors_NHEP_Refined_All(DS ds,PetscBool left)
      93             : {
      94           1 :   PetscInt       i;
      95             : 
      96           1 :   PetscFunctionBegin;
      97           2 :   for (i=0;i<ds->n;i++) PetscCall(DSVectors_NHEP_Refined_Some(ds,&i,NULL,left));
      98           1 :   PetscFunctionReturn(PETSC_SUCCESS);
      99             : }
     100             : 
     101        6373 : static PetscErrorCode DSVectors_NHEP_Eigen_Some(DS ds,PetscInt *k,PetscReal *rnorm,PetscBool left)
     102             : {
     103        6373 :   PetscInt          i;
     104        6373 :   PetscBLASInt      mm=1,mout,info,ld,n,*select,inc=1,cols=1,zero=0;
     105        6373 :   PetscScalar       sone=1.0,szero=0.0;
     106        6373 :   PetscReal         norm,done=1.0;
     107        6373 :   PetscBool         iscomplex = PETSC_FALSE;
     108        6373 :   PetscScalar       *X,*Y;
     109        6373 :   const PetscScalar *A,*Q;
     110             : 
     111        6373 :   PetscFunctionBegin;
     112        6373 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     113        6373 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     114        6373 :   PetscCall(DSAllocateWork_Private(ds,0,0,ld));
     115        6373 :   select = ds->iwork;
     116      149244 :   for (i=0;i<n;i++) select[i] = (PetscBLASInt)PETSC_FALSE;
     117             : 
     118             :   /* compute k-th eigenvector Y of A */
     119        6373 :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
     120       12406 :   PetscCall(MatDenseGetArray(ds->omat[left?DS_MAT_Y:DS_MAT_X],&X));
     121        6373 :   Y = X+(*k)*ld;
     122        6373 :   select[*k] = (PetscBLASInt)PETSC_TRUE;
     123             : #if !defined(PETSC_USE_COMPLEX)
     124             :   if ((*k)<n-1 && A[(*k)+1+(*k)*ld]!=0.0) iscomplex = PETSC_TRUE;
     125             :   mm = iscomplex? 2: 1;
     126             :   if (iscomplex) select[(*k)+1] = (PetscBLASInt)PETSC_TRUE;
     127             :   PetscCall(DSAllocateWork_Private(ds,3*ld,0,0));
     128             :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_(left?"L":"R","S",select,&n,(PetscScalar*)A,&ld,Y,&ld,Y,&ld,&mm,&mout,ds->work,&info));
     129             : #else
     130        6373 :   PetscCall(DSAllocateWork_Private(ds,2*ld,ld,0));
     131       12406 :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_(left?"L":"R","S",select,&n,(PetscScalar*)A,&ld,Y,&ld,Y,&ld,&mm,&mout,ds->work,ds->rwork,&info));
     132             : #endif
     133        6373 :   SlepcCheckLapackInfo("trevc",info);
     134        6373 :   PetscCheck(mout==mm,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Inconsistent arguments");
     135        6373 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
     136             : 
     137             :   /* accumulate and normalize eigenvectors */
     138        6373 :   if (ds->state>=DS_STATE_CONDENSED) {
     139        6373 :     PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
     140        6373 :     PetscCall(PetscArraycpy(ds->work,Y,mout*ld));
     141        6373 :     PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,Q,&ld,ds->work,&inc,&szero,Y,&inc));
     142             : #if !defined(PETSC_USE_COMPLEX)
     143             :     if (iscomplex) PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,Q,&ld,ds->work+ld,&inc,&szero,Y+ld,&inc));
     144             : #endif
     145        6373 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
     146        6373 :     cols = 1;
     147        6373 :     norm = BLASnrm2_(&n,Y,&inc);
     148             : #if !defined(PETSC_USE_COMPLEX)
     149             :     if (iscomplex) {
     150             :       norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,Y+ld,&inc));
     151             :       cols = 2;
     152             :     }
     153             : #endif
     154        6373 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Y,&ld,&info));
     155        6373 :     SlepcCheckLapackInfo("lascl",info);
     156             :   }
     157             : 
     158             :   /* set output arguments */
     159        6373 :   if (iscomplex) (*k)++;
     160        6373 :   if (rnorm) {
     161        5263 :     if (iscomplex) *rnorm = SlepcAbsEigenvalue(Y[n-1],Y[n-1+ld]);
     162        5263 :     else *rnorm = PetscAbsScalar(Y[n-1]);
     163             :   }
     164        6373 :   PetscCall(MatDenseRestoreArray(ds->omat[left?DS_MAT_Y:DS_MAT_X],&X));
     165        6373 :   PetscFunctionReturn(PETSC_SUCCESS);
     166             : }
     167             : 
     168         485 : static PetscErrorCode DSVectors_NHEP_Eigen_All(DS ds,PetscBool left)
     169             : {
     170         485 :   PetscInt          i;
     171         485 :   PetscBLASInt      n,ld,mout,info,inc=1,cols,zero=0;
     172         485 :   PetscBool         iscomplex;
     173         485 :   PetscScalar       *X,*Y,*Z;
     174         485 :   const PetscScalar *A,*Q;
     175         485 :   PetscReal         norm,done=1.0;
     176         485 :   const char        *side,*back;
     177             : 
     178         485 :   PetscFunctionBegin;
     179         485 :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
     180         485 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     181         485 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     182         485 :   if (left) {
     183           1 :     X = NULL;
     184           1 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
     185             :     side = "L";
     186             :   } else {
     187         484 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
     188         484 :     Y = NULL;
     189         484 :     side = "R";
     190             :   }
     191         485 :   Z = left? Y: X;
     192         485 :   if (ds->state>=DS_STATE_CONDENSED) {
     193             :     /* DSSolve() has been called, backtransform with matrix Q */
     194         121 :     back = "B";
     195         121 :     PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
     196         121 :     PetscCall(PetscArraycpy(Z,Q,ld*ld));
     197         121 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
     198             :   } else back = "A";
     199             : #if !defined(PETSC_USE_COMPLEX)
     200             :   PetscCall(DSAllocateWork_Private(ds,3*ld,0,0));
     201             :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_(side,back,NULL,&n,(PetscScalar*)A,&ld,Y,&ld,X,&ld,&n,&mout,ds->work,&info));
     202             : #else
     203         485 :   PetscCall(DSAllocateWork_Private(ds,2*ld,ld,0));
     204         485 :   PetscCallBLAS("LAPACKtrevc",LAPACKtrevc_(side,back,NULL,&n,(PetscScalar*)A,&ld,Y,&ld,X,&ld,&n,&mout,ds->work,ds->rwork,&info));
     205             : #endif
     206         485 :   SlepcCheckLapackInfo("trevc",info);
     207             : 
     208             :   /* normalize eigenvectors */
     209        5110 :   for (i=0;i<n;i++) {
     210        4625 :     iscomplex = (i<n-1 && A[i+1+i*ld]!=0.0)? PETSC_TRUE: PETSC_FALSE;
     211        4625 :     cols = 1;
     212        4625 :     norm = BLASnrm2_(&n,Z+i*ld,&inc);
     213             : #if !defined(PETSC_USE_COMPLEX)
     214             :     if (iscomplex) {
     215             :       norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,Z+(i+1)*ld,&inc));
     216             :       cols = 2;
     217             :     }
     218             : #endif
     219        4625 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Z+i*ld,&ld,&info));
     220        4625 :     SlepcCheckLapackInfo("lascl",info);
     221        4625 :     if (iscomplex) i++;
     222             :   }
     223         485 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
     224         969 :   PetscCall(MatDenseRestoreArray(ds->omat[left?DS_MAT_Y:DS_MAT_X],&Z));
     225         485 :   PetscFunctionReturn(PETSC_SUCCESS);
     226             : }
     227             : 
     228        6909 : static PetscErrorCode DSVectors_NHEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
     229             : {
     230        6909 :   PetscFunctionBegin;
     231        6909 :   switch (mat) {
     232        6568 :     case DS_MAT_X:
     233        6568 :       if (ds->refined) {
     234          51 :         PetscCheck(ds->extrarow,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Refined vectors require activating the extra row");
     235          51 :         if (j) PetscCall(DSVectors_NHEP_Refined_Some(ds,j,rnorm,PETSC_FALSE));
     236           1 :         else PetscCall(DSVectors_NHEP_Refined_All(ds,PETSC_FALSE));
     237             :       } else {
     238        6517 :         if (j) PetscCall(DSVectors_NHEP_Eigen_Some(ds,j,rnorm,PETSC_FALSE));
     239         484 :         else PetscCall(DSVectors_NHEP_Eigen_All(ds,PETSC_FALSE));
     240             :       }
     241             :       break;
     242         341 :     case DS_MAT_Y:
     243         341 :       PetscCheck(!ds->refined,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
     244         341 :       if (j) PetscCall(DSVectors_NHEP_Eigen_Some(ds,j,rnorm,PETSC_TRUE));
     245           1 :       else PetscCall(DSVectors_NHEP_Eigen_All(ds,PETSC_TRUE));
     246             :       break;
     247           0 :     case DS_MAT_U:
     248             :     case DS_MAT_V:
     249           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
     250           0 :     default:
     251           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
     252             :   }
     253        6909 :   PetscFunctionReturn(PETSC_SUCCESS);
     254             : }
     255             : 
     256          31 : static PetscErrorCode DSSort_NHEP_Arbitrary(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
     257             : {
     258          31 :   PetscInt       i;
     259          31 :   PetscBLASInt   info,n,ld,mout,lwork,*selection;
     260          31 :   PetscScalar    *T,*Q,*work;
     261          31 :   PetscReal      dummy;
     262             : #if !defined(PETSC_USE_COMPLEX)
     263             :   PetscBLASInt   *iwork,liwork;
     264             : #endif
     265             : 
     266          31 :   PetscFunctionBegin;
     267          31 :   PetscCheck(k,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Must supply argument k");
     268          31 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&T));
     269          31 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
     270          31 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     271          31 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     272             : #if !defined(PETSC_USE_COMPLEX)
     273             :   lwork = n;
     274             :   liwork = 1;
     275             :   PetscCall(DSAllocateWork_Private(ds,lwork,0,liwork+n));
     276             :   work = ds->work;
     277             :   PetscCall(PetscBLASIntCast(ds->lwork,&lwork));
     278             :   selection = ds->iwork;
     279             :   iwork = ds->iwork + n;
     280             :   PetscCall(PetscBLASIntCast(ds->liwork-n,&liwork));
     281             : #else
     282          31 :   lwork = 1;
     283          31 :   PetscCall(DSAllocateWork_Private(ds,lwork,0,n));
     284          31 :   work = ds->work;
     285          31 :   selection = ds->iwork;
     286             : #endif
     287             :   /* Compute the selected eigenvalue to be in the leading position */
     288          31 :   PetscCall(DSSortEigenvalues_Private(ds,rr,ri,ds->perm,PETSC_FALSE));
     289          31 :   PetscCall(PetscArrayzero(selection,n));
     290         138 :   for (i=0;i<*k;i++) selection[ds->perm[i]] = 1;
     291             : #if !defined(PETSC_USE_COMPLEX)
     292             :   PetscCallBLAS("LAPACKtrsen",LAPACKtrsen_("N","V",selection,&n,T,&ld,Q,&ld,wr,wi,&mout,&dummy,&dummy,work,&lwork,iwork,&liwork,&info));
     293             : #else
     294          31 :   PetscCallBLAS("LAPACKtrsen",LAPACKtrsen_("N","V",selection,&n,T,&ld,Q,&ld,wr,&mout,&dummy,&dummy,work,&lwork,&info));
     295             : #endif
     296          31 :   SlepcCheckLapackInfo("trsen",info);
     297          31 :   *k = mout;
     298          31 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&T));
     299          31 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
     300          31 :   PetscFunctionReturn(PETSC_SUCCESS);
     301             : }
     302             : 
     303        4004 : static PetscErrorCode DSSort_NHEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
     304             : {
     305        4004 :   PetscFunctionBegin;
     306        4004 :   if (!rr || wr == rr) PetscCall(DSSort_NHEP_Total(ds,DS_MAT_A,DS_MAT_Q,wr,wi));
     307          31 :   else PetscCall(DSSort_NHEP_Arbitrary(ds,wr,wi,rr,ri,k));
     308        4004 :   PetscFunctionReturn(PETSC_SUCCESS);
     309             : }
     310             : 
     311           1 : static PetscErrorCode DSSortWithPermutation_NHEP(DS ds,PetscInt *perm,PetscScalar *wr,PetscScalar *wi)
     312             : {
     313           1 :   PetscFunctionBegin;
     314           1 :   PetscCall(DSSortWithPermutation_NHEP_Private(ds,perm,DS_MAT_A,DS_MAT_Q,wr,wi));
     315           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     316             : }
     317             : 
     318        2892 : static PetscErrorCode DSUpdateExtraRow_NHEP(DS ds)
     319             : {
     320        2892 :   PetscInt          i;
     321        2892 :   PetscBLASInt      n,ld,incx=1;
     322        2892 :   PetscScalar       *A,*x,*y,one=1.0,zero=0.0;
     323        2892 :   const PetscScalar *Q;
     324             : 
     325        2892 :   PetscFunctionBegin;
     326        2892 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     327        2892 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     328        2892 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     329        2892 :   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
     330        2892 :   PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
     331        2892 :   x = ds->work;
     332        2892 :   y = ds->work+ld;
     333       59092 :   for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
     334        2892 :   PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
     335       59092 :   for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
     336        2892 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     337        2892 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
     338        2892 :   ds->k = n;
     339        2892 :   PetscFunctionReturn(PETSC_SUCCESS);
     340             : }
     341             : 
     342        3563 : static PetscErrorCode DSSolve_NHEP(DS ds,PetscScalar *wr,PetscScalar *wi)
     343             : {
     344        3563 :   PetscFunctionBegin;
     345             : #if !defined(PETSC_USE_COMPLEX)
     346             :   PetscAssertPointer(wi,3);
     347             : #endif
     348        3563 :   PetscCall(DSSolve_NHEP_Private(ds,DS_MAT_A,DS_MAT_Q,wr,wi));
     349        3563 :   PetscFunctionReturn(PETSC_SUCCESS);
     350             : }
     351             : 
     352             : #if !defined(PETSC_HAVE_MPIUNI)
     353          31 : static PetscErrorCode DSSynchronize_NHEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
     354             : {
     355          31 :   PetscInt       ld=ds->ld,l=ds->l,k;
     356          31 :   PetscMPIInt    n,rank,off=0,size,ldn;
     357          31 :   PetscScalar    *A,*Q;
     358             : 
     359          31 :   PetscFunctionBegin;
     360          31 :   k = (ds->n-l)*ld;
     361          31 :   if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
     362          31 :   if (eigr) k += ds->n-l;
     363          31 :   if (eigi) k += ds->n-l;
     364          31 :   PetscCall(DSAllocateWork_Private(ds,k,0,0));
     365          31 :   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
     366          31 :   PetscCall(PetscMPIIntCast(ds->n-l,&n));
     367          31 :   PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
     368          31 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     369          31 :   if (ds->state>DS_STATE_RAW) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
     370          31 :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
     371          31 :   if (!rank) {
     372          15 :     PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     373          15 :     if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Pack(Q+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     374          15 :     if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     375             : #if !defined(PETSC_USE_COMPLEX)
     376             :     if (eigi) PetscCallMPI(MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     377             : #endif
     378             :   }
     379          62 :   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
     380          31 :   if (rank) {
     381          16 :     PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     382          16 :     if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Unpack(ds->work,size,&off,Q+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     383          16 :     if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     384             : #if !defined(PETSC_USE_COMPLEX)
     385             :     if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     386             : #endif
     387             :   }
     388          31 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     389          31 :   if (ds->state>DS_STATE_RAW) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
     390          31 :   PetscFunctionReturn(PETSC_SUCCESS);
     391             : }
     392             : #endif
     393             : 
     394        2026 : static PetscErrorCode DSTruncate_NHEP(DS ds,PetscInt n,PetscBool trim)
     395             : {
     396        2026 :   PetscInt    i,ld=ds->ld,l=ds->l;
     397        2026 :   PetscScalar *A;
     398             : 
     399        2026 :   PetscFunctionBegin;
     400        2026 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     401             : #if defined(PETSC_USE_DEBUG)
     402             :   /* make sure diagonal 2x2 block is not broken */
     403        2026 :   PetscCheck(ds->state<DS_STATE_CONDENSED || n==0 || n==ds->n || A[n+(n-1)*ld]==0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"The given size would break a 2x2 block, call DSGetTruncateSize() first");
     404             : #endif
     405        2026 :   if (trim) {
     406         308 :     if (ds->extrarow) {   /* clean extra row */
     407        5761 :       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
     408             :     }
     409         308 :     ds->l = 0;
     410         308 :     ds->k = 0;
     411         308 :     ds->n = n;
     412         308 :     ds->t = ds->n;   /* truncated length equal to the new dimension */
     413             :   } else {
     414        1718 :     if (ds->extrarow && ds->k==ds->n) {
     415             :       /* copy entries of extra row to the new position, then clean last row */
     416       18080 :       for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
     417       35187 :       for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
     418             :     }
     419        1718 :     ds->k = (ds->extrarow)? n: 0;
     420        1718 :     ds->t = ds->n;   /* truncated length equal to previous dimension */
     421        1718 :     ds->n = n;
     422             :   }
     423        2026 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     424        2026 :   PetscFunctionReturn(PETSC_SUCCESS);
     425             : }
     426             : 
     427          41 : static PetscErrorCode DSCond_NHEP(DS ds,PetscReal *cond)
     428             : {
     429          41 :   PetscScalar    *work;
     430          41 :   PetscReal      *rwork;
     431          41 :   PetscBLASInt   *ipiv;
     432          41 :   PetscBLASInt   lwork,info,n,ld;
     433          41 :   PetscReal      hn,hin;
     434          41 :   PetscScalar    *A;
     435             : 
     436          41 :   PetscFunctionBegin;
     437          41 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     438          41 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     439          41 :   lwork = 8*ld;
     440          41 :   PetscCall(DSAllocateWork_Private(ds,lwork,ld,ld));
     441          41 :   work  = ds->work;
     442          41 :   rwork = ds->rwork;
     443          41 :   ipiv  = ds->iwork;
     444             : 
     445             :   /* use workspace matrix W to avoid overwriting A */
     446          41 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
     447          41 :   PetscCall(MatCopy(ds->omat[DS_MAT_A],ds->omat[DS_MAT_W],SAME_NONZERO_PATTERN));
     448          41 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&A));
     449             : 
     450             :   /* norm of A */
     451          41 :   if (ds->state<DS_STATE_INTERMEDIATE) hn = LAPACKlange_("I",&n,&n,A,&ld,rwork);
     452          41 :   else hn = LAPACKlanhs_("I",&n,A,&ld,rwork);
     453             : 
     454             :   /* norm of inv(A) */
     455          41 :   PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n,&n,A,&ld,ipiv,&info));
     456          41 :   SlepcCheckLapackInfo("getrf",info);
     457          41 :   PetscCallBLAS("LAPACKgetri",LAPACKgetri_(&n,A,&ld,ipiv,work,&lwork,&info));
     458          41 :   SlepcCheckLapackInfo("getri",info);
     459          41 :   hin = LAPACKlange_("I",&n,&n,A,&ld,rwork);
     460          41 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&A));
     461             : 
     462          41 :   *cond = hn*hin;
     463          41 :   PetscFunctionReturn(PETSC_SUCCESS);
     464             : }
     465             : 
     466          80 : static PetscErrorCode DSTranslateHarmonic_NHEP(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *gin,PetscReal *gammaout)
     467             : {
     468          80 :   PetscInt          i,j;
     469          80 :   PetscBLASInt      *ipiv,info,n,ld,one=1,ncol;
     470          80 :   PetscScalar       *A,*B,*g=gin,*ghat,done=1.0,dmone=-1.0,dzero=0.0;
     471          80 :   const PetscScalar *Q;
     472          80 :   PetscReal         gamma=1.0;
     473             : 
     474          80 :   PetscFunctionBegin;
     475          80 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     476          80 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     477          80 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     478             : 
     479          80 :   if (!recover) {
     480             : 
     481          56 :     PetscCall(DSAllocateWork_Private(ds,0,0,ld));
     482          56 :     ipiv = ds->iwork;
     483          56 :     if (!g) {
     484          27 :       PetscCall(DSAllocateWork_Private(ds,ld,0,0));
     485          27 :       g = ds->work;
     486             :     }
     487             :     /* use workspace matrix W to factor A-tau*eye(n) */
     488          56 :     PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
     489          56 :     PetscCall(MatCopy(ds->omat[DS_MAT_A],ds->omat[DS_MAT_W],SAME_NONZERO_PATTERN));
     490          56 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&B));
     491             : 
     492             :     /* Vector g initially stores b = beta*e_n^T */
     493          56 :     PetscCall(PetscArrayzero(g,n));
     494          56 :     g[n-1] = beta;
     495             : 
     496             :     /* g = (A-tau*eye(n))'\b */
     497        1105 :     for (i=0;i<n;i++) B[i+i*ld] -= tau;
     498          56 :     PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n,&n,B,&ld,ipiv,&info));
     499          56 :     SlepcCheckLapackInfo("getrf",info);
     500          56 :     PetscCall(PetscLogFlops(2.0*n*n*n/3.0));
     501          56 :     PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("C",&n,&one,B,&ld,ipiv,g,&ld,&info));
     502          56 :     SlepcCheckLapackInfo("getrs",info);
     503          56 :     PetscCall(PetscLogFlops(2.0*n*n-n));
     504          56 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&B));
     505             : 
     506             :     /* A = A + g*b' */
     507        1105 :     for (i=0;i<n;i++) A[i+(n-1)*ld] += g[i]*beta;
     508             : 
     509             :   } else { /* recover */
     510             : 
     511          24 :     PetscCall(DSAllocateWork_Private(ds,ld,0,0));
     512          24 :     ghat = ds->work;
     513          24 :     PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
     514             : 
     515             :     /* g^ = -Q(:,idx)'*g */
     516          24 :     PetscCall(PetscBLASIntCast(ds->l+ds->k,&ncol));
     517          24 :     PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&ncol,&dmone,Q,&ld,g,&one,&dzero,ghat,&one));
     518             : 
     519             :     /* A = A + g^*b' */
     520         236 :     for (i=0;i<ds->l+ds->k;i++)
     521        2092 :       for (j=ds->l;j<ds->l+ds->k;j++)
     522        1880 :         A[i+j*ld] += ghat[i]*Q[n-1+j*ld]*beta;
     523             : 
     524             :     /* g~ = (I-Q(:,idx)*Q(:,idx)')*g = g+Q(:,idx)*g^ */
     525          24 :     PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&ncol,&done,Q,&ld,ghat,&one,&done,g,&one));
     526          24 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
     527             :   }
     528             : 
     529             :   /* Compute gamma factor */
     530          80 :   if (gammaout || (recover && ds->extrarow)) gamma = SlepcAbs(1.0,BLASnrm2_(&n,g,&one));
     531          80 :   if (gammaout) *gammaout = gamma;
     532          80 :   if (recover && ds->extrarow) {
     533         236 :     for (j=ds->l;j<ds->l+ds->k;j++) A[ds->n+j*ld] *= gamma;
     534             :   }
     535          80 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     536          80 :   PetscFunctionReturn(PETSC_SUCCESS);
     537             : }
     538             : 
     539             : /*MC
     540             :    DSNHEP - Dense Non-Hermitian Eigenvalue Problem.
     541             : 
     542             :    Level: beginner
     543             : 
     544             :    Notes:
     545             :    The problem is expressed as A*X = X*Lambda, where A is the input matrix.
     546             :    Lambda is a diagonal matrix whose diagonal elements are the arguments of
     547             :    DSSolve(). After solve, A is overwritten with the upper quasi-triangular
     548             :    matrix T of the (real) Schur form, A*Q = Q*T.
     549             : 
     550             :    In the intermediate state A is reduced to upper Hessenberg form.
     551             : 
     552             :    Computation of left eigenvectors is supported, but two-sided Krylov solvers
     553             :    usually rely on the related DSNHEPTS.
     554             : 
     555             :    Used DS matrices:
     556             : +  DS_MAT_A - problem matrix
     557             : -  DS_MAT_Q - orthogonal/unitary transformation that reduces to Hessenberg form
     558             :    (intermediate step) or matrix of orthogonal Schur vectors
     559             : 
     560             :    Implemented methods:
     561             : .  0 - Implicit QR (_hseqr)
     562             : 
     563             : .seealso: DSCreate(), DSSetType(), DSType
     564             : M*/
     565         871 : SLEPC_EXTERN PetscErrorCode DSCreate_NHEP(DS ds)
     566             : {
     567         871 :   PetscFunctionBegin;
     568         871 :   ds->ops->allocate        = DSAllocate_NHEP;
     569         871 :   ds->ops->view            = DSView_NHEP;
     570         871 :   ds->ops->vectors         = DSVectors_NHEP;
     571         871 :   ds->ops->solve[0]        = DSSolve_NHEP;
     572         871 :   ds->ops->sort            = DSSort_NHEP;
     573         871 :   ds->ops->sortperm        = DSSortWithPermutation_NHEP;
     574             : #if !defined(PETSC_HAVE_MPIUNI)
     575         871 :   ds->ops->synchronize     = DSSynchronize_NHEP;
     576             : #endif
     577         871 :   ds->ops->gettruncatesize = DSGetTruncateSize_Default;
     578         871 :   ds->ops->truncate        = DSTruncate_NHEP;
     579         871 :   ds->ops->update          = DSUpdateExtraRow_NHEP;
     580         871 :   ds->ops->cond            = DSCond_NHEP;
     581         871 :   ds->ops->transharm       = DSTranslateHarmonic_NHEP;
     582         871 :   PetscFunctionReturn(PETSC_SUCCESS);
     583             : }