All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sparse.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
//
#include <iostream>
#include <vector>
#include <map>
#include <cmath>
//
// *** ViennaCL
//
#include "viennacl/scalar.hpp"
#include "viennacl/compressed_matrix.hpp"
#include "viennacl/compressed_compressed_matrix.hpp"
#include "viennacl/coordinate_matrix.hpp"
#include "viennacl/ell_matrix.hpp"
#include "viennacl/sliced_ell_matrix.hpp"
#include "viennacl/hyb_matrix.hpp"
#include "viennacl/vector.hpp"
#include "viennacl/vector_proxy.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/linalg/norm_2.hpp"
#include "viennacl/linalg/ilu.hpp"
#include "viennacl/linalg/detail/ilu/common.hpp"
#include "viennacl/io/matrix_market.hpp"
#include "viennacl/tools/random.hpp"
//
// -------------------------------------------------------------
//
template<typename ScalarType>
ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
{
if (s1 != s2)
return (s1 - s2) / std::max(fabs(s1), std::fabs(s2));
return 0;
}
template<typename ScalarType>
ScalarType diff(std::vector<ScalarType> & v1, viennacl::vector<ScalarType> & v2)
{
std::vector<ScalarType> v2_cpu(v2.size());
viennacl::backend::finish();
viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
for (unsigned int i=0;i<v1.size(); ++i)
{
if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
{
//if (std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) < 1e-10 ) //absolute tolerance (avoid round-off issues)
// v2_cpu[i] = 0;
//else
v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
}
else
v2_cpu[i] = 0.0;
if (v2_cpu[i] > 0.0001)
{
//std::cout << "Neighbor: " << i-1 << ": " << v1[i-1] << " vs. " << v2_cpu[i-1] << std::endl;
std::cout << "Error at entry " << i << ": Should: " << v1[i] << " vs. Is: " << v2[i] << std::endl;
//std::cout << "Neighbor: " << i+1 << ": " << v1[i+1] << " vs. " << v2_cpu[i+1] << std::endl;
exit(EXIT_FAILURE);
}
}
ScalarType norm_inf = 0;
for (std::size_t i=0; i<v2_cpu.size(); ++i)
norm_inf = std::max<ScalarType>(norm_inf, std::fabs(v2_cpu[i]));
return norm_inf;
}
template<typename IndexT, typename NumericT, typename SparseMatrixT>
NumericT diff(std::vector<std::map<IndexT, NumericT> > & cpu_A, SparseMatrixT & vcl_A)
{
typedef typename std::map<IndexT, NumericT>::const_iterator RowIterator;
std::vector<std::map<IndexT, NumericT> > from_gpu(vcl_A.size1());
viennacl::backend::finish();
viennacl::copy(vcl_A, from_gpu);
NumericT error = 0;
//step 1: compare all entries from cpu_A with vcl_A:
for (std::size_t i=0; i<cpu_A.size(); ++i)
{
//std::cout << "Row " << row_it.index1() << ": " << std::endl;
for (RowIterator it = cpu_A[i].begin(); it != cpu_A[i].end(); ++it)
{
//std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
NumericT current_error = 0;
NumericT val_cpu_A = it->second;
NumericT val_gpu_A = from_gpu[i][it->first];
NumericT max_val = std::max(std::fabs(val_cpu_A), std::fabs(val_gpu_A));
if (max_val > 0)
current_error = std::fabs(val_cpu_A - val_gpu_A) / max_val;
if (current_error > error)
error = current_error;
}
}
//step 2: compare all entries from gpu_matrix with cpu_matrix (sparsity pattern might differ):
//std::cout << "ViennaCL matrix: " << std::endl;
for (std::size_t i=0; i<from_gpu.size(); ++i)
{
//std::cout << "Row " << row_it.index1() << ": " << std::endl;
for (RowIterator it = from_gpu[i].begin(); it != from_gpu[i].end(); ++it)
{
//std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
NumericT current_error = 0;
NumericT val_gpu_A = it->second;
NumericT val_cpu_A = cpu_A[i][it->first];
NumericT max_val = std::max(std::fabs(val_cpu_A), std::fabs(val_gpu_A));
if (max_val > 0)
current_error = std::fabs(val_cpu_A - val_gpu_A) / max_val;
if (current_error > error)
error = current_error;
}
}
return error;
}
template<typename NumericT, typename VCL_MatrixT, typename Epsilon, typename STLVectorT, typename VCLVectorT>
int strided_matrix_vector_product_test(Epsilon epsilon,
STLVectorT & result, STLVectorT const & rhs,
VCLVectorT & vcl_result, VCLVectorT & vcl_rhs)
{
typedef typename std::map<unsigned int, NumericT>::const_iterator RowIterator;
int retval = EXIT_SUCCESS;
std::vector<std::map<unsigned int, NumericT> > std_A(5);
std_A[0][0] = NumericT(2.0); std_A[0][2] = NumericT(-1.0);
std_A[1][0] = NumericT(3.0); std_A[1][2] = NumericT(-5.0);
std_A[2][1] = NumericT(5.0); std_A[2][2] = NumericT(-2.0);
std_A[3][2] = NumericT(1.0); std_A[3][3] = NumericT(-6.0);
std_A[4][1] = NumericT(7.0); std_A[4][2] = NumericT(-5.0);
//the following computes project(result, slice(1, 3, 5)) = prod(std_A, project(rhs, slice(3, 2, 4)));
for (std::size_t i=0; i<5; ++i)
{
NumericT val = 0;
for (RowIterator it = std_A[i].begin(); it != std_A[i].end(); ++it)
val += it->second * rhs[3 + 2*it->first];
result[1 + 3*i] = val;
}
VCL_MatrixT vcl_sparse_matrix2;
viennacl::copy(std_A, vcl_sparse_matrix2);
viennacl::vector<NumericT> vec(4);
vec(0) = rhs[3];
vec(1) = rhs[5];
vec(2) = rhs[7];
vec(3) = rhs[9];
viennacl::project(vcl_result, viennacl::slice(1, 3, 5)) = viennacl::linalg::prod(vcl_sparse_matrix2, viennacl::project(vcl_rhs, viennacl::slice(3, 2, 4)));
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with strided vectors, part 1" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
retval = EXIT_FAILURE;
}
vcl_result(1) = NumericT(1.0);
vcl_result(4) = NumericT(1.0);
vcl_result(7) = NumericT(1.0);
vcl_result(10) = NumericT(1.0);
vcl_result(13) = NumericT(1.0);
viennacl::project(vcl_result, viennacl::slice(1, 3, 5)) = viennacl::linalg::prod(vcl_sparse_matrix2, vec);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with strided vectors, part 2" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
retval = EXIT_FAILURE;
}
return retval;
}
template< typename NumericT, typename VCL_MATRIX, typename Epsilon >
int resize_test(Epsilon const& epsilon)
{
int retval = EXIT_SUCCESS;
std::vector<std::map<unsigned int, NumericT> > std_A(5);
VCL_MATRIX vcl_matrix;
std_A[0][0] = NumericT(10.0); std_A[0][1] = NumericT(0.1); std_A[0][2] = NumericT(0.2); std_A[0][3] = NumericT(0.3); std_A[0][4] = NumericT(0.4);
std_A[1][0] = NumericT(1.0); std_A[1][1] = NumericT(1.1); std_A[1][2] = NumericT(1.2); std_A[1][3] = NumericT(1.3); std_A[1][4] = NumericT(1.4);
std_A[2][0] = NumericT(2.0); std_A[2][1] = NumericT(2.1); std_A[2][2] = NumericT(2.2); std_A[2][3] = NumericT(2.3); std_A[2][4] = NumericT(2.4);
std_A[3][0] = NumericT(3.0); std_A[3][1] = NumericT(3.1); std_A[3][2] = NumericT(3.2); std_A[3][3] = NumericT(3.3); std_A[3][4] = NumericT(3.4);
std_A[4][0] = NumericT(4.0); std_A[4][1] = NumericT(4.1); std_A[4][2] = NumericT(4.2); std_A[4][3] = NumericT(4.3); std_A[4][4] = NumericT(4.4);
viennacl::copy(std_A, vcl_matrix);
std::vector<std::map<unsigned int, NumericT> > std_B(std_A.size());
viennacl::copy(vcl_matrix, std_B);
std::cout << "Checking for equality after copy..." << std::endl;
if ( std::fabs(diff(std_A, vcl_matrix)) > epsilon )
{
std::cout << "# Error at operation: equality after copy with sparse matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(std_A, vcl_matrix)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing resize to larger..." << std::endl;
std_A.resize(10);
std_A[0][0] = NumericT(10.0); std_A[0][1] = NumericT(0.1); std_A[0][2] = NumericT(0.2); std_A[0][3] = NumericT(0.3); std_A[0][4] = NumericT(0.4);
std_A[1][0] = NumericT( 1.0); std_A[1][1] = NumericT(1.1); std_A[1][2] = NumericT(1.2); std_A[1][3] = NumericT(1.3); std_A[1][4] = NumericT(1.4);
std_A[2][0] = NumericT( 2.0); std_A[2][1] = NumericT(2.1); std_A[2][2] = NumericT(2.2); std_A[2][3] = NumericT(2.3); std_A[2][4] = NumericT(2.4);
std_A[3][0] = NumericT( 3.0); std_A[3][1] = NumericT(3.1); std_A[3][2] = NumericT(3.2); std_A[3][3] = NumericT(3.3); std_A[3][4] = NumericT(3.4);
std_A[4][0] = NumericT( 4.0); std_A[4][1] = NumericT(4.1); std_A[4][2] = NumericT(4.2); std_A[4][3] = NumericT(4.3); std_A[4][4] = NumericT(4.4);
vcl_matrix.resize(10, 10, true);
if ( std::fabs(diff(std_A, vcl_matrix)) > epsilon )
{
std::cout << "# Error at operation: resize (to larger) with sparse matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(std_A, vcl_matrix)) << std::endl;
return EXIT_FAILURE;
}
std_A[5][5] = NumericT(5.5); std_A[5][6] = NumericT(5.6); std_A[5][7] = NumericT(5.7); std_A[5][8] = NumericT(5.8); std_A[5][9] = NumericT(5.9);
std_A[6][5] = NumericT(6.5); std_A[6][6] = NumericT(6.6); std_A[6][7] = NumericT(6.7); std_A[6][8] = NumericT(6.8); std_A[6][9] = NumericT(6.9);
std_A[7][5] = NumericT(7.5); std_A[7][6] = NumericT(7.6); std_A[7][7] = NumericT(7.7); std_A[7][8] = NumericT(7.8); std_A[7][9] = NumericT(7.9);
std_A[8][5] = NumericT(8.5); std_A[8][6] = NumericT(8.6); std_A[8][7] = NumericT(8.7); std_A[8][8] = NumericT(8.8); std_A[8][9] = NumericT(8.9);
std_A[9][5] = NumericT(9.5); std_A[9][6] = NumericT(9.6); std_A[9][7] = NumericT(9.7); std_A[9][8] = NumericT(9.8); std_A[9][9] = NumericT(9.9);
viennacl::copy(std_A, vcl_matrix);
std::cout << "Testing resize to smaller..." << std::endl;
std_A.clear();
std_A.resize(7);
std_A[0][0] = NumericT(10.0); std_A[0][1] = NumericT(0.1); std_A[0][2] = NumericT(0.2); std_A[0][3] = NumericT(0.3); std_A[0][4] = NumericT(0.4);
std_A[1][0] = NumericT( 1.0); std_A[1][1] = NumericT(1.1); std_A[1][2] = NumericT(1.2); std_A[1][3] = NumericT(1.3); std_A[1][4] = NumericT(1.4);
std_A[2][0] = NumericT( 2.0); std_A[2][1] = NumericT(2.1); std_A[2][2] = NumericT(2.2); std_A[2][3] = NumericT(2.3); std_A[2][4] = NumericT(2.4);
std_A[3][0] = NumericT( 3.0); std_A[3][1] = NumericT(3.1); std_A[3][2] = NumericT(3.2); std_A[3][3] = NumericT(3.3); std_A[3][4] = NumericT(3.4);
std_A[4][0] = NumericT( 4.0); std_A[4][1] = NumericT(4.1); std_A[4][2] = NumericT(4.2); std_A[4][3] = NumericT(4.3); std_A[4][4] = NumericT(4.4);
std_A[5][5] = NumericT( 5.5); std_A[5][6] = NumericT(5.6); //std_A[5][7] = NumericT(5.7); std_A[5][8] = NumericT(5.8); std_A[5][9] = NumericT(5.9);
std_A[6][5] = NumericT( 6.5); std_A[6][6] = NumericT(6.6); //std_A[6][7] = NumericT(6.7); std_A[6][8] = NumericT(6.8); std_A[6][9] = NumericT(6.9);
vcl_matrix.resize(7, 7);
//std::cout << std_A << std::endl;
if ( std::fabs(diff(std_A, vcl_matrix)) > epsilon )
{
std::cout << "# Error at operation: resize (to smaller) with sparse matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(std_A, vcl_matrix)) << std::endl;
retval = EXIT_FAILURE;
}
std::vector<NumericT> std_vec(std_A.size(), NumericT(3.1415));
viennacl::vector<NumericT> vcl_vec(std_A.size());
std::cout << "Testing unit lower triangular solve: compressed_matrix" << std::endl;
viennacl::copy(std_vec, vcl_vec);
std::cout << "STL..." << std::endl;
//boost::numeric::ublas::inplace_solve((ublas_matrix), ublas_vec, boost::numeric::ublas::unit_lower_tag());
for (std::size_t i=1; i<std_A.size(); ++i)
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_A[i].begin(); it != std_A[i].end(); ++it)
{
if (it->first < static_cast<unsigned int>(i))
std_vec[i] -= it->second * std_vec[it->first];
else
continue;
}
std::cout << "ViennaCL..." << std::endl;
viennacl::linalg::inplace_solve((vcl_matrix), vcl_vec, viennacl::linalg::unit_lower_tag());
if ( std::fabs(diff(std_vec, vcl_vec)) > epsilon )
{
std::cout << "# Error at operation: unit lower triangular solve" << std::endl;
std::cout << " diff: " << std::fabs(diff(std_vec, vcl_vec)) << std::endl;
retval = EXIT_FAILURE;
}
return retval;
}
//
// -------------------------------------------------------------
//
template< typename NumericT, typename Epsilon >
int test(Epsilon const& epsilon)
{
viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
std::cout << "Testing resizing of compressed_matrix..." << std::endl;
int retval = resize_test<NumericT, viennacl::compressed_matrix<NumericT> >(epsilon);
if (retval != EXIT_SUCCESS)
return retval;
// --------------------------------------------------------------------------
std::vector<NumericT> rhs;
std::vector<NumericT> result;
std::vector<std::map<unsigned int, NumericT> > std_matrix;
if (viennacl::io::read_matrix_market_file(std_matrix, "../examples/testdata/mat65k.mtx") == EXIT_FAILURE)
{
std::cout << "Error reading Matrix file" << std::endl;
return EXIT_FAILURE;
}
//unsigned int cg_mat_size = cg_mat.size();
std::cout << "done reading matrix" << std::endl;
rhs.resize(std_matrix.size());
for (std::size_t i=0; i<rhs.size(); ++i)
{
std_matrix[i][static_cast<unsigned int>(i)] = NumericT(0.5); // Get rid of round-off errors by making row-sums unequal to zero:
rhs[i] = NumericT(1) + randomNumber();
}
// add some random numbers to the double-compressed matrix:
std::vector<std::map<unsigned int, NumericT> > std_cc_matrix(std_matrix.size());
std_cc_matrix[42][199] = NumericT(3.1415);
std_cc_matrix[31][69] = NumericT(2.71);
std_cc_matrix[23][32] = NumericT(6);
std_cc_matrix[177][57] = NumericT(4);
std_cc_matrix[21][97] = NumericT(-4);
std_cc_matrix[92][25] = NumericT(2);
std_cc_matrix[89][62] = NumericT(11);
std_cc_matrix[ 1][ 7] = NumericT(8);
std_cc_matrix[85][41] = NumericT(13);
std_cc_matrix[66][28] = NumericT(8);
std_cc_matrix[21][74] = NumericT(-2);
viennacl::tools::sparse_matrix_adapter<NumericT> adapted_std_cc_matrix(std_cc_matrix, std_matrix.size(), std_matrix.size());
result = rhs;
viennacl::vector<NumericT> vcl_rhs(rhs.size());
viennacl::vector<NumericT> vcl_result(result.size());
viennacl::vector<NumericT> vcl_result2(result.size());
viennacl::compressed_matrix<NumericT> vcl_compressed_matrix(rhs.size(), rhs.size());
viennacl::compressed_compressed_matrix<NumericT> vcl_compressed_compressed_matrix(rhs.size(), rhs.size());
viennacl::coordinate_matrix<NumericT> vcl_coordinate_matrix(rhs.size(), rhs.size());
viennacl::ell_matrix<NumericT> vcl_ell_matrix;
viennacl::sliced_ell_matrix<NumericT> vcl_sliced_ell_matrix;
viennacl::hyb_matrix<NumericT> vcl_hyb_matrix;
viennacl::copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
viennacl::copy(std_matrix, vcl_compressed_matrix);
viennacl::copy(adapted_std_cc_matrix, vcl_compressed_compressed_matrix);
viennacl::copy(std_matrix, vcl_coordinate_matrix);
// --------------------------------------------------------------------------
std::cout << "Testing products: STL" << std::endl;
result = viennacl::linalg::prod(std_matrix, rhs);
std::cout << "Testing products: compressed_matrix" << std::endl;
vcl_result = viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products: compressed_matrix, strided vectors" << std::endl;
retval = strided_matrix_vector_product_test<NumericT, viennacl::compressed_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
if (retval != EXIT_SUCCESS)
return retval;
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
vcl_result = vcl_rhs;
vcl_result += viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with compressed_matrix (+=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
vcl_result = vcl_rhs;
vcl_result -= viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with compressed_matrix (-=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
//
// Triangular solvers for A \ b:
//
std::cout << "Testing unit upper triangular solve: compressed_matrix" << std::endl;
result = rhs;
viennacl::copy(result, vcl_result);
//boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::unit_upper_tag());
for (std::size_t i2=0; i2<std_matrix.size(); ++i2)
{
std::size_t row = std_matrix.size() - i2 - 1;
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[row].begin(); it != std_matrix[row].end(); ++it)
{
if (it->first > static_cast<unsigned int>(row))
result[row] -= it->second * result[it->first];
else
continue;
}
}
viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::unit_upper_tag());
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: unit upper triangular solve with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing upper triangular solve: compressed_matrix" << std::endl;
result = rhs;
viennacl::copy(result, vcl_result);
//boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::upper_tag());
for (std::size_t i2=0; i2<std_matrix.size(); ++i2)
{
std::size_t row = std_matrix.size() - i2 - 1;
NumericT diag = 0;
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[row].begin(); it != std_matrix[row].end(); ++it)
{
if (it->first > static_cast<unsigned int>(row))
result[row] -= it->second * result[it->first];
else if (it->first == static_cast<unsigned int>(row))
diag = it->second;
else
continue;
}
result[row] /= diag;
}
viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::upper_tag());
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: upper triangular solve with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing unit lower triangular solve: compressed_matrix" << std::endl;
result = rhs;
viennacl::copy(result, vcl_result);
//boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::unit_lower_tag());
for (std::size_t i=1; i<std_matrix.size(); ++i)
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[i].begin(); it != std_matrix[i].end(); ++it)
{
if (it->first < static_cast<unsigned int>(i))
result[i] -= it->second * result[it->first];
else
continue;
}
viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::unit_lower_tag());
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: unit lower triangular solve with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing lower triangular solve: compressed_matrix" << std::endl;
result = rhs;
viennacl::copy(result, vcl_result);
//boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::lower_tag());
for (std::size_t i=0; i<std_matrix.size(); ++i)
{
NumericT diag = 0;
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[i].begin(); it != std_matrix[i].end(); ++it)
{
if (it->first < static_cast<unsigned int>(i))
result[i] -= it->second * result[it->first];
else if (it->first == static_cast<unsigned int>(i))
diag = it->second;
else
continue;
}
result[i] /= diag;
}
viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::lower_tag());
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: lower triangular solve with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
//
// Triangular solvers for A^T \ b
//
std::vector<std::map<unsigned int, NumericT> > std_matrix_trans(std_matrix.size());
// compute transpose:
for (std::size_t i=0; i<std_matrix.size(); ++i)
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[i].begin(); it != std_matrix[i].end(); ++it)
std_matrix_trans[i][it->first] = it->second;
std::cout << "Testing transposed unit upper triangular solve: compressed_matrix" << std::endl;
result = rhs;
viennacl::copy(result, vcl_result);
//boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::unit_upper_tag());
for (std::size_t i2=0; i2<std_matrix_trans.size(); ++i2)
{
std::size_t row = std_matrix_trans.size() - i2 - 1;
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[row].begin(); it != std_matrix_trans[row].end(); ++it)
{
if (it->first > static_cast<unsigned int>(row))
result[row] -= it->second * result[it->first];
else
continue;
}
}
viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::unit_upper_tag());
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: unit upper triangular solve with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing transposed upper triangular solve: compressed_matrix" << std::endl;
result = rhs;
viennacl::copy(result, vcl_result);
//boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::upper_tag());
for (std::size_t i2=0; i2<std_matrix_trans.size(); ++i2)
{
std::size_t row = std_matrix_trans.size() - i2 - 1;
NumericT diag = 0;
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[row].begin(); it != std_matrix_trans[row].end(); ++it)
{
if (it->first > static_cast<unsigned int>(row))
result[row] -= it->second * result[it->first];
else if (it->first == static_cast<unsigned int>(row))
diag = it->second;
else
continue;
}
result[row] /= diag;
}
viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::upper_tag());
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: upper triangular solve with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing transposed unit lower triangular solve: compressed_matrix" << std::endl;
result = rhs;
viennacl::copy(result, vcl_result);
//boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::unit_lower_tag());
for (std::size_t i=1; i<std_matrix_trans.size(); ++i)
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[i].begin(); it != std_matrix_trans[i].end(); ++it)
{
if (it->first < static_cast<unsigned int>(i))
result[i] -= it->second * result[it->first];
else
continue;
}
viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::unit_lower_tag());
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: unit lower triangular solve with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing transposed lower triangular solve: compressed_matrix" << std::endl;
result = rhs;
viennacl::copy(result, vcl_result);
//boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::lower_tag());
for (std::size_t i=0; i<std_matrix_trans.size(); ++i)
{
NumericT diag = 0;
for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[i].begin(); it != std_matrix_trans[i].end(); ++it)
{
if (it->first < static_cast<unsigned int>(i))
result[i] -= it->second * result[it->first];
else if (it->first == static_cast<unsigned int>(i))
diag = it->second;
else
continue;
}
result[i] /= diag;
}
viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::lower_tag());
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: lower triangular solve with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
//
//
std::cout << "Testing products: compressed_compressed_matrix" << std::endl;
result = viennacl::linalg::prod(std_cc_matrix, rhs);
vcl_result = viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
{
std::vector<std::map<unsigned int, NumericT> > temp(vcl_compressed_compressed_matrix.size1());
viennacl::copy(vcl_compressed_compressed_matrix, temp);
// check that entries are correct by computing the product again:
result = viennacl::linalg::prod(temp, rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (after copy back)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
}
result = rhs;
result = viennacl::linalg::prod(std_cc_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
vcl_result = vcl_rhs;
vcl_result += viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (+=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
result = rhs;
result = viennacl::linalg::prod(std_cc_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
vcl_result = vcl_rhs;
vcl_result -= viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (-=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
//
//
std::cout << "Testing products: coordinate_matrix" << std::endl;
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result = viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products: coordinate_matrix, strided vectors" << std::endl;
//std::cout << " --> SKIPPING <--" << std::endl;
retval = strided_matrix_vector_product_test<NumericT, viennacl::coordinate_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
if (retval != EXIT_SUCCESS)
return retval;
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
vcl_result = vcl_rhs;
vcl_result += viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with coordinate_matrix (+=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
vcl_result = vcl_rhs;
vcl_result -= viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with coordinate_matrix (-=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
//
//
//std::cout << "Copying ell_matrix" << std::endl;
viennacl::copy(std_matrix, vcl_ell_matrix);
std_matrix.clear();
viennacl::copy(vcl_ell_matrix, std_matrix);// just to check that it works
std::cout << "Testing products: ell_matrix" << std::endl;
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result.clear();
vcl_result = viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
//viennacl::linalg::prod_impl(vcl_ell_matrix, vcl_rhs, vcl_result);
//std::cout << vcl_result << "\n";
//std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
//std::cout << "First entry of result vector: " << vcl_result[0] << std::endl;
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with ell_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products: ell_matrix, strided vectors" << std::endl;
retval = strided_matrix_vector_product_test<NumericT, viennacl::ell_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
if (retval != EXIT_SUCCESS)
return retval;
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
vcl_result = vcl_rhs;
vcl_result += viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with ell_matrix (+=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
vcl_result = vcl_rhs;
vcl_result -= viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with ell_matrix (-=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
//
//
//std::cout << "Copying sliced_ell_matrix" << std::endl;
viennacl::copy(std_matrix, vcl_sliced_ell_matrix);
std::cout << "Testing products: sliced_ell_matrix" << std::endl;
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result.clear();
vcl_result = viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products: sliced_ell_matrix, strided vectors" << std::endl;
retval = strided_matrix_vector_product_test<NumericT, viennacl::sliced_ell_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
if (retval != EXIT_SUCCESS)
return retval;
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
vcl_result = vcl_rhs;
vcl_result += viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix (+=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
vcl_result = vcl_rhs;
vcl_result -= viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix (-=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
//
//
//std::cout << "Copying hyb_matrix" << std::endl;
viennacl::copy(std_matrix, vcl_hyb_matrix);
std_matrix.clear();
viennacl::copy(vcl_hyb_matrix, std_matrix);// just to check that it works
viennacl::copy(std_matrix, vcl_hyb_matrix);
std::cout << "Testing products: hyb_matrix" << std::endl;
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result.clear();
vcl_result = viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
//viennacl::linalg::prod_impl(vcl_hyb_matrix, vcl_rhs, vcl_result);
//std::cout << vcl_result << "\n";
//std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
//std::cout << "First entry of result vector: " << vcl_result[0] << std::endl;
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products: hyb_matrix, strided vectors" << std::endl;
retval = strided_matrix_vector_product_test<NumericT, viennacl::hyb_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
if (retval != EXIT_SUCCESS)
return retval;
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
vcl_result = vcl_rhs;
vcl_result += viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with hyb_matrix (+=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
result = rhs;
result = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
vcl_result = vcl_rhs;
vcl_result -= viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with hyb_matrix (-=)" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
// --------------------------------------------------------------------------
// --------------------------------------------------------------------------
NumericT alpha = static_cast<NumericT>(2.786);
NumericT beta = static_cast<NumericT>(1.432);
copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
copy(result.begin(), result.end(), vcl_result.begin());
copy(result.begin(), result.end(), vcl_result2.begin());
std::cout << "Testing scaled additions of products and vectors" << std::endl;
std::vector<NumericT> result2(result);
result2 = viennacl::linalg::prod(std_matrix, rhs);
for (std::size_t i=0; i<result.size(); ++i)
result[i] = alpha * result2[i] + beta * result[i];
vcl_result2 = alpha * viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs) + beta * vcl_result;
if ( std::fabs(diff(result, vcl_result2)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product (compressed_matrix) with scaled additions" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
return EXIT_FAILURE;
}
vcl_result2.clear();
vcl_result2 = alpha * viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs) + beta * vcl_result;
if ( std::fabs(diff(result, vcl_result2)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product (coordinate_matrix) with scaled additions" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
return EXIT_FAILURE;
}
vcl_result2.clear();
vcl_result2 = alpha * viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs) + beta * vcl_result;
if ( std::fabs(diff(result, vcl_result2)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product (ell_matrix) with scaled additions" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
return EXIT_FAILURE;
}
vcl_result2.clear();
vcl_result2 = alpha * viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs) + beta * vcl_result;
if ( std::fabs(diff(result, vcl_result2)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product (hyb_matrix) with scaled additions" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
return EXIT_FAILURE;
}
std_matrix.clear();
std::cout << "Testing products after clear(): compressed_matrix" << std::endl;
vcl_compressed_matrix.clear();
result = std::vector<NumericT>(result.size(), NumericT(1));
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result = viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products after clear(): compressed_compressed_matrix" << std::endl;
vcl_compressed_compressed_matrix.clear();
result = std::vector<NumericT>(result.size(), NumericT(1));
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result = viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products after clear(): coordinate_matrix" << std::endl;
vcl_coordinate_matrix.clear();
result = std::vector<NumericT>(result.size(), NumericT(1));
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result = viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products after clear(): ell_matrix" << std::endl;
vcl_ell_matrix.clear();
result = std::vector<NumericT>(result.size(), NumericT(1));
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result = viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with ell_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products after clear(): hyb_matrix" << std::endl;
vcl_hyb_matrix.clear();
result = std::vector<NumericT>(result.size(), NumericT(1));
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result = viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
std::cout << "Testing products after clear(): sliced_ell_matrix" << std::endl;
vcl_sliced_ell_matrix.clear();
result = std::vector<NumericT>(result.size(), NumericT(1));
result = viennacl::linalg::prod(std_matrix, rhs);
vcl_result = viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
if ( std::fabs(diff(result, vcl_result)) > epsilon )
{
std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix" << std::endl;
std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
return EXIT_FAILURE;
}
// --------------------------------------------------------------------------
return retval;
}
//
// -------------------------------------------------------------
//
int main()
{
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << "## Test :: Sparse Matrices" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << std::endl;
int retval = EXIT_SUCCESS;
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << std::endl;
{
typedef float NumericT;
NumericT epsilon = static_cast<NumericT>(1E-4);
std::cout << "# Testing setup:" << std::endl;
std::cout << " eps: " << epsilon << std::endl;
std::cout << " numeric: float" << std::endl;
retval = test<NumericT>(epsilon);
if ( retval == EXIT_SUCCESS )
std::cout << "# Test passed" << std::endl;
else
return retval;
}
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << std::endl;
#ifdef VIENNACL_WITH_OPENCL
if ( viennacl::ocl::current_device().double_support() )
#endif
{
{
typedef double NumericT;
NumericT epsilon = 1.0E-12;
std::cout << "# Testing setup:" << std::endl;
std::cout << " eps: " << epsilon << std::endl;
std::cout << " numeric: double" << std::endl;
retval = test<NumericT>(epsilon);
if ( retval == EXIT_SUCCESS )
std::cout << "# Test passed" << std::endl;
else
return retval;
}
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << std::endl;
}
#ifdef VIENNACL_WITH_OPENCL
else
std::cout << "No double precision support, skipping test..." << std::endl;
#endif
std::cout << std::endl;
std::cout << "------- Test completed --------" << std::endl;
std::cout << std::endl;
return retval;
}
viennacl::hyb_matrix
Sparse matrix class using a hybrid format composed of the ELL and CSR format for storing the nonzeros...
Definition: forwards.h:406
viennacl::sliced_ell_matrix::clear
void clear()
Resets all entries in the matrix back to zero without changing the matrix size. Resets the sparsity p...
Definition: sliced_ell_matrix.hpp:85
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
matrix_market.hpp
A reader and writer for the matrix market format is implemented here.
strided_matrix_vector_product_test
int strided_matrix_vector_product_test(Epsilon epsilon, STLVectorT &result, STLVectorT const &rhs, VCLVectorT &vcl_result, VCLVectorT &vcl_rhs)
Definition: sparse.cpp:156
diff
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
Definition: sparse.cpp:57
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...
test
int test(Epsilon const &epsilon)
Definition: sparse.cpp:318
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...
main
int main()
Definition: sparse.cpp:1076
viennacl::hyb_matrix::clear
void clear()
Resets all entries in the matrix back to zero without changing the matrix size. Resets the sparsity p...
Definition: hyb_matrix.hpp:69
viennacl::linalg::lower_tag
A tag class representing a lower triangular matrix.
Definition: forwards.h:849
viennacl::backend::finish
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
s2
viennacl::scalar< int > s2
Definition: global_variables.cpp:58
s1
viennacl::scalar< float > s1
Definition: global_variables.cpp:57
viennacl::linalg::detail::max
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
viennacl::ocl::current_device
viennacl::ocl::device const & current_device()
Convenience function for returning the active device in the current context.
Definition: backend.hpp:351
coordinate_matrix.hpp
Implementation of the coordinate_matrix class.
NumericT
float NumericT
Definition: bisect.cpp:40
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
hyb_matrix.hpp
Implementation of the hyb_matrix class.
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::ell_matrix
Sparse matrix class using the ELLPACK format for storing the nonzeros.
Definition: ell_matrix.hpp:53
viennacl::vector_base::begin
iterator begin()
Returns an iterator pointing to the beginning of the vector (STL like)
Definition: vector.hpp:841
ilu.hpp
Implementations of incomplete factorization preconditioners. Convenience header file.
viennacl::linalg::upper_tag
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
viennacl::sliced_ell_matrix
Sparse matrix class using the sliced ELLPACK with parameters C, .
Definition: forwards.h:403
compressed_matrix.hpp
Implementation of the compressed_matrix class.
viennacl::compressed_compressed_matrix
A sparse square matrix in compressed sparse rows format optimized for the case that only a few rows c...
Definition: compressed_compressed_matrix.hpp:265
sliced_ell_matrix.hpp
Implementation of the sliced_ell_matrix class.
viennacl::ocl::device::double_support
bool double_support() const
ViennaCL convenience function: Returns true if the device supports double precision.
Definition: device.hpp:956
viennacl::project
matrix_range< MatrixType > project(MatrixType const &A, viennacl::range const &r1, viennacl::range const &r2)
Definition: matrix_proxy.hpp:326
ell_matrix.hpp
Implementation of the ell_matrix class.
viennacl::vector< ScalarType >
viennacl::diag
vector_expression< const matrix_base< NumericT >, const int, op_matrix_diag > diag(const matrix_base< NumericT > &A, int k=0)
Definition: matrix.hpp:895
vector_proxy.hpp
Proxy classes for vectors.
viennacl::row
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
Definition: matrix.hpp:910
compressed_compressed_matrix.hpp
Implementation of the compressed_compressed_matrix class (CSR format with a relatively small number o...
viennacl::vector_base< NumericT >::clear
void clear()
Resets all entries to zero. Does not change the size of the vector.
v2
viennacl::vector< int > v2
Definition: global_variables.cpp:61
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::linalg::norm_inf
T norm_inf(std::vector< T, A > const &v1)
Definition: norm_inf.hpp:60
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.
viennacl::vector_base::size
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
viennacl::linalg::unit_lower_tag
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:859
viennacl::compressed_matrix< NumericT >
viennacl::tools::sparse_matrix_adapter
Adapts a non-const sparse matrix type made up from std::vector<std::map<SizeT, NumericT> > to basic u...
Definition: adapter.hpp:357
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
viennacl::vector_base::end
iterator end()
Returns an iterator pointing to the end of the vector (STL like)
Definition: vector.hpp:848
common.hpp
Common routines used within ILU-type preconditioners.
viennacl::basic_slice
A slice class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Definition: forwards.h:429
scalar.hpp
Implementation of the ViennaCL scalar class.
resize_test
int resize_test(Epsilon const &epsilon)
Definition: sparse.cpp:213
viennacl::linalg::unit_upper_tag
A tag class representing an upper triangular matrix with unit diagonal.
Definition: forwards.h:864
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