All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iterative.cpp
Go to the documentation of this file.
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the PDF manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
//
// include necessary system headers
//
#include <iostream>
//
// Necessary to obtain a suitable performance in ublas
#ifndef NDEBUG
#define BOOST_UBLAS_NDEBUG
#endif
//
// ublas includes
//
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/operation_sparse.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
// Must be set if you want to use ViennaCL algorithms on ublas objects
#define VIENNACL_WITH_UBLAS 1
//
// ViennaCL includes
//
#include "viennacl/scalar.hpp"
#include "viennacl/vector.hpp"
#include "viennacl/compressed_matrix.hpp"
#include "viennacl/coordinate_matrix.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/linalg/ilu.hpp"
#include "viennacl/linalg/jacobi_precond.hpp"
#include "viennacl/linalg/cg.hpp"
#include "viennacl/linalg/bicgstab.hpp"
#include "viennacl/linalg/gmres.hpp"
#include "viennacl/io/matrix_market.hpp"
// Some helper functions for this tutorial:
#include "vector-io.hpp"
using namespace boost::numeric;
int main()
{
typedef float ScalarType;
ublas::vector<ScalarType> rhs;
ublas::vector<ScalarType> ref_result;
ublas::vector<ScalarType> result;
ublas::compressed_matrix<ScalarType> ublas_matrix;
if (!viennacl::io::read_matrix_market_file(ublas_matrix, "../examples/testdata/mat65k.mtx"))
{
std::cout << "Error reading Matrix file" << std::endl;
return EXIT_FAILURE;
}
if (!readVectorFromFile("../examples/testdata/rhs65025.txt", rhs))
{
std::cout << "Error reading RHS file" << std::endl;
return EXIT_FAILURE;
}
if (!readVectorFromFile("../examples/testdata/result65025.txt", ref_result))
{
std::cout << "Error reading Result file" << std::endl;
return EXIT_FAILURE;
}
std::size_t vcl_size = rhs.size();
viennacl::compressed_matrix<ScalarType> vcl_compressed_matrix;
viennacl::coordinate_matrix<ScalarType> vcl_coordinate_matrix;
viennacl::vector<ScalarType> vcl_rhs(vcl_size);
viennacl::vector<ScalarType> vcl_result(vcl_size);
viennacl::vector<ScalarType> vcl_ref_result(vcl_size);
viennacl::copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
viennacl::copy(ref_result.begin(), ref_result.end(), vcl_ref_result.begin());
viennacl::copy(ublas_matrix, vcl_compressed_matrix);
std::vector< std::map< unsigned int, ScalarType> > stl_matrix(rhs.size());
for (ublas::compressed_matrix<ScalarType>::iterator1 iter1 = ublas_matrix.begin1();
iter1 != ublas_matrix.end1();
++iter1)
{
for (ublas::compressed_matrix<ScalarType>::iterator2 iter2 = iter1.begin();
iter2 != iter1.end();
++iter2)
stl_matrix[iter2.index1()][static_cast<unsigned int>(iter2.index2())] = *iter2;
}
std::vector<ScalarType> stl_rhs(rhs.size()), stl_result(result.size());
std::copy(rhs.begin(), rhs.end(), stl_rhs.begin());
viennacl::copy(stl_matrix, vcl_coordinate_matrix);
viennacl::copy(vcl_coordinate_matrix, stl_matrix);
std::cout << "Setting up preconditioners for uBLAS-matrix..." << std::endl;
viennacl::linalg::ilut_precond< ublas::compressed_matrix<ScalarType> > ublas_ilut(ublas_matrix, viennacl::linalg::ilut_tag());
viennacl::linalg::ilu0_precond< ublas::compressed_matrix<ScalarType> > ublas_ilu0(ublas_matrix, viennacl::linalg::ilu0_tag());
viennacl::linalg::block_ilu_precond< ublas::compressed_matrix<ScalarType>,
viennacl::linalg::ilu0_tag> ublas_block_ilu0(ublas_matrix, viennacl::linalg::ilu0_tag());
std::cout << "Setting up preconditioners for ViennaCL-matrix..." << std::endl;
viennacl::linalg::ilut_precond< viennacl::compressed_matrix<ScalarType> > vcl_ilut(vcl_compressed_matrix, viennacl::linalg::ilut_tag());
viennacl::linalg::ilu0_precond< viennacl::compressed_matrix<ScalarType> > vcl_ilu0(vcl_compressed_matrix, viennacl::linalg::ilu0_tag());
viennacl::linalg::block_ilu_precond< viennacl::compressed_matrix<ScalarType>,
viennacl::linalg::ilu0_tag> vcl_block_ilu0(vcl_compressed_matrix, viennacl::linalg::ilu0_tag());
viennacl::linalg::jacobi_precond< ublas::compressed_matrix<ScalarType> > ublas_jacobi(ublas_matrix, viennacl::linalg::jacobi_tag());
viennacl::linalg::jacobi_precond< viennacl::compressed_matrix<ScalarType> > vcl_jacobi(vcl_compressed_matrix, viennacl::linalg::jacobi_tag());
std::cout << "----- CG Method -----" << std::endl;
result = viennacl::linalg::solve(ublas_matrix, rhs, viennacl::linalg::cg_tag());
result = viennacl::linalg::solve(ublas_matrix, rhs, viennacl::linalg::cg_tag(1e-6, 20), ublas_ilut);
result = viennacl::linalg::solve(ublas_matrix, rhs, viennacl::linalg::cg_tag(1e-6, 20), ublas_jacobi);
vcl_result = viennacl::linalg::solve(vcl_compressed_matrix, vcl_rhs, viennacl::linalg::cg_tag());
vcl_result = viennacl::linalg::solve(vcl_compressed_matrix, vcl_rhs, viennacl::linalg::cg_tag(1e-6, 20), vcl_ilut);
vcl_result = viennacl::linalg::solve(vcl_compressed_matrix, vcl_rhs, viennacl::linalg::cg_tag(1e-6, 20), vcl_jacobi);
stl_result = viennacl::linalg::solve(stl_matrix, stl_rhs, viennacl::linalg::cg_tag());
stl_result = viennacl::linalg::solve(stl_matrix, stl_rhs, viennacl::linalg::cg_tag(1e-6, 20), vcl_ilut);
stl_result = viennacl::linalg::solve(stl_matrix, stl_rhs, viennacl::linalg::cg_tag(1e-6, 20), vcl_jacobi);
std::cout << "----- BiCGStab Method -----" << std::endl;
result = viennacl::linalg::solve(ublas_matrix, rhs, viennacl::linalg::bicgstab_tag()); //without preconditioner
result = viennacl::linalg::solve(ublas_matrix, rhs, viennacl::linalg::bicgstab_tag(1e-6, 20), ublas_ilut); //with preconditioner
result = viennacl::linalg::solve(ublas_matrix, rhs, viennacl::linalg::bicgstab_tag(1e-6, 20), ublas_jacobi); //with preconditioner
vcl_result = viennacl::linalg::solve(vcl_compressed_matrix, vcl_rhs, viennacl::linalg::bicgstab_tag()); //without preconditioner
vcl_result = viennacl::linalg::solve(vcl_compressed_matrix, vcl_rhs, viennacl::linalg::bicgstab_tag(1e-6, 20), vcl_ilut); //with preconditioner
vcl_result = viennacl::linalg::solve(vcl_compressed_matrix, vcl_rhs, viennacl::linalg::bicgstab_tag(1e-6, 20), vcl_jacobi); //with preconditioner
stl_result = viennacl::linalg::solve(stl_matrix, stl_rhs, viennacl::linalg::bicgstab_tag());
stl_result = viennacl::linalg::solve(stl_matrix, stl_rhs, viennacl::linalg::bicgstab_tag(1e-6, 20), vcl_ilut);
stl_result = viennacl::linalg::solve(stl_matrix, stl_rhs, viennacl::linalg::bicgstab_tag(1e-6, 20), vcl_jacobi);
std::cout << "----- GMRES Method -----" << std::endl;
result = viennacl::linalg::solve(ublas_matrix, rhs, viennacl::linalg::gmres_tag()); //without preconditioner
result = viennacl::linalg::solve(ublas_matrix, rhs, viennacl::linalg::gmres_tag(1e-6, 20), ublas_ilut);//with preconditioner
result = viennacl::linalg::solve(ublas_matrix, rhs, viennacl::linalg::gmres_tag(1e-6, 20), ublas_jacobi);//with preconditioner
vcl_result = viennacl::linalg::solve(vcl_compressed_matrix, vcl_rhs, viennacl::linalg::gmres_tag()); //without preconditioner
vcl_result = viennacl::linalg::solve(vcl_compressed_matrix, vcl_rhs, viennacl::linalg::gmres_tag(1e-6, 20), vcl_ilut);//with preconditioner
vcl_result = viennacl::linalg::solve(vcl_compressed_matrix, vcl_rhs, viennacl::linalg::gmres_tag(1e-6, 20), vcl_jacobi);//with preconditioner
stl_result = viennacl::linalg::solve(stl_matrix, stl_rhs, viennacl::linalg::gmres_tag());
stl_result = viennacl::linalg::solve(stl_matrix, stl_rhs, viennacl::linalg::gmres_tag(1e-6, 20), vcl_ilut);
stl_result = viennacl::linalg::solve(stl_matrix, stl_rhs, viennacl::linalg::gmres_tag(1e-6, 20), vcl_jacobi);
std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
return EXIT_SUCCESS;
}
matrix_market.hpp
A reader and writer for the matrix market format is implemented here.
viennacl::linalg::ilu0_precond
ILU0 preconditioner class, can be supplied to solve()-routines.
Definition: ilu0.hpp:146
viennacl::linalg::jacobi_precond
Jacobi preconditioner class, can be supplied to solve()-routines. Generic version for non-ViennaCL ma...
Definition: jacobi_precond.hpp:50
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
main
int main()
Definition: bisect.cpp:91
viennacl::linalg::ilu0_tag
A tag for incomplete LU factorization with static pattern (ILU0)
Definition: ilu0.hpp:58
bicgstab.hpp
The stabilized bi-conjugate gradient method is implemented here.
vector-io.hpp
viennacl::linalg::jacobi_tag
A tag for a jacobi preconditioner.
Definition: jacobi_precond.hpp:43
viennacl::linalg::block_ilu_precond
A block ILU preconditioner class, can be supplied to solve()-routines.
Definition: block_ilu.hpp:132
jacobi_precond.hpp
Implementation of a simple Jacobi preconditioner.
viennacl::linalg::gmres_tag
A tag for the solver GMRES. Used for supplying solver parameters and for dispatching the solve() func...
Definition: gmres.hpp:49
coordinate_matrix.hpp
Implementation of the coordinate_matrix class.
viennacl::linalg::solve
VectorT solve(MatrixT const &matrix, VectorT const &rhs, bicgstab_tag const &tag, PreconditionerT const &precond)
Definition: bicgstab.hpp:496
gmres.hpp
Implementations of the generalized minimum residual method are in this file.
ilu.hpp
Implementations of incomplete factorization preconditioners. Convenience header file.
viennacl::linalg::ilut_tag
A tag for incomplete LU factorization with threshold (ILUT)
Definition: ilut.hpp:45
compressed_matrix.hpp
Implementation of the compressed_matrix class.
readVectorFromFile
bool readVectorFromFile(const std::string &filename, VectorType &vec)
Definition: vector-io.hpp:104
viennacl::linalg::ilut_precond
ILUT preconditioner class, can be supplied to solve()-routines.
Definition: ilut.hpp:352
cg.hpp
The conjugate gradient method is implemented here.
viennacl::vector< ScalarType >
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::copy
void copy(std::vector< NumericT > &cpu_vec, circulant_matrix< NumericT, AlignmentV > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
viennacl::linalg::cg_tag
A tag for the conjugate gradient Used for supplying solver parameters and for dispatching the solve()...
Definition: cg.hpp:48
viennacl::compressed_matrix
A sparse square matrix in compressed sparse rows format.
Definition: compressed_matrix.hpp:559
viennacl::linalg::bicgstab_tag
A tag for the stabilized Bi-conjugate gradient solver. Used for supplying solver parameters and for d...
Definition: bicgstab.hpp:47
viennacl::io::read_matrix_market_file
long read_matrix_market_file(MatrixT &mat, const char *file, long index_base=1)
Reads a sparse matrix from a file (MatrixMarket format)
Definition: matrix_market.hpp:339
scalar.hpp
Implementation of the ViennaCL scalar class.
viennacl::coordinate_matrix
A sparse square matrix, where entries are stored as triplets (i,j, val), where i and j are the row an...
Definition: coordinate_matrix.hpp:186