matrix_product_test.cpp

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

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

#ifdef MTL_HAS_BLAS
    std::cout << "\n" << name << "  --- calling blas mult:\n"; std::cout.flush(); 
    gen_blas_dmat_dmat_mult_t<>  blas_mult;
    blas_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols()); 
#endif    

#ifdef MTL_USE_OPTERON_OPTIMIZATION
    std::cout << "\n" << name << "  --- calling platform specific mult (empty):\n"; std::cout.flush(); 
    gen_platform_dmat_dmat_mult_t<>  platform_mult;
    platform_mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols());

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

    platform_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_platform_dmat_dmat_mult_t<minus_sum>  platform_minus_mult;

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

#endif

    std::cout << "\n" << name << "  --- using mult(a, b, c) :\n"; std::cout.flush();
    mtl::mult(a, b, c);
    check_hessian_matrix_product(c, a.num_cols(), 1.0);

    std::cout << "\n" << name << "  --- called as c= a * b:\n"; std::cout.flush();
    c= a * b;
    check_hessian_matrix_product(c, a.num_cols());

    std::cout << "\n" << name << "  --- check c+= a * b:\n"; std::cout.flush();
    c+= a * b;
    check_hessian_matrix_product(c, a.num_cols(), 2.0);

    std::cout << "\n" << name << "  --- check c-= a * b:\n"; std::cout.flush();
    c-= a * b;
    check_hessian_matrix_product(c, a.num_cols(), 1.0);

    if (a.num_cols() <= 10) 
	std::cout << a << "\n" << b << "\n" << with_format(c, 4, 4) << "\n";

    if (0) {
	print_matrix_row_cursor(a); std::cout << "\n"; print_matrix_row_cursor(b); std::cout << "\n"; 
	print_matrix_row_cursor(c); std::cout << "\n"; } 
}
 



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

    std::cout << "\n\n before matrix multiplication:\n";
    std::cout << "A:\n"; print_matrix_row_cursor(a); 
    std::cout << "B:\n"; print_matrix_row_cursor(b); 
    std::cout << "C:\n"; print_matrix_row_cursor(c); std::cout << "\n"; 

    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); 
    
    std::cout << "\n\n after matrix multiplication:\n";
    std::cout << "A:\n"; print_matrix_row_cursor(a); 
    std::cout << "B:\n"; print_matrix_row_cursor(b); 
    std::cout << "C:\n"; print_matrix_row_cursor(c); std::cout << "\n"; 
}

#ifdef MTL_HAS_BLAS
extern "C" {
void dgemm_(const char* transa, const char* transb, 
	    const int* m, const int* n, const int* k,
	    const double* alpha,  const double *da,  const int* lda,
	    const double *db, const int* ldb, const double* dbeta,
	    double *dc, const int* ldc);
}

typedef dense2D<double, matrix::parameters<col_major> >        dc_t;

struct dgemm_t
{
    void operator()(const dc_t& a, const dc_t& b, dc_t& c)
    {
	int size= a.num_rows();
	double alpha= 1.0, beta= 0.0;
	dgemm_("N", "N", &size, &size, &size, &alpha, 
	       const_cast<double*>(&a[0][0]), &size, const_cast<double*>(&b[0][0]), 
	       &size, &beta, &c[0][0], &size);

    }
};


void test_blas()
{
    dense2D<double, matrix::parameters<col_major> > a(7, 7), b(7, 7), c(7, 7);
    hessian_setup(a, 1.0);
    hessian_setup(b, 2.0);
    dgemm_t()(a, b, c);

    print_matrix_row_cursor(c);
    check_hessian_matrix_product(c, a.num_cols());
    
}

#endif // MTL_HAS_BLAS

