/*
SLEPc eigensolver: "davidson"
Step: init subspace V
*/
#include "davidsones.h"
PetscInt dvd_initV_classic_0(dvdDashboard *d);
PetscErrorCode dvd_initV_classic_d(dvdDashboard *d);
PetscInt dvd_initV_krylov_0(dvdDashboard *d);
PetscErrorCode dvd_initV_krylov_d(dvdDashboard *d);
/*
Fill V with a random subspace
*/
typedef struct {
PetscInt k; /* number of vectors initialized */
void *old_initV_data; /* old initV data */
} dvdInitV_Classic;
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_initV_classic"
PetscInt dvd_initV_classic(dvdDashboard *d, dvdBlackboard *b, PetscInt k)
{
PetscErrorCode ierr;
dvdInitV_Classic
*data;
PetscFunctionBegin;
/* Setting configuration constrains */
b->max_size_V = PetscMax(b->max_size_V, k);
/* Setup the step */
if (b->state >= DVD_STATE_CONF) {
ierr = PetscMalloc(sizeof(dvdInitV_Classic), &data); CHKERRQ(ierr);
data->k = k;
data->old_initV_data = d->initV_data;
d->initV_data = data;
d->initV = dvd_initV_classic_0;
DVD_FL_ADD(d->destroyList, dvd_initV_classic_d);
}
PetscFunctionReturn(0);
}
EXTERN_C_END
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_initV_classic_0"
PetscInt dvd_initV_classic_0(dvdDashboard *d)
{
PetscErrorCode ierr;
dvdInitV_Classic
*data = (dvdInitV_Classic*)d->initV_data;
PetscInt i;
PetscFunctionBegin;
/* Generate a set of random initial vectors and orthonormalize them */
for (i=0; i<data->k; i++) {
ierr = SlepcVecSetRandom(d->V[i]); CHKERRQ(ierr);
}
d->size_V = data->k;
(*d->e_Vchanged)(d, 0, 0, 0, d->size_V);
CHKERRQ(ierr);
PetscFunctionReturn(0);
}
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_initV_classic_d"
PetscErrorCode dvd_initV_classic_d(dvdDashboard *d)
{
PetscErrorCode ierr;
dvdInitV_Classic
*data = (dvdInitV_Classic*)d->initV_data;
PetscFunctionBegin;
/* Restore changes in dvdDashboard */
d->initV_data = data->old_initV_data;
/* Free local data */
ierr = PetscFree(data); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
EXTERN_C_END
/*
Start with a krylov subspace with the matrix A
*/
#undef __FUNCT__
#define __FUNCT__ "dvd_BasicArnoldi"
/*
EPSBasicArnoldi - Computes an m-step Arnoldi factorization. The first k
columns are assumed to be locked and therefore they are not modified. On
exit, the following relation is satisfied:
OP * V - V * H = f * e_m^T
where the columns of V are the Arnoldi vectors (which are B-orthonormal),
H is an upper Hessenberg matrix, f is the residual vector and e_m is
the m-th vector of the canonical basis. The vector f is B-orthogonal to
the columns of V. On exit, beta contains the B-norm of f and the next
Arnoldi vector can be computed as v_{m+1} = f / beta.
auxS is an auxiliary vector of size *M
*/
PetscErrorCode dvd_BasicArnoldi(Mat OP, IP ip, PetscScalar *H, PetscInt ldh,
Vec *V, PetscInt k, PetscInt *M, Vec f, PetscReal *beta,
PetscTruth *breakdown, Vec auxV, PetscScalar *auxS)
{
PetscErrorCode ierr;
PetscInt j,m = *M;
PetscReal norm;
PetscFunctionBegin;
for (j=k; j<m-1; j++) {
ierr = MatMult(OP, V[j], V[j+1]); CHKERRQ(ierr);
ierr = IPOrthogonalize(ip, j+1, PETSC_NULL, V, V[j+1], H+ldh*j, &norm,
breakdown, auxV, auxS); CHKERRQ(ierr);
H[j+1+ldh*j] = norm;
if (*breakdown) {
*M = j+1;
*beta = norm;
PetscFunctionReturn(0);
} else {
ierr = VecScale(V[j+1], 1.0/norm); CHKERRQ(ierr);
}
}
ierr = MatMult(OP, V[m-1], f); CHKERRQ(ierr);
ierr = IPOrthogonalize(ip, m, PETSC_NULL, V, f, H+ldh*(m-1), beta, PETSC_NULL,
auxV, auxS); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
typedef struct {
PetscInt k; /* number of steps of arnoldi */
Mat OP; /* operator of the krylov subspace */
IP ip; /* orthogonalization routine */
void *old_initV_data; /* old initV data */
} dvdInitV_Krylov;
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_initV_krylov"
PetscInt dvd_initV_krylov(dvdDashboard *d, dvdBlackboard *b, PetscInt k,
Mat OP, IP ip)
{
PetscErrorCode ierr;
dvdInitV_Krylov *data;
PetscFunctionBegin;
/* Setting configuration constrains */
b->max_size_V = PetscMax(b->max_size_V, k);
b->max_size_auxV = PetscMax(b->max_size_auxV, 1);
b->max_size_auxS = PetscMax(b->max_size_auxS, k+1);
/* Setup the step */
if (b->state >= DVD_STATE_CONF) {
ierr = PetscMalloc(sizeof(dvdInitV_Krylov), &data); CHKERRQ(ierr);
data->k = k;
data->OP = OP;
data->ip = ip;
data->old_initV_data = d->initV_data;
d->initV_data = data;
d->initV = dvd_initV_krylov_0;
DVD_FL_ADD(d->destroyList, dvd_initV_krylov_d);
}
PetscFunctionReturn(0);
}
EXTERN_C_END
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_initV_krylov_0"
PetscInt dvd_initV_krylov_0(dvdDashboard *d)
{
PetscErrorCode ierr;
dvdInitV_Krylov *data = (dvdInitV_Krylov*)d->initV_data;
PetscScalar beta;
PetscReal norm;
PetscTruth breakdown;
PetscInt i;
PetscFunctionBegin;
/* Generate a random vector for starting the arnoldi method */
ierr = SlepcVecSetRandom(d->V[0]); CHKERRQ(ierr);
ierr = IPNorm(data->ip, d->V[0], &norm); CHKERRQ(ierr);
ierr = VecScale(d->V[0], 1.0/norm); CHKERRQ(ierr);
/* Perform k Arnoldi steps. After that,
V(:,0:k-1) will be filled with the subspace K(A,V[0]), and
H = V' * A * V
*/
i = data->k+1;
ierr = dvd_BasicArnoldi(d->A, data->ip, d->H, d->ldH, d->V, 0, &i,
d->W[data->k-1]/* f */, &beta, &breakdown,
d->auxV[0], d->auxS); CHKERRQ(ierr);
d->size_V = data->k;
d->size_H = 0/* data->k */;
/* Calc W <- V * H + f * e_m^T */
/*for (i=0; i<data->k; i++) {
if (i != data->k - 1) {
ierr = VecSet(d->W[i], 0.0); CHKERRQ(ierr);
}
ierr = VecMAXPY(d->W[i], data->k, d->H+d->ldH*i, d->V); CHKERRQ(ierr);
}
d->size_W = data->k;*/
d->size_W = 0;
(*d->e_Vchanged)(d, 0, 0, 0, d->size_V);
PetscFunctionReturn(0);
}
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_initV_krylov_d"
PetscErrorCode dvd_initV_krylov_d(dvdDashboard *d)
{
PetscErrorCode ierr;
dvdInitV_Krylov *data = (dvdInitV_Krylov*)d->initV_data;
PetscFunctionBegin;
/* Restore changes in dvdDashboard */
d->initV_data = data->old_initV_data;
/* Free local data */
ierr = PetscFree(data); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
EXTERN_C_END