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             :    SLEPc eigensolver: "davidson"
      12             : 
      13             :    Step: initialize subspace V
      14             : */
      15             : 
      16             : #include "davidson.h"
      17             : 
      18             : typedef struct {
      19             :   PetscInt k;                 /* desired initial subspace size */
      20             :   PetscInt user;              /* number of user initial vectors */
      21             :   void     *old_initV_data;   /* old initV data */
      22             : } dvdInitV;
      23             : 
      24          87 : static PetscErrorCode dvd_initV_classic_0(dvdDashboard *d)
      25             : {
      26          87 :   dvdInitV       *data = (dvdInitV*)d->initV_data;
      27          87 :   PetscInt       i,user = PetscMin(data->user,d->eps->mpd), l,k;
      28             : 
      29          87 :   PetscFunctionBegin;
      30          87 :   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
      31             :   /* User vectors are added at the beginning, so no active column should be in V */
      32          87 :   PetscAssert(data->user==0 || l==0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
      33             :   /* Generate a set of random initial vectors and orthonormalize them */
      34         590 :   for (i=l+user;i<l+data->k && i<d->eps->ncv && i-l<d->eps->mpd;i++) PetscCall(BVSetRandomColumn(d->eps->V,i));
      35          87 :   d->V_tra_s = 0; d->V_tra_e = 0;
      36          87 :   d->V_new_s = 0; d->V_new_e = i-l;
      37             : 
      38             :   /* After that the user vectors will be destroyed */
      39          87 :   data->user = 0;
      40          87 :   PetscFunctionReturn(PETSC_SUCCESS);
      41             : }
      42             : 
      43          10 : static PetscErrorCode dvd_initV_krylov_0(dvdDashboard *d)
      44             : {
      45          10 :   dvdInitV       *data = (dvdInitV*)d->initV_data;
      46          10 :   PetscInt       i,user = PetscMin(data->user,d->eps->mpd),l,k;
      47          10 :   Vec            av,v1,v2;
      48             : 
      49          10 :   PetscFunctionBegin;
      50          10 :   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
      51             :   /* User vectors are added at the beginning, so no active column should be in V */
      52          10 :   PetscAssert(data->user==0 || l==0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
      53             : 
      54             :   /* If needed, generate a random vector for starting the arnoldi method */
      55          10 :   if (user == 0) {
      56           9 :     PetscCall(BVSetRandomColumn(d->eps->V,l));
      57             :     user = 1;
      58             :   }
      59             : 
      60             :   /* Perform k steps of Arnoldi with the operator K^{-1}*(t[1]*A-t[2]*B) */
      61          10 :   PetscCall(dvd_orthV(d->eps->V,l,l+user));
      62          60 :   for (i=l+user;i<l+data->k && i<d->eps->ncv && i-l<d->eps->mpd;i++) {
      63             :     /* aux <- theta[1]A*in - theta[0]*B*in */
      64          50 :     PetscCall(BVGetColumn(d->eps->V,i,&v1));
      65          50 :     PetscCall(BVGetColumn(d->eps->V,i-user,&v2));
      66          50 :     PetscCall(BVGetColumn(d->auxBV,0,&av));
      67          50 :     if (d->B) {
      68           0 :       PetscCall(MatMult(d->A,v2,v1));
      69           0 :       PetscCall(MatMult(d->B,v2,av));
      70           0 :       PetscCall(VecAXPBY(av,d->target[1],-d->target[0],v1));
      71             :     } else {
      72          50 :       PetscCall(MatMult(d->A,v2,av));
      73          50 :       PetscCall(VecAXPBY(av,-d->target[0],d->target[1],v2));
      74             :     }
      75          50 :     PetscCall(d->improvex_precond(d,0,av,v1));
      76          50 :     PetscCall(BVRestoreColumn(d->eps->V,i,&v1));
      77          50 :     PetscCall(BVRestoreColumn(d->eps->V,i-user,&v2));
      78          50 :     PetscCall(BVRestoreColumn(d->auxBV,0,&av));
      79          50 :     PetscCall(dvd_orthV(d->eps->V,i,i+1));
      80             :   }
      81             : 
      82          10 :   d->V_tra_s = 0; d->V_tra_e = 0;
      83          10 :   d->V_new_s = 0; d->V_new_e = i-l;
      84             : 
      85             :   /* After that the user vectors will be destroyed */
      86          10 :   data->user = 0;
      87          10 :   PetscFunctionReturn(PETSC_SUCCESS);
      88             : }
      89             : 
      90          97 : static PetscErrorCode dvd_initV_d(dvdDashboard *d)
      91             : {
      92          97 :   dvdInitV       *data = (dvdInitV*)d->initV_data;
      93             : 
      94          97 :   PetscFunctionBegin;
      95             :   /* Restore changes in dvdDashboard */
      96          97 :   d->initV_data = data->old_initV_data;
      97             : 
      98             :   /* Free local data */
      99          97 :   PetscCall(PetscFree(data));
     100          97 :   PetscFunctionReturn(PETSC_SUCCESS);
     101             : }
     102             : 
     103         291 : PetscErrorCode dvd_initV(dvdDashboard *d, dvdBlackboard *b, PetscInt k,PetscInt user, PetscBool krylov)
     104             : {
     105         291 :   dvdInitV       *data;
     106             : 
     107         291 :   PetscFunctionBegin;
     108             :   /* Setting configuration constrains */
     109         291 :   b->max_size_V = PetscMax(b->max_size_V, k);
     110             : 
     111             :   /* Setup the step */
     112         291 :   if (b->state >= DVD_STATE_CONF) {
     113          97 :     PetscCall(PetscNew(&data));
     114          97 :     data->k = k;
     115          97 :     data->user = user;
     116          97 :     data->old_initV_data = d->initV_data;
     117          97 :     d->initV_data = data;
     118          97 :     if (krylov) d->initV = dvd_initV_krylov_0;
     119          87 :     else d->initV = dvd_initV_classic_0;
     120          97 :     PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_initV_d));
     121             :   }
     122         291 :   PetscFunctionReturn(PETSC_SUCCESS);
     123             : }
     124             : 
     125        8732 : PetscErrorCode dvd_orthV(BV V,PetscInt V_new_s,PetscInt V_new_e)
     126             : {
     127        8732 :   PetscInt       i;
     128             : 
     129        8732 :   PetscFunctionBegin;
     130       21194 :   for (i=V_new_s;i<V_new_e;i++) PetscCall(BVOrthonormalizeColumn(V,i,PETSC_TRUE,NULL,NULL));
     131        8732 :   PetscFunctionReturn(PETSC_SUCCESS);
     132             : }