/*
SLEPc eigensolver: "davidson"
Step: calc the best eigenpairs in the subspace V.
For that, performs these steps:
1) Update W <- A * V
2) Update H <- V' * W
3) Obtain eigenpairs of H
4) Select some eigenpairs
5) Compute the Ritz pairs of the selected ones
*/
#include "slepc.h"
#include "private/epsimpl.h"
#include "davidsones.h"
PetscInt dvd_eVchanged_rr(dvdDashboard *d, PetscInt s_imm, PetscInt e_imm,
PetscInt s_new, PetscInt e_new);
PetscInt dvd_calcpairs_rr_0(dvdDashboard *d);
PetscInt dvd_calcpairs_residual_rr(dvdDashboard *d, PetscInt s, PetscInt e,
Vec *R);
PetscInt dvd_calcpairs_rr_d(dvdDashboard *d);
PetscInt dvd_calcpairs_rr_sym(dvdDashboard *d);
typedef struct {
PetscInt s_imm, /* start of immutable vectors in V */
e_imm, /* end of immutable vectors in V */
s_new, /* start of added vectors in V */
e_new, /* end of added vectors in V */
imm_size_V, /* size of V unchanged */
old_size_V; /* size of last iteration V */
Vec *real_W; /* real start of W */
PetscScalar *U, /* accumulated transformations */
*real_H; /* real start of H */
void
*old_calcPairs_data; /* old calcPairs data */
e_Vchanged_type
old_e_Vchanged; /* old e_Vchanged */
void
*old_calcpairs_residual;
/* old calcpairs_residual */
} dvdCalcpairs_rr;
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_rr"
PetscInt dvd_calcpairs_rr(dvdDashboard *d, dvdBlackboard *b)
{
PetscErrorCode ierr;
dvdCalcpairs_rr
*data;
PetscFunctionBegin;
/* Setting configuration constrains */
b->own_vecs+= b->max_size_V /* W */;
b->own_scalars+= b->max_size_V*b->max_size_V /* H */ +
b->max_size_V*b->max_size_V /* U */;
b->max_size_auxS = PetscMax(b->max_size_auxS,
b->max_size_V*b->max_size_V /* T */ +
b->max_size_V*b->max_size_oldX /* oldU */);
/* Setup the step */
if (b->state >= DVD_STATE_CONF) {
ierr = PetscMalloc(sizeof(dvdCalcpairs_rr), &data); CHKERRQ(ierr);
d->size_W = 0;
data->real_W = d->W = b->free_vecs; b->free_vecs+= b->max_size_V;
d->H = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_V;
data->real_H = d->H;
data->U = b->free_scalars; b->free_scalars+= b->max_size_V*b->max_size_V;
data->imm_size_V = 0;
data->old_size_V = 0;
d->ldH = b->max_size_V;
data->old_calcPairs_data = d->calcPairs_data;
d->calcPairs_data = data;
d->calcPairs = dvd_calcpairs_rr_sym;
data->old_e_Vchanged = d->e_Vchanged;
d->e_Vchanged = dvd_eVchanged_rr;
d->calcpairs_residual = dvd_calcpairs_residual_rr;
DVD_FL_ADD(d->destroyList, dvd_calcpairs_rr_d);
}
PetscFunctionReturn(0);
}
EXTERN_C_END
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_rr"
PetscInt dvd_eVchanged_rr(dvdDashboard *d, PetscInt s_imm, PetscInt e_imm,
PetscInt s_new, PetscInt e_new)
{
dvdCalcpairs_rr *data = (dvdCalcpairs_rr*)d->calcPairs_data;
PetscFunctionBegin;
data->s_imm = s_imm;
data->e_imm = e_imm;
data->s_new = s_new;
data->e_new = e_new;
/* Callback old e_Vchanged */
if (data->old_e_Vchanged)
(*data->old_e_Vchanged)(d, s_imm, e_imm, s_new, e_new);
PetscFunctionReturn(0);
}
EXTERN_C_END
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_rr"
PetscInt dvd_calcpairs_rr_0(dvdDashboard *d)
{
PetscErrorCode ierr;
dvdCalcpairs_rr
*data = (dvdCalcpairs_rr*)d->calcPairs_data;
PetscInt i,j,addToU=0;
PetscScalar *T = d->auxS, *oldU = d->auxS + d->ldH*d->size_V;
PetscFunctionBegin;
/* Compute a Rayleigh-Ritz projection step */
/* If the dimension of V is cut, then perform the cutting */
if (data->s_imm != 0) {
if (d->size_W - data->s_imm >= d->size_V) {
d->W+= data->s_imm;
d->size_W-= data->s_imm;
d->H+= d->ldH*data->s_imm + data->s_imm;
} else {
for(i=data->s_imm,j=0; i<data->e_imm; i++,j++) {
ierr = VecCopy(d->W[i], d->W[j]); CHKERRQ(ierr);
ierr = PetscMemcpy(d->H + d->ldH*j, d->H + d->ldH*i,
sizeof(PetscScalar)*d->ldH); CHKERRQ(ierr);
}
d->size_W = 0;
}
data->e_imm-= data->s_imm;
data->s_new-= data->s_imm;
data->e_new-= data->s_imm;
data->s_imm = 0;
}
/* W <- A * V */
for (i=0; i<d->size_V; i++) {
ierr = MatMult(d->A, d->V[i], d->W[i]); CHKERRQ(ierr);
}
/*for (i=data->s_new; i<data->e_new; i++) {
ierr = MatMult(d->A, d->V[i], d->W[i]); CHKERRQ(ierr);
}*/
/* H <- V' * W */
for (i=0; i<d->size_V; i++) {
ierr = VecMDot(d->W[i], d->size_V, d->V, d->H+i*d->ldH);
CHKERRQ(ierr);
}
/*for (i=data->s_new; i<data->e_new; i++) {
ierr = VecMDot(d->W[i], d->size_V, d->V, d->H+i*d->ldH);
CHKERRQ(ierr);
}*/
/* H(s_new:e_new-1,0:s_new-1) <- H(0:s_new-1,s_new:e_new-1)'
H is read contiguously and written with a stroke of ld */
/*for(j=data->s_new; j<data->e_new; j++)
for (i=0; i<data->s_new; i++)
d->H[d->ldH*i+j] = d->H[d->ldH*j+i];*/
/* Save the U if needed and possible */
if (data->imm_size_V > 0) addToU = PetscMin(d->n_oldX, d->max_size_oldX);
if (addToU > 0) {
ierr = PetscMemcpy(oldU, data->U,
sizeof(PetscScalar)*data->old_size_V*addToU);
CHKERRQ(ierr);
}
/* Reduce projected matrix to Hessenberg form: [U,T] = hess(H) */
ierr = PetscMemcpy(T, d->H, sizeof(PetscScalar)*d->ldH*d->size_V);
CHKERRQ(ierr);
ierr = EPSDenseHessenberg(d->size_V, 0, T, d->ldH, data->U);
CHKERRQ(ierr);
/* Reduce T to quasi-triangular (Schur) form */
/* NOTE: U has ld as d->size_V */
ierr = EPSDenseSchur(d->size_V, 0, T, d->ldH, data->U, d->eigr,
d->eigi); CHKERRQ(ierr);
/* Sort diagonal elements in T and accumulate rotations on U */
if (d->withTarget == PETSC_TRUE) {
ierr = EPSSortDenseSchurTarget(d->size_V, 0, T, d->ldH, data->U,
d->eigr, d->eigi, d->target, d->which);
CHKERRQ(ierr);
} else {
ierr = EPSSortDenseSchur(d->size_V, 0, T, d->ldH, data->U,
d->eigr, d->eigi, d->which); CHKERRQ(ierr);
}
/* X <- V * U(:,0:bs-1); oldX <- V * oldU */
if (d->X != d->V) {
for (i=0; i<d->n_calcPairs; i++) {
ierr = VecSet(d->X[i], 0.0); CHKERRQ(ierr);
ierr = VecMAXPY(d->X[i], d->size_V, data->U+d->size_V*i, d->V);
CHKERRQ(ierr);
}
for (i=0; i<addToU; i++) {
ierr = VecSet(d->oldX[i], 0.0); CHKERRQ(ierr);
ierr = VecMAXPY(d->oldX[i], data->imm_size_V,
oldU+data->old_size_V*i, d->V); CHKERRQ(ierr);
}
data->imm_size_V = data->old_size_V = d->size_V;
/* If V starts at X, and X and oldX are contiguous */
} else if ((addToU == 0) || (d->X+d->n_calcPairs == d->oldX)) {
/* U <- [ U(:,0:bs-1) oldU(0:(size_V-1),:) ] */
for (i=0; i<addToU; i++) {
ierr = PetscMemcpy(data->U+d->size_V*(i+d->n_calcPairs),
oldU+data->old_size_V*i,
sizeof(PetscScalar)*PetscMin(data->old_size_V, d->size_V));
CHKERRQ(ierr);
for (j=data->old_size_V; j<d->size_V; j++)
data->U[d->size_V*(i+d->n_calcPairs)+j] = 0.0;
}
ierr = SlepcUpdateVectors(d->size_V, d->V, 0, d->n_calcPairs+addToU,
data->U, d->size_V, PETSC_FALSE); CHKERRQ(ierr);
data->imm_size_V = data->old_size_V = 0;
} else {
SETERRQ(1, "This implementation of calcPairs not support this state of "
"Dashboard!");
}
d->size_X = d->n_calcPairs;
d->size_oldX = addToU;
PetscFunctionReturn(0);
}
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_rr_sym"
PetscInt dvd_calcpairs_rr_sym(dvdDashboard *d)
{
PetscErrorCode ierr;
dvdCalcpairs_rr
*data = (dvdCalcpairs_rr*)d->calcPairs_data;
PetscInt i,j,addToU=0;
PetscScalar *T = d->auxS, *oldU = d->auxS + d->ldH*d->size_V;
PetscFunctionBegin;
/* Compute a Rayleigh-Ritz projection step */
/* If the dimension of V is cut, then perform the cutting */
if (data->s_imm != 0) {
if (d->size_W - data->s_imm >= d->size_V) {
d->W+= data->s_imm;
d->size_W-= data->s_imm;
d->H+= d->ldH*data->s_imm + data->s_imm;
} else {
for(i=data->s_imm,j=0; i<data->e_imm; i++,j++) {
ierr = VecCopy(d->W[i], d->W[j]); CHKERRQ(ierr);
ierr = PetscMemcpy(d->H + d->ldH*j, d->H + d->ldH*i,
sizeof(PetscScalar)*d->ldH); CHKERRQ(ierr);
}
d->size_W = 0;
}
data->e_imm-= data->s_imm;
data->s_new-= data->s_imm;
data->e_new-= data->s_imm;
data->s_imm = 0;
}
/* W <- A * V */
if (d->size_W < data->s_new) {
SETERRQ(1, "This implementation of calcPairs not support this state of "
"Dashboard!");
}
/*for (i=0; i<d->size_V; i++) {
ierr = MatMult(d->A, d->V[i], d->W[i]); CHKERRQ(ierr);
}*/
for (i=d->size_W; i<data->e_new; i++) {
ierr = MatMult(d->A, d->V[i], d->W[i]); CHKERRQ(ierr);
}
d->size_W = i;
/* H <- V' * W */
/*for (i=0; i<d->size_V; i++) {
ierr = VecMDot(d->W[i], d->size_V, d->V, d->H+i*d->ldH);
CHKERRQ(ierr);
}
d->size_H = i;*/
if (d->size_H < data->s_new) {
SETERRQ(1, "This implementation of calcPairs not support this state of "
"Dashboard!");
}
for (i=d->size_H; i<d->size_V; i++) {
ierr = VecMDot(d->W[i], d->size_V, d->V, d->H+i*d->ldH);
CHKERRQ(ierr);
}
/* H(size_H:e_new-1,0:size_H-1) <- H(0:size_H-1,size_H:e_new-1)'
H is read contiguously and written with a stroke of ld */
for(j=d->size_H; j<d->size_V; j++)
for (i=0; i<d->size_H; i++)
d->H[d->ldH*i+j] = d->H[d->ldH*j+i];
d->size_H = j;
/* Save the U if needed and possible */
if (data->imm_size_V > 0) addToU = PetscMin(d->n_oldX, d->max_size_oldX);
if (addToU > 0) {
ierr = PetscMemcpy(oldU, data->U,
sizeof(PetscScalar)*data->old_size_V*addToU);
CHKERRQ(ierr);
}
/* Reduce projected matrix to Hessenberg form: [U,T] = hess(H) */
ierr = PetscMemcpy(T, d->H, sizeof(PetscScalar)*d->ldH*d->size_V);
CHKERRQ(ierr);
ierr = EPSDenseHessenberg(d->size_V, 0, T, d->ldH, data->U);
CHKERRQ(ierr);
/* Reduce T to quasi-triangular (Schur) form */
/* NOTE: U has ld as d->size_V */
ierr = EPSDenseSchur(d->size_V, 0, T, d->ldH, data->U, d->eigr,
d->eigi); CHKERRQ(ierr);
/* Sort diagonal elements in T and accumulate rotations on U */
if (d->withTarget == PETSC_TRUE) {
ierr = EPSSortDenseSchurTarget(d->size_V, 0, T, d->ldH, data->U,
d->eigr, d->eigi, d->target, d->which);
CHKERRQ(ierr);
} else {
ierr = EPSSortDenseSchur(d->size_V, 0, T, d->ldH, data->U,
d->eigr, d->eigi, d->which); CHKERRQ(ierr);
}
/* X <- V * U(:,0:bs-1); oldX <- V * oldU,
W[pos(X)] <- V * U(:,0:bs-1), W[pos(oldX)] <- V * oldU */
if (d->X != d->V) {
/*for (i=0; i<d->n_calcPairs; i++) {
ierr = VecSet(d->X[i], 0.0); CHKERRQ(ierr);
ierr = VecMAXPY(d->X[i], d->size_V, data->U+d->size_V*i, d->V);
CHKERRQ(ierr);
}*/
ierr = SlepcUpdateVectorsZ(d->X, d->V, d->size_V, data->U, d->size_V,
d->size_V, d->n_calcPairs); CHKERRQ(ierr);
/*for (i=0; i<addToU; i++) {
ierr = VecSet(d->oldX[i], 0.0); CHKERRQ(ierr);
ierr = VecMAXPY(d->oldX[i], data->old_size_V,
oldU+data->old_size_V*i, d->V); CHKERRQ(ierr);
}*/
ierr = SlepcUpdateVectorsZ(d->oldX, d->V, data->old_size_V, oldU,
data->old_size_V, data->old_size_V, addToU);
CHKERRQ(ierr);
/*for (i=0; i<d->n_calcPairs; i++) {
ierr = VecSet(d->W[(d->X-d->V)+i], 0.0); CHKERRQ(ierr);
ierr = VecMAXPY(d->W[(d->X-d->V)+i], d->size_V, data->U+d->size_V*i, d->W);
CHKERRQ(ierr);
}*/
ierr = SlepcUpdateVectorsZ(d->W+(d->X-d->V), d->W, d->size_V, data->U,
d->size_V, d->size_V, d->n_calcPairs);
CHKERRQ(ierr);
for (i=0; i<addToU; i++) {
ierr = VecSet(d->W[(d->oldX-d->V)+i], 0.0); CHKERRQ(ierr);
ierr = VecMAXPY(d->W[(d->oldX-d->V)+i], data->imm_size_V,
oldU+data->old_size_V*i, d->V); CHKERRQ(ierr);
}
ierr = SlepcUpdateVectorsZ(d->W+(d->oldX-d->V), d->V, data->old_size_V,
oldU, data->old_size_V, data->old_size_V, addToU);
CHKERRQ(ierr);
/* If X will be in V, the update in W will be permanent */
if (d->X != d->D) {
d->size_W = d->n_calcPairs + (d->X - d->V);
if (d->X + d->n_calcPairs == d->oldX)
d->size_W+= addToU;
}
data->imm_size_V = data->old_size_V = d->size_V;
/* If V starts at X, and X and oldX are contiguous */
} else if ((addToU == 0) || (d->X+d->n_calcPairs == d->oldX)) {
/* U <- [ U(:,0:bs-1) oldU(0:(size_V-1),:) ] */
for (i=0; i<addToU; i++) {
ierr = PetscMemcpy(data->U+d->size_V*(i+d->n_calcPairs),
oldU+data->old_size_V*i,
sizeof(PetscScalar)*PetscMin(data->old_size_V, d->size_V));
CHKERRQ(ierr);
for (j=data->old_size_V; j<d->size_V; j++)
data->U[d->size_V*(i+d->n_calcPairs)+j] = 0.0;
}
if (d->size_V != d->size_W) {
SETERRQ(1, "This implementation of calcPairs not support this state of "
"Dashboard!");
}
// TODO: funny bug: I cannot update the oldX columns of W
ierr = SlepcUpdateVectors(d->size_V, d->V, 0, d->n_calcPairs+addToU,
data->U, d->size_V, PETSC_FALSE); CHKERRQ(ierr);
ierr = SlepcUpdateVectors(d->size_V, d->W, 0, d->n_calcPairs/*+addToU*/,
data->U, d->size_V, PETSC_FALSE); CHKERRQ(ierr);
data->imm_size_V = data->old_size_V = 0;
d->size_W = d->n_calcPairs/*+addToU*/;
/* H <- U^T * (H * U) */
ierr = SlepcDenseMatProd(T, d->ldH, 0.0, 1.0,
d->H, d->ldH, d->size_H, d->size_H, PETSC_FALSE,
data->U, d->size_V, d->size_V, d->n_calcPairs/*+addToU*/, PETSC_FALSE);
CHKERRQ(ierr);
ierr = SlepcDenseMatProd(d->H, d->ldH, 0.0, 1.0,
data->U, d->size_V, d->size_V, d->n_calcPairs/*+addToU*/, PETSC_TRUE,
T, d->ldH, d->size_V, d->n_calcPairs/*+addToU*/, PETSC_FALSE);
CHKERRQ(ierr);
d->size_H = d->n_calcPairs/*+addToU*/;
//d->size_H = 0;
} else {
SETERRQ(1, "This implementation of calcPairs not support this state of "
"Dashboard!");
}
d->size_X = d->n_calcPairs;
d->size_oldX = addToU;
PetscFunctionReturn(0);
}
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_residual_rr"
PetscInt dvd_calcpairs_residual_rr(dvdDashboard *d, PetscInt s, PetscInt e,
Vec *R)
{
PetscErrorCode ierr;
dvdCalcpairs_rr
*data = (dvdCalcpairs_rr*)d->calcPairs_data;
PetscInt i,j;
PetscFunctionBegin;
/* R <- W * U(:,s:e-1) - X(:,s:e-1) * L(s:e-1) */
for (i=s,j=0; i<e; i++,j++) {
if (d->B) {
/* R[j] <- B * X[i] */
ierr = MatMult(d->B, d->X[i], R[j]); CHKERRQ(ierr);
/* R[j] <- R[j] * (- eigv) + W * U[:,i] */
ierr = VecAYPX(R[j], -d->eigr[i], d->W[(d->X-d->V)+i]);
CHKERRQ(ierr);
} else {
/* R[j] <- X[j] * (- eigv) + W * U[:,i] */
ierr = VecWAXPY(R[j], -d->eigr[i], d->X[i], d->W[(d->X-d->V)+i]);
CHKERRQ(ierr);
}
}
PetscFunctionReturn(0);
}
EXTERN_C_END
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "dvd_calcpairs_rr_d"
PetscInt dvd_calcpairs_rr_d(dvdDashboard *d)
{
PetscErrorCode ierr;
dvdCalcpairs_rr
*data = (dvdCalcpairs_rr*)d->calcPairs_data;
PetscFunctionBegin;
/* Restore changes in dvdDashboard */
d->calcPairs_data = data->old_calcPairs_data;
/* Free local data */
ierr = PetscFree(data); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
EXTERN_C_END