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-sparse.C
/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
/* Copyright (C) 2014 FFLAS-FFPACK
* Written by : Bastien Vialla <bastien.vialla@lirmm.fr>
* 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========
*/
#include "fflas-ffpack/fflas/fflas.h"
#include "fflas-ffpack/fflas/fflas_sparse.h"
#include "fflas-ffpack/utils/args-parser.h"
#include "givaro/modular-double.h"
#include "givaro/zring.h"
#include <vector>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <string>
#include <iterator>
#include <cstdlib>
#include <cstdio>
// #include <stdlib.h>
#include <sstream>
using namespace FFLAS;
using namespace FFPACK;
using namespace std;
using namespace Givaro;
template <typename T> T from_string(std::string const & s) {
std::stringstream ss(s);
T result;
ss >> result; // TODO handle errors
return result;
}
template <class PtrT> void testEq(PtrT y1, PtrT y2, uint64_t n) {
for (uint64_t i = 0; i < n; ++i) {
if (y1[i] != y2[i]) {
cout << "Error " << i << endl;
cout << y1[i] << " != " << y2[i] << endl;
break;
}
}
}
template <class MatT, class Field, class IndexT>
void test_spmv(const Field &F, IndexT *row, IndexT *col,
typename Field::Element_ptr dat, index_t rowdim, index_t coldim,
uint64_t nnz, typename Field::Element_ptr x,
typename Field::Element_ptr y, typename Field::Element beta) {
MatT matrix;
sparse_init(F, matrix, row, col, dat, rowdim, coldim, nnz);
fspmv(F, matrix, x, 1, y);
sparse_delete(matrix);
}
template <class Field, class IndexT>
void
test_spmv_sell(const Field &F, IndexT *row, IndexT *col,
typename Field::Element_ptr dat, index_t rowdim, index_t coldim,
uint64_t nnz, Sparse<Field, SparseMatrix_t::SELL> &matrix,
typename Field::Element_ptr x, typename Field::Element_ptr y,
typename Field::Element beta) {
sparse_init(F, matrix, row, col, dat, rowdim, coldim, nnz);
fspmv(F, matrix, x, 1, y);
auto tmp = fflas_new(F, rowdim);
for (size_t i = 0; i < rowdim; ++i) {
tmp[i] = y[matrix.perm[i]];
}
for (size_t i = 0; i < rowdim; ++i) {
y[i] = tmp[i];
}
sparse_delete(matrix);
fflas_delete(tmp);
}
template <class MatT, class Field, class IndexT>
void test_spmm(const Field &F, IndexT *row, IndexT *col,
typename Field::Element_ptr dat, index_t rowdim, index_t coldim,
uint64_t nnz, int blockSize, typename Field::Element_ptr x,
int ldx, typename Field::Element_ptr y, int ldy,
typename Field::Element beta) {
MatT matrix;
sparse_init(F, matrix, row, col, dat, rowdim, coldim, nnz);
fspmm(F, matrix, blockSize, x, ldx, beta, y, ldy);
sparse_delete(matrix);
}
#if 0
template <class MatT, class Field, class IndexT>
void test_pspmm(const Field &F, IndexT *row, IndexT *col,
typename Field::Element_ptr dat, index_t rowdim, index_t coldim,
uint64_t nnz, int blockSize, typename Field::Element_ptr x,
int ldx, typename Field::Element_ptr y, int ldy,
typename Field::Element beta) {
MatT matrix;
sparse_init(F, matrix, row, col, dat, rowdim, coldim, nnz);
pfspmm(F, matrix, blockSize, x, ldx, beta, y, ldy);
sparse_delete(matrix);
}
template <class MatT, class Field, class IndexT>
void test_pspmv(const Field &F, IndexT *row, IndexT *col,
typename Field::Element_ptr dat, index_t rowdim, index_t coldim,
uint64_t nnz, typename Field::Element_ptr x,
typename Field::Element_ptr y, typename Field::Element beta) {
MatT matrix;
sparse_init(F, matrix, row, col, dat, rowdim, coldim, nnz);
pfspmv(F, matrix, x, 1, y);
sparse_delete(matrix);
}
template <class Field, class IndexT>
void
test_pspmv_sell(const Field &F, IndexT *row, IndexT *col,
typename Field::Element_ptr dat, index_t rowdim, index_t coldim,
uint64_t nnz, Sparse<Field, SparseMatrix_t::SELL> &matrix,
typename Field::Element_ptr x, typename Field::Element_ptr y,
typename Field::Element beta) {
sparse_init(F, matrix, row, col, dat, rowdim, coldim, nnz);
pfspmv(F, matrix, x, 1, y);
auto tmp = fflas_new(F, rowdim, 1);
for (size_t i = 0; i < rowdim; ++i) {
tmp[i] = y[matrix.perm[i]];
}
for (size_t i = 0; i < rowdim; ++i) {
y[i] = tmp[i];
}
sparse_delete(matrix);
fflas_delete(tmp);
}
#endif
int main(int argc, char **argv) {
using Field = Modular<float>;
Field F(101);
std::string path;
index_t *row = nullptr, *col = nullptr;
typename Field::Element_ptr dat;
index_t rowdim, coldim;
uint64_t nnz;
if(argc > 1)
path = argv[1];
// path = "data/mat11.sms";
index_t * st = nullptr ;
readSmsFormat(path, F, row, col, dat, rowdim, coldim, nnz);
row = fflas_new<index_t>(nnz);
for (index_t j = 0 ; j < rowdim ; ++j) {
for (index_t k = row[j] ; k < row[j+1] ; ++k)
row[k] = j ;
}
auto x = fflas_new(F, coldim, Alignment::CACHE_LINE);
auto y = fflas_new(F, rowdim, Alignment::CACHE_LINE);
auto y1 = fflas_new(F, rowdim, Alignment::CACHE_LINE);
for (size_t i = 0; i < coldim; ++i) {
x[i] = 1;
}
for (size_t i = 0; i < rowdim; ++i) {
y[i] = 0;
y1[i] = 0;
}
/************************************************************************************
*
* SPMV
*
*************************************************************************************/
cout << "=== spmv ===" << endl;
test_spmv<Sparse<Field, SparseMatrix_t::CSR>>(F, row, col, dat, rowdim,
coldim, nnz, x, y, 1);
cout << "CSR: OK" << endl;
test_spmv<Sparse<Field, SparseMatrix_t::CSR_ZO>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
// for(size_t i = 0 ; i < 10 ; ++i)
// {
// cout << y[i] << " ";
// }
// cout << endl;
// for(size_t i = 0 ; i < 10 ; ++i)
// {
// cout << y1[i] << " ";
// }
// cout << endl;
cout << "CSR_ZO: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
test_spmv<Sparse<Field, SparseMatrix_t::COO>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
cout << "COO: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
test_spmv<Sparse<Field, SparseMatrix_t::ELL>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
cout << "ELL: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
test_spmv<Sparse<Field, SparseMatrix_t::ELL_simd>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
cout << "ELL_simd: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
test_spmv<Sparse<Field, SparseMatrix_t::CSR_HYB>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
cout << "CSR_HYB: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
test_spmv<Sparse<Field, SparseMatrix_t::HYB_ZO>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
cout << "HYB_ZO: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
Sparse<Field, SparseMatrix_t::SELL> A;
test_spmv_sell(F, row, col, dat, rowdim, coldim, nnz, A, x, y1, 1);
cout << "SELL: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
/************************************************************************************
*
* pSPMV
*
*************************************************************************************/
#if 0
cout << "=== pspmv ===" << endl;
test_pspmv<Sparse<Field, SparseMatrix_t::CSR>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
cout << "CSR: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
test_pspmv<Sparse<Field, SparseMatrix_t::ELL>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
cout << "ELL: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
test_pspmv<Sparse<Field, SparseMatrix_t::ELL_simd>>(
F, row, col, dat, rowdim, coldim, nnz, x, y1, 1);
cout << "ELL_simd: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
test_pspmv<Sparse<Field, SparseMatrix_t::CSR_HYB>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
cout << "CSR_HYB: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
test_pspmv<Sparse<Field, SparseMatrix_t::HYB_ZO>>(F, row, col, dat, rowdim,
coldim, nnz, x, y1, 1);
cout << "HYB_ZO: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
Sparse<Field, SparseMatrix_t::SELL> A1;
test_pspmv_sell(F, row, col, dat, rowdim, coldim, nnz, A1, x, y1, 1);
cout << "SELL: " << ((std::equal(y, y + rowdim, y1)) ? "OK" : "ERROR")
<< endl;
for (size_t i = 0; i < rowdim; ++i) {
y1[i] = 0;
}
#endif
// // test_spmm<Sparse<Field, SparseMatrix_t::CSR>>(F, row, col, dat,
// rowdim,
// coldim, nnz, 1, x, 1, y, 1, 1);
// // test_pspmm<Sparse<Field, SparseMatrix_t::CSR>>(F, row, col, dat,
// rowdim,
// coldim, nnz, 1, x, 1, y, 1, 1);
// // test_spmm<Sparse<Field, SparseMatrix_t::COO_ZO>>(F, row, col, dat,
// rowdim, coldim, nnz, 1, x, 1, y, 1, 1);
// // test_spmm<Sparse<Field, SparseMatrix_t::CSR>>(F, row, col, dat,
// rowdim,
// coldim, nnz, 1, x, 1, y, 1, 1);
// // test_spmv<Sparse<Field, SparseMatrix_t::ELL_ZO>>(F, row, col, dat,
// rowdim, coldim, nnz, x, y1, 1);
// for(size_t i = 0 ; i < 11 ; ++i)
// {
// cout << y[i] << " ";
// }
// cout << endl;
// for(size_t i = 0 ; i < 11 ; ++i)
// {
// cout << y1[i] << " ";
// }
// cout << endl;
// auto bb = std::equal(y, y+rowdim, y1);
// cout << ((bb) ? "CORRECT" : "ERROR") << endl;
// if(!bb)
// testEq(y, y1, rowdim);
fflas_delete(x);
fflas_delete(y);
fflas_delete(y1);
return 0;
}