LCOV - code coverage report
Current view: top level - src/sys/classes/ds/impls/nhep - dsnhep.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 378 409 92.4 %
Date: 2019-03-05 01:06:09 Functions: 17 18 94.4 %

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-2018, 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         253 : PetscErrorCode DSAllocate_NHEP(DS ds,PetscInt ld)
      15             : {
      16             :   PetscErrorCode ierr;
      17             : 
      18         253 :   PetscFunctionBegin;
      19         253 :   ierr = DSAllocateMat_Private(ds,DS_MAT_A);CHKERRQ(ierr);
      20         253 :   ierr = DSAllocateMat_Private(ds,DS_MAT_Q);CHKERRQ(ierr);
      21         253 :   ierr = PetscFree(ds->perm);CHKERRQ(ierr);
      22         253 :   ierr = PetscMalloc1(ld,&ds->perm);CHKERRQ(ierr);
      23         253 :   ierr = PetscLogObjectMemory((PetscObject)ds,ld*sizeof(PetscInt));CHKERRQ(ierr);
      24         253 :   PetscFunctionReturn(0);
      25             : }
      26             : 
      27           4 : PetscErrorCode DSView_NHEP(DS ds,PetscViewer viewer)
      28             : {
      29             :   PetscErrorCode ierr;
      30             : 
      31           4 :   PetscFunctionBegin;
      32           4 :   ierr = DSViewMat(ds,viewer,DS_MAT_A);CHKERRQ(ierr);
      33           4 :   if (ds->state>DS_STATE_INTERMEDIATE) {
      34           1 :     ierr = DSViewMat(ds,viewer,DS_MAT_Q);CHKERRQ(ierr);
      35             :   }
      36           4 :   if (ds->mat[DS_MAT_X]) {
      37           3 :     ierr = DSViewMat(ds,viewer,DS_MAT_X);CHKERRQ(ierr);
      38             :   }
      39           4 :   if (ds->mat[DS_MAT_Y]) {
      40           0 :     ierr = DSViewMat(ds,viewer,DS_MAT_Y);CHKERRQ(ierr);
      41             :   }
      42           4 :   PetscFunctionReturn(0);
      43             : }
      44             : 
      45          78 : static PetscErrorCode DSVectors_NHEP_Refined_Some(DS ds,PetscInt *k,PetscReal *rnorm,PetscBool left)
      46             : {
      47             : #if defined(PETSC_MISSING_LAPACK_GESVD)
      48             :   PetscFunctionBegin;
      49             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable");
      50             : #else
      51             :   PetscErrorCode ierr;
      52             :   PetscInt       i,j;
      53          78 :   PetscBLASInt   info,ld,n,n1,lwork,inc=1;
      54          78 :   PetscScalar    sdummy,done=1.0,zero=0.0;
      55             :   PetscReal      *sigma;
      56          78 :   PetscBool      iscomplex = PETSC_FALSE;
      57          78 :   PetscScalar    *A = ds->mat[DS_MAT_A];
      58          78 :   PetscScalar    *Q = ds->mat[DS_MAT_Q];
      59          78 :   PetscScalar    *X = ds->mat[left?DS_MAT_Y:DS_MAT_X];
      60             :   PetscScalar    *W;
      61             : 
      62          78 :   PetscFunctionBegin;
      63          78 :   if (left) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for left vectors");
      64          78 :   ierr = PetscBLASIntCast(ds->n,&n);CHKERRQ(ierr);
      65          78 :   ierr = PetscBLASIntCast(ds->ld,&ld);CHKERRQ(ierr);
      66          78 :   n1 = n+1;
      67          78 :   if ((*k)<n-1 && A[(*k)+1+(*k)*ld]!=0.0) iscomplex = PETSC_TRUE;
      68          78 :   if (iscomplex) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for complex eigenvalues yet");
      69          78 :   ierr = DSAllocateWork_Private(ds,5*ld,6*ld,0);CHKERRQ(ierr);
      70          78 :   ierr = DSAllocateMat_Private(ds,DS_MAT_W);CHKERRQ(ierr);
      71          78 :   W = ds->mat[DS_MAT_W];
      72          78 :   lwork = 5*ld;
      73          78 :   sigma = ds->rwork+5*ld;
      74             : 
      75             :   /* build A-w*I in W */
      76        1311 :   for (j=0;j<n;j++)
      77       22179 :     for (i=0;i<=n;i++)
      78       20946 :       W[i+j*ld] = A[i+j*ld];
      79        1311 :   for (i=0;i<n;i++)
      80        1233 :     W[i+i*ld] -= A[(*k)+(*k)*ld];
      81             : 
      82             :   /* compute SVD of W */
      83             : #if !defined(PETSC_USE_COMPLEX)
      84          78 :   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&n1,&n,W,&ld,sigma,&sdummy,&ld,&sdummy,&ld,ds->work,&lwork,&info));
      85             : #else
      86             :   PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","O",&n1,&n,W,&ld,sigma,&sdummy,&ld,&sdummy,&ld,ds->work,&lwork,ds->rwork,&info));
      87             : #endif
      88          78 :   SlepcCheckLapackInfo("gesvd",info);
      89             : 
      90             :   /* the smallest singular value is the new error estimate */
      91          78 :   if (rnorm) *rnorm = sigma[n-1];
      92             : 
      93             :   /* update vector with right singular vector associated to smallest singular value,
      94             :      accumulating the transformation matrix Q */
      95          78 :   PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&done,Q,&ld,W+n-1,&ld,&zero,X+(*k)*ld,&inc));
      96          78 :   PetscFunctionReturn(0);
      97             : #endif
      98             : }
      99             : 
     100           1 : static PetscErrorCode DSVectors_NHEP_Refined_All(DS ds,PetscBool left)
     101             : {
     102             :   PetscErrorCode ierr;
     103             :   PetscInt       i;
     104             : 
     105           1 :   PetscFunctionBegin;
     106           2 :   for (i=0;i<ds->n;i++) {
     107           1 :     ierr = DSVectors_NHEP_Refined_Some(ds,&i,NULL,left);CHKERRQ(ierr);
     108             :   }
     109           1 :   PetscFunctionReturn(0);
     110             : }
     111             : 
     112        6488 : static PetscErrorCode DSVectors_NHEP_Eigen_Some(DS ds,PetscInt *k,PetscReal *rnorm,PetscBool left)
     113             : {
     114             : #if defined(SLEPC_MISSING_LAPACK_TREVC)
     115             :   PetscFunctionBegin;
     116             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TREVC - Lapack routine is unavailable");
     117             : #else
     118             :   PetscErrorCode ierr;
     119             :   PetscInt       i;
     120        6488 :   PetscBLASInt   mm=1,mout,info,ld,n,*select,inc = 1;
     121        6488 :   PetscScalar    tmp,done=1.0,zero=0.0;
     122             :   PetscReal      norm;
     123        6488 :   PetscBool      iscomplex = PETSC_FALSE;
     124        6488 :   PetscScalar    *A = ds->mat[DS_MAT_A];
     125        6488 :   PetscScalar    *Q = ds->mat[DS_MAT_Q];
     126        6488 :   PetscScalar    *X = ds->mat[left?DS_MAT_Y:DS_MAT_X];
     127             :   PetscScalar    *Y;
     128             : 
     129        6488 :   PetscFunctionBegin;
     130        6488 :   ierr = PetscBLASIntCast(ds->n,&n);CHKERRQ(ierr);
     131        6488 :   ierr = PetscBLASIntCast(ds->ld,&ld);CHKERRQ(ierr);
     132        6488 :   ierr = DSAllocateWork_Private(ds,0,0,ld);CHKERRQ(ierr);
     133        6488 :   select = ds->iwork;
     134        6488 :   for (i=0;i<n;i++) select[i] = (PetscBLASInt)PETSC_FALSE;
     135             : 
     136             :   /* compute k-th eigenvector Y of A */
     137        6488 :   Y = X+(*k)*ld;
     138        6488 :   select[*k] = (PetscBLASInt)PETSC_TRUE;
     139             : #if !defined(PETSC_USE_COMPLEX)
     140        6488 :   if ((*k)<n-1 && A[(*k)+1+(*k)*ld]!=0.0) iscomplex = PETSC_TRUE;
     141        6488 :   mm = iscomplex? 2: 1;
     142        6488 :   if (iscomplex) select[(*k)+1] = (PetscBLASInt)PETSC_TRUE;
     143        6488 :   ierr = DSAllocateWork_Private(ds,3*ld,0,0);CHKERRQ(ierr);
     144        6488 :   PetscStackCallBLAS("LAPACKtrevc",LAPACKtrevc_(left?"L":"R","S",select,&n,A,&ld,Y,&ld,Y,&ld,&mm,&mout,ds->work,&info));
     145             : #else
     146             :   ierr = DSAllocateWork_Private(ds,2*ld,ld,0);CHKERRQ(ierr);
     147             :   PetscStackCallBLAS("LAPACKtrevc",LAPACKtrevc_(left?"L":"R","S",select,&n,A,&ld,Y,&ld,Y,&ld,&mm,&mout,ds->work,ds->rwork,&info));
     148             : #endif
     149        6488 :   SlepcCheckLapackInfo("trevc",info);
     150        6488 :   if (mout != mm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Inconsistent arguments");
     151             : 
     152             :   /* accumulate and normalize eigenvectors */
     153        6488 :   if (ds->state>=DS_STATE_CONDENSED) {
     154        6488 :     ierr = PetscMemcpy(ds->work,Y,mout*ld*sizeof(PetscScalar));CHKERRQ(ierr);
     155        6488 :     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&done,Q,&ld,ds->work,&inc,&zero,Y,&inc));
     156             : #if !defined(PETSC_USE_COMPLEX)
     157        6488 :     if (iscomplex) PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&done,Q,&ld,ds->work+ld,&inc,&zero,Y+ld,&inc));
     158             : #endif
     159        6488 :     norm = BLASnrm2_(&n,Y,&inc);
     160             : #if !defined(PETSC_USE_COMPLEX)
     161        6488 :     if (iscomplex) {
     162         612 :       tmp = BLASnrm2_(&n,Y+ld,&inc);
     163         612 :       norm = SlepcAbsEigenvalue(norm,tmp);
     164             :     }
     165             : #endif
     166        6488 :     tmp = 1.0 / norm;
     167        6488 :     PetscStackCallBLAS("BLASscal",BLASscal_(&n,&tmp,Y,&inc));
     168             : #if !defined(PETSC_USE_COMPLEX)
     169        6488 :     if (iscomplex) PetscStackCallBLAS("BLASscal",BLASscal_(&n,&tmp,Y+ld,&inc));
     170             : #endif
     171             :   }
     172             : 
     173             :   /* set output arguments */
     174        6488 :   if (iscomplex) (*k)++;
     175        6488 :   if (rnorm) {
     176        3592 :     if (iscomplex) *rnorm = SlepcAbsEigenvalue(Y[n-1],Y[n-1+ld]);
     177        3163 :     else *rnorm = PetscAbsScalar(Y[n-1]);
     178             :   }
     179        6488 :   PetscFunctionReturn(0);
     180             : #endif
     181             : }
     182             : 
     183         256 : static PetscErrorCode DSVectors_NHEP_Eigen_All(DS ds,PetscBool left)
     184             : {
     185             : #if defined(SLEPC_MISSING_LAPACK_TREVC)
     186             :   PetscFunctionBegin;
     187             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TREVC - Lapack routine is unavailable");
     188             : #else
     189             :   PetscErrorCode ierr;
     190             :   PetscInt       i;
     191         256 :   PetscBLASInt   n,ld,mout,info,inc = 1;
     192             :   PetscBool      iscomplex;
     193         256 :   PetscScalar    *X,*Y,*Z,*A = ds->mat[DS_MAT_A],tmp;
     194             :   PetscReal      norm;
     195             :   const char     *side,*back;
     196             : 
     197         256 :   PetscFunctionBegin;
     198         256 :   ierr = PetscBLASIntCast(ds->n,&n);CHKERRQ(ierr);
     199         256 :   ierr = PetscBLASIntCast(ds->ld,&ld);CHKERRQ(ierr);
     200         256 :   if (left) {
     201           0 :     X = NULL;
     202           0 :     Y = ds->mat[DS_MAT_Y];
     203           0 :     side = "L";
     204             :   } else {
     205         256 :     X = ds->mat[DS_MAT_X];
     206         256 :     Y = NULL;
     207         256 :     side = "R";
     208             :   }
     209         256 :   Z = left? Y: X;
     210         256 :   if (ds->state>=DS_STATE_CONDENSED) {
     211             :     /* DSSolve() has been called, backtransform with matrix Q */
     212          55 :     back = "B";
     213          55 :     ierr = PetscMemcpy(Z,ds->mat[DS_MAT_Q],ld*ld*sizeof(PetscScalar));CHKERRQ(ierr);
     214         201 :   } else back = "A";
     215             : #if !defined(PETSC_USE_COMPLEX)
     216         256 :   ierr = DSAllocateWork_Private(ds,3*ld,0,0);CHKERRQ(ierr);
     217         256 :   PetscStackCallBLAS("LAPACKtrevc",LAPACKtrevc_(side,back,NULL,&n,A,&ld,Y,&ld,X,&ld,&n,&mout,ds->work,&info));
     218             : #else
     219             :   ierr = DSAllocateWork_Private(ds,2*ld,ld,0);CHKERRQ(ierr);
     220             :   PetscStackCallBLAS("LAPACKtrevc",LAPACKtrevc_(side,back,NULL,&n,A,&ld,Y,&ld,X,&ld,&n,&mout,ds->work,ds->rwork,&info));
     221             : #endif
     222         256 :   SlepcCheckLapackInfo("trevc",info);
     223             : 
     224             :   /* normalize eigenvectors */
     225        2055 :   for (i=0;i<n;i++) {
     226        1799 :     iscomplex = (i<n-1 && A[i+1+i*ld]!=0.0)? PETSC_TRUE: PETSC_FALSE;
     227        1799 :     norm = BLASnrm2_(&n,Z+i*ld,&inc);
     228             : #if !defined(PETSC_USE_COMPLEX)
     229        1799 :     if (iscomplex) {
     230         111 :       tmp = BLASnrm2_(&n,Z+(i+1)*ld,&inc);
     231         111 :       norm = SlepcAbsEigenvalue(norm,tmp);
     232             :     }
     233             : #endif
     234        1799 :     tmp = 1.0 / norm;
     235        1799 :     PetscStackCallBLAS("BLASscal",BLASscal_(&n,&tmp,Z+i*ld,&inc));
     236             : #if !defined(PETSC_USE_COMPLEX)
     237        1799 :     if (iscomplex) PetscStackCallBLAS("BLASscal",BLASscal_(&n,&tmp,Z+(i+1)*ld,&inc));
     238             : #endif
     239        1799 :     if (iscomplex) i++;
     240             :   }
     241         256 :   PetscFunctionReturn(0);
     242             : #endif
     243             : }
     244             : 
     245        6822 : PetscErrorCode DSVectors_NHEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
     246             : {
     247             :   PetscErrorCode ierr;
     248             : 
     249        6822 :   PetscFunctionBegin;
     250        6822 :   switch (mat) {
     251             :     case DS_MAT_X:
     252        5581 :       if (ds->refined) {
     253          78 :         if (!ds->extrarow) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Refined vectors require activating the extra row");
     254          78 :         if (j) {
     255          77 :           ierr = DSVectors_NHEP_Refined_Some(ds,j,rnorm,PETSC_FALSE);CHKERRQ(ierr);
     256             :         } else {
     257           1 :           ierr = DSVectors_NHEP_Refined_All(ds,PETSC_FALSE);CHKERRQ(ierr);
     258             :         }
     259             :       } else {
     260        5503 :         if (j) {
     261        5247 :           ierr = DSVectors_NHEP_Eigen_Some(ds,j,rnorm,PETSC_FALSE);CHKERRQ(ierr);
     262             :         } else {
     263         256 :           ierr = DSVectors_NHEP_Eigen_All(ds,PETSC_FALSE);CHKERRQ(ierr);
     264             :         }
     265             :       }
     266        5581 :       break;
     267             :     case DS_MAT_Y:
     268        1241 :       if (ds->refined) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
     269        1241 :       if (j) {
     270        1241 :         ierr = DSVectors_NHEP_Eigen_Some(ds,j,rnorm,PETSC_TRUE);CHKERRQ(ierr);
     271             :       } else {
     272           0 :         ierr = DSVectors_NHEP_Eigen_All(ds,PETSC_TRUE);CHKERRQ(ierr);
     273             :       }
     274        1241 :       break;
     275             :     case DS_MAT_U:
     276             :     case DS_MAT_VT:
     277           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
     278             :       break;
     279             :     default:
     280           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
     281             :   }
     282        6822 :   if (ds->state < DS_STATE_CONDENSED) {
     283         202 :     ierr = DSSetState(ds,DS_STATE_CONDENSED);CHKERRQ(ierr);
     284             :   }
     285        6822 :   PetscFunctionReturn(0);
     286             : }
     287             : 
     288          15 : static PetscErrorCode DSSort_NHEP_Arbitrary(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
     289             : {
     290             : #if defined(PETSC_MISSING_LAPACK_TRSEN)
     291             :   PetscFunctionBegin;
     292             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TRSEN - Lapack routine is unavailable");
     293             : #else
     294             :   PetscErrorCode ierr;
     295             :   PetscInt       i;
     296             :   PetscBLASInt   info,n,ld,mout,lwork,*selection;
     297          15 :   PetscScalar    *T = ds->mat[DS_MAT_A],*Q = ds->mat[DS_MAT_Q],*work;
     298             : #if !defined(PETSC_USE_COMPLEX)
     299             :   PetscBLASInt   *iwork,liwork;
     300             : #endif
     301             : 
     302          15 :   PetscFunctionBegin;
     303          15 :   if (!k) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Must supply argument k");
     304          15 :   ierr = PetscBLASIntCast(ds->n,&n);CHKERRQ(ierr);
     305          15 :   ierr = PetscBLASIntCast(ds->ld,&ld);CHKERRQ(ierr);
     306             : #if !defined(PETSC_USE_COMPLEX)
     307          15 :   lwork = n;
     308          15 :   liwork = 1;
     309          15 :   ierr = DSAllocateWork_Private(ds,lwork,0,liwork+n);CHKERRQ(ierr);
     310          15 :   work = ds->work;
     311          15 :   lwork = ds->lwork;
     312          15 :   selection = ds->iwork;
     313          15 :   iwork = ds->iwork + n;
     314          15 :   liwork = ds->liwork - n;
     315             : #else
     316             :   lwork = 1;
     317             :   ierr = DSAllocateWork_Private(ds,lwork,0,n);CHKERRQ(ierr);
     318             :   work = ds->work;
     319             :   selection = ds->iwork;
     320             : #endif
     321             :   /* Compute the selected eigenvalue to be in the leading position */
     322          15 :   ierr = DSSortEigenvalues_Private(ds,rr,ri,ds->perm,PETSC_FALSE);CHKERRQ(ierr);
     323          15 :   ierr = PetscMemzero(selection,n*sizeof(PetscBLASInt));CHKERRQ(ierr);
     324          15 :   for (i=0;i<*k;i++) selection[ds->perm[i]] = 1;
     325             : #if !defined(PETSC_USE_COMPLEX)
     326          15 :   PetscStackCallBLAS("LAPACKtrsen",LAPACKtrsen_("N","V",selection,&n,T,&ld,Q,&ld,wr,wi,&mout,NULL,NULL,work,&lwork,iwork,&liwork,&info));
     327             : #else
     328             :   PetscStackCallBLAS("LAPACKtrsen",LAPACKtrsen_("N","V",selection,&n,T,&ld,Q,&ld,wr,&mout,NULL,NULL,work,&lwork,&info));
     329             : #endif
     330          15 :   SlepcCheckLapackInfo("trsen",info);
     331          15 :   *k = mout;
     332          15 :   PetscFunctionReturn(0);
     333             : #endif
     334             : }
     335             : 
     336        5838 : static PetscErrorCode DSSort_NHEP_Total(DS ds,PetscScalar *wr,PetscScalar *wi)
     337             : {
     338             : #if defined(SLEPC_MISSING_LAPACK_TREXC)
     339             :   PetscFunctionBegin;
     340             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"TREXC - Lapack routine is unavailable");
     341             : #else
     342             :   PetscErrorCode ierr;
     343             :   PetscScalar    re;
     344             :   PetscInt       i,j,pos,result;
     345             :   PetscBLASInt   ifst,ilst,info,n,ld;
     346        5838 :   PetscScalar    *T = ds->mat[DS_MAT_A];
     347        5838 :   PetscScalar    *Q = ds->mat[DS_MAT_Q];
     348             : #if !defined(PETSC_USE_COMPLEX)
     349             :   PetscScalar    *work,im;
     350             : #endif
     351             : 
     352        5838 :   PetscFunctionBegin;
     353        5838 :   ierr = PetscBLASIntCast(ds->n,&n);CHKERRQ(ierr);
     354        5838 :   ierr = PetscBLASIntCast(ds->ld,&ld);CHKERRQ(ierr);
     355             : #if !defined(PETSC_USE_COMPLEX)
     356        5838 :   ierr = DSAllocateWork_Private(ds,ld,0,0);CHKERRQ(ierr);
     357        5838 :   work = ds->work;
     358             : #endif
     359             :   /* selection sort */
     360       72003 :   for (i=ds->l;i<n-1;i++) {
     361       66165 :     re = wr[i];
     362             : #if !defined(PETSC_USE_COMPLEX)
     363       66165 :     im = wi[i];
     364             : #endif
     365       66165 :     pos = 0;
     366       66165 :     j=i+1; /* j points to the next eigenvalue */
     367             : #if !defined(PETSC_USE_COMPLEX)
     368       66165 :     if (im != 0) j=i+2;
     369             : #endif
     370             :     /* find minimum eigenvalue */
     371      616050 :     for (;j<n;j++) {
     372             : #if !defined(PETSC_USE_COMPLEX)
     373      549885 :       ierr = SlepcSCCompare(ds->sc,re,im,wr[j],wi[j],&result);CHKERRQ(ierr);
     374             : #else
     375             :       ierr = SlepcSCCompare(ds->sc,re,0.0,wr[j],0.0,&result);CHKERRQ(ierr);
     376             : #endif
     377      549885 :       if (result > 0) {
     378      206573 :         re = wr[j];
     379             : #if !defined(PETSC_USE_COMPLEX)
     380      206573 :         im = wi[j];
     381             : #endif
     382      206573 :         pos = j;
     383             :       }
     384             : #if !defined(PETSC_USE_COMPLEX)
     385      549885 :       if (wi[j] != 0) j++;
     386             : #endif
     387             :     }
     388       66165 :     if (pos) {
     389             :       /* interchange blocks */
     390       38460 :       ierr = PetscBLASIntCast(pos+1,&ifst);CHKERRQ(ierr);
     391       38460 :       ierr = PetscBLASIntCast(i+1,&ilst);CHKERRQ(ierr);
     392             : #if !defined(PETSC_USE_COMPLEX)
     393       38460 :       PetscStackCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,work,&info));
     394             : #else
     395             :       PetscStackCallBLAS("LAPACKtrexc",LAPACKtrexc_("V",&n,T,&ld,Q,&ld,&ifst,&ilst,&info));
     396             : #endif
     397       38460 :       SlepcCheckLapackInfo("trexc",info);
     398             :       /* recover original eigenvalues from T matrix */
     399      446094 :       for (j=i;j<n;j++) {
     400      407634 :         wr[j] = T[j+j*ld];
     401             : #if !defined(PETSC_USE_COMPLEX)
     402      407634 :         if (j<n-1 && T[j+1+j*ld] != 0.0) {
     403             :           /* complex conjugate eigenvalue */
     404       62622 :           wi[j] = PetscSqrtReal(PetscAbsReal(T[j+1+j*ld])) *
     405       31311 :                   PetscSqrtReal(PetscAbsReal(T[j+(j+1)*ld]));
     406       31311 :           wr[j+1] = wr[j];
     407       31311 :           wi[j+1] = -wi[j];
     408       31311 :           j++;
     409             :         } else {
     410      376323 :           wi[j] = 0.0;
     411             :         }
     412             : #endif
     413             :       }
     414             :     }
     415             : #if !defined(PETSC_USE_COMPLEX)
     416       66165 :     if (wi[i] != 0) i++;
     417             : #endif
     418             :   }
     419        5838 :   PetscFunctionReturn(0);
     420             : #endif
     421             : }
     422             : 
     423        5853 : PetscErrorCode DSSort_NHEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
     424             : {
     425             :   PetscErrorCode ierr;
     426             : 
     427        5853 :   PetscFunctionBegin;
     428        5853 :   if (!rr || wr == rr) {
     429        5838 :     ierr = DSSort_NHEP_Total(ds,wr,wi);CHKERRQ(ierr);
     430             :   } else {
     431          15 :     ierr = DSSort_NHEP_Arbitrary(ds,wr,wi,rr,ri,k);CHKERRQ(ierr);
     432             :   }
     433        5853 :   PetscFunctionReturn(0);
     434             : }
     435             : 
     436        1728 : PetscErrorCode DSUpdateExtraRow_NHEP(DS ds)
     437             : {
     438             :   PetscErrorCode ierr;
     439             :   PetscInt       i;
     440        1728 :   PetscBLASInt   n,ld,incx=1;
     441        1728 :   PetscScalar    *A,*Q,*x,*y,one=1.0,zero=0.0;
     442             : 
     443        1728 :   PetscFunctionBegin;
     444        1728 :   ierr = PetscBLASIntCast(ds->n,&n);CHKERRQ(ierr);
     445        1728 :   ierr = PetscBLASIntCast(ds->ld,&ld);CHKERRQ(ierr);
     446        1728 :   A  = ds->mat[DS_MAT_A];
     447        1728 :   Q  = ds->mat[DS_MAT_Q];
     448        1728 :   ierr = DSAllocateWork_Private(ds,2*ld,0,0);CHKERRQ(ierr);
     449        1728 :   x = ds->work;
     450        1728 :   y = ds->work+ld;
     451        1728 :   for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
     452        1728 :   PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
     453        1728 :   for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
     454        1728 :   ds->k = n;
     455        1728 :   PetscFunctionReturn(0);
     456             : }
     457             : 
     458             : /*
     459             :    Reduce a matrix A to upper Hessenberg form Q'*A*Q, where Q is an orthogonal matrix.
     460             :    The result overwrites A. Matrix A has the form
     461             : 
     462             :          [ S | * ]
     463             :      A = [-------]
     464             :          [ R | H ]
     465             : 
     466             :   where S is an upper (quasi-)triangular matrix of order k, H is an upper Hessenberg
     467             :   matrix of order n-k, and R is all zeros except the first row (the arrow).
     468             :   The algorithm uses elementary reflectors to annihilate entries in the arrow, and
     469             :   then proceeds upwards. 
     470             :   If ilo>1, then it is assumed that the first ilo-1 entries of the arrow are zero, and
     471             :   hence the first ilo-1 rows and columns of Q are set to the identity matrix.
     472             : 
     473             :   Required workspace is 2*n.
     474             : */
     475           0 : static PetscErrorCode ArrowHessenberg(PetscBLASInt n,PetscBLASInt k,PetscBLASInt ilo,PetscScalar *A,PetscBLASInt lda,PetscScalar *Q,PetscBLASInt ldq,PetscScalar *work)
     476             : {
     477             : #if defined(SLEPC_MISSING_LAPACK_LARFG) || defined(SLEPC_MISSING_LAPACK_LARF)
     478             :   PetscFunctionBegin;
     479             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"LARFG/LARF - Lapack routines are unavailable");
     480             : #else
     481           0 :   PetscBLASInt i,j,n0,m,inc=1,incn=-1;
     482           0 :   PetscScalar  t,*v=work,*w=work+n,tau,tauc;
     483             : 
     484           0 :   PetscFunctionBegin;
     485           0 :   m = n-ilo+1;
     486           0 :   for (i=k;i>ilo;i--) {
     487           0 :     for (j=0;j<=i-ilo;j++) v[j] = A[i+(i-j-1)*lda];  /* _larfg does not allow negative inc, so use buffer */
     488           0 :     n0 = i-ilo+1;
     489           0 :     PetscStackCallBLAS("LAPACKlarfg",LAPACKlarfg_(&n0,v,v+1,&inc,&tau));
     490           0 :     for (j=1;j<=i-ilo;j++) v[j] = PetscConj(v[j]);
     491           0 :     tauc = PetscConj(tau);
     492           0 :     A[i+(i-1)*lda] = v[0];
     493           0 :     v[0] = 1.0;
     494           0 :     PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("R",&i,&n0,v,&incn,&tauc,A+(ilo-1)*lda,&lda,w));
     495           0 :     for (j=1;j<=i-ilo;j++) A[i+(i-j-1)*lda] = 0.0;
     496           0 :     PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n0,&m,v,&incn,&tau,A+ilo-1+(ilo-1)*lda,&lda,w));
     497           0 :     PetscStackCallBLAS("LAPACKlarf",LAPACKlarf_("L",&n0,&m,v,&incn,&tau,Q+ilo-1+(ilo-1)*ldq,&ldq,w));
     498             :   }
     499             :   /* trivial in-place transposition of Q */
     500           0 :   for (j=ilo-1;j<k;j++) {
     501           0 :     for (i=j;i<k;i++) {
     502           0 :       t = Q[i+j*ldq];
     503           0 :       if (i!=j) Q[i+j*ldq] = PetscConj(Q[j+i*ldq]);
     504           0 :       Q[j+i*ldq] = PetscConj(t);
     505             :     }
     506             :   }
     507           0 :   PetscFunctionReturn(0);
     508             : #endif
     509             : }
     510             : 
     511        4246 : PetscErrorCode DSSolve_NHEP(DS ds,PetscScalar *wr,PetscScalar *wi)
     512             : {
     513             : #if defined(SLEPC_MISSING_LAPACK_GEHRD) || defined(SLEPC_MISSING_LAPACK_ORGHR) || defined(PETSC_MISSING_LAPACK_HSEQR)
     514             :   PetscFunctionBegin;
     515             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEHRD/ORGHR/HSEQR - Lapack routines are unavailable");
     516             : #else
     517             :   PetscErrorCode ierr;
     518             :   PetscScalar    *work,*tau;
     519             :   PetscInt       i,j;
     520             :   PetscBLASInt   ilo,lwork,info,n,k,ld;
     521        4246 :   PetscScalar    *A = ds->mat[DS_MAT_A];
     522        4246 :   PetscScalar    *Q = ds->mat[DS_MAT_Q];
     523             : 
     524        4246 :   PetscFunctionBegin;
     525             : #if !defined(PETSC_USE_COMPLEX)
     526        4246 :   PetscValidPointer(wi,3);
     527             : #endif
     528        4246 :   ierr = PetscBLASIntCast(ds->n,&n);CHKERRQ(ierr);
     529        4246 :   ierr = PetscBLASIntCast(ds->ld,&ld);CHKERRQ(ierr);
     530        4246 :   ierr = PetscBLASIntCast(ds->l+1,&ilo);CHKERRQ(ierr);
     531        4246 :   ierr = PetscBLASIntCast(ds->k,&k);CHKERRQ(ierr);
     532        4246 :   ierr = DSAllocateWork_Private(ds,ld+6*ld,0,0);CHKERRQ(ierr);
     533        4246 :   tau  = ds->work;
     534        4246 :   work = ds->work+ld;
     535        4246 :   lwork = 6*ld;
     536             : 
     537             :   /* initialize orthogonal matrix */
     538        4246 :   ierr = PetscMemzero(Q,ld*ld*sizeof(PetscScalar));CHKERRQ(ierr);
     539        4246 :   for (i=0;i<n;i++) Q[i+i*ld] = 1.0;
     540        4246 :   if (n==1) { /* quick return */
     541          22 :     wr[0] = A[0];
     542          22 :     wi[0] = 0.0;
     543          22 :     PetscFunctionReturn(0);
     544             :   }
     545             : 
     546             :   /* reduce to upper Hessenberg form */
     547        4224 :   if (ds->state<DS_STATE_INTERMEDIATE) {
     548             :     if (PETSC_FALSE && k>0) {
     549             :       ierr = ArrowHessenberg(n,k,ilo,A,ld,Q,ld,work);CHKERRQ(ierr);
     550             :     } else {
     551        3106 :       PetscStackCallBLAS("LAPACKgehrd",LAPACKgehrd_(&n,&ilo,&n,A,&ld,tau,work,&lwork,&info));
     552        3106 :       SlepcCheckLapackInfo("gehrd",info);
     553       47162 :       for (j=0;j<n-1;j++) {
     554      477992 :         for (i=j+2;i<n;i++) {
     555      433936 :           Q[i+j*ld] = A[i+j*ld];
     556      433936 :           A[i+j*ld] = 0.0;
     557             :         }
     558             :       }
     559        3106 :       PetscStackCallBLAS("LAPACKorghr",LAPACKorghr_(&n,&ilo,&n,Q,&ld,tau,work,&lwork,&info));
     560        3106 :       SlepcCheckLapackInfo("orghr",info);
     561             :     }
     562             :   }
     563             : 
     564             :   /* compute the (real) Schur form */
     565             : #if !defined(PETSC_USE_COMPLEX)
     566        4224 :   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,wi,Q,&ld,work,&lwork,&info));
     567        8501 :   for (j=0;j<ds->l;j++) {
     568        4277 :     if (j==n-1 || A[j+1+j*ld] == 0.0) {
     569             :       /* real eigenvalue */
     570        4101 :       wr[j] = A[j+j*ld];
     571        4101 :       wi[j] = 0.0;
     572             :     } else {
     573             :       /* complex eigenvalue */
     574         176 :       wr[j] = A[j+j*ld];
     575         176 :       wr[j+1] = A[j+j*ld];
     576         352 :       wi[j] = PetscSqrtReal(PetscAbsReal(A[j+1+j*ld])) *
     577         176 :               PetscSqrtReal(PetscAbsReal(A[j+(j+1)*ld]));
     578         176 :       wi[j+1] = -wi[j];
     579         176 :       j++;
     580             :     }
     581             :   }
     582             : #else
     583             :   PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("S","V",&n,&ilo,&n,A,&ld,wr,Q,&ld,work,&lwork,&info));
     584             :   if (wi) for (i=ds->l;i<n;i++) wi[i] = 0.0;
     585             : #endif
     586        4224 :   SlepcCheckLapackInfo("hseqr",info);
     587        4224 :   PetscFunctionReturn(0);
     588             : #endif
     589             : }
     590             : 
     591           5 : PetscErrorCode DSSynchronize_NHEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
     592             : {
     593             :   PetscErrorCode ierr;
     594           5 :   PetscInt       ld=ds->ld,l=ds->l,k;
     595           5 :   PetscMPIInt    n,rank,off=0,size,ldn;
     596             : 
     597           5 :   PetscFunctionBegin;
     598           5 :   k = (ds->n-l)*ld;
     599           5 :   if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
     600           5 :   if (eigr) k += ds->n-l; 
     601           5 :   if (eigi) k += ds->n-l; 
     602           5 :   ierr = DSAllocateWork_Private(ds,k,0,0);CHKERRQ(ierr);
     603           5 :   ierr = PetscMPIIntCast(k*sizeof(PetscScalar),&size);CHKERRQ(ierr);
     604           5 :   ierr = PetscMPIIntCast(ds->n-l,&n);CHKERRQ(ierr);
     605           5 :   ierr = PetscMPIIntCast(ld*(ds->n-l),&ldn);CHKERRQ(ierr);
     606           5 :   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);CHKERRQ(ierr);
     607           5 :   if (!rank) {
     608           2 :     ierr = MPI_Pack(ds->mat[DS_MAT_A]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRQ(ierr);
     609           2 :     if (ds->state>DS_STATE_RAW) {
     610           2 :       ierr = MPI_Pack(ds->mat[DS_MAT_Q]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRQ(ierr);
     611             :     }
     612           2 :     if (eigr) {
     613           2 :       ierr = MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRQ(ierr);
     614             :     }
     615           2 :     if (eigi) {
     616           2 :       ierr = MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));CHKERRQ(ierr);
     617             :     }
     618             :   }
     619           5 :   ierr = MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));CHKERRQ(ierr);
     620           5 :   if (rank) {
     621           3 :     ierr = MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_A]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));CHKERRQ(ierr);
     622           3 :     if (ds->state>DS_STATE_RAW) {
     623           3 :       ierr = MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_Q]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));CHKERRQ(ierr);
     624             :     }
     625           3 :     if (eigr) {
     626           3 :       ierr = MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));CHKERRQ(ierr);
     627             :     }
     628           3 :     if (eigi) {
     629           3 :       ierr = MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));CHKERRQ(ierr);
     630             :     }
     631             :   }
     632           5 :   PetscFunctionReturn(0);
     633             : }
     634             : 
     635         688 : PetscErrorCode DSTruncate_NHEP(DS ds,PetscInt n)
     636             : {
     637         688 :   PetscInt    i,newn,ld=ds->ld,l=ds->l;
     638             :   PetscScalar *A;
     639             : 
     640         688 :   PetscFunctionBegin;
     641         688 :   if (ds->state==DS_STATE_CONDENSED) ds->t = ds->n;
     642         688 :   A = ds->mat[DS_MAT_A];
     643             :   /* be careful not to break a diagonal 2x2 block */
     644         688 :   if (A[n+(n-1)*ld]==0.0) newn = n;
     645             :   else {
     646         141 :     if (n<ds->n-1) newn = n+1;
     647           0 :     else newn = n-1;
     648             :   }
     649         688 :   if (ds->extrarow && ds->k==ds->n) {
     650             :     /* copy entries of extra row to the new position, then clean last row */
     651         688 :     for (i=l;i<newn;i++) A[newn+i*ld] = A[ds->n+i*ld];
     652         688 :     for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
     653             :   }
     654         688 :   ds->k = 0;
     655         688 :   ds->n = newn;
     656         688 :   PetscFunctionReturn(0);
     657             : }
     658             : 
     659           3 : PetscErrorCode DSCond_NHEP(DS ds,PetscReal *cond)
     660             : {
     661             : #if defined(PETSC_MISSING_LAPACK_GETRF) || defined(PETSC_MISSING_LAPACK_GETRI) || defined(SLEPC_MISSING_LAPACK_LANGE) || defined(SLEPC_MISSING_LAPACK_LANHS)
     662             :   PetscFunctionBegin;
     663             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF/GETRI/LANGE/LANHS - Lapack routines are unavailable");
     664             : #else
     665             :   PetscErrorCode ierr;
     666             :   PetscScalar    *work;
     667             :   PetscReal      *rwork;
     668             :   PetscBLASInt   *ipiv;
     669             :   PetscBLASInt   lwork,info,n,ld;
     670             :   PetscReal      hn,hin;
     671             :   PetscScalar    *A;
     672             : 
     673           3 :   PetscFunctionBegin;
     674           3 :   ierr = PetscBLASIntCast(ds->n,&n);CHKERRQ(ierr);
     675           3 :   ierr = PetscBLASIntCast(ds->ld,&ld);CHKERRQ(ierr);
     676           3 :   lwork = 8*ld;
     677           3 :   ierr = DSAllocateWork_Private(ds,lwork,ld,ld);CHKERRQ(ierr);
     678           3 :   work  = ds->work;
     679           3 :   rwork = ds->rwork;
     680           3 :   ipiv  = ds->iwork;
     681             : 
     682             :   /* use workspace matrix W to avoid overwriting A */
     683           3 :   ierr = DSAllocateMat_Private(ds,DS_MAT_W);CHKERRQ(ierr);
     684           3 :   A = ds->mat[DS_MAT_W];
     685           3 :   ierr = PetscMemcpy(A,ds->mat[DS_MAT_A],sizeof(PetscScalar)*ds->ld*ds->ld);CHKERRQ(ierr);
     686             : 
     687             :   /* norm of A */
     688           3 :   if (ds->state<DS_STATE_INTERMEDIATE) hn = LAPACKlange_("I",&n,&n,A,&ld,rwork);
     689           3 :   else hn = LAPACKlanhs_("I",&n,A,&ld,rwork);
     690             : 
     691             :   /* norm of inv(A) */
     692           3 :   PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n,&n,A,&ld,ipiv,&info));
     693           3 :   SlepcCheckLapackInfo("getrf",info);
     694           3 :   PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,A,&ld,ipiv,work,&lwork,&info));
     695           3 :   SlepcCheckLapackInfo("getri",info);
     696           3 :   hin = LAPACKlange_("I",&n,&n,A,&ld,rwork);
     697             : 
     698           3 :   *cond = hn*hin;
     699           3 :   PetscFunctionReturn(0);
     700             : #endif
     701             : }
     702             : 
     703          70 : PetscErrorCode DSTranslateHarmonic_NHEP(DS ds,PetscScalar tau,PetscReal beta,PetscBool recover,PetscScalar *gin,PetscReal *gamma)
     704             : {
     705             : #if defined(PETSC_MISSING_LAPACK_GETRF) || defined(PETSC_MISSING_LAPACK_GETRS)
     706             :   PetscFunctionBegin;
     707             :   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF/GETRS - Lapack routines are unavailable");
     708             : #else
     709             :   PetscErrorCode ierr;
     710             :   PetscInt       i,j;
     711          70 :   PetscBLASInt   *ipiv,info,n,ld,one=1,ncol;
     712          70 :   PetscScalar    *A,*B,*Q,*g=gin,*ghat;
     713          70 :   PetscScalar    done=1.0,dmone=-1.0,dzero=0.0;
     714             :   PetscReal      gnorm;
     715             : 
     716          70 :   PetscFunctionBegin;
     717          70 :   ierr = PetscBLASIntCast(ds->n,&n);CHKERRQ(ierr);
     718          70 :   ierr = PetscBLASIntCast(ds->ld,&ld);CHKERRQ(ierr);
     719          70 :   A  = ds->mat[DS_MAT_A];
     720             : 
     721          70 :   if (!recover) {
     722             : 
     723          51 :     ierr = DSAllocateWork_Private(ds,0,0,ld);CHKERRQ(ierr);
     724          51 :     ipiv = ds->iwork;
     725          51 :     if (!g) {
     726          28 :       ierr = DSAllocateWork_Private(ds,ld,0,0);CHKERRQ(ierr);
     727          28 :       g = ds->work;
     728             :     }
     729             :     /* use workspace matrix W to factor A-tau*eye(n) */
     730          51 :     ierr = DSAllocateMat_Private(ds,DS_MAT_W);CHKERRQ(ierr);
     731          51 :     B = ds->mat[DS_MAT_W];
     732          51 :     ierr = PetscMemcpy(B,A,sizeof(PetscScalar)*ld*ld);CHKERRQ(ierr);
     733             : 
     734             :     /* Vector g initialy stores b = beta*e_n^T */
     735          51 :     ierr = PetscMemzero(g,n*sizeof(PetscScalar));CHKERRQ(ierr);
     736          51 :     g[n-1] = beta;
     737             : 
     738             :     /* g = (A-tau*eye(n))'\b */
     739          51 :     for (i=0;i<n;i++) B[i+i*ld] -= tau;
     740          51 :     PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n,&n,B,&ld,ipiv,&info));
     741          51 :     SlepcCheckLapackInfo("getrf",info);
     742          51 :     ierr = PetscLogFlops(2.0*n*n*n/3.0);CHKERRQ(ierr);
     743          51 :     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("C",&n,&one,B,&ld,ipiv,g,&ld,&info));
     744          51 :     SlepcCheckLapackInfo("getrs",info);
     745          51 :     ierr = PetscLogFlops(2.0*n*n-n);CHKERRQ(ierr);
     746             : 
     747             :     /* A = A + g*b' */
     748          51 :     for (i=0;i<n;i++) A[i+(n-1)*ld] += g[i]*beta;
     749             : 
     750             :   } else { /* recover */
     751             : 
     752          19 :     PetscValidPointer(g,6);
     753          19 :     ierr = DSAllocateWork_Private(ds,ld,0,0);CHKERRQ(ierr);
     754          19 :     ghat = ds->work;
     755          19 :     Q    = ds->mat[DS_MAT_Q];
     756             : 
     757             :     /* g^ = -Q(:,idx)'*g */
     758          19 :     ierr = PetscBLASIntCast(ds->l+ds->k,&ncol);CHKERRQ(ierr);
     759          19 :     PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&n,&ncol,&dmone,Q,&ld,g,&one,&dzero,ghat,&one));
     760             : 
     761             :     /* A = A + g^*b' */
     762         188 :     for (i=0;i<ds->l+ds->k;i++)
     763        1682 :       for (j=ds->l;j<ds->l+ds->k;j++)
     764        1513 :         A[i+j*ld] += ghat[i]*Q[n-1+j*ld]*beta;
     765             : 
     766             :     /* g~ = (I-Q(:,idx)*Q(:,idx)')*g = g+Q(:,idx)*g^ */
     767          19 :     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&ncol,&done,Q,&ld,ghat,&one,&done,g,&one));
     768             :   }
     769             : 
     770             :   /* Compute gamma factor */
     771          70 :   if (gamma) {
     772          70 :     gnorm = 0.0;
     773          70 :     for (i=0;i<n;i++) gnorm = gnorm + PetscRealPart(g[i]*PetscConj(g[i]));
     774          70 :     *gamma = PetscSqrtReal(1.0+gnorm);
     775             :   }
     776          70 :   PetscFunctionReturn(0);
     777             : #endif
     778             : }
     779             : 
     780         556 : SLEPC_EXTERN PetscErrorCode DSCreate_NHEP(DS ds)
     781             : {
     782         556 :   PetscFunctionBegin;
     783         556 :   ds->ops->allocate      = DSAllocate_NHEP;
     784         556 :   ds->ops->view          = DSView_NHEP;
     785         556 :   ds->ops->vectors       = DSVectors_NHEP;
     786         556 :   ds->ops->solve[0]      = DSSolve_NHEP;
     787         556 :   ds->ops->sort          = DSSort_NHEP;
     788         556 :   ds->ops->synchronize   = DSSynchronize_NHEP;
     789         556 :   ds->ops->truncate      = DSTruncate_NHEP;
     790         556 :   ds->ops->update        = DSUpdateExtraRow_NHEP;
     791         556 :   ds->ops->cond          = DSCond_NHEP;
     792         556 :   ds->ops->transharm     = DSTranslateHarmonic_NHEP;
     793         556 :   PetscFunctionReturn(0);
     794             : }
     795             : 

Generated by: LCOV version 1.13