itl_bicgstab.cpp

Go to the documentation of this file.

/*
    Sheep - A Rigid Body Dynamics Engine
    Copyright (C) 2001-2004 Francois Beaune
    Contact: http://toxicengine.sourceforge.net/

    This file is part of Sheep.

    Sheep is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    Sheep 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 General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sheep; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

// Disable warnings about identifier truncation
// with Microsoft Visual C++ (R).
#pragma warning(disable : 4786)

#include "itl_bicgstab.h"   // include first
#include "itl/interface/mtl.h"
#include "itl/krylov/bicgstab.h"
#include "itl/preconditioner/diagonal.h"
#include "itl/preconditioner/ssor.h"

#include "mtl/matrix.h"
#include "mtl/mtl.h"
#include "mtl/utils.h"

#include "real.h"

#include <cassert>
#include <iostream>

using namespace itl;
using namespace mtl;

using sheep::Real;

typedef matrix<Real, rectangle<>, array<compressed<> >, row_major>::type mtl_matrix;

int sheep::itl_bicgstab(const sheep::Matrix &a, const sheep::Vector &b, sheep::Vector *x) {
    assert(b.GetDimension() == a.GetColumns());
    assert(x);
    assert(x->GetDimension() == a.GetColumns());

    // Build the MTL version of the 'a' matrix.

    mtl_matrix mtl_a(a.GetRows(), a.GetColumns());

    for(int i = 0; i < mtl_a.nrows(); i++) {
        for(int j = 0; j < mtl_a.ncols(); j++) {
            mtl_a(i, j) = a(i, j);
        }
    }

#ifdef VERBOSE
    std::cout << "A:" << std::endl;
    for(int i = 0; i < mtl_a.nrows(); i++) {
        for(int j = 0; j < mtl_a.ncols(); j++) {
            std::cout << mtl_a(i, j) << "\t";
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;
#endif  // VERBOSE

    // Build the MTL version of the 'b' vector.

    dense1D<Real> mtl_b(mtl_a.ncols());

    int idx = 0;

    for(dense1D<Real>::iterator i = mtl_b.begin(); i != mtl_b.end(); i++) {
        *i = b[idx++];
    }

#ifdef VERBOSE
    std::cout << "B:" << std::endl;
    for(dense1D<Real>::iterator i = mtl_b.begin(); i != mtl_b.end(); i++) {
        std::cout << *i << "\t";
    }
    std::cout << std::endl;
#endif  // VERBOSE

    // Run the Stabilized Biconjugate Gradient Solver.

    dense1D<Real> mtl_x(mtl_a.nrows(), Real(0));

//  SSOR<mtl_matrix> precond(mtl_a);
    diagonal_precond<mtl_matrix> precond(mtl_a);
//  identity_preconditioner precond;

    basic_iteration<Real> iter(mtl_b, 50, 1e-6);

    bicgstab(mtl_a, mtl_x, mtl_b, precond(), iter);

#ifdef VERBOSE
    std::cout << "X:" << std::endl;
    for(dense1D<Real>::iterator i = mtl_x.begin(); i != mtl_x.end(); i++) {
        std::cout << *i << "\t";
    }
    std::cout << std::endl;
#endif  // VERBOSE

    // Copy back the solution vector 'mtl_x' to 'x'.

    idx = 0;

    for(dense1D<Real>::iterator i = mtl_x.begin(); i != mtl_x.end(); i++) {
        (*x)[idx++] = *i;
    }

    return 0;
}

---