ViennaCL - The Vienna Computing Library: /data/development/ViennaCL/dev/viennacl/linalg/bicgstab.hpp Source File

Go to the documentation of this file.
00001 #ifndef VIENNACL_BICGSTAB_HPP_
00002 #define VIENNACL_BICGSTAB_HPP_
00003 
00004 /* =========================================================================
00005    Copyright (c) 2010-2011, Institute for Microelectronics,
00006                             Institute for Analysis and Scientific Computing,
00007                             TU Wien.
00008 
00009                             -----------------
00010                   ViennaCL - The Vienna Computing Library
00011                             -----------------
00012 
00013    Project Head:    Karl Rupp                   rupp@iue.tuwien.ac.at
00014                
00015    (A list of authors and contributors can be found in the PDF manual)
00016 
00017    License:         MIT (X11), see file LICENSE in the base directory
00018 ============================================================================= */
00019 
00024 #include <vector>
00025 #include <cmath>
00026 #include "viennacl/forwards.h"
00027 #include "viennacl/tools/tools.hpp"
00028 #include "viennacl/linalg/prod.hpp"
00029 #include "viennacl/linalg/inner_prod.hpp"
00030 #include "viennacl/traits/clear.hpp"
00031 #include "viennacl/traits/size.hpp"
00032 #include "viennacl/meta/result_of.hpp"
00033 
00034 namespace viennacl
00035 {
00036   namespace linalg
00037   {
00038     
00041     class bicgstab_tag
00042     {
00043       public:
00049         bicgstab_tag(double tol = 1e-8, unsigned int max_iterations = 300) : _tol(tol), _iterations(max_iterations) {};
00050       
00052         double tolerance() const { return _tol; }
00054         unsigned int max_iterations() const { return _iterations; }
00055         
00057         unsigned int iters() const { return iters_taken_; }
00058         void iters(unsigned int i) const { iters_taken_ = i; }
00059         
00061         double error() const { return last_error_; }
00063         void error(double e) const { last_error_ = e; }
00064         
00065       private:
00066         double _tol;
00067         unsigned int _iterations;
00068 
00069         //return values from solver
00070         mutable unsigned int iters_taken_;
00071         mutable double last_error_;
00072     };
00073     
00074 
00084     template <typename MatrixType, typename VectorType>
00085     VectorType solve(const MatrixType & matrix, VectorType const & rhs, bicgstab_tag const & tag)
00086     {
00087       typedef typename viennacl::result_of::value_type<VectorType>::type        ScalarType;
00088       typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type    CPU_ScalarType;
00089       unsigned int problem_size = viennacl::traits::size(rhs);
00090       VectorType result(problem_size);
00091       viennacl::traits::clear(result);
00092 
00093       VectorType residual = rhs;
00094       VectorType p = rhs;
00095       VectorType r0star = rhs;
00096       VectorType tmp0(problem_size);
00097       VectorType tmp1(problem_size);
00098       VectorType s(problem_size);
00099 
00100       CPU_ScalarType ip_rr0star = viennacl::linalg::inner_prod(rhs,r0star);
00101       CPU_ScalarType norm_rhs_host = ip_rr0star;
00102       CPU_ScalarType beta;
00103       CPU_ScalarType alpha;
00104       CPU_ScalarType omega;
00105       ScalarType inner_prod_temp; //temporary variable for inner product computation
00106       ScalarType new_ip_rr0star = 0;
00107       
00108       for (unsigned int i = 0; i < tag.max_iterations(); ++i)
00109       {
00110         tag.iters(i+1);
00111         tmp0 = viennacl::linalg::prod(matrix, p);
00112         //alpha = ip_rr0star / viennacl::linalg::inner_prod(tmp0, r0star);
00113         inner_prod_temp = viennacl::linalg::inner_prod(tmp0, r0star);
00114         alpha = ip_rr0star / static_cast<CPU_ScalarType>(inner_prod_temp);
00115 
00116         //s = residual - alpha*tmp0;
00117         s = residual;
00118         s -= alpha*tmp0;
00119         
00120         tmp1 = viennacl::linalg::prod(matrix, s);
00121         //omega = viennacl::linalg::inner_prod(tmp1, s) / viennacl::linalg::inner_prod(tmp1, tmp1);
00122         inner_prod_temp = viennacl::linalg::inner_prod(tmp1, s);
00123         omega = inner_prod_temp;
00124         inner_prod_temp = viennacl::linalg::inner_prod(tmp1, tmp1);
00125         omega /= inner_prod_temp;
00126         
00127         //result += alpha * p + omega * s;
00128         result += alpha * p;
00129         result += omega * s;
00130         
00131         //residual = s - omega * tmp1;
00132         residual = s;
00133         residual -= omega*tmp1;
00134         
00135         new_ip_rr0star = viennacl::linalg::inner_prod(residual,r0star);
00136         if (fabs(CPU_ScalarType(viennacl::linalg::inner_prod(residual, residual)) / norm_rhs_host) < tag.tolerance() * tag.tolerance())
00137           break;
00138         
00139         //beta = new_ip_rr0star / ip_rr0star * alpha/omega;
00140         CPU_ScalarType cpu_temp = new_ip_rr0star; //read from device only once
00141         beta = cpu_temp / ip_rr0star * alpha/omega;
00142         ip_rr0star = cpu_temp;
00143 
00144         // Execution of
00145         //  p = residual + beta * (p - omega*tmp0);
00146         // without introducing temporary vectors:
00147         p -= omega * tmp0;
00148         p *= beta;
00149         p += residual;
00150       }
00151       
00152       //store last error estimate:
00153       tag.error(std::sqrt(fabs(CPU_ScalarType(viennacl::linalg::inner_prod(residual, residual)) / norm_rhs_host)));
00154       
00155       return result;
00156     }
00157 
00158     template <typename MatrixType, typename VectorType>
00159     VectorType solve(const MatrixType & matrix, VectorType const & rhs, bicgstab_tag const & tag, viennacl::linalg::no_precond)
00160     {
00161       return solve(matrix, rhs, tag);
00162     }
00163 
00174     template <typename MatrixType, typename VectorType, typename PreconditionerType>
00175     VectorType solve(const MatrixType & matrix, VectorType const & rhs, bicgstab_tag const & tag, PreconditionerType const & precond)
00176     {
00177       typedef typename viennacl::result_of::value_type<VectorType>::type        ScalarType;
00178       typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type    CPU_ScalarType;
00179       unsigned int problem_size = viennacl::traits::size(rhs);
00180       VectorType result(problem_size);
00181       result.clear();
00182 
00183       VectorType residual = rhs;
00184       precond.apply(residual);
00185       VectorType r0star = residual;  //can be chosen arbitrarily in fact
00186       VectorType tmp0(problem_size);
00187       VectorType tmp1(problem_size);
00188       VectorType s(problem_size);
00189       
00190       VectorType p = residual;
00191 
00192       CPU_ScalarType ip_rr0star = viennacl::linalg::inner_prod(residual,r0star);
00193       CPU_ScalarType norm_rhs_host = ip_rr0star;
00194       CPU_ScalarType beta;
00195       CPU_ScalarType alpha;
00196       CPU_ScalarType omega;
00197       ScalarType new_ip_rr0star = 0;
00198       ScalarType inner_prod_temp; //temporary variable for inner product
00199       
00200       for (unsigned int i = 0; i < tag.max_iterations(); ++i)
00201       {
00202         tag.iters(i+1);
00203         tmp0 = viennacl::linalg::prod(matrix, p);
00204         precond.apply(tmp0);
00205         //alpha = ip_rr0star / viennacl::linalg::inner_prod(tmp0, r0star);
00206         inner_prod_temp = viennacl::linalg::inner_prod(tmp0, r0star);
00207         alpha = ip_rr0star / static_cast<CPU_ScalarType>(inner_prod_temp);
00208 
00209         //s = residual - alpha*tmp0;
00210         s = residual;
00211         s -= alpha*tmp0;
00212 
00213         tmp1 = viennacl::linalg::prod(matrix, s);
00214         precond.apply(tmp1);
00215         //omega = viennacl::linalg::inner_prod(tmp1, s) / viennacl::linalg::inner_prod(tmp1, tmp1);
00216         inner_prod_temp = viennacl::linalg::inner_prod(tmp1, s);
00217         omega = inner_prod_temp;
00218         inner_prod_temp = viennacl::linalg::inner_prod(tmp1, tmp1);
00219         omega /= inner_prod_temp;
00220         
00221         //result += alpha * p + omega * s;
00222         result += alpha * p;
00223         result += omega * s;
00224         //residual = s - omega * tmp1;
00225         residual = s;
00226         residual -= omega*tmp1;
00227         
00228         new_ip_rr0star = viennacl::linalg::inner_prod(residual,r0star);
00229         if (fabs(CPU_ScalarType(viennacl::linalg::inner_prod(residual, residual) / norm_rhs_host)) < tag.tolerance() * tag.tolerance() )
00230           break;
00231         
00232         //beta = new_ip_rr0star / ip_rr0star * alpha/omega;
00233         CPU_ScalarType cpu_temp = new_ip_rr0star; //read from device only once
00234         beta = cpu_temp / ip_rr0star * alpha/omega;
00235         ip_rr0star = cpu_temp;
00236 
00237         // Execution of
00238         //  p = residual + beta * (p - omega*tmp0);
00239         // without introducing temporary vectors:
00240         p -= omega * tmp0;
00241         p *= beta;
00242         p += residual;
00243         
00244         //std::cout << "Rel. Residual in current step: " << std::sqrt(std::fabs(viennacl::linalg::inner_prod(residual, residual) / norm_rhs_host)) << std::endl;
00245       }
00246       
00247       //store last error estimate:
00248       tag.error(std::sqrt(fabs(CPU_ScalarType(viennacl::linalg::inner_prod(residual, residual)) / norm_rhs_host)));
00249       
00250       return result;
00251     }
00252 
00253   }
00254 }
00255 
00256 #endif