matrix_product_test.cpp

来自「矩阵运算源码最新版本」· C++ 代码 · 共 394 行 · 第 1/2 页

// Software License for MTL
// 
// Copyright (c) 2007 The Trustees of Indiana University. All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
// 
// This file is part of the Matrix Template Library
// 
// See also license.mtl.txt in the distribution.

#include <iostream>
#include <cmath>
#include <boost/test/minimal.hpp>

// #define MTL_HAS_BLAS
// #define MTL_USE_OPTERON_OPTIMIZATION

#include <boost/numeric/mtl/utility/glas_tag.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
#include <boost/numeric/mtl/matrix/morton_dense.hpp> 
#include <boost/numeric/mtl/matrix/transposed_view.hpp>
#include <boost/numeric/mtl/recursion/bit_masking.hpp>
#include <boost/numeric/mtl/operation/print.hpp>
#include <boost/numeric/mtl/operation/dmat_dmat_mult.hpp>
#include <boost/numeric/mtl/operation/mult.hpp>
#include <boost/numeric/mtl/operation/operators.hpp>
#include <boost/numeric/mtl/matrix/hessian_setup.hpp>
#include <boost/numeric/mtl/operation/assign_mode.hpp>
#include <boost/numeric/mtl/operation/mult_assign_mode.hpp>
#include <boost/numeric/mtl/recursion/base_case_test.hpp>

using namespace mtl;
using namespace std;  


template <typename MatrixA, typename MatrixB, typename MatrixC>
void test(MatrixA& a, MatrixB& b, MatrixC& c, const char* name)
{
    using assign::assign_sum; using assign::plus_sum; 
    using assign::minus_sum; using assign::mult_assign_mode; 
    using recursion::bound_test_static;

    hessian_setup(a, 1.0);
    hessian_setup(b, 2.0);

    std::cout << "\n" << name << "  --- calling simple mult:\n"; std::cout.flush();
    typedef gen_dmat_dmat_mult_t<>  mult_t;
    mult_t                              mult;

    mult(a, b, c);
    // cout << "correct result is:\n" << with_format(c, 5, 3);
    check_hessian_matrix_product(c, a.num_cols());

    std::cout << "\n" << name << "  --- check += :\n"; std::cout.flush();
    typedef gen_dmat_dmat_mult_t<plus_sum>  add_mult_t;
    add_mult_t add_mult;

    add_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 2.0);
    
    std::cout << "\n" << name << "  --- check -= :\n"; std::cout.flush();
    typedef gen_dmat_dmat_mult_t<minus_sum>  minus_mult_t;
    minus_mult_t minus_mult;

    minus_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 1.0);

#if 0
    std::cout << "\n" << name << "  --- calling mult with cursors and property maps:\n"; std::cout.flush();
    gen_cursor_dmat_dmat_mult_t<>  cursor_mult;

    cursor_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols());

    std::cout << "\n" << name << "  --- check += :\n"; std::cout.flush();
    gen_cursor_dmat_dmat_mult_t<plus_sum>  cursor_add_mult;

    cursor_add_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 2.0);
    
    std::cout << "\n" << name << "  --- check -= :\n"; std::cout.flush();
    gen_cursor_dmat_dmat_mult_t<minus_sum>  cursor_minus_mult; 

    cursor_minus_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 1.0);
