jacobi.cc

// Copyright (C) 2011-2014 Vincent Heuveline
//
// HiFlow3 is free software: you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License as published by the Free
// Software Foundation, either version 3 of the License, or (at your option) any
// later version.
//  
// HiFlow3 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 Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with HiFlow3.  If not, see <http://www.gnu.org/licenses/>.




/// @author Chandramowli Subramanian

#include <cassert>
#include <cmath>

#include "jacobi.h"
#include "linear_algebra/la_descriptor.h"
#include "common/log.h"

namespace hiflow {
namespace la {

/// standard constructor
template<class LAD>
Jacobi<LAD>::Jacobi() {
  this->op_ = NULL;
  LOG_INFO("Linear solver", "Jacobi");
  //LOG_INFO("Preconditioning", this->method());
}

/// destructor
template<class LAD>
Jacobi<LAD>::~Jacobi() {
  this->op_ = NULL;
}

/// Sets the operator of the linear system.
template<class LAD>
void Jacobi<LAD>::SetupOperator(const OperatorType& op) {
  this->op_ = &op;
}


/// Solves the linear system.
/// @param [in] b right hand side vector
/// @param [in,out] x start and solution vector
template<class LAD>
LinearSolverState Jacobi<LAD>::Solve(const VectorType& b, VectorType* x) {
	IterateControl::State conv = IterateControl::kIterate;
	VectorType out;
	out.CloneFromWithoutContent(*x);
	VectorType res_h;
	res_h.CloneFromWithoutContent(*x);
	this->op_->VectorMult(*x, &res_h); // = Ax
	res_h.Axpy(b, -1.); // = b - Ax

	this->res_ = res_h.Norm2(); // r = b - Ax
	int iter = 0;
	conv = this->control().Check(iter, this->res());
	// op get diag elems once
	//this->op_->PrepareForJacobi();
	
	while (conv == IterateControl::kIterate){
		++iter;
		this->op_->JacobiIter(b, *x, out);
		x->CopyFrom(out);

		res_h.CloneFromWithoutContent(*x);
		this->op_->VectorMult(*x, &res_h);
		res_h.Axpy(b, -1.);

		this->res_ = res_h.Norm2();
		std::cout < "Residuum: " << this->res() << "\n";
		conv = this->control().Check(iter, this->res());
		if (conv != IterateControl::kIterate) {
			break;
		}
	}
	std::cout << "Iterations done: " << iter << "\n";

	if (conv == IterateControl::kFailureDivergenceTol ||
		conv == IterateControl::kFailureMaxitsExceeded)
		return kSolverExceeded;
	else
		return kSolverSuccess;
}

/// template instantiation
template class Jacobi<LADescriptorCoupledD>;
template class Jacobi<LADescriptorCoupledS>;

} // namespace la
} // namespace hiflow