diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index d36ba1f15b63244666ee1fdc9ae93a904dae4073..c9163e823dfe730ce1833efb2339bfe330a70750 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -20,6 +20,8 @@ SET( HEADERS
LinAlg/Solvers/TriangularSolve.h
LinAlg/Solvers/IterativeLinearSolver.h
LinAlg/Solvers/solver.h
+ LinAlg/Solvers/CG.h
+ LinAlg/Solvers/BiCGStab.h
LinAlg/Sparse/amuxCRS.h
LinAlg/Sparse/CRSMatrix.h
LinAlg/Sparse/CRSMatrixPThreads.h
@@ -36,6 +38,7 @@ SET( SOURCES
LinkedTriangle.cpp
LinAlg/Sparse/amuxCRS.cpp
LinAlg/Solvers/CG.cpp
+ LinAlg/Solvers/BiCGStab.cpp
LinAlg/Solvers/GaussAlgorithm.cpp
LinAlg/Solvers/TriangularSolve.cpp
LinAlg/Preconditioner/generateDiagPrecond.cpp
diff --git a/MathLib/LinAlg/Solvers/BiCGStab.cpp b/MathLib/LinAlg/Solvers/BiCGStab.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bfd57b11377e20433343510affc282ab4e533ad3
--- /dev/null
+++ b/MathLib/LinAlg/Solvers/BiCGStab.cpp
@@ -0,0 +1,135 @@
+/*
+ * BiCGStab.cpp
+ *
+ * Created on: Oct 4, 2011
+ * Author: TF
+ */
+
+#include "BiCGStab.h"
+
+#include "MathTools.h"
+#include "blas.h"
+
+namespace MathLib {
+
+unsigned BiCGStab(CRSMatrix<double> const& A, double* const b, double* const x,
+ double& eps, unsigned& nsteps)
+{
+ unsigned N(A.getNRows());
+ double *v (new double[8* N]);
+ double *p (v + N);
+ double *phat (p + N);
+ double *s (phat + N);
+ double *shat (s + N);
+ double *t (shat + N);
+ double *r (t + N);
+ double *r0 (r + N);
+ double resid;
+
+ // normb = |b|
+ double nrmb = blas::nrm2(N, b);
+ if (nrmb < D_PREC) nrmb = D_ONE;
+
+ // r = r0 = b - A x0
+ blas::copy(N, b, r0);
+ A.amux(D_MONE, x, r0);
+ blas::copy(N, r0, r);
+
+ resid = blas::nrm2(N, r) / nrmb;
+
+ if (resid < eps) {
+ eps = resid;
+ nsteps = 0;
+ delete[] v;
+ return 0;
+ }
+
+ double alpha = D_ZERO, omega = D_ZERO, rho2 = D_ZERO;
+
+ for (unsigned l = 1; l <= nsteps; ++l) {
+ // rho1 = r0 * r
+ const double rho1 = blas::scpr(N, r0, r);
+ if (fabs(rho1) < D_PREC) {
+ eps = blas::nrm2(N, r) / nrmb;
+ delete[] v;
+ return 2;
+ }
+
+ if (l == 1)
+ blas::copy(N, r, p); // p = r
+ else {
+ blas::axpy(N, -omega, v, p); // p = (p-omega v)*beta+r
+ const double beta = rho1 * alpha / (rho2 * omega);
+ blas::scal(N, beta, p);
+ blas::axpy(N, D_ONE, r, p);
+ }
+
+ // p^ = C p
+ blas::copy(N, p, phat);
+ A.precondApply(phat);
+ // v = A p^
+ blas::setzero(N, v);
+ A.amux(D_ONE, phat, v);
+
+ alpha = rho1 / blas::scpr(N, r0, v);
+
+ // s = r - alpha v
+ blas::copy(N, r, s);
+ blas::axpy(N, -alpha, v, s);
+
+ resid = blas::nrm2(N, s) / nrmb;
+#ifndef NDEBUG
+ std::cout << "Step " << l << ", resid=" << resid << std::endl;
+#endif
+ if (resid < eps) {
+ // x += alpha p^
+ blas::axpy(N, alpha, phat, x);
+ eps = resid;
+ nsteps = l;
+ delete[] v;
+ return 0;
+ }
+
+ // s^ = C s
+ blas::copy(N, s, shat);
+ A.precondApply(shat);
+
+ // t = A s^
+ blas::setzero(N, t);
+ A.amux(D_ONE, shat, t);
+
+ // omega = t*s / t*t
+ omega = blas::scpr(N, t, s) / blas::scpr(N, t, t);
+
+ // x += alpha p^ + omega s^
+ blas::axpy(N, alpha, phat, x);
+ blas::axpy(N, omega, shat, x);
+
+ // r = s - omega t
+ blas::copy(N, s, r);
+ blas::axpy(N, -omega, t, r);
+
+ rho2 = rho1;
+
+ resid = blas::nrm2(N, r) / nrmb;
+
+ if (resid < eps) {
+ eps = resid;
+ nsteps = l;
+ delete[] v;
+ return 0;
+ }
+
+ if (fabs(omega) < D_PREC) {
+ eps = resid;
+ delete[] v;
+ return 3;
+ }
+ }
+
+ eps = resid;
+ delete[] v;
+ return 1;
+}
+
+} // end namespace MathLib
diff --git a/MathLib/LinAlg/Solvers/BiCGStab.h b/MathLib/LinAlg/Solvers/BiCGStab.h
new file mode 100644
index 0000000000000000000000000000000000000000..1f7821b7f29494955e229bfa7d0d9923013ce725
--- /dev/null
+++ b/MathLib/LinAlg/Solvers/BiCGStab.h
@@ -0,0 +1,22 @@
+/*
+ * BiCGStab.h
+ *
+ * Created on: Oct 4, 2011
+ * Author: TF
+ */
+
+#ifndef BICGSTAB_H_
+#define BICGSTAB_H_
+
+#include "blas.h"
+#include "../Sparse/CRSMatrix.h"
+#include "../Sparse/CRSMatrixDiagPrecond.h"
+
+namespace MathLib {
+
+unsigned BiCGStab(CRSMatrix<double> const& A, double* const b, double* const x,
+ double& eps, unsigned& nsteps);
+
+} // end namespace MathLib
+
+#endif /* BICGSTAB_H_ */
diff --git a/MathLib/LinAlg/Solvers/solver.h b/MathLib/LinAlg/Solvers/solver.h
index d89a01b897ba2a07e99a3539b26318241ef468d7..280a62ac9afdea49467446ab5a61417d06b748a2 100644
--- a/MathLib/LinAlg/Solvers/solver.h
+++ b/MathLib/LinAlg/Solvers/solver.h
@@ -8,11 +8,7 @@
#ifndef SOLVER_H_
#define SOLVER_H_
-namespace MathLib {
-
-unsigned CG(CRSMatrix<double> const * mat, double const * const b,
- double* const x, double& eps, unsigned& nsteps, unsigned num_threads = 1);
-
-} // end namespace MathLib
+#include "CG.h"
+#include "BiCGStab.h"
#endif /* SOLVER_H_ */
diff --git a/SimpleTests/SolverTests/BiCGStabDiagPrecond.cpp b/SimpleTests/SolverTests/BiCGStabDiagPrecond.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..764ad1b963e5de459d7bf975c5f84a9fd611dac0
--- /dev/null
+++ b/SimpleTests/SolverTests/BiCGStabDiagPrecond.cpp
@@ -0,0 +1,82 @@
+#include <fstream>
+#include <iostream>
+#include <cmath>
+#include <cstdlib>
+#include <omp.h>
+#include "LinAlg/Solvers/BiCGStab.h"
+#include "LinAlg/Sparse/CRSMatrixDiagPrecond.h"
+#include "sparse.h"
+#include "vector_io.h"
+#include "RunTimeTimer.h"
+#include "CPUTimeTimer.h"
+
+int main(int argc, char *argv[])
+{
+ if (argc != 4) {
+ std::cout << "Usage: " << argv[0] << " matrix rhs number-of-threads" << std::endl;
+ return -1;
+ }
+
+ // read number of threads
+ unsigned num_omp_threads (1);
+ num_omp_threads = atoi (argv[3]);
+
+ // *** reading matrix in crs format from file
+ std::string fname(argv[1]);
+ MathLib::CRSMatrixDiagPrecond *mat (new MathLib::CRSMatrixDiagPrecond(fname));
+
+ unsigned n (mat->getNRows());
+ bool verbose (true);
+ if (verbose)
+ std::cout << "Parameters read: n=" << n << std::endl;
+
+ double *x(new double[n]);
+ double *b(new double[n]);
+
+ // *** init start vector x
+ for (size_t k(0); k<n; k++) {
+ x[k] = 0.0;
+ }
+ // *** read rhs
+ fname = argv[2];
+ std::ifstream in(fname.c_str());
+ if (in) {
+ read (in, n, b);
+ in.close();
+ } else {
+ std::cout << "problem reading rhs - initializing b with 1.0" << std::endl;
+ for (size_t k(0); k<n; k++) {
+ b[k] = 1.0;
+ }
+ }
+
+
+ if (verbose)
+ std::cout << "solving system with PCG method (diagonal preconditioner) ... " << std::flush;
+
+ double eps (1.0e-6);
+ unsigned steps (4000);
+ RunTimeTimer run_timer;
+ CPUTimeTimer cpu_timer;
+ run_timer.start();
+ cpu_timer.start();
+
+ MathLib::BiCGStab ((*mat), b, x, eps, steps);
+
+ cpu_timer.stop();
+ run_timer.stop();
+
+ if (verbose) {
+ std::cout << " in " << steps << " iterations" << std::endl;
+ std::cout << "\t(residuum is " << eps << ") took " << cpu_timer.elapsed() << " sec time and " << run_timer.elapsed() << " sec" << std::endl;
+ } else {
+ std::cout << cpu_timer.elapsed() << std::endl;
+ }
+
+ delete mat;
+ delete [] x;
+ delete [] b;
+
+ return 0;
+}
+
diff --git a/SimpleTests/SolverTests/CMakeLists.txt b/SimpleTests/SolverTests/CMakeLists.txt
index 439de4a1d68c60c31c5d62213b7168bcf08e945f..952fb044734203e65de086dbc7d61f869b9613dd 100644
--- a/SimpleTests/SolverTests/CMakeLists.txt
+++ b/SimpleTests/SolverTests/CMakeLists.txt
@@ -31,6 +31,12 @@ ADD_EXECUTABLE( ConjugateGradientDiagPrecond
${HEADERS}
)
+ADD_EXECUTABLE( BiCGStabDiagPrecond
+ BiCGStabDiagPrecond.cpp
+ ${SOURCES}
+ ${HEADERS}
+)
+
TARGET_LINK_LIBRARIES ( ConjugateGradientUnpreconditioned
${BLAS_LIBRARIES}
${LAPACK_LIBRARIES}
@@ -45,3 +51,10 @@ TARGET_LINK_LIBRARIES ( ConjugateGradientDiagPrecond
Base
)
+TARGET_LINK_LIBRARIES( BiCGStabDiagPrecond
+ ${BLAS_LIBRARIES}
+ ${LAPACK_LIBRARIES}
+ MathLib
+ Base
+)
+