Re: [vsipl++] SAL Solvers

[Jules Bergmann <jules@xxxxxxxxxxxxxxxx>]

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@xxxxxxxxxxxxxxxx
(650) 331-3385 x705