/* MLBiCGSTABL.c */
#include "../Iter.h"
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------------------------
purpose -- to solve a unsymmetric real matrix equation
Ax=b
using left preconditioned ML(k)BiCGSTAB method
by M-C Yeung and T. Chan
x -- Initial guess as zeros
A -- Input matrix
Precond -- Front Matrix as the preconditioner
Q -- Starting vectors
b -- Right-hand side
tol -- Convergence tolerance
type -- type of entries
SPOOLES_REAL or SPOOLES_COMPLEX
symmetryflag -- symmetry of the matrix
SPOOLES_SYMMETRIC, SPOOLES_HERMITIAN or SPOOLES_NONSYMMETRIC
nrhs -- number of right hand sides
msglvl -- message level
msgFile -- message file
return value -- error flag
created -- Nov. 28, 1998 Wei-Pai Tang
---------------------------------------------------------------------
*/
int
mlbicgstabl (
int n_matrixSize,
int type,
int symmetryflag,
InpMtx *mtxA,
FrontMtx *Precond,
DenseMtx *mtxQ,
DenseMtx *mtxX,
DenseMtx *mtxB,
int itermax,
double convergetol,
int msglvl,
FILE *msgFile
)
{
Chv *chv, *rootchv ;
ChvManager *chvmanager ;
DenseMtx *mtxD, *mtxG, *mtxU, *mtxW;
DenseMtx *vecGt, *vecR, *vecT, *vecT1;
DenseMtx *vecX, *vecZD, *vecZG, *vecZW ;
double Alpha, Beta, Rho ;
double c[100] ;
double Init_norm, ratio, Res_norm;
double error_trol, m, Rtmp;
double t1, t2, cpus[9] ;
double one[2] = {1.0, 0.0} ;
double zero[2] = {0.0, 0.0} ;
double minusone[2] = {-1.0, 0.0};
double Tiny = 0.1e-28;
int Iter, Imv, neqns, Ik, ii, is;
int stats[6] ;
int return_flag;
neqns = n_matrixSize;
Ik = mtxQ->ncol;
if (Ik > 100){
fprintf(msgFile, "\n\n Fatal Error, \n"
" Too many starting vectors in Q !!") ;
return_flag = -1;
return (return_flag);
};
return_flag = 1;
/*
--------------------
init the vectors in MLBiCGSTABL
--------------------
*/
mtxD = DenseMtx_new() ;
DenseMtx_init(mtxD, type, 0, 0, neqns, Ik, 1, neqns) ;
mtxG = DenseMtx_new() ;
DenseMtx_init(mtxG, type, 0, 0, neqns, Ik, 1, neqns) ;
mtxU = DenseMtx_new() ;
DenseMtx_init(mtxU, type, 0, 0, neqns, 2, 1, neqns) ;
mtxW = DenseMtx_new() ;
DenseMtx_init(mtxW, type, 0, 0, neqns, Ik, 1, neqns) ;
vecGt = DenseMtx_new() ;
DenseMtx_init(vecGt, type, 0, 0, neqns, 1, 1, neqns) ;
vecR = DenseMtx_new() ;
DenseMtx_init(vecR, type, 0, 0, neqns, 1, 1, neqns) ;
vecT = DenseMtx_new() ;
DenseMtx_init(vecT, type, 0, 0, neqns, 1, 1, neqns) ;
vecT1 = DenseMtx_new() ;
DenseMtx_init(vecT1, type, 0, 0, neqns, 1, 1, neqns) ;
vecX = DenseMtx_new() ;
DenseMtx_init(vecX, type, 0, 0, neqns, 1, 1, neqns) ;
vecZD = DenseMtx_new() ;
DenseMtx_init(vecZD, type, 0, 0, neqns, 1, 1, neqns) ;
vecZG = DenseMtx_new() ;
DenseMtx_init(vecZG, type, 0, 0, neqns, 1, 1, neqns) ;
vecZW = DenseMtx_new() ;
DenseMtx_init(vecZW, type, 0, 0, neqns, 1, 1, neqns) ;
/*
c = DV_new();
DV_init(c, Ik, NULL);
*/
for ( ii = 0; ii <100; ii++){
c[ii] = 0;
}
/*
--------------------------
Initialize the iterations
--------------------------
*/
/* ---- Set initial guess as zero ---- */
DenseMtx_zero(vecX) ;
/* ---- If x_0 is not zero ---- */
/*
DenseMtx_colCopy(vecX, 0, mtxX, 0);
*/
/*
InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecX) ;
*/
switch ( symmetryflag ) {
case SPOOLES_SYMMETRIC :
InpMtx_sym_gmmm(mtxA, zero, vecT, one, vecX) ;
break ;
case SPOOLES_HERMITIAN :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
case SPOOLES_NONSYMMETRIC :
InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecX) ;
break ;
default :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
}
DenseMtx_colCopy(vecR, 0, mtxB, 0);
DenseMtx_sub(vecR, vecT) ;
FrontMtx_solve(Precond, vecT, vecR, Precond->manager,
cpus, msglvl, msgFile) ;
DenseMtx_colCopy(vecR, 0, vecT, 0);
Init_norm = DenseMtx_twoNormOfColumn(vecR, 0);
if ( Init_norm == 0.0 ){
Init_norm = 1.0;
};
error_trol = Init_norm * convergetol ;
fprintf(msgFile, "\n MLBiCGSTABL Initial norml: %6.2e ", Init_norm ) ;
fprintf(msgFile, "\n MLBiCGSTABL Conveg. Control: %7.3e ", convergetol ) ;
fprintf(msgFile, "\n MLBiCGSTABL Convergen Control: %7.3e ",error_trol ) ;
DenseMtx_zero(mtxG) ;
DenseMtx_zero(mtxD) ;
DenseMtx_zero(mtxW) ;
Iter = 0;
Imv = 0;
DenseMtx_colCopy (mtxG, Ik-1, vecR, 0);
/*
------------------------------
MLBiCGSTABL Iteration start
------------------------------
*/
MARKTIME(t1) ;
while ( Iter <= itermax ){
Iter++;
/*
g_tld = G(:,k);
W(:,k) = U\(L\( A*g_tld));
*/
DenseMtx_colCopy (vecGt, 0, mtxG, Ik-1);
/*
InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecGt) ;
*/
switch ( symmetryflag ) {
case SPOOLES_SYMMETRIC :
InpMtx_sym_gmmm(mtxA, zero, vecT, one, vecGt) ;
break ;
case SPOOLES_HERMITIAN :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
case SPOOLES_NONSYMMETRIC :
InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecGt) ;
break ;
default :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
}
FrontMtx_solveOneColumn(Precond, mtxW, Ik-1, vecT, 0,
Precond->manager, cpus, msglvl, msgFile) ;
Imv++;
/* c[Ik] = *DenseMtx_colDotProduct (mtxQ, 0, mtxG, Ik-1); */
DenseMtx_colDotProduct (mtxQ, 0, mtxG, Ik-1, c+Ik);
if (c[Ik] == 0){
fprintf(msgFile, "\n\n Fatal Error, \n"
" MLBiCGSTABL Breakdown, c[k] = 0 !!") ;
return_flag = -1;
goto end;
};
/* Alpha = Q(:,1)'*r/c(k); ; */
DenseMtx_colDotProduct (mtxQ,0, vecR, 0, &Rtmp);
Alpha = Rtmp/c[Ik];
DenseMtx_colCopy (mtxU, 0, vecR, 0);
Rtmp = -Alpha;
DenseMtx_colGenAxpy (one, mtxU, 0, &Rtmp, mtxW, Ik-1);
/*
u = r - alpha*W(:,k);
u_tld = u;
temp = U\(L\(A*u_tld));
rho =temp'*temp;
*/
DenseMtx_colCopy (mtxU, 1, mtxU, 0);
DenseMtx_colCopy (vecT, 0, mtxU, 0);
/*
InpMtx_nonsym_gmmm(mtxA, zero, vecT1, one, vecT) ;
*/
switch ( symmetryflag ) {
case SPOOLES_SYMMETRIC :
InpMtx_sym_gmmm(mtxA, zero, vecT1, one, vecT) ;
break ;
case SPOOLES_HERMITIAN :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
case SPOOLES_NONSYMMETRIC :
InpMtx_nonsym_gmmm(mtxA, zero, vecT1, one, vecT) ;
break ;
default :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
}
FrontMtx_solveOneColumn(Precond, vecT, 0, vecT1, 0,
Precond->manager, cpus, msglvl, msgFile) ;
Imv++;
/* Rho = temp'*temp; */
DenseMtx_colDotProduct (vecT, 0, vecT, 0, &Rho);
if (Rho == 0){
fprintf(msgFile, "\n\n Fatal Error, \n"
" MLBiCGSTABL Breakdown, Rho = 0 !!") ;
return_flag = -1;
goto end;
};
/* Rho = -(*DenseMtx_colDotProduct (mtxU, 0, vecT, 0))/Rho; */
DenseMtx_colDotProduct (mtxU, 0, vecT, 0, &Rtmp);
Rho = -Rtmp/Rho;
DenseMtx_colGenAxpy(one, vecX, 0, &Alpha, vecGt, 0);
Rtmp = -Rho;
DenseMtx_colGenAxpy (one, vecX, 0, &Rtmp, mtxU, 1) ;
DenseMtx_colCopy (vecR, 0, mtxU, 0);
DenseMtx_colGenAxpy (one, vecR, 0, &Rho, vecT, 0) ;
/* Iter++; */
/*
----------------
Convergence Test
---------------
*/
Rtmp = DenseMtx_twoNormOfColumn ( vecR, 0);
ratio = Rtmp/Init_norm;
fprintf(msgFile, "\n\n At iteration %d"
" the convergence ratio is %12.4e"
"\n Residual norm is %6.2ee",
Imv, ratio, Rtmp) ;
fflush(msgFile) ;
if ( Rtmp <= error_trol ) {
fprintf(msgFile, "\n # iterations = %d", Imv) ;
return_flag = Imv;
goto end;
};
for (ii = 1; ii < Ik+1; ii++ ){
if (Iter > itermax ){
fprintf(msgFile, "\n # iterations = %d", Imv) ;
fprintf(msgFile, "\n\n MLBiCGSTABL did not Converge !") ;
return_flag = Imv;
goto end;
};
DenseMtx_colCopy (vecZD,0, mtxU, 0);
DenseMtx_colCopy (vecZG,0, vecR, 0);
DenseMtx_zero(vecZW) ;
if ( Iter > 1 ){
for ( is = ii ; is < Ik ; is++ ){
/* Beta =-( *(DenseMtx_colDotProduct(mtxQ, is, vecZD, 0)))/c[is]; */
DenseMtx_colDotProduct(mtxQ, is, vecZD, 0, &Rtmp);
Beta = -Rtmp/c[is];
DenseMtx_colGenAxpy (one, vecZD, 0, &Beta, mtxD, is-1);
DenseMtx_colGenAxpy (one, vecZG, 0, &Beta, mtxG, is-1);
DenseMtx_colGenAxpy (one, vecZW, 0, &Beta, mtxW, is-1);
};
};
Beta = Rho * c[Ik];
if (Beta == 0){
fprintf(msgFile, "\n\n Fatal Error, \n"
" MLBiCGSTABL Breakdown, Beta = 0 !!") ;
return_flag = -1;
goto end;
};
/*
Beta = - Q(:,1)'* (r + Rho* zw )/ Beta;
*/
DenseMtx_colCopy (vecT, 0, vecR, 0);
DenseMtx_colGenAxpy (one, vecT, 0, &Rho, vecZW, 0);
/* Beta = - (*DenseMtx_colDotProduct(mtxQ, 0, vecT, 0))/Beta; */
DenseMtx_colDotProduct(mtxQ, 0, vecT, 0, &Rtmp);
Beta = - Rtmp/Beta;
/*
zg = zg + beta*G(:,k);
zw = rho*(zw + beta*W(:,k));
zd = r + zw;
*/
DenseMtx_colGenAxpy (one, vecZG, 0, &Beta, mtxG, Ik-1);
Rtmp = Rho*Beta;
DenseMtx_colGenAxpy (&Rho, vecZW, 0, &Rtmp, mtxW, Ik-1);
DenseMtx_colCopy (vecZD, 0, vecR, 0);
DenseMtx_colGenAxpy (one, vecZD, 0, one, vecZW, 0);
/*
for s = 1:i-1
beta = -Q(:,s+1)'*zd/c(s);
zd = zd + beta*D(:,s);
zg = zg + beta*G(:,s);
end
*/
for ( is = 1; is < ii - 1; is ++){
/* Beta = -(*DenseMtx_colDotProduct(mtxQ, is, vecZD, 0))/c[is] ; */
DenseMtx_colDotProduct(mtxQ, is, vecZD, 0, &Rtmp);
Beta = - Rtmp/c[is];
DenseMtx_colGenAxpy (one, vecZD, 0, &Beta, mtxD, is-1);
DenseMtx_colGenAxpy (one, vecZG, 0, &Beta, mtxG, is-1);
};
/*
D(:,i) = zd - u;
G(:,i) = zg + zw;
*/
DenseMtx_colCopy (mtxD, ii-1, vecZD, 0);
DenseMtx_colGenAxpy (one, mtxD, ii-1, minusone, mtxU, 0);
DenseMtx_colCopy (mtxG, ii-1, vecZG, 0);
DenseMtx_colGenAxpy (one, mtxG, ii-1, one, vecZW, 0);
/*
if i < k
c(i) = Q(:,i+1)'*D(:,i);
*/
if ( ii < Ik ){
/* c[ii] = *DenseMtx_colDotProduct(mtxQ, ii, mtxD, ii-1); */
DenseMtx_colDotProduct(mtxQ, ii, mtxD, ii-1, c+ii);
/*
If breakdown ?
*/
if (c[ii] == 0){
fprintf(msgFile, "\n\n Fatal Error, \n"
" MLBiCGSTABL Breakdown, c[ii] = 0 !!") ;
return_flag = -1;
goto end;
};
/*
alpha = Q(:,i+1)'*u/c(i);
u = u - alpha*D(:,i);
g_tld = G(:,i);
*/
DenseMtx_colDotProduct(mtxQ, ii, mtxU, 0, &Rtmp);
Alpha = Rtmp/c[ii];
Rtmp = -Alpha;
DenseMtx_colGenAxpy (one, mtxU, 0, &Rtmp, mtxD, ii-1);
DenseMtx_colCopy (vecGt, 0, mtxG, ii-1);
/*
x = x + rho*alpha*g_tld;
W(:,i) = U\(L\(A*g_tld));
r = r - rho*alpha*W(:,i);
*/
Rtmp = Rho * Alpha;
DenseMtx_colGenAxpy (one, vecX, 0, &Rtmp, vecGt, 0);
/*
InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecGt) ;
*/
switch ( symmetryflag ) {
case SPOOLES_SYMMETRIC :
InpMtx_sym_gmmm(mtxA, zero, vecT, one, vecGt) ;
break ;
case SPOOLES_HERMITIAN :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
case SPOOLES_NONSYMMETRIC :
InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, vecGt) ;
break ;
default :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
}
FrontMtx_solveOneColumn(Precond, mtxW, ii-1, vecT, 0,
Precond->manager, cpus, msglvl, msgFile) ;
Imv++;
Rtmp = -Rtmp;
DenseMtx_colGenAxpy (one, vecR, 0, &Rtmp, mtxW, ii-1);
/*
----------------
Convergence Test
---------------
*/
Rtmp = DenseMtx_twoNormOfColumn ( vecR, 0);
if ( Rtmp <= error_trol ) {
fprintf(msgFile, "\n # iterations = %d", Imv) ;
return_flag = Imv;
goto end;
};
};
};
}
/* End of while loop */
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU : Total iteration time is : %8.3f ", t2 - t1) ;
fprintf(msgFile, "\n # iterations = %d", Imv) ;
fprintf(msgFile, "\n\n MLBiCGSTABL did not Converge !") ;
DenseMtx_colCopy(mtxX, 0, vecX, 0);
/*
------------------------
free the working storage
------------------------
*/
end:
MARKTIME(t2) ;
fprintf(msgFile, "\n CPU : Total iteration time is : %8.3f ", t2 - t1) ;
DenseMtx_colCopy(mtxX, 0, vecX, 0);
/*
InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, mtxX) ;
*/
switch ( symmetryflag ) {
case SPOOLES_SYMMETRIC :
InpMtx_sym_gmmm(mtxA, zero, vecT, one, mtxX) ;
break ;
case SPOOLES_HERMITIAN :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
case SPOOLES_NONSYMMETRIC :
InpMtx_nonsym_gmmm(mtxA, zero, vecT, one, mtxX) ;
break ;
default :
fprintf(msgFile, "\n MLBiCGSTABL Matrix type wrong");
fprintf(msgFile, "\n Fatal error");
return_flag = -1;
goto end;
}
DenseMtx_sub(vecT, mtxB) ;
Rtmp = DenseMtx_twoNormOfColumn(vecT, 0);
fprintf(msgFile, "\n MLBiCGSTABL True Residual norm: %6.2e ",
Rtmp) ;
fprintf(msgFile, "\n\n after MLBiCGSTABL") ;
DenseMtx_free(mtxD) ;
DenseMtx_free(mtxG) ;
DenseMtx_free(mtxU) ;
DenseMtx_free(mtxW) ;
DenseMtx_free(vecGt) ;
DenseMtx_free(vecR) ;
DenseMtx_free(vecT) ;
DenseMtx_free(vecT1) ;
DenseMtx_free(vecX) ;
DenseMtx_free(vecZD) ;
DenseMtx_free(vecZG) ;
DenseMtx_free(vecZW) ;
fprintf(msgFile, "\n") ;
return(return_flag) ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1