#endif 
    std::cout << "\n" << name << "  --- calling mult with tiling:\n"; std::cout.flush();
    typedef gen_tiling_dmat_dmat_mult_t<2, 2>  tiling_mult_t;
    tiling_mult_t tiling_mult;

    tiling_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols()); 

    std::cout << "\n" << name << "  --- check += :\n"; std::cout.flush();
    typedef gen_tiling_dmat_dmat_mult_t<2, 2, plus_sum>  tiling_add_mult_t;
    tiling_add_mult_t tiling_add_mult;

    tiling_add_mult(a, b, c); 
    check_hessian_matrix_product(c, a.num_cols(), 2.0);
    
    std::cout << "\n" << name << "  --- check -= :\n"; std::cout.flush();
    typedef gen_tiling_dmat_dmat_mult_t<2, 2, minus_sum>  tiling_minus_mult_t;
    tiling_minus_mult_t tiling_minus_mult;

    tiling_minus_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 1.0);

 
    std::cout << "\n" << name << "  --- calling mult with tiling 2x2:\n"; std::cout.flush();
    typedef gen_tiling_22_dmat_dmat_mult_t<>  tiling_22_mult_t;
    tiling_22_mult_t tiling_22_mult;

    tiling_22_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols()); 

    std::cout << "\n" << name << "  --- check += :\n"; std::cout.flush();
    // typedef gen_tiling_22_dmat_dmat_mult_t<plus_sum>  tiling_22_add_mult_t;
    typedef typename mult_assign_mode<tiling_22_mult_t, plus_sum>::type   tiling_22_add_mult_t;
    tiling_22_add_mult_t tiling_22_add_mult;

    tiling_22_add_mult(a, b, c); 
    check_hessian_matrix_product(c, a.num_cols(), 2.0);
    
    std::cout << "\n" << name << "  --- check -= :\n"; std::cout.flush();
    typedef gen_tiling_22_dmat_dmat_mult_t<minus_sum>  tiling_22_minus_mult_t;
    tiling_22_minus_mult_t tiling_22_minus_mult;

    tiling_22_minus_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 1.0);


    std::cout << "\n" << name << "  --- calling mult with tiling 4x4:\n"; std::cout.flush();
    typedef gen_tiling_44_dmat_dmat_mult_t<>  tiling_44_mult_t;
    tiling_44_mult_t tiling_44_mult;

    tiling_44_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols()); 

    std::cout << "\n" << name << "  --- check += :\n"; std::cout.flush();
    typedef gen_tiling_44_dmat_dmat_mult_t<plus_sum>  tiling_44_add_mult_t;
    tiling_44_add_mult_t tiling_44_add_mult;

    tiling_44_add_mult(a, b, c); 
    check_hessian_matrix_product(c, a.num_cols(), 2.0);
    
    std::cout << "\n" << name << "  --- check -= :\n"; std::cout.flush();
    typedef gen_tiling_44_dmat_dmat_mult_t<minus_sum>  tiling_44_minus_mult_t;
    tiling_44_minus_mult_t tiling_44_minus_mult;

    tiling_44_minus_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 1.0);

    typedef gen_recursive_dmat_dmat_mult_t<add_mult_t, bound_test_static<2>, plus_sum>  recursive_add_mult_t;

    std::cout << "\n" << name << "  --- calling mult recursively:\n"; std::cout.flush();
    // The recursive functor is C= A*B but the base case must be C+= A*B !!!!!!
    // gen_recursive_dmat_dmat_mult_t<add_mult_t, bound_test_static<32> >  recursive_mult;
    typename mult_assign_mode<recursive_add_mult_t, assign_sum>::type	recursive_mult;

    recursive_mult(a, b, c);
    // cout << "recursive result is:\n" << with_format(c, 5, 3);
    check_hessian_matrix_product(c, a.num_cols()); 

    std::cout << "\n" << name << "  --- check += :\n"; std::cout.flush();
    recursive_add_mult_t   recursive_add_mult;

    recursive_add_mult(a, b, c); 
    check_hessian_matrix_product(c, a.num_cols(), 2.0);
    
    std::cout << "\n" << name << "  --- check -= :\n"; std::cout.flush();
    // gen_recursive_dmat_dmat_mult_t<minus_mult_t, bound_test_static<32>, minus_sum>  recursive_minus_mult; 
    // check assign mode substitution both on matrix and on block level
    typename mult_assign_mode<recursive_add_mult_t, minus_sum>::type	recursive_minus_mult;
    recursive_minus_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 1.0);

    std::cout << "\n" << name << "  --- calling mult recursively with tiling:\n"; std::cout.flush();
    // The recursive functor is C= A*B but the base case must be C+= A*B !!!!!!
    gen_recursive_dmat_dmat_mult_t<tiling_add_mult_t, bound_test_static<32> >  recursive_tiling_mult;

    recursive_tiling_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols()); 
 
    std::cout << "\n" << name << "  --- check += :\n"; std::cout.flush();
    gen_recursive_dmat_dmat_mult_t<tiling_add_mult_t, bound_test_static<32>, plus_sum>  recursive_tiling_add_mult;

    recursive_tiling_add_mult(a, b, c); 
    check_hessian_matrix_product(c, a.num_cols(), 2.0);
    
    std::cout << "\n" << name << "  --- check -= :\n"; std::cout.flush();
    gen_recursive_dmat_dmat_mult_t<tiling_minus_mult_t, bound_test_static<32>, minus_sum>  recursive_tiling_minus_mult; 

    recursive_tiling_minus_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 1.0);

    std::cout << "\n" << name << "  --- calling mult recursively platform specific plus tiling:\n"; std::cout.flush();
    typedef gen_platform_dmat_dmat_mult_t<plus_sum, tiling_add_mult_t> platform_tiling_add_mult_t;
    gen_recursive_dmat_dmat_mult_t<platform_tiling_add_mult_t, bound_test_static<32> >  recursive_platform_tiling_mult;