LCOV - code coverage report
Current view: top level - eps/impls/davidson - dvdcalcpairs.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 312 349 89.4 %
Date: 2019-07-19 01:26:51 Functions: 13 13 100.0 %

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-2019, 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             :    SLEPc eigensolver: "davidson"
      12             : 
      13             :    Step: calculate the best eigenpairs in the subspace V
      14             : 
      15             :    For that, performs these steps:
      16             :      1) Update W <- A * V
      17             :      2) Update H <- V' * W
      18             :      3) Obtain eigenpairs of H
      19             :      4) Select some eigenpairs
      20             :      5) Compute the Ritz pairs of the selected ones
      21             : */
      22             : 
      23             : #include "davidson.h"
      24             : #include <slepcblaslapack.h>
      25             : 
      26          90 : static PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d)
      27             : {
      28             :   PetscErrorCode ierr;
      29             : 
      30          90 :   PetscFunctionBegin;
      31          90 :   ierr = BVSetActiveColumns(d->eps->V,0,0);CHKERRQ(ierr);
      32          90 :   if (d->W) { ierr = BVSetActiveColumns(d->W,0,0);CHKERRQ(ierr); }
      33          90 :   ierr = BVSetActiveColumns(d->AX,0,0);CHKERRQ(ierr);
      34          90 :   if (d->BX) { ierr = BVSetActiveColumns(d->BX,0,0);CHKERRQ(ierr); }
      35          90 :   PetscFunctionReturn(0);
      36             : }
      37             : 
      38          90 : static PetscErrorCode dvd_calcpairs_qz_d(dvdDashboard *d)
      39             : {
      40             :   PetscErrorCode  ierr;
      41             : 
      42          90 :   PetscFunctionBegin;
      43          90 :   ierr = BVDestroy(&d->W);CHKERRQ(ierr);
      44          90 :   ierr = BVDestroy(&d->AX);CHKERRQ(ierr);
      45          90 :   ierr = BVDestroy(&d->BX);CHKERRQ(ierr);
      46          90 :   ierr = BVDestroy(&d->auxBV);CHKERRQ(ierr);
      47          90 :   ierr = MatDestroy(&d->H);CHKERRQ(ierr);
      48          90 :   ierr = MatDestroy(&d->G);CHKERRQ(ierr);
      49          90 :   ierr = MatDestroy(&d->auxM);CHKERRQ(ierr);
      50          90 :   ierr = SlepcVecPoolDestroy(&d->auxV);CHKERRQ(ierr);
      51          90 :   ierr = PetscFree(d->nBds);CHKERRQ(ierr);
      52          90 :   PetscFunctionReturn(0);
      53             : }
      54             : 
      55             : /* in complex, d->size_H real auxiliary values are needed */
      56        9949 : static PetscErrorCode dvd_calcpairs_projeig_solve(dvdDashboard *d)
      57             : {
      58             :   PetscErrorCode    ierr;
      59             :   Vec               v;
      60             :   PetscScalar       *pA;
      61             :   const PetscScalar *pv;
      62             :   PetscInt          i,lV,kV,n,ld;
      63             : 
      64        9949 :   PetscFunctionBegin;
      65        9949 :   ierr = BVGetActiveColumns(d->eps->V,&lV,&kV);CHKERRQ(ierr);
      66        9949 :   n = kV-lV;
      67        9949 :   ierr = DSSetDimensions(d->eps->ds,n,0,0,0);CHKERRQ(ierr);
      68        9949 :   ierr = DSCopyMat(d->eps->ds,DS_MAT_A,0,0,d->H,lV,lV,n,n,PETSC_FALSE);CHKERRQ(ierr);
      69        9949 :   if (d->G) {
      70        1316 :     ierr = DSCopyMat(d->eps->ds,DS_MAT_B,0,0,d->G,lV,lV,n,n,PETSC_FALSE);CHKERRQ(ierr);
      71             :   }
      72             :   /* Set the signature on projected matrix B */
      73        9949 :   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
      74           0 :     ierr = DSGetLeadingDimension(d->eps->ds,&ld);CHKERRQ(ierr);
      75           0 :     ierr = DSGetArray(d->eps->ds,DS_MAT_B,&pA);CHKERRQ(ierr);
      76           0 :     ierr = PetscArrayzero(pA,n*ld);CHKERRQ(ierr);
      77           0 :     ierr = VecCreateSeq(PETSC_COMM_SELF,kV,&v);CHKERRQ(ierr);
      78           0 :     ierr = BVGetSignature(d->eps->V,v);CHKERRQ(ierr);
      79           0 :     ierr = VecGetArrayRead(v,&pv);CHKERRQ(ierr);
      80           0 :     for (i=0;i<n;i++) {
      81           0 :       pA[i+ld*i] = d->nBds[i] = PetscRealPart(pv[lV+i]);
      82             :     }
      83           0 :     ierr = VecRestoreArrayRead(v,&pv);CHKERRQ(ierr);
      84           0 :     ierr = VecDestroy(&v);CHKERRQ(ierr);
      85           0 :     ierr = DSRestoreArray(d->eps->ds,DS_MAT_B,&pA);CHKERRQ(ierr);
      86             :   }
      87        9949 :   ierr = DSSetState(d->eps->ds,DS_STATE_RAW);CHKERRQ(ierr);
      88        9949 :   ierr = DSSolve(d->eps->ds,d->eigr,d->eigi);CHKERRQ(ierr);
      89        9949 :   PetscFunctionReturn(0);
      90             : }
      91             : 
      92             : /*
      93             :    A(lA:kA-1,lA:kA-1) <- Z(l:k-1)'*A(l:k-1,l:k-1)*Q(l,k-1), where k=l+kA-lA
      94             :  */
      95        1946 : static PetscErrorCode EPSXDUpdateProj(Mat Q,Mat Z,PetscInt l,Mat A,PetscInt lA,PetscInt kA,Mat aux)
      96             : {
      97             :   PetscErrorCode ierr;
      98        1946 :   PetscScalar    one=1.0,zero=0.0;
      99        1946 :   PetscInt       i,j,dA_=kA-lA,m0,n0,ldA_,ldQ_,ldZ_,nQ_;
     100             :   PetscBLASInt   dA,nQ,ldA,ldQ,ldZ;
     101             :   PetscScalar    *pA,*pQ,*pZ,*pW;
     102        1946 :   PetscBool      symm=PETSC_FALSE,set,flg;
     103             : 
     104        1946 :   PetscFunctionBegin;
     105        1946 :   ierr = MatGetSize(A,&m0,&n0);CHKERRQ(ierr); ldA_=m0;
     106        1946 :   if (m0!=n0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"A should be square");
     107        1946 :   if (lA<0 || lA>m0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial row, column in A");
     108        1946 :   if (kA<0 || kA<lA || kA>m0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid final row, column in A");
     109        1946 :   ierr = MatIsHermitianKnown(A,&set,&flg);CHKERRQ(ierr);
     110        1946 :   symm = set? flg: PETSC_FALSE;
     111        1946 :   ierr = MatGetSize(Q,&m0,&n0);CHKERRQ(ierr); ldQ_=nQ_=m0;
     112        1946 :   if (l<0 || l>n0 || l+dA_>n0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Q");
     113        1946 :   ierr = MatGetSize(Z,&m0,&n0);CHKERRQ(ierr); ldZ_=m0;
     114        1946 :   if (l<0 || l>n0 || l+dA_>n0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Z");
     115        1946 :   ierr = MatGetSize(aux,&m0,&n0);CHKERRQ(ierr);
     116        1946 :   if (m0*n0<nQ_*dA_) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"aux should be larger");
     117        1946 :   ierr = PetscBLASIntCast(dA_,&dA);CHKERRQ(ierr);
     118        1946 :   ierr = PetscBLASIntCast(nQ_,&nQ);CHKERRQ(ierr);
     119        1946 :   ierr = PetscBLASIntCast(ldA_,&ldA);CHKERRQ(ierr);
     120        1946 :   ierr = PetscBLASIntCast(ldQ_,&ldQ);CHKERRQ(ierr);
     121        1946 :   ierr = PetscBLASIntCast(ldZ_,&ldZ);CHKERRQ(ierr);
     122        1946 :   ierr = MatDenseGetArray(A,&pA);CHKERRQ(ierr);
     123        1946 :   ierr = MatDenseGetArray(Q,&pQ);CHKERRQ(ierr);
     124        1946 :   if (Q!=Z) { ierr = MatDenseGetArray(Z,&pZ);CHKERRQ(ierr); }
     125        1476 :   else pZ = pQ;
     126             : #if PETSC_USE_DEBUG
     127             :   /* Avoid valgrind warning in xgemm and xsymm */
     128        1946 :   ierr = MatZeroEntries(aux);CHKERRQ(ierr);
     129             : #endif
     130        1946 :   ierr = MatDenseGetArray(aux,&pW);CHKERRQ(ierr);
     131             :   /* W = A*Q */
     132        1946 :   if (symm) {
     133             :     /* symmetrize before multiplying */
     134       17262 :     for (i=lA+1;i<lA+nQ;i++) {
     135      117670 :       for (j=lA;j<i;j++) pA[i+j*ldA] = PetscConj(pA[j+i*ldA]);
     136             :     }
     137             :   }
     138        1946 :   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&nQ,&dA,&nQ,&one,&pA[ldA*lA+lA],&ldA,&pQ[ldQ*l+l],&ldQ,&zero,pW,&nQ));
     139             :   /* A = Q'*W */
     140        1946 :   PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&dA,&dA,&nQ,&one,&pZ[ldZ*l+l],&ldZ,pW,&nQ,&zero,&pA[ldA*lA+lA],&ldA));
     141        1946 :   ierr = MatDenseRestoreArray(A,&pA);CHKERRQ(ierr);
     142        1946 :   ierr = MatDenseRestoreArray(Q,&pQ);CHKERRQ(ierr);
     143        1946 :   if (Q!=Z) { ierr = MatDenseGetArray(Z,&pZ);CHKERRQ(ierr); }
     144        1476 :   else pZ = pQ;
     145        1946 :   ierr = MatDenseRestoreArray(aux,&pW);CHKERRQ(ierr);
     146        1946 :   PetscFunctionReturn(0);
     147             : }
     148             : 
     149        1711 : static PetscErrorCode dvd_calcpairs_updateproj(dvdDashboard *d)
     150             : {
     151             :   PetscErrorCode ierr;
     152             :   Mat            Q,Z;
     153             :   PetscInt       lV,kV;
     154             :   PetscBool      symm;
     155             : 
     156        1711 :   PetscFunctionBegin;
     157        1711 :   ierr = DSGetMat(d->eps->ds,DS_MAT_Q,&Q);CHKERRQ(ierr);
     158        1711 :   if (d->W) { ierr = DSGetMat(d->eps->ds,DS_MAT_Z,&Z);CHKERRQ(ierr); }
     159        1476 :   else Z = Q;
     160        1711 :   ierr = BVGetActiveColumns(d->eps->V,&lV,&kV);CHKERRQ(ierr);
     161        1711 :   ierr = EPSXDUpdateProj(Q,Z,0,d->H,lV,lV+d->V_tra_e,d->auxM);CHKERRQ(ierr);
     162        1711 :   if (d->G) { ierr = EPSXDUpdateProj(Q,Z,0,d->G,lV,lV+d->V_tra_e,d->auxM);CHKERRQ(ierr); }
     163        1711 :   ierr = MatDestroy(&Q);CHKERRQ(ierr);
     164        1711 :   if (d->W) { ierr = MatDestroy(&Z);CHKERRQ(ierr); }
     165             : 
     166        1711 :   ierr = PetscObjectTypeCompareAny((PetscObject)d->eps->ds,&symm,DSHEP,DSGHIEP,DSGHEP,"");CHKERRQ(ierr);
     167        1711 :   if (d->V_tra_s==0 || symm) PetscFunctionReturn(0);
     168             :   /* Compute upper part of H (and G): H(0:l-1,l:k-1) <- W(0:l-1)' * AV(l:k-1), where
     169             :      k=l+d->V_tra_s */
     170          69 :   ierr = BVSetActiveColumns(d->W?d->W:d->eps->V,0,lV);CHKERRQ(ierr);
     171          69 :   ierr = BVSetActiveColumns(d->AX,lV,lV+d->V_tra_s);CHKERRQ(ierr);
     172          69 :   ierr = BVDot(d->AX,d->W?d->W:d->eps->V,d->H);CHKERRQ(ierr);
     173          69 :   if (d->G) {
     174          57 :     ierr = BVSetActiveColumns(d->BX?d->BX:d->eps->V,lV,lV+d->V_tra_s);CHKERRQ(ierr);
     175          57 :     ierr = BVDot(d->BX?d->BX:d->eps->V,d->W?d->W:d->eps->V,d->G);CHKERRQ(ierr);
     176             :   }
     177          69 :   ierr = PetscObjectTypeCompareAny((PetscObject)d->eps->ds,&symm,DSGHEP,"");CHKERRQ(ierr);
     178          69 :   if (!symm) {
     179          69 :     ierr = BVSetActiveColumns(d->W?d->W:d->eps->V,lV,lV+d->V_tra_s);CHKERRQ(ierr);
     180          69 :     ierr = BVSetActiveColumns(d->AX,0,lV);CHKERRQ(ierr);
     181          69 :     ierr = BVDot(d->AX,d->W?d->W:d->eps->V,d->H);CHKERRQ(ierr);
     182          69 :     if (d->G) {
     183          57 :       ierr = BVSetActiveColumns(d->BX?d->BX:d->eps->V,0,lV);CHKERRQ(ierr);
     184          57 :       ierr = BVDot(d->BX?d->BX:d->eps->V,d->W?d->W:d->eps->V,d->G);CHKERRQ(ierr);
     185             :     }
     186             :   }
     187          69 :   ierr = BVSetActiveColumns(d->eps->V,lV,kV);CHKERRQ(ierr);
     188          69 :   ierr = BVSetActiveColumns(d->AX,lV,kV);CHKERRQ(ierr);
     189          69 :   if (d->BX) { ierr = BVSetActiveColumns(d->BX,lV,kV);CHKERRQ(ierr); }
     190          69 :   if (d->W) { ierr = BVSetActiveColumns(d->W,lV,kV);CHKERRQ(ierr); }
     191          69 :   if (d->W) { ierr = dvd_harm_updateproj(d);CHKERRQ(ierr); }
     192          69 :   PetscFunctionReturn(0);
     193             : }
     194             : 
     195             : /*
     196             :    BV <- BV*MT
     197             :  */
     198        3957 : PETSC_STATIC_INLINE PetscErrorCode dvd_calcpairs_updateBV0_gen(dvdDashboard *d,BV bv,DSMatType mat)
     199             : {
     200             :   PetscErrorCode ierr;
     201             :   PetscInt       l,k,n;
     202             :   Mat            auxM;
     203             : 
     204        3957 :   PetscFunctionBegin;
     205        3957 :   ierr = BVGetActiveColumns(d->eps->V,&l,&k);CHKERRQ(ierr);
     206        3957 :   ierr = MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&auxM);CHKERRQ(ierr);
     207        3957 :   ierr = MatZeroEntries(auxM);CHKERRQ(ierr);
     208        3957 :   ierr = DSGetDimensions(d->eps->ds,&n,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
     209        3957 :   if (k-l!=n) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
     210        3957 :   ierr = DSCopyMat(d->eps->ds,mat,0,0,auxM,l,l,n,d->V_tra_e,PETSC_TRUE);CHKERRQ(ierr);
     211        3957 :   ierr = BVMultInPlace(bv,auxM,l,l+d->V_tra_e);CHKERRQ(ierr);
     212        3957 :   ierr = MatDestroy(&auxM);CHKERRQ(ierr);
     213        3957 :   PetscFunctionReturn(0);
     214             : }
     215             : 
     216        9949 : static PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d)
     217             : {
     218             :   PetscErrorCode ierr;
     219             :   PetscInt       i,l,k;
     220             :   Vec            v1,v2;
     221             :   PetscScalar    *pv;
     222             : 
     223        9949 :   PetscFunctionBegin;
     224        9949 :   ierr = BVGetActiveColumns(d->eps->V,&l,&k);CHKERRQ(ierr);
     225             :   /* Update AV, BV, W and the projected matrices */
     226             :   /* 1. S <- S*MT */
     227        9949 :   if (d->V_tra_s != d->V_tra_e || d->V_tra_e > 0) {
     228        1711 :     ierr = dvd_calcpairs_updateBV0_gen(d,d->eps->V,DS_MAT_Q);CHKERRQ(ierr);
     229        1711 :     if (d->W) { ierr = dvd_calcpairs_updateBV0_gen(d,d->W,DS_MAT_Z);CHKERRQ(ierr); }
     230        1711 :     ierr = dvd_calcpairs_updateBV0_gen(d,d->AX,DS_MAT_Q);CHKERRQ(ierr);
     231        1711 :     if (d->BX) { ierr = dvd_calcpairs_updateBV0_gen(d,d->BX,DS_MAT_Q);CHKERRQ(ierr); }
     232        1711 :     ierr = dvd_calcpairs_updateproj(d);CHKERRQ(ierr);
     233             :     /* Update signature */
     234        1711 :     if (d->nBds) {
     235           0 :       ierr = VecCreateSeq(PETSC_COMM_SELF,l+d->V_tra_e,&v1);CHKERRQ(ierr);
     236           0 :       ierr = BVSetActiveColumns(d->eps->V,0,l+d->V_tra_e);CHKERRQ(ierr);
     237           0 :       ierr = BVGetSignature(d->eps->V,v1);CHKERRQ(ierr);
     238           0 :       ierr = VecGetArray(v1,&pv);CHKERRQ(ierr);
     239           0 :       for (i=0;i<d->V_tra_e;i++) pv[l+i] = d->nBds[i];
     240           0 :       ierr = VecRestoreArray(v1,&pv);CHKERRQ(ierr);
     241           0 :       ierr = BVSetSignature(d->eps->V,v1);CHKERRQ(ierr);
     242           0 :       ierr = BVSetActiveColumns(d->eps->V,l,k);CHKERRQ(ierr);
     243           0 :       ierr = VecDestroy(&v1);CHKERRQ(ierr);
     244             :     }
     245        1711 :     k = l+d->V_tra_e;
     246        1711 :     l+= d->V_tra_s;
     247             :   } else {
     248             :     /* 2. V <- orth(V, V_new) */
     249        8238 :     ierr = dvd_orthV(d->eps->V,l+d->V_new_s,l+d->V_new_e);CHKERRQ(ierr);
     250             :     /* 3. AV <- [AV A * V(V_new_s:V_new_e-1)] */
     251             :     /* Check consistency */
     252        8238 :     if (k-l != d->V_new_s) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
     253       19826 :     for (i=l+d->V_new_s;i<l+d->V_new_e;i++) {
     254       11588 :       ierr = BVGetColumn(d->eps->V,i,&v1);CHKERRQ(ierr);
     255       11588 :       ierr = BVGetColumn(d->AX,i,&v2);CHKERRQ(ierr);
     256       11588 :       ierr = MatMult(d->A,v1,v2);CHKERRQ(ierr);
     257       11588 :       ierr = BVRestoreColumn(d->eps->V,i,&v1);CHKERRQ(ierr);
     258       11588 :       ierr = BVRestoreColumn(d->AX,i,&v2);CHKERRQ(ierr);
     259             :     }
     260             :     /* 4. BV <- [BV B * V(V_new_s:V_new_e-1)] */
     261        8238 :     if (d->BX) {
     262             :       /* Check consistency */
     263        1514 :       if (k-l != d->V_new_s) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
     264        3315 :       for (i=l+d->V_new_s;i<l+d->V_new_e;i++) {
     265        1801 :         ierr = BVGetColumn(d->eps->V,i,&v1);CHKERRQ(ierr);
     266        1801 :         ierr = BVGetColumn(d->BX,i,&v2);CHKERRQ(ierr);
     267        1801 :         ierr = MatMult(d->B,v1,v2);CHKERRQ(ierr);
     268        1801 :         ierr = BVRestoreColumn(d->eps->V,i,&v1);CHKERRQ(ierr);
     269        1801 :         ierr = BVRestoreColumn(d->BX,i,&v2);CHKERRQ(ierr);
     270             :       }
     271             :     }
     272             :     /* 5. W <- [W f(AV,BV)] */
     273        8238 :     if (d->W) {
     274        1081 :       ierr = d->calcpairs_W(d);CHKERRQ(ierr);
     275        1081 :       ierr = dvd_orthV(d->W,l+d->V_new_s,l+d->V_new_e);CHKERRQ(ierr);
     276             :     }
     277             :     /* 6. H <- W' * AX; G <- W' * BX */
     278        8238 :     ierr = BVSetActiveColumns(d->eps->V,l+d->V_new_s,l+d->V_new_e);CHKERRQ(ierr);
     279        8238 :     ierr = BVSetActiveColumns(d->AX,l+d->V_new_s,l+d->V_new_e);CHKERRQ(ierr);
     280        8238 :     if (d->BX) { ierr = BVSetActiveColumns(d->BX,l+d->V_new_s,l+d->V_new_e);CHKERRQ(ierr); }
     281        8238 :     if (d->W) { ierr = BVSetActiveColumns(d->W,l+d->V_new_s,l+d->V_new_e);CHKERRQ(ierr); }
     282        8238 :     ierr = BVMatProject(d->AX,NULL,d->W?d->W:d->eps->V,d->H);CHKERRQ(ierr);
     283        8238 :     if (d->G) { ierr = BVMatProject(d->BX?d->BX:d->eps->V,NULL,d->W?d->W:d->eps->V,d->G);CHKERRQ(ierr); }
     284        8238 :     ierr = BVSetActiveColumns(d->eps->V,l,k);CHKERRQ(ierr);
     285        8238 :     ierr = BVSetActiveColumns(d->AX,l,k);CHKERRQ(ierr);
     286        8238 :     if (d->BX) { ierr = BVSetActiveColumns(d->BX,l,k);CHKERRQ(ierr); }
     287        8238 :     if (d->W) { ierr = BVSetActiveColumns(d->W,l,k);CHKERRQ(ierr); }
     288             : 
     289             :     /* Perform the transformation on the projected problem */
     290        8238 :     if (d->W) {
     291        1081 :       ierr = d->calcpairs_proj_trans(d);CHKERRQ(ierr);
     292             :     }
     293        8238 :     k = l+d->V_new_e;
     294             :   }
     295        9949 :   ierr = BVSetActiveColumns(d->eps->V,l,k);CHKERRQ(ierr);
     296        9949 :   ierr = BVSetActiveColumns(d->AX,l,k);CHKERRQ(ierr);
     297        9949 :   if (d->BX) { ierr = BVSetActiveColumns(d->BX,l,k);CHKERRQ(ierr); }
     298        9949 :   if (d->W) { ierr = BVSetActiveColumns(d->W,l,k);CHKERRQ(ierr); }
     299             : 
     300             :   /* Solve the projected problem */
     301        9949 :   ierr = dvd_calcpairs_projeig_solve(d);CHKERRQ(ierr);
     302             : 
     303        9949 :   d->V_tra_s = d->V_tra_e = 0;
     304        9949 :   d->V_new_s = d->V_new_e;
     305        9949 :   PetscFunctionReturn(0);
     306             : }
     307             : 
     308       11322 : static PetscErrorCode dvd_calcpairs_apply_arbitrary(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscScalar *rr,PetscScalar *ri)
     309             : {
     310             :   PetscInt       i,k,ld;
     311             :   PetscScalar    *pX;
     312             :   Vec            *X,xr,xi;
     313             :   PetscErrorCode ierr;
     314             : #if defined(PETSC_USE_COMPLEX)
     315             :   PetscInt       N=1;
     316             : #else
     317       11322 :   PetscInt       N=2,j;
     318             : #endif
     319             : 
     320       11322 :   PetscFunctionBegin;
     321             :   /* Quick exit without neither arbitrary selection nor harmonic extraction */
     322       11322 :   if (!d->eps->arbitrary && !d->calcpairs_eig_backtrans) PetscFunctionReturn(0);
     323             : 
     324             :   /* Quick exit without arbitrary selection, but with harmonic extraction */
     325        2926 :   if (d->calcpairs_eig_backtrans) {
     326       28110 :     for (i=r_s; i<r_e; i++) {
     327       28110 :       ierr = d->calcpairs_eig_backtrans(d,d->eigr[i],d->eigi[i],&rr[i-r_s],&ri[i-r_s]);CHKERRQ(ierr);
     328             :     }
     329             :   }
     330        2926 :   if (!d->eps->arbitrary) PetscFunctionReturn(0);
     331             : 
     332        1168 :   ierr = SlepcVecPoolGetVecs(d->auxV,N,&X);CHKERRQ(ierr);
     333        1168 :   ierr = DSGetLeadingDimension(d->eps->ds,&ld);CHKERRQ(ierr);
     334       13644 :   for (i=r_s;i<r_e;i++) {
     335       13644 :     k = i;
     336       13644 :     ierr = DSVectors(d->eps->ds,DS_MAT_X,&k,NULL);CHKERRQ(ierr);
     337       13644 :     ierr = DSGetArray(d->eps->ds,DS_MAT_X,&pX);CHKERRQ(ierr);
     338       13644 :     ierr = dvd_improvex_compute_X(d,i,k+1,X,pX,ld);CHKERRQ(ierr);
     339       13644 :     ierr = DSRestoreArray(d->eps->ds,DS_MAT_X,&pX);CHKERRQ(ierr);
     340             : #if !defined(PETSC_USE_COMPLEX)
     341       13644 :     if (d->nX[i] != 1.0) {
     342           0 :       for (j=i;j<k+1;j++) {
     343           0 :         ierr = VecScale(X[j-i],1.0/d->nX[i]);CHKERRQ(ierr);
     344             :       }
     345             :     }
     346       13644 :     xr = X[0];
     347       13644 :     xi = X[1];
     348       13644 :     if (i == k) {
     349       13644 :       ierr = VecSet(xi,0.0);CHKERRQ(ierr);
     350             :     }
     351             : #else
     352             :     xr = X[0];
     353             :     xi = NULL;
     354             :     if (d->nX[i] != 1.0) {
     355             :       ierr = VecScale(xr,1.0/d->nX[i]);CHKERRQ(ierr);
     356             :     }
     357             : #endif
     358       13644 :     ierr = (d->eps->arbitrary)(rr[i-r_s],ri[i-r_s],xr,xi,&rr[i-r_s],&ri[i-r_s],d->eps->arbitraryctx);CHKERRQ(ierr);
     359             : #if !defined(PETSC_USE_COMPLEX)
     360       13644 :     if (i != k) {
     361           0 :       rr[i+1-r_s] = rr[i-r_s];
     362           0 :       ri[i+1-r_s] = ri[i-r_s];
     363           0 :       i++;
     364             :     }
     365             : #endif
     366             :   }
     367        1168 :   ierr = SlepcVecPoolRestoreVecs(d->auxV,N,&X);CHKERRQ(ierr);
     368        1168 :   PetscFunctionReturn(0);
     369             : }
     370             : 
     371        9859 : static PetscErrorCode dvd_calcpairs_selectPairs(dvdDashboard *d,PetscInt n)
     372             : {
     373             :   PetscInt       k,lV,kV,nV;
     374             :   PetscScalar    *rr,*ri;
     375             :   PetscErrorCode ierr;
     376             : 
     377        9859 :   PetscFunctionBegin;
     378        9859 :   ierr = BVGetActiveColumns(d->eps->V,&lV,&kV);CHKERRQ(ierr);
     379        9859 :   nV = kV - lV;
     380        9859 :   n = PetscMin(n,nV);
     381        9859 :   if (n <= 0) PetscFunctionReturn(0);
     382             :   /* Put the best n pairs at the beginning. Useful for restarting */
     383        9859 :   if (d->eps->arbitrary || d->calcpairs_eig_backtrans) {
     384        1463 :     ierr = PetscMalloc1(nV,&rr);CHKERRQ(ierr);
     385        1463 :     ierr = PetscMalloc1(nV,&ri);CHKERRQ(ierr);
     386        1463 :     ierr = dvd_calcpairs_apply_arbitrary(d,0,nV,rr,ri);CHKERRQ(ierr);
     387             :   } else {
     388        8396 :     rr = d->eigr;
     389        8396 :     ri = d->eigi;
     390             :   }
     391        9859 :   k = n;
     392        9859 :   ierr = DSSort(d->eps->ds,d->eigr,d->eigi,rr,ri,&k);CHKERRQ(ierr);
     393             :   /* Put the best pair at the beginning. Useful to check its residual */
     394             : #if !defined(PETSC_USE_COMPLEX)
     395        9859 :   if (n != 1 && (n != 2 || d->eigi[0] == 0.0))
     396             : #else
     397             :   if (n != 1)
     398             : #endif
     399             :   {
     400        9859 :     ierr = dvd_calcpairs_apply_arbitrary(d,0,nV,rr,ri);CHKERRQ(ierr);
     401        9859 :     k = 1;
     402        9859 :     ierr = DSSort(d->eps->ds,d->eigr,d->eigi,rr,ri,&k);CHKERRQ(ierr);
     403             :   }
     404        9859 :   ierr = DSSynchronize(d->eps->ds,d->eigr,d->eigi);CHKERRQ(ierr);
     405             : 
     406        9859 :   if (d->calcpairs_eigs_trans) {
     407        1295 :     ierr = d->calcpairs_eigs_trans(d);CHKERRQ(ierr);
     408             :   }
     409        9859 :   if (d->eps->arbitrary || d->calcpairs_eig_backtrans) {
     410        1463 :     ierr = PetscFree(rr);CHKERRQ(ierr);
     411        1463 :     ierr = PetscFree(ri);CHKERRQ(ierr);
     412             :   }
     413        9859 :   PetscFunctionReturn(0);
     414             : }
     415             : 
     416          90 : static PetscErrorCode EPSXDComputeDSConv(dvdDashboard *d)
     417             : {
     418             :   PetscErrorCode    ierr;
     419             :   PetscInt          i,ld;
     420             :   Vec               v;
     421             :   PetscScalar       *pA;
     422             :   const PetscScalar *pv;
     423             :   PetscBool         symm;
     424             : 
     425          90 :   PetscFunctionBegin;
     426          90 :   ierr = BVSetActiveColumns(d->eps->V,0,d->eps->nconv);CHKERRQ(ierr);
     427          90 :   ierr = PetscObjectTypeCompareAny((PetscObject)d->eps->ds,&symm,DSHEP,"");CHKERRQ(ierr);
     428          90 :   if (symm) PetscFunctionReturn(0);
     429          24 :   ierr = DSSetDimensions(d->eps->ds,d->eps->nconv,0,0,0);CHKERRQ(ierr);
     430          24 :   ierr = DSCopyMat(d->eps->ds,DS_MAT_A,0,0,d->H,0,0,d->eps->nconv,d->eps->nconv,PETSC_FALSE);CHKERRQ(ierr);
     431          24 :   if (d->G) {
     432          21 :     ierr = DSCopyMat(d->eps->ds,DS_MAT_B,0,0,d->G,0,0,d->eps->nconv,d->eps->nconv,PETSC_FALSE);CHKERRQ(ierr);
     433             :   }
     434             :   /* Set the signature on projected matrix B */
     435          24 :   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
     436           0 :     ierr = DSGetLeadingDimension(d->eps->ds,&ld);CHKERRQ(ierr);
     437           0 :     ierr = DSGetArray(d->eps->ds,DS_MAT_B,&pA);CHKERRQ(ierr);
     438           0 :     ierr = PetscArrayzero(pA,d->eps->nconv*ld);CHKERRQ(ierr);
     439           0 :     ierr = VecCreateSeq(PETSC_COMM_SELF,d->eps->nconv,&v);CHKERRQ(ierr);
     440           0 :     ierr = BVGetSignature(d->eps->V,v);CHKERRQ(ierr);
     441           0 :     ierr = VecGetArrayRead(v,&pv);CHKERRQ(ierr);
     442           0 :     for (i=0;i<d->eps->nconv;i++) pA[i+ld*i] = pv[i];
     443           0 :     ierr = VecRestoreArrayRead(v,&pv);CHKERRQ(ierr);
     444           0 :     ierr = VecDestroy(&v);CHKERRQ(ierr);
     445           0 :     ierr = DSRestoreArray(d->eps->ds,DS_MAT_B,&pA);CHKERRQ(ierr);
     446             :   }
     447          24 :   ierr = DSSetState(d->eps->ds,DS_STATE_RAW);CHKERRQ(ierr);
     448          24 :   ierr = DSSolve(d->eps->ds,d->eps->eigr,d->eps->eigi);CHKERRQ(ierr);
     449          24 :   ierr = DSSynchronize(d->eps->ds,d->eps->eigr,d->eps->eigi);CHKERRQ(ierr);
     450          24 :   if (d->W) {
     451          59 :     for (i=0; i<d->eps->nconv; i++) {
     452          59 :       ierr = d->calcpairs_eig_backtrans(d,d->eps->eigr[i],d->eps->eigi[i],&d->eps->eigr[i],&d->eps->eigi[i]);CHKERRQ(ierr);
     453             :     }
     454             :   }
     455          24 :   PetscFunctionReturn(0);
     456             : }
     457             : 
     458             : /*
     459             :    Compute the residual vectors R(i) <- (AV - BV*eigr(i))*pX(i), and also
     460             :    the norm associated to the Schur pair, where i = r_s..r_e
     461             : */
     462        2644 : static PetscErrorCode dvd_calcpairs_res_0(dvdDashboard *d,PetscInt r_s,PetscInt r_e)
     463             : {
     464             :   PetscInt       i,ldpX;
     465             :   PetscScalar    *pX;
     466             :   PetscErrorCode ierr;
     467        2644 :   BV             BX = d->BX?d->BX:d->eps->V;
     468             :   Vec            *R;
     469             : 
     470        2644 :   PetscFunctionBegin;
     471        2644 :   ierr = DSGetLeadingDimension(d->eps->ds,&ldpX);CHKERRQ(ierr);
     472        2644 :   ierr = DSGetArray(d->eps->ds,DS_MAT_Q,&pX);CHKERRQ(ierr);
     473             :   /* nX(i) <- ||X(i)|| */
     474        2644 :   ierr = dvd_improvex_compute_X(d,r_s,r_e,NULL,pX,ldpX);CHKERRQ(ierr);
     475        2644 :   ierr = SlepcVecPoolGetVecs(d->auxV,r_e-r_s,&R);CHKERRQ(ierr);
     476        2667 :   for (i=r_s;i<r_e;i++) {
     477             :     /* R(i-r_s) <- AV*pX(i) */
     478        2667 :     ierr = BVMultVec(d->AX,1.0,0.0,R[i-r_s],&pX[ldpX*i]);CHKERRQ(ierr);
     479             :     /* R(i-r_s) <- R(i-r_s) - eigr(i)*BX*pX(i) */
     480        2667 :     ierr = BVMultVec(BX,-d->eigr[i],1.0,R[i-r_s],&pX[ldpX*i]);CHKERRQ(ierr);
     481             :   }
     482        2644 :   ierr = DSRestoreArray(d->eps->ds,DS_MAT_Q,&pX);CHKERRQ(ierr);
     483        2644 :   ierr = d->calcpairs_proj_res(d,r_s,r_e,R);CHKERRQ(ierr);
     484        2644 :   ierr = SlepcVecPoolRestoreVecs(d->auxV,r_e-r_s,&R);CHKERRQ(ierr);
     485        2644 :   PetscFunctionReturn(0);
     486             : }
     487             : 
     488       11393 : static PetscErrorCode dvd_calcpairs_proj_res(dvdDashboard *d,PetscInt r_s,PetscInt r_e,Vec *R)
     489             : {
     490             :   PetscInt       i,l,k;
     491             :   PetscErrorCode ierr;
     492       11393 :   PetscBool      lindep=PETSC_FALSE;
     493             :   BV             cX;
     494             : 
     495       11393 :   PetscFunctionBegin;
     496       11393 :   if (d->W) cX = d->W; /* If left subspace exists, R <- orth(cY, R), nR[i] <- ||R[i]|| */
     497        9924 :   else if (!(DVD_IS(d->sEP, DVD_EP_STD) && DVD_IS(d->sEP, DVD_EP_HERMITIAN))) cX = d->eps->V; /* If not HEP, R <- orth(cX, R), nR[i] <- ||R[i]|| */
     498             :   else cX = NULL; /* Otherwise, nR[i] <- ||R[i]|| */
     499             : 
     500       11393 :   if (cX) {
     501        3182 :     ierr = BVGetActiveColumns(cX,&l,&k);CHKERRQ(ierr);
     502        3182 :     ierr = BVSetActiveColumns(cX,0,l);CHKERRQ(ierr);
     503        3341 :     for (i=0;i<r_e-r_s;i++) {
     504        3341 :       ierr = BVOrthogonalizeVec(cX,R[i],NULL,&d->nR[r_s+i],&lindep);CHKERRQ(ierr);
     505             :     }
     506        3182 :     ierr = BVSetActiveColumns(cX,l,k);CHKERRQ(ierr);
     507        3182 :     if (lindep || (PetscAbs(d->nR[r_s+i]) < PETSC_MACHINE_EPSILON)) {
     508           0 :       ierr = PetscInfo2(d->eps,"The computed eigenvector residual %D is too low, %g!\n",r_s+i,(double)(d->nR[r_s+i]));CHKERRQ(ierr);
     509             :     }
     510             :   } else {
     511        8211 :     for (i=0;i<r_e-r_s;i++) {
     512        8211 :       ierr = VecNormBegin(R[i],NORM_2,&d->nR[r_s+i]);CHKERRQ(ierr);
     513             :     }
     514        8211 :     for (i=0;i<r_e-r_s;i++) {
     515        8211 :       ierr = VecNormEnd(R[i],NORM_2,&d->nR[r_s+i]);CHKERRQ(ierr);
     516             :     }
     517             :   }
     518       11393 :   PetscFunctionReturn(0);
     519             : }
     520             : 
     521         270 : PetscErrorCode dvd_calcpairs_qz(dvdDashboard *d,dvdBlackboard *b,PetscBool borth,PetscBool harm)
     522             : {
     523             :   PetscErrorCode ierr;
     524             :   PetscBool      std_probl,her_probl,ind_probl;
     525             :   DSType         dstype;
     526             :   Vec            v1;
     527             : 
     528         270 :   PetscFunctionBegin;
     529         270 :   std_probl = DVD_IS(d->sEP,DVD_EP_STD)? PETSC_TRUE: PETSC_FALSE;
     530         270 :   her_probl = DVD_IS(d->sEP,DVD_EP_HERMITIAN)? PETSC_TRUE: PETSC_FALSE;
     531         270 :   ind_probl = DVD_IS(d->sEP,DVD_EP_INDEFINITE)? PETSC_TRUE: PETSC_FALSE;
     532             : 
     533             :   /* Setting configuration constrains */
     534         270 :   b->max_size_proj = PetscMax(b->max_size_proj,b->max_size_V);
     535         270 :   d->W_shift = d->B? PETSC_TRUE: PETSC_FALSE;
     536             : 
     537             :   /* Setup the step */
     538         270 :   if (b->state >= DVD_STATE_CONF) {
     539          90 :     d->max_size_P = b->max_size_P;
     540          90 :     d->max_size_proj = b->max_size_proj;
     541             :     /* Create a DS if the method works with Schur decompositions */
     542          90 :     d->calcPairs = dvd_calcpairs_proj;
     543          90 :     d->calcpairs_residual = dvd_calcpairs_res_0;
     544          90 :     d->calcpairs_proj_res = dvd_calcpairs_proj_res;
     545          90 :     d->calcpairs_selectPairs = dvd_calcpairs_selectPairs;
     546             :     /* Create and configure a DS for solving the projected problems */
     547          90 :     if (d->W) dstype = DSGNHEP;    /* If we use harmonics */
     548             :     else {
     549          90 :       if (ind_probl) dstype = DSGHIEP;
     550          90 :       else if (std_probl) dstype = her_probl? DSHEP : DSNHEP;
     551          21 :       else dstype = her_probl? DSGHEP : DSGNHEP;
     552             :     }
     553          90 :     ierr = DSSetType(d->eps->ds,dstype);CHKERRQ(ierr);
     554          90 :     ierr = DSAllocate(d->eps->ds,d->eps->ncv);CHKERRQ(ierr);
     555             :     /* Create various vector basis */
     556          90 :     if (harm) {
     557          21 :       ierr = BVDuplicateResize(d->eps->V,d->eps->ncv,&d->W);CHKERRQ(ierr);
     558          21 :       ierr = BVSetMatrix(d->W,NULL,PETSC_FALSE);CHKERRQ(ierr);
     559          69 :     } else d->W = NULL;
     560          90 :     ierr = BVDuplicateResize(d->eps->V,d->eps->ncv,&d->AX);CHKERRQ(ierr);
     561          90 :     ierr = BVSetMatrix(d->AX,NULL,PETSC_FALSE);CHKERRQ(ierr);
     562          90 :     ierr = BVDuplicateResize(d->eps->V,d->eps->ncv,&d->auxBV);CHKERRQ(ierr);
     563          90 :     ierr = BVSetMatrix(d->auxBV,NULL,PETSC_FALSE);CHKERRQ(ierr);
     564          90 :     if (d->B) {
     565          12 :       ierr = BVDuplicateResize(d->eps->V,d->eps->ncv,&d->BX);CHKERRQ(ierr);
     566          12 :       ierr = BVSetMatrix(d->BX,NULL,PETSC_FALSE);CHKERRQ(ierr);
     567          78 :     } else d->BX = NULL;
     568          90 :     ierr = MatCreateVecsEmpty(d->A,&v1,NULL);CHKERRQ(ierr);
     569          90 :     ierr = SlepcVecPoolCreate(v1,0,&d->auxV);CHKERRQ(ierr);
     570          90 :     ierr = VecDestroy(&v1);CHKERRQ(ierr);
     571             :     /* Create projected problem matrices */
     572          90 :     ierr = MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->H);CHKERRQ(ierr);
     573          90 :     if (!std_probl) {
     574          21 :       ierr = MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->G);CHKERRQ(ierr);
     575          69 :     } else d->G = NULL;
     576          90 :     if (her_probl) {
     577          66 :       ierr = MatSetOption(d->H,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr);
     578          66 :       if (d->G) { ierr = MatSetOption(d->G,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr); }
     579             :     }
     580             : 
     581          90 :     if (ind_probl) {
     582           0 :       ierr = PetscMalloc1(d->eps->ncv,&d->nBds);CHKERRQ(ierr);
     583          90 :     } else d->nBds = NULL;
     584          90 :     ierr = MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->auxM);CHKERRQ(ierr);
     585             : 
     586          90 :     ierr = EPSDavidsonFLAdd(&d->startList,dvd_calcpairs_qz_start);CHKERRQ(ierr);
     587          90 :     ierr = EPSDavidsonFLAdd(&d->endList,EPSXDComputeDSConv);CHKERRQ(ierr);
     588          90 :     ierr = EPSDavidsonFLAdd(&d->destroyList,dvd_calcpairs_qz_d);CHKERRQ(ierr);
     589             :   }
     590         270 :   PetscFunctionReturn(0);
     591             : }
     592             : 

Generated by: LCOV version 1.13