Actual source code: dvdinitv.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: {
27: dvdInitV *data = (dvdInitV*)d->initV_data;
28: PetscInt i,user = PetscMin(data->user,d->eps->mpd), l,k;
31: BVGetActiveColumns(d->eps->V,&l,&k);
32: /* User vectors are added at the beginning, so no active column should be in V */
33: if (data->user>0&&l>0) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
34: /* Generate a set of random initial vectors and orthonormalize them */
35: for (i=l+user;i<l+data->k && i<d->eps->ncv && i-l<d->eps->mpd;i++) {
36: BVSetRandomColumn(d->eps->V,i);
37: }
38: d->V_tra_s = 0; d->V_tra_e = 0;
39: d->V_new_s = 0; d->V_new_e = i-l;
41: /* After that the user vectors will be destroyed */
42: data->user = 0;
43: return(0);
44: }
46: static PetscErrorCode dvd_initV_krylov_0(dvdDashboard *d)
47: {
49: dvdInitV *data = (dvdInitV*)d->initV_data;
50: PetscInt i,user = PetscMin(data->user,d->eps->mpd),l,k;
51: Vec av,v1,v2;
54: BVGetActiveColumns(d->eps->V,&l,&k);
55: /* User vectors are added at the beginning, so no active column should be in V */
56: if (data->user>0&&l>0) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
58: /* If needed, generate a random vector for starting the arnoldi method */
59: if (user == 0) {
60: BVSetRandomColumn(d->eps->V,l);
61: user = 1;
62: }
64: /* Perform k steps of Arnoldi with the operator K^{-1}*(t[1]*A-t[2]*B) */
65: dvd_orthV(d->eps->V,l,l+user);
66: for (i=l+user;i<l+data->k && i<d->eps->ncv && i-l<d->eps->mpd;i++) {
67: /* aux <- theta[1]A*in - theta[0]*B*in */
68: BVGetColumn(d->eps->V,i,&v1);
69: BVGetColumn(d->eps->V,i-user,&v2);
70: BVGetColumn(d->auxBV,0,&av);
71: if (d->B) {
72: MatMult(d->A,v2,v1);
73: MatMult(d->B,v2,av);
74: VecAXPBY(av,d->target[1],-d->target[0],v1);
75: } else {
76: MatMult(d->A,v2,av);
77: VecAXPBY(av,-d->target[0],d->target[1],v2);
78: }
79: d->improvex_precond(d,0,av,v1);
80: BVRestoreColumn(d->eps->V,i,&v1);
81: BVRestoreColumn(d->eps->V,i-user,&v2);
82: BVRestoreColumn(d->auxBV,0,&av);
83: dvd_orthV(d->eps->V,i,i+1);
84: }
86: d->V_tra_s = 0; d->V_tra_e = 0;
87: d->V_new_s = 0; d->V_new_e = i-l;
89: /* After that the user vectors will be destroyed */
90: data->user = 0;
91: return(0);
92: }
94: static PetscErrorCode dvd_initV_d(dvdDashboard *d)
95: {
97: dvdInitV *data = (dvdInitV*)d->initV_data;
100: /* Restore changes in dvdDashboard */
101: d->initV_data = data->old_initV_data;
103: /* Free local data */
104: PetscFree(data);
105: return(0);
106: }
108: PetscErrorCode dvd_initV(dvdDashboard *d, dvdBlackboard *b, PetscInt k,PetscInt user, PetscBool krylov)
109: {
111: dvdInitV *data;
114: /* Setting configuration constrains */
115: b->max_size_V = PetscMax(b->max_size_V, k);
117: /* Setup the step */
118: if (b->state >= DVD_STATE_CONF) {
119: PetscNewLog(d->eps,&data);
120: data->k = k;
121: data->user = user;
122: data->old_initV_data = d->initV_data;
123: d->initV_data = data;
124: if (krylov) d->initV = dvd_initV_krylov_0;
125: else d->initV = dvd_initV_classic_0;
126: EPSDavidsonFLAdd(&d->destroyList,dvd_initV_d);
127: }
128: return(0);
129: }
131: PetscErrorCode dvd_orthV(BV V,PetscInt V_new_s,PetscInt V_new_e)
132: {
134: PetscInt i;
137: for (i=V_new_s;i<V_new_e;i++) {
138: BVOrthonormalizeColumn(V,i,PETSC_TRUE,NULL,NULL);
139: }
140: return(0);
141: }