/* Jacobi (diagonal) preconditioner */
/* Copyright John B. Pormann, 21 July 2000, all rights reserved */
#include "CardioWave.h"
static char* RCSID = "$Id: PrecondJac.c,v 1.4 2005/05/18 21:29:22 john Exp $";
static vector InvDiag? = _INIT_VECTOR;
static real* Stn_InvDiag = NULL;
int InitPreconditioner( char** res ) {
if( (DebugLevel>0) and (SelfPE==0) ) {
printf("Preconditioner: Jacobi\n");
if( DebugLevel > 1 ) {
printf("PrecondJac?: RCSID: %s\n",RCSID);
}
}
return( 0 );
}
void ExitPreconditioner( void ) {
vecfree( &InvDiag? );
return;
}
void ComputeFactors( bigmatrix A ) {
int i,b,ofs;
int nblks = A.NumBlocks?;
real* av;
real* xv;
sparse* S;
/* A.Int may be NULL, so we'll look at A.Ext for matrix type */
if( A.MatrixInfo?[0].type == _STENCIL ) {
Stn_InvDiag = (real*)malloc( nblks*sizeof(real) );
for(b=0;b<nblks;b++) {
/* find dimensionality of system */
S = &A.MatrixInfo?[b+b*nblks];
switch( S->jcoef[0] ) {
case 0: Stn_InvDiag[b] = 1.0/S->coef[0];
break;
case 1: Stn_InvDiag[b] = 1.0/S->coef[1];
break;
case 2: Stn_InvDiag[b] = 1.0/S->coef[4];
break;
case 3: Stn_InvDiag[b] = 1.0/S->coef[13];
break;
}
}
} else {
vecalloc( &InvDiag?, Solver );
/* go thru each block in turn */
/* : but put all 1/Aii terms into 1 big vector */
/* for better performance during Prec.Solve */
xv = InvDiag?.data;
for(b=0;b<nblks;b++) {
av = A.MatrixInfo?[b+b*nblks].coef;
ofs = A.BlockSplit?[b];
for(i=0;i<A.MatrixInfo?[b+b*nblks].rows;i++) {
xv[i+ofs] = 1.0/av[i];
}
}
}
return;
}
int PreconditionSystem( bigmatrix A, vector x, vector b ) {
int i,j,k,nd,ne;
int ofsI,ofsR;
int nblks = A.NumBlocks?;
real val;
vector xtmp,btmp;
/* x = b/diag(A) */
if( A.MatrixInfo?[0].type == _STENCIL ) {
for(i=0;i<nblks;i++) {
xtmp.size = A.BlockSplit?[i+1] - A.BlockSplit?[i];
xtmp.dtype = A.MatrixInfo?[i+i*nblks].dtype;
xtmp.data = x.data + A.BlockSplit?[i];
btmp.size = A.BlockSplit?[i+1] - A.BlockSplit?[i];
btmp.dtype = A.MatrixInfo?[i+i*nblks].dtype;
btmp.data = b.data + A.BlockSplit?[i];
vvaddb( Stn_InvDiag[i], btmp, 0.0, xtmp );
/* take care of the DBCs? */
nd = A.MatrixInfo?[i+i*nblks].jcoef[0];
switch( nd ) {
case 3: ofsI = 5; /* where do DBCs? start? */
ofsR = 27; /* (I=jcoef, R=coef offsets) */
break;
case 2: ofsI = 4;
ofsR = 9;
break;
case 1: ofsI = 3;
ofsR = 3;
break;
case 0: ofsI = 2;
ofsR = 1;
break;
}
ne = A.MatrixInfo?[i+i*nblks].jcoef[ofsI];
for(j=0;j<ne;j++) {
k = A.MatrixInfo?[i+i*nblks].jcoef[ofsI+1+j];
/* xtmp.data[k] = A.MatrixInfo?[i+i*nblks].coef[ofsR+j]; */
xtmp.data[k] = btmp.data[k];
}
}
} else {
vvmpyb( 1.0, b, InvDiag?, 0.0, x );
}
return( 0 );
}
