/*
SLEPc eigensolver: "davidson"
Some utils
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SLEPc - Scalable Library for Eigenvalue Problem Computations
Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
This file is part of SLEPc.
SLEPc is free software: you can redistribute it and/or modify it under the
terms of version 3 of the GNU Lesser General Public License as published by
the Free Software Foundation.
SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
more details.
You should have received a copy of the GNU Lesser General Public License
along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
#include "davidson.h"
PetscErrorCode dvd_static_precond_PC_0(dvdDashboard *d, PetscInt i, Vec x,
Vec Px);
PetscErrorCode dvd_jacobi_precond_0(dvdDashboard *d, PetscInt i, Vec x, Vec Px);
PetscErrorCode dvd_precond_none(dvdDashboard *d, PetscInt i, Vec x, Vec Px);
PetscErrorCode dvd_improvex_precond_d(dvdDashboard *d);
typedef struct {
PC pc;
} dvdPCWrapper;
/*
Create a static preconditioner from a PC
*/
#undef __FUNCT__
#define __FUNCT__ "dvd_static_precond_PC"
PetscErrorCode dvd_static_precond_PC(dvdDashboard *d, dvdBlackboard *b, PC pc)
{
PetscErrorCode ierr;
dvdPCWrapper *dvdpc;
Mat P;
MatStructure str;
PetscFunctionBegin;
/* Setup the step */
if (b->state >= DVD_STATE_CONF) {
/* If the preconditioner is valid */
if (pc) {
ierr = PetscMalloc(sizeof(dvdPCWrapper), &dvdpc); CHKERRQ(ierr);
dvdpc->pc = pc;
ierr = PetscObjectReference((PetscObject)pc); CHKERRQ(ierr);
d->improvex_precond_data = dvdpc;
d->improvex_precond = dvd_static_precond_PC_0;
/* PC saves the matrix associated with the linear system, and it has to
be initialize to a valid matrix */
ierr = PCGetOperators(pc, PETSC_NULL, &P, &str); CHKERRQ(ierr);
ierr = PCSetOperators(pc, P, P, str); CHKERRQ(ierr);
ierr = PCSetUp(pc); CHKERRQ(ierr);
DVD_FL_ADD(d->destroyList, dvd_improvex_precond_d);
/* Else, use no preconditioner */
} else
d->improvex_precond = dvd_precond_none;
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "dvd_improvex_precond_d"
PetscErrorCode dvd_improvex_precond_d(dvdDashboard *d)
{
PetscErrorCode ierr;
dvdPCWrapper *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;
PetscFunctionBegin;
/* Free local data */
if (dvdpc->pc) { ierr = PCDestroy(dvdpc->pc); CHKERRQ(ierr); }
ierr = PetscFree(d->improvex_precond_data); CHKERRQ(ierr);
d->improvex_precond_data = PETSC_NULL;
PetscFunctionReturn(0);
}
PetscErrorCode dvd_static_precond_PC_0(dvdDashboard *d, PetscInt i, Vec x,
Vec Px)
{
PetscErrorCode ierr;
dvdPCWrapper *dvdpc = (dvdPCWrapper*)d->improvex_precond_data;
PetscFunctionBegin;
ierr = PCApply(dvdpc->pc, x, Px); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
typedef struct {
Vec diagA, diagB;
} dvdJacobiPrecond;
/*
Create the Jacobi preconditioner for Generalized Eigenproblems
*/
PetscErrorCode dvd_jacobi_precond(dvdDashboard *d, dvdBlackboard *b)
{
PetscErrorCode ierr;
dvdJacobiPrecond
*dvdjp;
PetscTruth t;
PetscFunctionBegin;
/* Check if the problem matrices support GetDiagonal */
ierr = MatHasOperation(d->A, MATOP_GET_DIAGONAL, &t); CHKERRQ(ierr);
if ((t == PETSC_TRUE) && d->B) {
ierr = MatHasOperation(d->B, MATOP_GET_DIAGONAL, &t); CHKERRQ(ierr);
}
/* Setting configuration constrains */
b->own_vecs+= t==PETSC_TRUE?( (d->B == 0)?1:2 ) : 0;
/* Setup the step */
if (b->state >= DVD_STATE_CONF) {
if (t == PETSC_TRUE) {
ierr = PetscMalloc(sizeof(dvdJacobiPrecond), &dvdjp); CHKERRQ(ierr);
dvdjp->diagA = *b->free_vecs; b->free_vecs++;
ierr = MatGetDiagonal(d->A,dvdjp->diagA); CHKERRQ(ierr);
if (d->B) {
dvdjp->diagB = *b->free_vecs; b->free_vecs++;
ierr = MatGetDiagonal(d->B, dvdjp->diagB); CHKERRQ(ierr);
} else
dvdjp->diagB = 0;
d->improvex_precond_data = dvdjp;
d->improvex_precond = dvd_jacobi_precond_0;
DVD_FL_ADD(d->destroyList, dvd_improvex_precond_d);
/* Else, use no preconditioner */
} else
d->improvex_precond = dvd_precond_none;
}
PetscFunctionReturn(0);
}
PetscErrorCode dvd_jacobi_precond_0(dvdDashboard *d, PetscInt i, Vec x, Vec Px)
{
PetscErrorCode ierr;
dvdJacobiPrecond
*dvdjp = (dvdJacobiPrecond*)d->improvex_precond_data;
PetscFunctionBegin;
/* Compute inv(D - eig)*x */
if (dvdjp->diagB == 0) {
/* Px <- diagA - l */
ierr = VecCopy(dvdjp->diagA, Px); CHKERRQ(ierr);
ierr = VecShift(Px, -d->eigr[i]); CHKERRQ(ierr);
} else {
/* Px <- diagA - l*diagB */
ierr = VecWAXPY(Px, -d->eigr[i], dvdjp->diagB, dvdjp->diagA);
CHKERRQ(ierr);
}
/* Px(i) <- x/Px(i) */
ierr = VecPointwiseDivide(Px, x, Px); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
/*
Create a trivial preconditioner
*/
PetscErrorCode dvd_precond_none(dvdDashboard *d, PetscInt i, Vec x, Vec Px)
{
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = VecCopy(x, Px); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
/*
Use of PETSc profiler functions
*/
/* Define stages */
#define DVD_STAGE_INITV 0
#define DVD_STAGE_NEWITER 1
#define DVD_STAGE_CALCPAIRS 2
#define DVD_STAGE_IMPROVEX 3
#define DVD_STAGE_UPDATEV 4
#define DVD_STAGE_ORTHV 5
PetscErrorCode dvd_profiler_d(dvdDashboard *d);
typedef struct {
PetscErrorCode (*old_initV)(struct _dvdDashboard*);
PetscErrorCode (*old_calcPairs)(struct _dvdDashboard*);
PetscErrorCode (*old_improveX)(struct _dvdDashboard*, Vec *D,
PetscInt max_size_D, PetscInt r_s,
PetscInt r_e, PetscInt *size_D);
PetscErrorCode (*old_updateV)(struct _dvdDashboard*);
PetscErrorCode (*old_orthV)(struct _dvdDashboard*);
} DvdProfiler;
PetscLogStage stages[6] = {0,0,0,0,0,0};
/*** Other things ****/
#undef __FUNCT__
#define __FUNCT__ "dvd_prof_init"
PetscErrorCode dvd_prof_init() {
PetscErrorCode ierr;
PetscFunctionBegin;
if (!stages[0]) {
ierr = PetscLogStageRegister("Dvd_step_initV", &stages[DVD_STAGE_INITV]);
CHKERRQ(ierr);
ierr = PetscLogStageRegister("Dvd_step_calcPairs",
&stages[DVD_STAGE_CALCPAIRS]); CHKERRQ(ierr);
ierr = PetscLogStageRegister("Dvd_step_improveX",
&stages[DVD_STAGE_IMPROVEX]); CHKERRQ(ierr);
ierr = PetscLogStageRegister("Dvd_step_updateV",
&stages[DVD_STAGE_UPDATEV]); CHKERRQ(ierr);
ierr = PetscLogStageRegister("Dvd_step_orthV",
&stages[DVD_STAGE_ORTHV]); CHKERRQ(ierr);
}
ierr = dvd_blas_prof_init(); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
PetscErrorCode dvd_initV_prof(dvdDashboard* d) {
DvdProfiler *p = (DvdProfiler*)d->prof_data;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscLogStagePush(stages[DVD_STAGE_INITV]);
ierr = p->old_initV(d); CHKERRQ(ierr);
PetscLogStagePop();
PetscFunctionReturn(0);
}
PetscErrorCode dvd_calcPairs_prof(dvdDashboard* d) {
DvdProfiler *p = (DvdProfiler*)d->prof_data;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscLogStagePush(stages[DVD_STAGE_CALCPAIRS]);
ierr = p->old_calcPairs(d); CHKERRQ(ierr);
PetscLogStagePop();
PetscFunctionReturn(0);
}
PetscErrorCode dvd_improveX_prof(dvdDashboard* d, Vec *D, PetscInt max_size_D,
PetscInt r_s, PetscInt r_e, PetscInt *size_D) {
DvdProfiler *p = (DvdProfiler*)d->prof_data;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscLogStagePush(stages[DVD_STAGE_IMPROVEX]);
ierr = p->old_improveX(d, D, max_size_D, r_s, r_e, size_D); CHKERRQ(ierr);
PetscLogStagePop();
PetscFunctionReturn(0);
}
PetscErrorCode dvd_updateV_prof(dvdDashboard *d) {
DvdProfiler *p = (DvdProfiler*)d->prof_data;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscLogStagePush(stages[DVD_STAGE_UPDATEV]);
ierr = p->old_updateV(d); CHKERRQ(ierr);
PetscLogStagePop();
PetscFunctionReturn(0);
}
PetscErrorCode dvd_orthV_prof(dvdDashboard *d) {
DvdProfiler *p = (DvdProfiler*)d->prof_data;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscLogStagePush(stages[DVD_STAGE_ORTHV]);
ierr = p->old_orthV(d); CHKERRQ(ierr);
PetscLogStagePop();
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "dvd_profiler"
PetscErrorCode dvd_profiler(dvdDashboard *d, dvdBlackboard *b)
{
PetscErrorCode ierr;
DvdProfiler *p;
PetscFunctionBegin;
/* Setup the step */
if (b->state >= DVD_STATE_CONF) {
if (d->prof_data) {
ierr = PetscFree(d->prof_data); CHKERRQ(ierr);
}
ierr = PetscMalloc(sizeof(DvdProfiler), &p); CHKERRQ(ierr);
d->prof_data = p;
p->old_initV = d->initV; d->initV = dvd_initV_prof;
p->old_calcPairs = d->calcPairs; d->calcPairs = dvd_calcPairs_prof;
p->old_improveX = d->improveX; d->improveX = dvd_improveX_prof;
p->old_updateV = d->updateV; d->updateV = dvd_updateV_prof;
DVD_FL_ADD(d->destroyList, dvd_profiler_d);
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "dvd_profiler_d"
PetscErrorCode dvd_profiler_d(dvdDashboard *d)
{
PetscErrorCode ierr;
DvdProfiler *p = (DvdProfiler*)d->prof_data;
PetscFunctionBegin;
/* Free local data */
ierr = PetscFree(p); CHKERRQ(ierr);
d->prof_data = PETSC_NULL;
PetscFunctionReturn(0);
}
/*
Configure the harmonics.
switch(mode) {
DVD_HARM_RR: harmonic RR
DVD_HARM_RRR: relative harmonic RR
DVD_HARM_REIGS: rightmost eigenvalues
DVD_HARM_LEIGS: largest eigenvalues
}
fixedTarged, if true use the target instead of the best eigenvalue
target, the fixed target to be used
*/
typedef struct {
PetscScalar
Wa, Wb, /* span{W} = span{Wa*AV - Wb*BV} */
Pa, Pb; /* H=W'*(Pa*AV - Pb*BV), G=W'*(Wa*AV - Wb*BV) */
PetscTruth
withTarget;
HarmType_t
mode;
/* old values of eps */
EPSWhich
old_which;
PetscErrorCode
(*old_which_func)(EPS,PetscScalar,PetscScalar,PetscScalar,PetscScalar,
PetscInt*,void*);
void
*old_which_ctx;
} dvdHarmonic;
PetscErrorCode dvd_harm_start(dvdDashboard *d);
PetscErrorCode dvd_harm_end(dvdDashboard *d);
PetscErrorCode dvd_harm_d(dvdDashboard *d);
PetscErrorCode dvd_harm_transf(dvdHarmonic *dvdh, PetscScalar t);
PetscErrorCode dvd_harm_updateW(dvdDashboard *d);
PetscErrorCode dvd_harm_proj(dvdDashboard *d);
PetscErrorCode dvd_harm_sort(EPS eps, PetscScalar ar, PetscScalar ai,
PetscScalar br, PetscScalar bi, PetscInt *r,
void *ctx);
PetscErrorCode dvd_harm_eigs_trans(dvdDashboard *d);
PetscErrorCode dvd_harm_conf(dvdDashboard *d, dvdBlackboard *b,
HarmType_t mode, PetscTruth fixedTarget,
PetscScalar t)
{
PetscErrorCode ierr;
dvdHarmonic *dvdh;
PetscFunctionBegin;
/* Set the problem to GNHEP */
// TODO: d->G maybe is upper triangular due to biorthogonality of V and W
d->sEP = d->sA = d->sB = 0;
/* Setup the step */
if (b->state >= DVD_STATE_CONF) {
ierr = PetscMalloc(sizeof(dvdHarmonic), &dvdh); CHKERRQ(ierr);
dvdh->withTarget = fixedTarget;
dvdh->mode = mode;
if (fixedTarget == PETSC_TRUE) dvd_harm_transf(dvdh, t);
d->calcpairs_W_data = dvdh;
d->calcpairs_W = dvd_harm_updateW;
d->calcpairs_proj_trans = dvd_harm_proj;
d->calcpairs_eigs_trans = dvd_harm_eigs_trans;
DVD_FL_ADD(d->startList, dvd_harm_start);
DVD_FL_ADD(d->endList, dvd_harm_end);
DVD_FL_ADD(d->destroyList, dvd_harm_d);
}
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "dvd_harm_d"
PetscErrorCode dvd_harm_d(dvdDashboard *d)
{
PetscErrorCode ierr;
PetscFunctionBegin;
/* Free local data */
ierr = PetscFree(d->calcpairs_W_data); CHKERRQ(ierr);
d->calcpairs_W_data = PETSC_NULL;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "dvd_harm_start"
PetscErrorCode dvd_harm_start(dvdDashboard *d)
{
dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
PetscFunctionBegin;
/* Overload the eigenpairs selection routine */
data->old_which = d->eps->which;
data->old_which_func = d->eps->which_func;
data->old_which_ctx = d->eps->which_ctx;
d->eps->which = EPS_WHICH_USER;
d->eps->which_func = dvd_harm_sort;
d->eps->which_ctx = data;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "dvd_harm_end"
PetscErrorCode dvd_harm_end(dvdDashboard *d)
{
dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
PetscFunctionBegin;
/* Restore the eigenpairs selection routine */
d->eps->which = data->old_which;
d->eps->which_func = data->old_which_func;
d->eps->which_ctx = data->old_which_ctx;
PetscFunctionReturn(0);
}
PetscErrorCode dvd_harm_transf(dvdHarmonic *dvdh, PetscScalar t)
{
PetscFunctionBegin;
switch(dvdh->mode) {
case DVD_HARM_RR: /* harmonic RR */
dvdh->Wa = 1.0; dvdh->Wb = t; dvdh->Pa = 0.0; dvdh->Pb = -1.0; break;
case DVD_HARM_RRR: /* relative harmonic RR */
dvdh->Wa = 1.0; dvdh->Wb = t; dvdh->Pa = 1.0; dvdh->Pb = 0.0; break;
case DVD_HARM_REIGS: /* rightmost eigenvalues */
dvdh->Wa = 1.0; dvdh->Wb = t; dvdh->Pa = 1.0; dvdh->Pb = -PetscConj(t);
break;
case DVD_HARM_LEIGS: /* largest eigenvalues */
dvdh->Wa = 0.0; dvdh->Wb = 1.0; dvdh->Pa = 1.0; dvdh->Pb = 0.0; break;
case DVD_HARM_NONE:
default:
SETERRQ(1, "The harmonic type is not supported!");
}
/* Check the transformation does not change the sign of the imaginary part */
#if !defined(PETSC_USE_COMPLEX)
if (dvdh->Pb*dvdh->Wa - dvdh->Wb*dvdh->Pa < 0.0)
dvdh->Pa*= -1.0, dvdh->Pb*= -1.0;
#endif
PetscFunctionReturn(0);
}
PetscErrorCode dvd_harm_updateW(dvdDashboard *d)
{
dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
PetscErrorCode ierr;
PetscInt i;
PetscFunctionBegin;
/* Update the target if it is necessary */
if (data->withTarget == PETSC_FALSE) dvd_harm_transf(data, d->eigr[0]);
for(i=d->V_new_s; i<d->V_new_e; i++) {
/* W(i) <- Wa*AV(i) - Wb*BV(i) */
ierr = VecCopy(d->AV[i], d->W[i]); CHKERRQ(ierr);
ierr = VecAXPBY(d->W[i], -data->Wb, data->Wa, (d->BV?d->BV:d->V)[i]);
CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
PetscErrorCode dvd_harm_proj(dvdDashboard *d)
{
dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
PetscInt i,j;
PetscFunctionBegin;
if (d->sH != d->sG) {
SETERRQ(1, "Error: Projected matrices H and G must have the same structure!");
PetscFunctionReturn(1);
}
/* [H G] <- [Pa*H - Pb*G, Wa*H - Wb*G] */
if (DVD_ISNOT(d->sH,DVD_MAT_LTRIANG)) /* Upper triangular part */
for(i=d->V_new_s; i<d->V_new_e; i++)
for(j=0; j<=i; j++) {
PetscScalar h = d->H[d->ldH*i+j], g = d->G[d->ldH*i+j];
d->H[d->ldH*i+j] = data->Pa*h - data->Pb*g;
d->G[d->ldH*i+j] = data->Wa*h - data->Wb*g;
}
if (DVD_ISNOT(d->sH,DVD_MAT_UTRIANG)) /* Lower triangular part */
for(i=0; i<d->V_new_e; i++)
for(j=PetscMax(d->V_new_s,i+(DVD_ISNOT(d->sH,DVD_MAT_LTRIANG)?1:0));
j<d->V_new_e; j++) {
PetscScalar h = d->H[d->ldH*i+j], g = d->G[d->ldH*i+j];
d->H[d->ldH*i+j] = data->Pa*h - data->Pb*g;
d->G[d->ldH*i+j] = data->Wa*h - data->Wb*g;
}
PetscFunctionReturn(0);
}
PetscErrorCode dvd_harm_backtrans(dvdHarmonic *data, PetscScalar *ar,
PetscScalar *ai)
{
PetscScalar xr;
#if !defined(PETSC_USE_COMPLEX)
PetscScalar xi, k;
#endif
PetscFunctionBegin;
if(!ar) SETERRQ(1, "The real part has to be present!");
xr = *ar;
#if !defined(PETSC_USE_COMPLEX)
if(!ai) SETERRQ(1, "The imaginary part has to be present!");
xi = *ai;
if (xi != 0.0) {
k = (data->Pa - data->Wa*xr)*(data->Pa - data->Wa*xr) +
data->Wa*data->Wa*xi*xi;
*ar = (data->Pb*data->Pa - (data->Pb*data->Wa + data->Wb*data->Pa)*xr +
data->Wb*data->Wa*(xr*xr + xi*xi))/k;
*ai = (data->Pb*data->Wa - data->Wb*data->Pa)*xi/k;
} else
#endif
*ar = (data->Pb - data->Wb*xr) / (data->Pa - data->Wa*xr);
PetscFunctionReturn(0);
}
PetscErrorCode dvd_harm_sort(EPS eps, PetscScalar ar, PetscScalar ai,
PetscScalar br, PetscScalar bi, PetscInt *r,
void *ctx)
{
dvdHarmonic *data = (dvdHarmonic*)ctx;
PetscErrorCode ierr;
PetscFunctionBegin;
/* Back-transform the harmonic values */
dvd_harm_backtrans(data, &ar, &ai);
dvd_harm_backtrans(data, &br, &bi);
/* Compare values using the user options for the eigenpairs selection */
eps->which = data->old_which;
eps->which_func = data->old_which_func;
eps->which_ctx = data->old_which_ctx;
ierr = EPSCompareEigenvalues(eps, ar, ai, br, bi, r); CHKERRQ(ierr);
/* Restore the eps values */
eps->which = EPS_WHICH_USER;
eps->which_func = dvd_harm_sort;
eps->which_ctx = data;
PetscFunctionReturn(0);
}
PetscErrorCode dvd_harm_eigs_trans(dvdDashboard *d)
{
dvdHarmonic *data = (dvdHarmonic*)d->calcpairs_W_data;
PetscInt i;
PetscFunctionBegin;
for(i=0; i<d->size_H; i++)
dvd_harm_backtrans(data, &d->eigr[i], &d->eigi[i]);
PetscFunctionReturn(0);
}