/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/* ***************************************************************************
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License version 2 as 
 *  published by the Free Software Foundation.
 *
 *  This program 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 *
 *  As a special exception, you may use this file as part of a free software
 *  library without restriction.  Specifically, if other files instantiate
 *  templates or use macros or inline functions from this file, or you compile
 *  this file and link it with other files to produce an executable, this
 *  file does not by itself cause the resulting executable to be covered by
 *  the GNU General Public License.  This exception does not however
 *  invalidate any other reasons why the executable file might be covered by
 *  the GNU General Public License.
 *
 ****************************************************************************
 */
 
/* 
 * NxN square matrix multiplication: C = A x B
 * Let's suppose B is memorized trasposed
 *
 * Author: Massimo Torquati <torquati@di.unipi.it>
 */
#include <vector>
#include <iostream>
#include <cmath>
#include <ff/mdf.hpp>
 
using namespace ff;
 
// it initialises one single row of the matrix A
void initA(const size_t N, size_t i, double *A) {
    for(size_t j=0;j<N;++j) 
        A[j] = (i+j); 
}
// it initialises one single column (row) of the matrix B
void initB(const size_t N, size_t i, double *B) {
    for(size_t j=0;j<N;++j) 
        B[j] = i*j; 
}
// vector multiplication
void mm(const size_t N, double *Cij, const double *A, const double *B) {
    double C=0;
    for(size_t k=0;k<N;++k) {
        C += A[k]*B[k];
    }
    *Cij = C;
}
void printarray(const double *A, const size_t N) {
    for (size_t i=0; i<N; i++) {
        for (size_t j=0; j<N; j++)
            printf("%g\t", A[i*N+j]);
        printf("\n");
    }
}
 
// used to pass all parameters needed by the taskGen function
template<typename T>
struct Parameters {
    double *A,*B,*C;
    size_t N;
    T* mdf;
};
 
void taskGen(Parameters<ff_mdf > *const P){
    std::vector<param_info> Param;
    auto mdf = P->mdf;
    double *A = P->A;
    double *B = P->B;
    double *C = P->C;
    const size_t N = P->N;
 
    for(size_t i=0;i<N; ++i) {
        double *Ai         = &A[i*N];
        const param_info _1={(uintptr_t)Ai,OUTPUT};
        Param.push_back(_1);  
        mdf->AddTask(Param, initA, N, i, Ai);
 
        Param.clear();
        double *Bi       = &B[i*N];
        const param_info _11={(uintptr_t)Bi,OUTPUT};
        Param.push_back(_11);  
        mdf->AddTask(Param, initB, N, i, Bi);
    }
 
    for(size_t i=0;i<N; ++i)
        for(size_t j=0;j<N;++j) {
            Param.clear();
            double *Cij       = &C[i*N+j];
            double *Ai        = &A[i*N];
            double *Bi        = &B[j*N];
            const param_info _1 ={(uintptr_t)Ai, INPUT};
            const param_info _2 ={(uintptr_t)Bi, INPUT};
            Param.push_back(_1); Param.push_back(_2); 
            mdf->AddTask(Param, mm, N, Cij, Ai, Bi);            
        }  
}
 
 
int main(int argc, char * argv[]) {
    if (argc<3) {
        std::cerr << "use: " << argv[0] << " nworkers size [check=0]\n";
        return -1;
    }
    bool   check    = false;  
    int    nworkers =atoi(argv[1]);
    size_t N        =atol(argv[2]);
    assert(N>0);
    if (argc==4)      check=true;  // checks result
 
    double *A = new double[N*N];
    double *B = new double[N*N];
    double *C = new double[N*N];
 
    assert(A && B && C);
 
    Parameters<ff_mdf > P;
    ff_mdf mm(taskGen, &P, 1024, nworkers);
    P.A=A,P.B=B,P.C=C,P.N=N;
    P.mdf = &mm;
 
    mm.run_and_wait_end(); 
 
    //printarray(A,N); printf("\n");
    //printarray(B,N); printf("\n");
    printarray(C,N);  printf("\n");
 
    // checking the result
    if (check) {
        double R=0;
        for(size_t i=0;i<N;++i)
            for(size_t j=0;j<N;++j) {
                for(size_t k=0;k<N;++k)
                    R += A[i*N+k]*B[j*N+k];
 
                if (abs(C[i*N+j]-R)>1e-06) {
                    std::cerr << "Wrong result\n";                    
                    return -1;
                }
                R=0;
            }
        std::cout << "OK\n";
    }
    // to check for memory leaks
    //delete [] A; delete [] B; delete [] C;
 
    return 0;
}