https://github.com/linbox-team/fflas-ffpack
Tip revision: a7801a65e9972b71558322e43812f5a7e08bbb4d authored by Clement Pernet on 14 November 2017, 16:52:10 UTC
fix parallel transpose
fix parallel transpose
Tip revision: a7801a6
test-pfgemv.C
/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
/*
* Copyright (C) the FFLAS-FFPACK group
* Written by Clément Pernet
* Brice Boyer (briceboyer) <boyer.brice@gmail.com>
* This file is Free Software and part of FFLAS-FFPACK.
*
* ========LICENCE========
* This file is part of the library FFLAS-FFPACK.
*
* FFLAS-FFPACK is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
* ========LICENCE========
*.
*/
// #ifndef NEWINO
// #define NEWWINO
// #endif
// #define WINOTHRESHOLD 100
// #define OLD_DYNAMIC_PEELING
#define ENABLE_CHECKER_fgemv 1
#include "fflas-ffpack/fflas-ffpack-config.h"
#include <iomanip>
#include <iostream>
#include <givaro/modular.h>
#include <givaro/udl.h>
#include <recint/rint.h>
#include <givaro/givintprime.h>
#include "fflas-ffpack/utils/timer.h"
#include "fflas-ffpack/fflas/fflas.h"
#include "fflas-ffpack/utils/args-parser.h"
#include "fflas-ffpack/utils/test-utils.h"
#include <random>
#include <chrono>
#ifdef __FFLASFFPACK_USE_OPENMP
#include "fflas-ffpack/paladin/pfgemv.inl"
#include "fflas-ffpack/utils/fflas_io.h"
#include <omp.h>
#endif
#include<vector>
using namespace std;
using namespace FFPACK;
using namespace FFLAS;
using Givaro::Modular;
using Givaro::ModularBalanced;
// checks that D = alpha . C + beta . A ^ta * B ^tb
template<class Field>
bool check_MV(const Field & F,
const typename Field::Element_ptr Cd, // c0
enum FFLAS::FFLAS_TRANSPOSE & ta,
const size_t m,
const size_t k,
const typename Field::Element & alpha,
const typename Field::Element_ptr A, size_t lda,
const typename Field::Element_ptr X, size_t incX,
const typename Field::Element & beta,
const typename Field::Element_ptr Y, size_t incY)
{
bool wrong = false;
typename Field::Element_ptr D;
if (ta == FflasNoTrans){
D = fflas_new(F,m);
fassign (F, m, Cd, 1, D, 1);
for (size_t i=0; i<m; i++){
F.mulin (D[i], beta);
typename Field::Element tmp;
F.init(tmp);
for (size_t j=0; j<k; j++){
F.axpyin (tmp, A[i*lda+j], X[j*incX]);
}
F.axpyin(D[i],alpha,tmp);
}
wrong = !fequal(F, m, D, 1, Y, incY);
} else {
D = fflas_new(F,k);
fassign (F, k, Cd, 1, D, 1);
for (size_t i=0; i<k; i++){
F.mulin (D[i], beta);
typename Field::Element tmp;
F.init(tmp);
for (size_t j=0; j<m; j++){
F.axpyin (tmp, A[i+j*lda], X[j*incX]);
}
F.axpyin(D[i],alpha,tmp);
}
wrong = !fequal(F, k, D, 1, Y, incY);
}
size_t Ydim = (ta==FflasNoTrans)? m : k;
if ( wrong ){
size_t canprint = 20 ;
std::cerr<<"FAIL"<<std::endl;
std::cerr << "alpha :" << alpha<<", beta : " << beta << std::endl;
std::cerr << "m :" << m << ", k : " << k << std::endl;
std::cerr << "ldA :" << lda << ", incX : " << incX << ", incY : " << incY << std::endl;
for (size_t i=0; i<Ydim && canprint; ++i){
if (!F.areEqual( Y[i*incY], D[i] ) ) {
std::cerr<<"Error Y["<<i<<"]="<<Y[i*incY]<<" D["<<i<<"]="<<D[i]<<std::endl;
canprint--;
}
}
if (Ydim<80) {
for (size_t i=0; i<Ydim ; ++i){
if (!F.areEqual( Y[i*incY], D[i] ) )
std::cout << 'X' ;
else
std::cout << '.' ;
}
std::cout << std::endl;
}
}
FFLAS::fflas_delete (D);
return !wrong ;
}
template<class Field, class RandIter>
bool launch_MV(const Field & F,
const size_t m,
const size_t k,
const typename Field::Element alpha,
const typename Field::Element beta,
const size_t lda,
enum FFLAS::FFLAS_TRANSPOSE ta,
const size_t incX,
const size_t incY,
size_t iters,
bool par,
RandIter& G)
{
bool ok = true;
typedef typename Field::Element_ptr Element_ptr;
Element_ptr A ;
for(size_t i = 0;i<iters;++i){
FFLASFFPACK_check(lda >= k);
A = FFLAS::fflas_new (F, m, lda);
FFLAS::fzero(F,m,lda,A,lda);
RandomMatrix(F, m, k, A, lda, G);
size_t Xdim = (ta == FFLAS::FflasNoTrans)? k : m;
size_t Ydim = (ta == FFLAS::FflasNoTrans)? m : k;
Element_ptr X = FFLAS::fflas_new (F, Xdim, incX);
Element_ptr Y = FFLAS::fflas_new (F, Ydim, incY);
FFLAS::fzero (F, Xdim, incX, X, incX);
FFLAS::fzero (F, Ydim, incY, Y, incY);
Element_ptr D = FFLAS::fflas_new (F, Ydim);
RandomMatrix (F, Xdim, 1, X, incX, G);
RandomMatrix (F, Ydim, 1, Y, incY, G);
FFLAS::fassign (F, Ydim, Y, incY, D, 1);
//Y will be modified so keep a copy of Y as Y2
Element_ptr Y2 = FFLAS::fflas_new (F, Ydim, incY);
FFLAS::fassign (F, Ydim, Y, incY, Y2, incY);
if (par){
{
FFLAS::MMHelper<Field, FFLAS::MMHelperAlgo::Classic, FFLAS::ModeTraits<Field>, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive, FFLAS::StrategyParameter::Threads> > WH;
PAR_BLOCK{
FFLAS::pfgemv(F, ta, m,k,alpha, A,lda, X, incX, beta, Y, incY, WH);
}
}
}else{
//FFLAS::MMHelper<Field,FFLAS::MMHelperAlgo::Auto,typename FFLAS::ModeTraits<Field>::value> WH(F,nbw,FFLAS::ParSeqHelper::Sequential());
FFLAS::fgemv(F, ta, m, k,alpha, A,lda, X, incX, beta, Y, incY);
}
ok = ok && check_MV(F, D, ta, m, k,alpha, A, lda, X, incX, beta, Y, incY);
if (!ok){
FFLAS::fflas_delete (A, X, Y, Y2, D);
break;
}
FFLAS::fassign (F, Ydim, Y2, incY, D, 1);
FFLAS::fassign (F, Ydim, Y2, incY, Y, incY);
if (par){
{
FFLAS::MMHelper<Field, FFLAS::MMHelperAlgo::Classic, FFLAS::ModeTraits<Field>, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Row, FFLAS::StrategyParameter::Threads> > WH;
PAR_BLOCK{
FFLAS::pfgemv(F, ta, m,k,alpha, A,lda, X, incX, beta, Y, incY, WH);
}
}
}else{
//FFLAS::MMHelper<Field,FFLAS::MMHelperAlgo::Auto,typename FFLAS::ModeTraits<Field>::value> WH(F,nbw,FFLAS::ParSeqHelper::Sequential());
FFLAS::fgemv(F, ta, m, k,alpha, A,lda, X, incX, beta, Y, incY);
}
ok = ok && check_MV(F, D, ta, m, k,alpha, A, lda, X, incX, beta, Y, incY);
if (!ok){
FFLAS::fflas_delete (A, X, Y, Y2, D);
break;
}
FFLAS::fassign (F, Ydim, Y2, incY, D, 1);
FFLAS::fassign (F, Ydim, Y2, incY, Y, incY);
if (par){
{
FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Row,FFLAS::StrategyParameter::Grain> WH(4);;
PAR_BLOCK{
FFLAS::pfgemv(F, ta, m,k,alpha, A,lda, X, incX, beta, Y, incY, WH);
}
}
}else{
//FFLAS::MMHelper<Field,FFLAS::MMHelperAlgo::Auto,typename FFLAS::ModeTraits<Field>::value> WH(F,nbw,FFLAS::ParSeqHelper::Sequential());
FFLAS::fgemv(F, ta, m, k,alpha, A,lda, X, incX, beta, Y, incY);
}
ok = ok && check_MV(F, D, ta, m, k,alpha, A, lda, X, incX, beta, Y, incY);
FFLAS::fflas_delete (A, X, Y, Y2, D);
if (!ok){
break;
}
}
return ok ;
}
template<class Field, class RandIter>
bool launch_MV_dispatch(const Field &F,
const int mm,
const int kk,
const typename Field::Element alpha,
const typename Field::Element beta,
const size_t iters,
const bool par,
RandIter& G)
{
bool ok = true;
size_t m,k;
size_t lda,incX, incY;
size_t ld = 13 ;
{
//FFLAS::FFLAS_TRANSPOSE ta = FFLAS::FflasNoTrans ;
//if (! par) {
//if (random()%2) ta = FFLAS::FflasTrans ;
//}
if (mm<0)
m = 1+(size_t)random() % -mm;
else m = mm;
if (kk<0)
k = 1+(size_t)random() % -kk;
else k = kk;
lda = k+(size_t)random()%ld;
incX = 1+(size_t)random()%ld;
incY = 1+(size_t)random()%ld;
ok = ok && launch_MV (F, m, k, alpha,beta, lda, FflasNoTrans, incX, incY, iters, par, G);
ok = ok && launch_MV (F, m, k, alpha,beta, lda, FflasTrans, incX, incY, iters, par, G);
}
return ok ;
}
template <class Field>
bool run_with_field (Givaro::Integer q, uint64_t b, int m, int k, size_t iters, bool par, size_t seed){
bool ok = true ;
int nbit=(int)iters;
while (ok && nbit){
typedef typename Field::Element Element ;
// choose Field
Field* F= chooseField<Field>(q,b,seed);
if (F==nullptr)
return true;
std::ostringstream oss;
F->write(oss);
std::cout.fill('.');
std::cout<<"Checking ";
std::cout.width(50);
std::cout<<oss.str();
std::cout<<" ... ";
#ifdef __FFLASFFPACK_DEBUG
F->write(std::cerr) << std::endl;
#endif
typedef typename Field::Element Element ;
typename Field::RandIter R(*F,b,seed++);
typename Field::NonZeroRandIter NZR(R);
//size_t k = 0 ;
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one,F->zero,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,F->zero,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,F->zero,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one ,F->one,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,F->one,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,F->one,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one ,F->mOne,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,F->mOne,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,F->mOne,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
Element alpha,beta ;
NZR.random(alpha);
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->one ,alpha,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->zero,alpha,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,F->mOne,alpha,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,F->one ,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,F->zero,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,F->mOne,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
for (size_t j = 0 ; j < 3 ; ++j) {
R.random(alpha);
R.random(beta);
ok = ok && launch_MV_dispatch<Field>(*F,m,k,alpha,beta,iters, par, R);
//std::cout << k << "/24" << std::endl; ++k;
}
//std::cout<<std::endl;
nbit--;
if ( !ok )
//std::cout << "\033[1;31mFAILED\033[0m "<<std::endl;
std::cout << "FAILED "<<std::endl;
else
//std::cout << "\033[1;32mPASSED\033[0m "<<std::endl;
std::cout << "PASSED "<<std::endl;
delete F;
}
return ok;
}
int main(int argc, char** argv)
{
std::cout<<setprecision(17);
std::cerr<<setprecision(17);
size_t seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
size_t iters = 3 ;
Givaro::Integer q = -1 ;
uint64_t b = 0 ;
int m = -50 ;
int k = -50 ;
int nbw = -1 ;
bool loop = false;
bool p = false;
Argument as[] = {
{ 'q', "-q Q", "Set the field characteristic (-1 for random).", TYPE_INTEGER , &q },
{ 'b', "-b B", "Set the bitsize of the random characteristic.", TYPE_INT , &b },
{ 'm', "-m M", "Set the dimension of the matrix (negative values, mean, any random value between 0 and |n|).", TYPE_INT , &m },
{ 'k', "-k K", "Set the dimension of the matrix (negative values, mean, any random value between 0 and |k|).", TYPE_INT , &k },
{ 'w', "-w N", "Set the number of winograd levels (-1 for random).", TYPE_INT , &nbw },
{ 'i', "-i R", "Set number of repetitions.", TYPE_INT , &iters },
{ 'l', "-l Y/N", "run the test in an infinte loop.", TYPE_BOOL , &loop },
{ 'p', "-p Y/N", "run the parallel fgemv.", TYPE_BOOL , &p },
{ 's', "-s seed", "Set seed for the random generator", TYPE_INT, &seed },
END_OF_ARGUMENTS
};
FFLAS::parseArguments(argc,argv,as);
bool ok = true;
srand(seed);
do{
ok = ok && run_with_field<Modular<double> >(q,b,m,k,iters,p, seed);
ok = ok && run_with_field<ModularBalanced<double> >(q,b,m,k,iters,p, seed);
ok = ok && run_with_field<Modular<float> >(q,b,m,k,iters,p, seed);
ok = ok && run_with_field<ModularBalanced<float> >(q,b,m,k,iters,p, seed);
ok = ok && run_with_field<Modular<int32_t> >(q,b,m,k,iters,p, seed);
ok = ok && run_with_field<ModularBalanced<int32_t> >(q,b,m,k,iters,p, seed);
ok = ok && run_with_field<Modular<int64_t> >(q,b,m,k,iters, p, seed);
ok = ok && run_with_field<ModularBalanced<int64_t> >(q,b,m,k,iters, p, seed);
ok = ok && run_with_field<Modular<Givaro::Integer> >(q,(b?b:512_ui64),m,k,iters,p, seed);
ok = ok && run_with_field<Givaro::ZRing<Givaro::Integer> >(0,(b?b:512_ui64),m,k,iters,p, seed);
} while (loop && ok);
return !ok ;
}