All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
blas2.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
============================================================================= */
// System headers
#include <iostream>
// uBLAS headers
#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/vector_proxy.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
// Must be set if you want to use ViennaCL algorithms on ublas objects
#define VIENNACL_WITH_UBLAS 1
// ViennaCL headers
#include "viennacl/scalar.hpp"
#include "viennacl/vector.hpp"
#include "viennacl/matrix.hpp"
#include "viennacl/linalg/direct_solve.hpp"
#include "viennacl/linalg/prod.hpp" //generic matrix-vector product
#include "viennacl/linalg/norm_2.hpp" //generic l2-norm for vectors
#include "viennacl/linalg/lu.hpp" //LU substitution routines
#include "viennacl/tools/random.hpp"
// Make `boost::numeric::ublas` available under the shortcut `ublas`:
using namespace boost::numeric;
int main()
{
typedef float ScalarType;
viennacl::tools::uniform_random_numbers<ScalarType> randomNumber;
ublas::vector<ScalarType> rhs(12);
for (unsigned int i = 0; i < rhs.size(); ++i)
rhs(i) = randomNumber();
ublas::vector<ScalarType> rhs2 = rhs;
ublas::vector<ScalarType> result = ublas::zero_vector<ScalarType>(10);
ublas::vector<ScalarType> result2 = result;
ublas::vector<ScalarType> rhs_trans = rhs;
rhs_trans.resize(result.size(), true);
ublas::vector<ScalarType> result_trans = ublas::zero_vector<ScalarType>(rhs.size());
ublas::matrix<ScalarType> matrix(result.size(),rhs.size());
for (unsigned int i = 0; i < matrix.size1(); ++i)
for (unsigned int j = 0; j < matrix.size2(); ++j)
matrix(i,j) = randomNumber();
std::vector< ScalarType > stl_result(result.size());
std::vector< ScalarType > stl_rhs(rhs.size());
std::vector< std::vector<ScalarType> > stl_matrix(result.size());
for (unsigned int i=0; i < result.size(); ++i)
{
stl_matrix[i].resize(rhs.size());
for (unsigned int j = 0; j < matrix.size2(); ++j)
{
stl_rhs[j] = rhs[j];
stl_matrix[i][j] = matrix(i,j);
}
}
viennacl::vector<ScalarType> vcl_rhs(rhs.size());
viennacl::vector<ScalarType> vcl_result(result.size());
viennacl::matrix<ScalarType> vcl_matrix(result.size(), rhs.size());
viennacl::matrix<ScalarType> vcl_matrix2(result.size(), rhs.size());
viennacl::copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
viennacl::copy(matrix, vcl_matrix); //copy from ublas dense matrix type to ViennaCL type
vcl_matrix2 = vcl_matrix;
vcl_matrix2 += vcl_matrix;
vcl_matrix2 -= vcl_matrix;
vcl_matrix2 = vcl_matrix2 + vcl_matrix;
vcl_matrix2 = vcl_matrix2 - vcl_matrix;
viennacl::scalar<ScalarType> vcl_3(3.0);
vcl_matrix2 *= ScalarType(2.0);
vcl_matrix2 /= ScalarType(2.0);
vcl_matrix2 *= vcl_3;
vcl_matrix2 /= vcl_3;
vcl_matrix.clear();
viennacl::copy(stl_matrix, vcl_matrix); //alternative: copy from STL vector< vector<> > type to ViennaCL type
//for demonstration purposes (no effect):
viennacl::copy(vcl_matrix, matrix); //copy back from ViennaCL to ublas type.
viennacl::copy(vcl_matrix, stl_matrix); //copy back from ViennaCL to STL type.
std::cout << "----- Matrix-Vector product -----" << std::endl;
result = ublas::prod(matrix, rhs); //the ublas way
stl_result = viennacl::linalg::prod(stl_matrix, stl_rhs); //using STL
vcl_result = viennacl::linalg::prod(vcl_matrix, vcl_rhs); //the ViennaCL way
std::cout << "----- Transposed Matrix-Vector product -----" << std::endl;
result_trans = prod(trans(matrix), rhs_trans);
viennacl::vector<ScalarType> vcl_rhs_trans(rhs_trans.size());
viennacl::vector<ScalarType> vcl_result_trans(result_trans.size());
viennacl::copy(rhs_trans.begin(), rhs_trans.end(), vcl_rhs_trans.begin());
vcl_result_trans = viennacl::linalg::prod(trans(vcl_matrix), vcl_rhs_trans);
ublas::matrix<ScalarType> tri_matrix(10,10);
for (std::size_t i=0; i<tri_matrix.size1(); ++i)
{
for (std::size_t j=0; j<i; ++j)
tri_matrix(i,j) = 0.0;
for (std::size_t j=i; j<tri_matrix.size2(); ++j)
tri_matrix(i,j) = matrix(i,j);
}
viennacl::matrix<ScalarType> vcl_tri_matrix = viennacl::identity_matrix<ScalarType>(tri_matrix.size1());
viennacl::copy(tri_matrix, vcl_tri_matrix);
// Bring vectors to correct size:
rhs.resize(tri_matrix.size1(), true);
rhs2.resize(tri_matrix.size1(), true);
vcl_rhs.resize(tri_matrix.size1(), true);
viennacl::copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
vcl_result.resize(10);
std::cout << "----- Upper Triangular solve -----" << std::endl;
result = ublas::solve(tri_matrix, rhs, ublas::upper_tag()); //ublas
vcl_result = viennacl::linalg::solve(vcl_tri_matrix, vcl_rhs, viennacl::linalg::upper_tag()); //ViennaCL
ublas::inplace_solve(tri_matrix, rhs, ublas::upper_tag()); //ublas
viennacl::linalg::inplace_solve(vcl_tri_matrix, vcl_rhs, viennacl::linalg::upper_tag()); //ViennaCL
std::cout << "----- LU factorization -----" << std::endl;
std::size_t lu_dim = 300;
ublas::matrix<ScalarType> square_matrix(lu_dim, lu_dim);
ublas::vector<ScalarType> lu_rhs(lu_dim);
viennacl::matrix<ScalarType> vcl_square_matrix(lu_dim, lu_dim);
viennacl::vector<ScalarType> vcl_lu_rhs(lu_dim);
for (std::size_t i=0; i<lu_dim; ++i)
for (std::size_t j=0; j<lu_dim; ++j)
square_matrix(i,j) = randomNumber();
//put some more weight on diagonal elements:
for (std::size_t j=0; j<lu_dim; ++j)
{
square_matrix(j,j) += ScalarType(10.0);
lu_rhs(j) = randomNumber();
}
viennacl::copy(square_matrix, vcl_square_matrix);
viennacl::copy(lu_rhs, vcl_lu_rhs);
viennacl::linalg::lu_factorize(vcl_square_matrix);
viennacl::linalg::lu_substitute(vcl_square_matrix, vcl_lu_rhs);
viennacl::copy(square_matrix, vcl_square_matrix);
viennacl::copy(lu_rhs, vcl_lu_rhs);
ublas::lu_factorize(square_matrix);
ublas::inplace_solve (square_matrix, lu_rhs, ublas::unit_lower_tag ());
ublas::inplace_solve (square_matrix, lu_rhs, ublas::upper_tag ());
viennacl::linalg::lu_factorize(vcl_square_matrix);
viennacl::linalg::lu_substitute(vcl_square_matrix, vcl_lu_rhs);
std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
return EXIT_SUCCESS;
}
viennacl::linalg::cuda::inplace_solve
void inplace_solve(matrix_base< NumericT > const &A, matrix_base< NumericT > &B, SolverTagT tag)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:253
viennacl::linalg::inplace_solve
void inplace_solve(const matrix_base< NumericT > &A, matrix_base< NumericT > &B, SolverTagT)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:217
viennacl::scalar
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
norm_2.hpp
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
trans
std::vector< std::vector< NumericT > > trans(std::vector< std::vector< NumericT > > const &A)
Definition: blas3_solve.cpp:195
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
matrix.hpp
Implementation of the dense matrix class.
main
int main()
Definition: bisect.cpp:91
viennacl::linalg::lu_substitute
void lu_substitute(matrix< NumericT, F1, AlignmentV1 > const &A, matrix< NumericT, F2, AlignmentV2 > &B)
LU substitution for the system LU = rhs.
Definition: lu.hpp:201
viennacl::matrix
A dense matrix class.
Definition: forwards.h:375
viennacl::identity_matrix
Represents a vector consisting of 1 at a given index and zeros otherwise. To be used as an initialize...
Definition: matrix_def.hpp:69
viennacl::linalg::solve
VectorT solve(MatrixT const &matrix, VectorT const &rhs, bicgstab_tag const &tag, PreconditionerT const &precond)
Definition: bicgstab.hpp:496
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
viennacl::tools::uniform_random_numbers
Random number generator for returning uniformly distributed values in the closed interval [0...
Definition: random.hpp:44
viennacl::linalg::upper_tag
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
viennacl::vector< ScalarType >
lu.hpp
Implementations of LU factorization for row-major and column-major dense matrices.
direct_solve.hpp
Implementations of dense direct solvers are found here.
prod
void prod(std::vector< std::map< IndexT, NumericT > > const &stl_A, std::vector< std::map< IndexT, NumericT > > const &stl_B, std::vector< std::map< IndexT, NumericT > > &stl_C)
Definition: sparse_prod.cpp:114
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
random.hpp
A small collection of sequential random number generators.
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
scalar.hpp
Implementation of the ViennaCL scalar class.
viennacl::linalg::lu_factorize
void lu_factorize(matrix< NumericT, viennacl::row_major > &A)
LU factorization of a row-major dense matrix.
Definition: lu.hpp:42