1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  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"

 13:    Step: initialize subspace V
 14: */

 16: #include "davidson.h"

 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;

 24: static PetscErrorCode dvd_initV_classic_0(dvdDashboard *d)
 25: {
 26:   dvdInitV       *data = (dvdInitV*)d->initV_data;
 27:   PetscInt       i,user = PetscMin(data->user,d->eps->mpd), l,k;

 29:   PetscFunctionBegin;
 30:   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
 31:   /* User vectors are added at the beginning, so no active column should be in V */
 32:   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:   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:   d->V_tra_s = 0; d->V_tra_e = 0;
 36:   d->V_new_s = 0; d->V_new_e = i-l;

 38:   /* After that the user vectors will be destroyed */
 39:   data->user = 0;
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: static PetscErrorCode dvd_initV_krylov_0(dvdDashboard *d)
 44: {
 45:   dvdInitV       *data = (dvdInitV*)d->initV_data;
 46:   PetscInt       i,user = PetscMin(data->user,d->eps->mpd),l,k;
 47:   Vec            av,v1,v2;

 49:   PetscFunctionBegin;
 50:   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
 51:   /* User vectors are added at the beginning, so no active column should be in V */
 52:   PetscAssert(data->user==0 || l==0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");

 54:   /* If needed, generate a random vector for starting the arnoldi method */
 55:   if (user == 0) {
 56:     PetscCall(BVSetRandomColumn(d->eps->V,l));
 57:     user = 1;
 58:   }

 60:   /* Perform k steps of Arnoldi with the operator K^{-1}*(t[1]*A-t[2]*B) */
 61:   PetscCall(dvd_orthV(d->eps->V,l,l+user));
 62:   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:     PetscCall(BVGetColumn(d->eps->V,i,&v1));
 65:     PetscCall(BVGetColumn(d->eps->V,i-user,&v2));
 66:     PetscCall(BVGetColumn(d->auxBV,0,&av));
 67:     if (d->B) {
 68:       PetscCall(MatMult(d->A,v2,v1));
 69:       PetscCall(MatMult(d->B,v2,av));
 70:       PetscCall(VecAXPBY(av,d->target[1],-d->target[0],v1));
 71:     } else {
 72:       PetscCall(MatMult(d->A,v2,av));
 73:       PetscCall(VecAXPBY(av,-d->target[0],d->target[1],v2));
 74:     }
 75:     PetscCall(d->improvex_precond(d,0,av,v1));
 76:     PetscCall(BVRestoreColumn(d->eps->V,i,&v1));
 77:     PetscCall(BVRestoreColumn(d->eps->V,i-user,&v2));
 78:     PetscCall(BVRestoreColumn(d->auxBV,0,&av));
 79:     PetscCall(dvd_orthV(d->eps->V,i,i+1));
 80:   }

 82:   d->V_tra_s = 0; d->V_tra_e = 0;
 83:   d->V_new_s = 0; d->V_new_e = i-l;

 85:   /* After that the user vectors will be destroyed */
 86:   data->user = 0;
 87:   PetscFunctionReturn(PETSC_SUCCESS);
 88: }

 90: static PetscErrorCode dvd_initV_d(dvdDashboard *d)
 91: {
 92:   dvdInitV       *data = (dvdInitV*)d->initV_data;

 94:   PetscFunctionBegin;
 95:   /* Restore changes in dvdDashboard */
 96:   d->initV_data = data->old_initV_data;

 98:   /* Free local data */
 99:   PetscCall(PetscFree(data));
100:   PetscFunctionReturn(PETSC_SUCCESS);
101: }

103: PetscErrorCode dvd_initV(dvdDashboard *d, dvdBlackboard *b, PetscInt k,PetscInt user, PetscBool krylov)
104: {
105:   dvdInitV       *data;

107:   PetscFunctionBegin;
108:   /* Setting configuration constrains */
109:   b->max_size_V = PetscMax(b->max_size_V, k);

111:   /* Setup the step */
112:   if (b->state >= DVD_STATE_CONF) {
113:     PetscCall(PetscNew(&data));
114:     data->k = k;
115:     data->user = user;
116:     data->old_initV_data = d->initV_data;
117:     d->initV_data = data;
118:     if (krylov) d->initV = dvd_initV_krylov_0;
119:     else d->initV = dvd_initV_classic_0;
120:     PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_initV_d));
121:   }
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: PetscErrorCode dvd_orthV(BV V,PetscInt V_new_s,PetscInt V_new_e)
126: {
127:   PetscInt       i;

129:   PetscFunctionBegin;
130:   for (i=V_new_s;i<V_new_e;i++) PetscCall(BVOrthonormalizeColumn(V,i,PETSC_TRUE,NULL,NULL));
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }