/* Preconditioned Conjugate-Gradient linear solver */
/* original paper?: MR Hestenes and E Stiefel, "Method of conjugate */
/* gradients for solving linear systems," J. of Res. of the Nat'l */
/* Bureau of Standards, v.49, p.409, 1952 */
/* Copyright John B. Pormann, 21 July 2000, all rights reserved */
#include "CardioWave.h"
static char* RCSID = "$Id: SolverPCG.c,v 1.3 2004/08/24 19:37:38 jpormann Exp $";
static int MaxIter? = 100;
static real Tolerance = 1.0e-6;
static real Tolerance2?;
static vector CG_r = _INIT_VECTOR;
static vector CG_v = _INIT_VECTOR;
static vector CG_z = _INIT_VECTOR;
static rword resources[] = {
{ "maxit", 5001 },
{ "maxiter", 5001 },
{ "solvertol", 5002 },
{ "stol", 5002 },
{ "stolerance", 5002 },
{ NULL, 0 }
};
int InitLinearSolver( char** res ) {
int i;
int cmd;
i = 0;
while( res[i] != NULL ) {
cmd = FindCommand( resources, res[i] );
switch( cmd ) {
case 5001:
MaxIter? = GetIntValue( res[i] );
break;
case 5002:
Tolerance = GetRealValue( res[i] );
break;
}
i++;
}
Tolerance2? = Tolerance*Tolerance;
if( vecalloc(&CG_r,Solver) < 0 ) {
return( -1 );
}
if( vecalloc(&CG_v,Solver) < 0 ) {
return( -2 );
}
if( vecalloc(&CG_z,Solver) < 0 ) {
return( -3 );
}
if( (DebugLevel>0) and (SelfPE==0) ) {
printf("SolverPCG?: maxit = %i, tol = [le]\n",
MaxIter?,Tolerance,Tolerance2?);
if( DebugLevel > 1 ) {
printf("SolverPCG?: RCSID: %s\n",RCSID);
}
}
return( 0 );
}
void ExitLinearSolver( void ) {
vecfree( &CG_r );
vecfree( &CG_v );
vecfree( &CG_z );
return;
}
int Solve( bigmatrix A, vector x, vector b ) {
int i;
real t,d,c,vz;
/* r = b - A*x */
vcopy( b, CG_r );
MatVec( -1.0, A, x, 1.0, CG_r );
/* solve Qv=r for v */
PreconditionSystem( A, CG_v, CG_r );
/* c = <v,r> */
c = innerprod( CG_v, CG_r );
/* if( c < Tolerance2? ) {
return( 0 );
}
*/
/* this loop guarantees that we eventually exit */
for(i=0;i<MaxIter?;i++) {
/* check for near-0 vectors */
d = innerprod( CG_v, CG_v );
if( d < Tolerance2? ) {
break;
}
/* z = A*v */
MatVec( 1.0, A, CG_v, 0.0, CG_z );
/* t = c/<v,z> */
vz = innerprod( CG_v, CG_z );
t = c / vz;
/* x = x + t*v */
vvadd( t, CG_v, x );
/* r = r - t*z */
vvadd( -t, CG_z, CG_r );
/* solve Qz=r for z */
PreconditionSystem( A, CG_z, CG_r );
/* d = <z,r> */
d = innerprod( CG_z, CG_r );
/* when d gets small enough, we're done */
/* : usually look at sqrt(d) */
if( d < Tolerance2? ) {
/* now check <r,r> */
t = innerprod( CG_r, CG_r );
if( t < Tolerance2? ) {
break;
}
}
/* v = z + (d/c)*v */
vvaddb( 1.0, CG_z, (d/c), CG_v );
c = d;
}
if( i>=MaxIter? ) {
if( DebugLevel > 1 ) {
fprintf(stderr,"Solver: maximum iterations exceeded\n");
}
return( -1 );
}
return( i );
}
