[vsipl++] SAL Solvers
Jules Bergmann
jules at codesourcery.com
Fri Apr 14 14:23:24 UTC 2006
Assem Salama wrote:
> Everyone,
> This patch adds the SAL LU and Cholesky solvers to the library. This
> has support for interleaved and split complex formats.
Assem,
This is looking good. I have a couple of minor comments below. Once
you address those, please check it in.
One note: the declaration of Is_impl_chold_avail is missing a template
parameters. For this code to compile, VSIP_IMPL_USE_SAL_SOL is not
defined. Before checking this code in, please double check that
VSIP_IMPL_USE_SAL_SOL is defined and that the SAL solvers are being
exercised.
thanks,
-- Jules
>
> Index: src/vsip/impl/solver_common.hpp
> ===================================================================
> RCS file: src/vsip/impl/solver_common.hpp
> diff -N src/vsip/impl/solver_common.hpp
> --- /dev/null 1 Jan 1970 00:00:00 -0000
> +++ src/vsip/impl/solver_common.hpp 14 Apr 2006 01:14:06 -0000
> @@ -0,0 +1,57 @@
> +/* Copyright (c) 2005, 2006 by CodeSourcery, LLC. All rights reserved. */
> +
> +/** @file vsip/impl/solver_common.hpp
> + @author Assem Salama
> + @date 2005-04-13
> + @brief VSIPL++ Library: Common stuff for linear system solvers.
> +
> +*/
> +
> +#ifndef VSIP_IMPL_SOLVER_COMMON_HPP
> +#define VSIP_IMPL_SOLVER_COMMON_HPP
> +
> +namespace vsip
> +{
> +namespace impl
> +{
> +
Important! the template <...> associates with the subsequent 'struct'
(or 'class' or function, etc). It is important that the template
declaration be contiguous with the struct it templates. Otherwise,
someone reading the code might think the struct is a regular struct and
not a template struct. The comment below "Structures for availability"
should go before the 'template' decl.
> +template <typename ImplTag,
> + typename T>
> +
> +// Structures for availability
> +struct Is_lud_impl_avail
> +{
> + static bool const value = false;
> +};
> +
The Is_chold_impl_avail struct is not templated, but it should be:
> +struct Is_chold_impl_avail
> +{
> + static bool const value = false;
> +};
> +
> +// LUD solver impl class
> +template <typename T,
> + typename ImplTag>
> +class Lud_impl;
> +
> +// CHOLESKY solver impl class
> +template <typename T,
> + typename ImplTag>
> +class Chold_impl;
> +
> +// Implementation tags
> +struct Lapack_tag;
> +
> +
> +} // namespace vsip::impl
> +
> +// Common enums
> +enum mat_uplo
> +{
> + lower,
> + upper
> +};
> +
> +} // namespace vsip
> +
> +#endif
> Index: src/vsip/impl/lapack/solver_cholesky.hpp
> ===================================================================
> RCS file: src/vsip/impl/lapack/solver_cholesky.hpp
> diff -N src/vsip/impl/lapack/solver_cholesky.hpp
> --- /dev/null 1 Jan 1970 00:00:00 -0000
> +++ src/vsip/impl/lapack/solver_cholesky.hpp 14 Apr 2006 01:14:06 -0000
> @@ -0,0 +1,192 @@
> +/* Copyright (c) 2005, 2006 by CodeSourcery, LLC. All rights reserved. */
> +
> +/** @file vsip/impl/lapack/solver_cholesky.hpp
> + @author Assem Salama
> + @date 2006-04-13
> + @brief VSIPL++ Library: Cholesky Linear system solver using LAPACK.
> +
> +*/
> +
> +#ifndef VSIP_IMPL_LAPACK_SOLVER_CHOLESKY_HPP
> +#define VSIP_IMPL_LAPACK_SOLVER_CHOLESKY_HPP
> +
> +/***********************************************************************
> + Included Files
> +***********************************************************************/
> +
> +#include <algorithm>
> +
> +#include <vsip/support.hpp>
> +#include <vsip/matrix.hpp>
> +#include <vsip/impl/math-enum.hpp>
> +#include <vsip/impl/lapack.hpp>
> +#include <vsip/impl/temp_buffer.hpp>
> +#include <vsip/impl/solver_common.hpp>
> +
> +
> +/***********************************************************************
> + Declarations
> +***********************************************************************/
> +
> +namespace vsip
> +{
> +
> +namespace impl
> +{
Need to specialize Is_chold_avail for types that Lapack supports
> +
> +/// Cholesky factorization implementation class. Common functionality
> +/// for chold by-value and by-reference classes.
> +
> +template <typename T>
> +class Chold_impl<T, Lapack_tag>
> + : Compile_time_assert<blas::Blas_traits<T>::valid>
> +{
> + // BLAS/LAPACK require complex data to be in interleaved format.
> + typedef Layout<2, col2_type, Stride_unit_dense, Cmplx_inter_fmt> data_LP;
> + typedef Fast_block<2, T, data_LP> data_block_type;
> +
> + // Constructors, copies, assignments, and destructors.
> +public:
> + Chold_impl(mat_uplo, length_type)
> + VSIP_THROW((std::bad_alloc));
> + Chold_impl(Chold_impl const&)
> + VSIP_THROW((std::bad_alloc));
> +
> + Chold_impl& operator=(Chold_impl const&) VSIP_NOTHROW;
> + ~Chold_impl() VSIP_NOTHROW;
> +
> + // Accessors.
> +public:
> + length_type length()const VSIP_NOTHROW { return length_; }
> + mat_uplo uplo() const VSIP_NOTHROW { return uplo_; }
> +
> + // Solve systems.
> +public:
> + template <typename Block>
> + bool decompose(Matrix<T, Block>) VSIP_NOTHROW;
> +
> +protected:
> + template <typename Block0,
> + typename Block1>
> + bool impl_solve(const_Matrix<T, Block0>, Matrix<T, Block1>)
> + VSIP_NOTHROW;
> +
> + // Member data.
> +private:
> + typedef std::vector<T, Aligned_allocator<T> > vector_type;
> +
> + mat_uplo uplo_; // A upper/lower triangular
> + length_type length_; // Order of A.
> +
> + Matrix<T, data_block_type> data_; // Factorized Cholesky matrix (A)
> +};
You should put a definitions comment block here before starting the
member function definitions. I.e.
/****************** ...
Definitions
***************** ... */
> +
> +
> +
> +template <typename T>
> +Chold_impl<T,Lapack_tag>::Chold_impl(
> + mat_uplo uplo,
> + length_type length
> + )
> +VSIP_THROW((std::bad_alloc))
> + : uplo_ (uplo),
> + length_ (length),
> + data_ (length_, length_)
> +{
> + assert(length_ > 0);
> + assert(uplo_ == upper || uplo_ == lower);
> +}
> +
> Index: src/vsip/impl/lapack/solver_lu.hpp
> ===================================================================
> RCS file: src/vsip/impl/lapack/solver_lu.hpp
> diff -N src/vsip/impl/lapack/solver_lu.hpp
> --- /dev/null 1 Jan 1970 00:00:00 -0000
> +++ src/vsip/impl/lapack/solver_lu.hpp 14 Apr 2006 01:14:06 -0000
> @@ -0,0 +1,225 @@
> +/* Copyright (c) 2005, 2006 by CodeSourcery, LLC. All rights reserved. */
> +
> +/** @file vsip/impl/lapack/solver_lu.hpp
> + @author Assem Salama
> + @date 2006-04-13
> + @brief VSIPL++ Library: LU linear system solver using lapack.
> +
> +*/
> +
> +#ifndef VSIP_IMPL_LAPACK_SOLVER_LU_HPP
> +#define VSIP_IMPL_LAPACK_SOLVER_LU_HPP
> +
> +/***********************************************************************
> + Included Files
> +***********************************************************************/
> +
> +#include <algorithm>
> +
> +#include <vsip/support.hpp>
> +#include <vsip/matrix.hpp>
> +#include <vsip/impl/math-enum.hpp>
> +#include <vsip/impl/lapack.hpp>
> +#include <vsip/impl/temp_buffer.hpp>
> +#include <vsip/impl/solver_common.hpp>
> +
> +
> +
> +/***********************************************************************
> + Declarations
> +***********************************************************************/
> +
> +namespace vsip
> +{
> +
> +namespace impl
> +{
Need to specialize Is_lud_impl_avail for types Lapack supports.
> Index: src/vsip/impl/sal/solver_cholesky.hpp
> ===================================================================
> RCS file: src/vsip/impl/sal/solver_cholesky.hpp
> diff -N src/vsip/impl/sal/solver_cholesky.hpp
> --- /dev/null 1 Jan 1970 00:00:00 -0000
> +++ src/vsip/impl/sal/solver_cholesky.hpp 14 Apr 2006 01:14:06 -0000
> @@ -0,0 +1,274 @@
> +/* Copyright (c) 2005, 2006 by CodeSourcery, LLC. All rights reserved. */
> +
> +/** @file vsip/impl/sal/solver_cholesky.hpp
> + @author Assem Salama
> + @date 2006-04-13
> + @brief VSIPL++ Library: Cholesky linear system solver using SAL.
> +
> +*/
> +
> +#ifndef VSIP_IMPL_SAL_SOLVER_CHOLESKY_HPP
> +#define VSIP_IMPL_SAL_SOLVER_CHOLESKY_HPP
> +
> +/***********************************************************************
> + Included Files
> +***********************************************************************/
> +
> +#include <algorithm>
> +
> +#include <vsip/support.hpp>
> +#include <vsip/matrix.hpp>
> +#include <vsip/impl/math-enum.hpp>
> +#include <vsip/impl/temp_buffer.hpp>
> +#include <vsip/impl/working-view.hpp>
> +#include <vsip/impl/fns_elementwise.hpp>
> +#include <vsip/impl/solver_common.hpp>
> +
> +
> +/***********************************************************************
> + Declarations
> +***********************************************************************/
> +
> +namespace vsip
> +{
> +namespace impl
> +{
Need to specialize Is_chold_impl_avail for types that Mercury_sal_impl
supports.
> +/// Cholesky factorization implementation class. Common functionality
> +/// for chold by-value and by-reference classes.
> +
> +template <typename T>
> +class Chold_impl<T,Mercury_sal_tag>
> + : impl::Compile_time_assert<blas::Blas_traits<T>::valid>
This isn't the right assertion (Blas_traits<double>::valid is true, but
the SAL Chold_impl doesn't support double ... yet!)
Let's drop this compile_time_assert. We can craft the right one, but
Compile_time_assert's do impact compile time, and we're already using
Choose_chold_impl to select a good ImplTag.
> +/// Form Cholesky factorization of matrix A
> +///
> +/// Requires
> +/// A to be a square matrix, either
> +/// symmetric positive definite (T real), or
> +/// hermitian positive definite (T complex).
> +///
> +/// FLOPS:
> +/// real : (1/3) n^3
> +/// complex: (4/3) n^3
> +
> +template <typename T>
> +template <typename Block>
> +bool
> +Chold_impl<T,Mercury_sal_tag>::decompose(Matrix<T, Block> m)
> + VSIP_NOTHROW
> +{
> + assert(m.size(0) == length_ && m.size(1) == length_);
> +
> + data_ = m;
> + Ext_data<data_block_type> ext(data_.block());
> +
> + if(length_ > 1)
> + sal_mat_chol_dec(
> + ext.data(), // matrix A, will also store output
> + &idv_[0], // diagnal vector
^^ diagonal (spelling)
> + length_); // order of matrix A
> + return true;
> +}
--
Jules Bergmann
CodeSourcery
jules at codesourcery.com
(650) 331-3385 x705
More information about the vsipl++
mailing list