All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
spai.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_SPAI_HPP
#define VIENNACL_LINALG_SPAI_HPP
/* =========================================================================
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 manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include <utility>
#include <iostream>
#include <fstream>
#include <string>
#include <algorithm>
#include <vector>
#include <math.h>
#include <map>
// ViennaCL includes
#include "viennacl/linalg/detail/spai/spai_tag.hpp"
#include "viennacl/linalg/qr.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/linalg/detail/spai/spai-dynamic.hpp"
#include "viennacl/linalg/detail/spai/spai-static.hpp"
#include "viennacl/linalg/detail/spai/sparse_vector.hpp"
#include "viennacl/linalg/detail/spai/block_matrix.hpp"
#include "viennacl/linalg/detail/spai/block_vector.hpp"
#include "viennacl/linalg/detail/spai/fspai.hpp"
#include "viennacl/linalg/detail/spai/spai.hpp"
//boost includes
#include "boost/numeric/ublas/vector.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/storage.hpp"
#include "boost/numeric/ublas/io.hpp"
#include "boost/numeric/ublas/lu.hpp"
#include "boost/numeric/ublas/triangular.hpp"
#include "boost/numeric/ublas/matrix_expression.hpp"
namespace viennacl
{
namespace linalg
{
typedef viennacl::linalg::detail::spai::spai_tag spai_tag;
typedef viennacl::linalg::detail::spai::fspai_tag fspai_tag;
//UBLAS version
template<typename MatrixType>
class spai_precond
{
public:
typedef typename MatrixType::value_type ScalarType;
typedef typename boost::numeric::ublas::vector<ScalarType> VectorType;
spai_precond(const MatrixType& A,
const spai_tag& tag): tag_(tag){
//VCLMatrixType vcl_Ap((unsigned int)A.size2(), (unsigned int)A.size1()), vcl_A((unsigned int)A.size1(), (unsigned int)A.size2()),
//vcl_At((unsigned int)A.size1(), (unsigned int)A.size2());
//UBLASDenseMatrixType dA = A;
MatrixType pA(A.size1(), A.size2());
MatrixType At;
//std::cout<<A<<std::endl;
if (!tag_.getIsRight()){
viennacl::linalg::detail::spai::sparse_transpose(A, At);
}else{
At = A;
}
pA = At;
viennacl::linalg::detail::spai::initPreconditioner(pA, spai_m_);
viennacl::linalg::detail::spai::computeSPAI(At, spai_m_, tag_);
//(At, pA, tag_.getIsRight(), tag_.getIsStatic(), (ScalarType)_tag.getResidualNormThreshold(), (unsigned int)_tag.getIterationLimit(),
//_spai_m);
}
void apply(VectorType& vec) const {
vec = viennacl::linalg::prod(spai_m_, vec);
}
private:
// variables
spai_tag tag_;
// result of SPAI
MatrixType spai_m_;
};
//VIENNACL version
template<typename ScalarType, unsigned int MAT_ALIGNMENT>
class spai_precond< viennacl::compressed_matrix<ScalarType, MAT_ALIGNMENT> >
{
typedef viennacl::compressed_matrix<ScalarType, MAT_ALIGNMENT> MatrixType;
typedef boost::numeric::ublas::compressed_matrix<ScalarType> UBLASSparseMatrixType;
typedef viennacl::vector<ScalarType> VectorType;
typedef viennacl::matrix<ScalarType> VCLDenseMatrixType;
typedef boost::numeric::ublas::vector<ScalarType> UBLASVectorType;
public:
spai_precond(const MatrixType& A,
const spai_tag& tag): tag_(tag), spai_m_(viennacl::traits::context(A))
{
viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(A).context());
viennacl::linalg::opencl::kernels::spai<ScalarType>::init(ctx);
MatrixType At(A.size1(), A.size2(), viennacl::context(ctx));
UBLASSparseMatrixType ubls_A(A.size1(), A.size2()), ubls_spai_m;
UBLASSparseMatrixType ubls_At;
viennacl::copy(A, ubls_A);
if (!tag_.getIsRight()){
viennacl::linalg::detail::spai::sparse_transpose(ubls_A, ubls_At);
}
else{
ubls_At = ubls_A;
}
//current pattern is A
//pA = ubls_At;
//execute SPAI with ublas matrix types
viennacl::linalg::detail::spai::initPreconditioner(ubls_At, ubls_spai_m);
viennacl::copy(ubls_At, At);
viennacl::linalg::detail::spai::computeSPAI(At, ubls_At, ubls_spai_m, spai_m_, tag_);
//viennacl::copy(ubls_spai_m, spai_m_);
tmp_.resize(A.size1(), viennacl::traits::context(A), false);
}
void apply(VectorType& vec) const {
tmp_ = viennacl::linalg::prod(spai_m_, vec);
vec = tmp_;
}
private:
// variables
spai_tag tag_;
// result of SPAI
MatrixType spai_m_;
mutable VectorType tmp_;
};
//
// FSPAI
//
//UBLAS version
template<typename MatrixType>
class fspai_precond
{
typedef typename MatrixType::value_type ScalarType;
typedef typename boost::numeric::ublas::vector<ScalarType> VectorType;
typedef typename boost::numeric::ublas::matrix<ScalarType> UBLASDenseMatrixType;
typedef typename viennacl::matrix<ScalarType> VCLMatrixType;
public:
fspai_precond(const MatrixType& A,
const fspai_tag& tag): tag_(tag)
{
MatrixType pA = A;
viennacl::linalg::detail::spai::computeFSPAI(A, pA, L, L_trans, tag_);
}
void apply(VectorType& vec) const
{
VectorType temp = viennacl::linalg::prod(L_trans, vec);
vec = viennacl::linalg::prod(L, temp);
}
private:
// variables
const fspai_tag & tag_;
// result of SPAI
MatrixType L;
MatrixType L_trans;
};
//
// ViennaCL version
//
template<typename ScalarType, unsigned int MAT_ALIGNMENT>
class fspai_precond< viennacl::compressed_matrix<ScalarType, MAT_ALIGNMENT> >
{
typedef viennacl::compressed_matrix<ScalarType, MAT_ALIGNMENT> MatrixType;
typedef viennacl::vector<ScalarType> VectorType;
typedef viennacl::matrix<ScalarType> VCLDenseMatrixType;
typedef boost::numeric::ublas::compressed_matrix<ScalarType> UBLASSparseMatrixType;
typedef boost::numeric::ublas::vector<ScalarType> UBLASVectorType;
public:
fspai_precond(const MatrixType & A,
const fspai_tag & tag) : tag_(tag), L(viennacl::traits::context(A)), L_trans(viennacl::traits::context(A)), temp_apply_vec_(A.size1(), viennacl::traits::context(A))
{
//UBLASSparseMatrixType ubls_A;
UBLASSparseMatrixType ublas_A(A.size1(), A.size2());
UBLASSparseMatrixType pA(A.size1(), A.size2());
UBLASSparseMatrixType ublas_L(A.size1(), A.size2());
UBLASSparseMatrixType ublas_L_trans(A.size1(), A.size2());
viennacl::copy(A, ublas_A);
//viennacl::copy(ubls_A, vcl_A);
//vcl_At = viennacl::linalg::prod(vcl_A, vcl_A);
//vcl_pA = viennacl::linalg::prod(vcl_A, vcl_At);
//viennacl::copy(vcl_pA, pA);
pA = ublas_A;
//execute SPAI with ublas matrix types
viennacl::linalg::detail::spai::computeFSPAI(ublas_A, pA, ublas_L, ublas_L_trans, tag_);
//copy back to GPU
viennacl::copy(ublas_L, L);
viennacl::copy(ublas_L_trans, L_trans);
}
void apply(VectorType& vec) const
{
temp_apply_vec_ = viennacl::linalg::prod(L_trans, vec);
vec = viennacl::linalg::prod(L, temp_apply_vec_);
}
private:
// variables
const fspai_tag & tag_;
MatrixType L;
MatrixType L_trans;
mutable VectorType temp_apply_vec_;
};
}
}
#endif
viennacl::compressed_matrix::size2
const vcl_size_t & size2() const
Returns the number of columns.
Definition: compressed_matrix.hpp:929
viennacl::compressed_matrix::size1
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:927
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
viennacl::linalg::detail::spai::sparse_transpose
void sparse_transpose(MatrixT const &A_in, MatrixT &A)
Transposition of sparse matrix.
Definition: spai.hpp:640
fspai.hpp
Implementation of FSPAI. Experimental.
viennacl::traits::size1
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:163
viennacl::linalg::detail::spai::computeSPAI
void computeSPAI(MatrixT const &A, MatrixT &M, spai_tag &tag)
Construction of SPAI preconditioner on CPU.
Definition: spai.hpp:686
viennacl::ocl::context
Manages an OpenCL context and provides the respective convenience functions for creating buffers...
Definition: context.hpp:55
viennacl::matrix
A dense matrix class.
Definition: forwards.h:375
spai.hpp
Main implementation of SPAI (not FSPAI). Experimental.
viennacl::linalg::detail::spai::spai_tag
A tag for SPAI.
Definition: spai_tag.hpp:64
sparse_vector.hpp
Implementation of a helper sparse vector class for SPAI. Experimental.
viennacl::linalg::spai_precond< viennacl::compressed_matrix< ScalarType, MAT_ALIGNMENT > >::apply
void apply(VectorType &vec) const
Application of current preconditioner, multiplication on the right-hand side vector.
Definition: spai.hpp:166
viennacl::context
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
viennacl::linalg::spai_precond::apply
void apply(VectorType &vec) const
Application of current preconditioner, multiplication on the right-hand side vector.
Definition: spai.hpp:108
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
viennacl::linalg::detail::spai::fspai_tag
A tag for FSPAI. Experimental.
Definition: fspai.hpp:71
block_matrix.hpp
Implementation of a bunch of (small) matrices on GPU. Experimental.
viennacl::linalg::fspai_precond
Implementation of the Factored SParse Approximate Inverse Algorithm for a generic, uBLAS-compatible matrix type.
Definition: spai.hpp:189
viennacl::linalg::spai_precond< viennacl::compressed_matrix< ScalarType, MAT_ALIGNMENT > >::spai_precond
spai_precond(const MatrixType &A, const spai_tag &tag)
Constructor.
Definition: spai.hpp:138
spai-static.hpp
Implementation of a static SPAI. Experimental.
viennacl::linalg::fspai_tag
viennacl::linalg::detail::spai::fspai_tag fspai_tag
Definition: spai.hpp:67
viennacl::linalg::detail::spai::spai_tag::getIsRight
bool getIsRight() const
Definition: spai_tag.hpp:94
block_vector.hpp
Implementation of a bunch of vectors on GPU. Experimental.
viennacl::linalg::spai_precond
Implementation of the SParse Approximate Inverse Algorithm for a generic, uBLAS-compatible matrix typ...
Definition: spai.hpp:75
qr.hpp
Provides a QR factorization using a block-based approach.
spai_tag.hpp
Implementation of the spai tag holding SPAI configuration parameters. Experimental.
viennacl::linalg::spai_precond::spai_precond
spai_precond(const MatrixType &A, const spai_tag &tag)
Constructor.
Definition: spai.hpp:84
viennacl::vector< ScalarType >
viennacl::linalg::spai_precond::ScalarType
MatrixType::value_type ScalarType
Definition: spai.hpp:78
viennacl::linalg::fspai_precond::fspai_precond
fspai_precond(const MatrixType &A, const fspai_tag &tag)
Constructor.
Definition: spai.hpp:201
viennacl::linalg::spai_precond::VectorType
boost::numeric::ublas::vector< ScalarType > VectorType
Definition: spai.hpp:79
viennacl::traits::context
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:40
viennacl::linalg::fspai_precond< viennacl::compressed_matrix< ScalarType, MAT_ALIGNMENT > >::fspai_precond
fspai_precond(const MatrixType &A, const fspai_tag &tag)
Constructor.
Definition: spai.hpp:250
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
viennacl::linalg::fspai_precond< viennacl::compressed_matrix< ScalarType, MAT_ALIGNMENT > >::apply
void apply(VectorType &vec) const
Application of current preconditioner, multiplication on the right-hand side vector.
Definition: spai.hpp:275
viennacl::linalg::opencl::kernels::spai::init
static void init(viennacl::ocl::context &ctx)
Definition: spai.hpp:594
viennacl::linalg::detail::spai::computeFSPAI
void computeFSPAI(MatrixT const &A, MatrixT const &PatternA, MatrixT &L, MatrixT &L_trans, fspai_tag)
Definition: fspai.hpp:313
viennacl::compressed_matrix
A sparse square matrix in compressed sparse rows format.
Definition: compressed_matrix.hpp:559
viennacl::linalg::fspai_precond::apply
void apply(VectorType &vec) const
Application of current preconditioner, multiplication on the right-hand side vector.
Definition: spai.hpp:211
viennacl::linalg::detail::spai::initPreconditioner
void initPreconditioner(SparseMatrixT const &A, SparseMatrixT &M)
Initialize preconditioner with sparcity pattern = p(A)
Definition: spai-static.hpp:154
spai-dynamic.hpp
Implementation of a dynamic SPAI. Provides the routines for automatic pattern updates Experimental...
viennacl::linalg::spai_tag
viennacl::linalg::detail::spai::spai_tag spai_tag
Definition: spai.hpp:66