int test_main(int argc, char* argv[])
{
    // Bitmasks:
    const unsigned long morton_mask= generate_mask<true, 0, row_major, 0>::value,
	morton_z_mask= generate_mask<false, 0, row_major, 0>::value,
	doppled_32_row_mask_no_shark= generate_mask<true, 5, row_major, 0>::value,
	doppled_32_col_mask_no_shark= generate_mask<true, 5, col_major, 0>::value,
	doppled_32_row_mask= generate_mask<true, 5, row_major, 1>::value,
	doppled_32_col_mask= generate_mask<true, 5, col_major, 1>::value,
	doppled_z_32_row_mask= generate_mask<false, 5, row_major, 1>::value,
	doppled_z_32_col_mask= generate_mask<false, 5, col_major, 1>::value;
 
    unsigned size= 5; 
    if (argc > 1) size= atoi(argv[1]); 
    if (size < 2) size= 2;

    dense2D<double>                                  da(size, size-1), db(size-1, size-2), dc(size, size-2); 
    dense2D<double, matrix::parameters<col_major> >  dca(size, size-1), dcb(size-1, size-2), dcc(size, size-2);
    dense2D<float>                                   fa(size, size-1), fb(size-1, size-2), fc(size, size-2);
    dense2D<float, matrix::parameters<col_major> >   fca(size, size-1), fcb(size-1, size-2), fcc(size, size-2);
    morton_dense<double,  morton_mask>               mda(size, size-1), mdb(size-1, size-2), mdc(size, size-2);

    typedef morton_dense<double, doppled_32_row_mask_no_shark>  morton_t;
    morton_dense<double, doppled_32_row_mask_no_shark>      mrans(size, size-1), mrbns(size-1, size-2), mrcns(size, size-2);;
    morton_dense<double, doppled_32_col_mask_no_shark>      mcans(size, size-1), mcbns(size-1, size-2), mccns(size, size-2); 
    morton_dense<double, doppled_32_col_mask>      mca(size, size-1), mcb(size-1, size-2), mcc(size, size-2);
    morton_dense<double, doppled_32_row_mask>      mra(size, size-1), mrb(size-1, size-2), mrc(size, size-2);
    morton_dense<double, doppled_z_32_col_mask>    mzca(size, size-1), mzcb(size-1, size-2), mzcc(size, size-2);
    morton_dense<double, doppled_z_32_row_mask>    mzra(size, size-1), mzrb(size-1, size-2), mzrc(size, size-2);
    morton_dense<float, doppled_32_col_mask>       mcaf(size, size-1), mcbf(size-1, size-2), mccf(size, size-2);
    morton_dense<float, doppled_32_row_mask>       mraf(size, size-1), mrbf(size-1, size-2), mrcf(size, size-2);

    transposed_view<dense2D<double> > trans_db(db); 
    transposed_view<morton_t >        trans_mrbns(mrbns); 


    std::cout << "Testing different products\n";

#if 0
    test(da, trans_db, dc, "dense2D and transposed dense2D");
    test(mrans, trans_mrbns, mrcns, "hybrid with transposed matrix");
#endif

    test(da, db, dc, "dense2D");
    test(dca, dcb, dcc, "dense2D col-major");
    test(da, dcb, dc, "dense2D row x column-major");
    test(fa, fcb, fc, "dense2D float, row x column-major");
    test(da, fcb, fc, "dense2D mixed, dense and float"); 
    test(mda, mdb, mdc, "pure Morton");
    test(mca, mcb, mcc, "Hybrid col-major");
    test(mra, mrb, mrc, "Hybrid row-major");
    test(mrans, mcbns, mrcns, "Hybrid col-major and row-major, no shark tooth");
    test(mrans, mrbns, mrcns, "Hybrid row-major, no shark tooth");
    test(mraf, mcbf, mrcf, "Hybrid col-major and row-major with float");
    test(mra, mcb, mrc, "Hybrid col-major and row-major");
    test(mzra, mzcb, mzrc, "Hybrid col-major and row-major, Z-order");
    test(mra, mzcb, mzrc, "Hybrid col-major and row-major, Z and E-order");
    test(mra, dcb, mzrc, "Hybrid col-major and row-major, Z and E-order mixed with dense2D");
    test(mra, db, mrcns, "Hybric matrix = Shark * dense2D");
    test(mrans, db, mccns, "Hybric matrix (col-major) = hybrid (row) * dense2D");

    return 0;
}