// $Id: matrix.cc 374 2007-11-10 11:34:27Z cactus $ -*- c++ -*-

#include "matrix.h"

#include <stdexcept>
#include <pvm3.h>

Matrix::Matrix (unsigned int cRow_, unsigned int cCol_):
    cRow (cRow_),
    cCol (cCol_),
    rgf (cCol * cRow)
{
}
    
void Matrix::pvm_pack () const
{
    pvm_pkuint (const_cast<unsigned int*>(&cRow), 1, 1);
    pvm_pkuint (const_cast<unsigned int*>(&cCol), 1, 1);
    pvm_pkdouble (const_cast<double*>(&rgf[0]), cRow * cCol, 1);
}

Matrix Matrix::pvm_unpack ()
{
    unsigned int cRow, cCol;
    pvm_upkuint (&cRow, 1, 1);
    pvm_upkuint (&cCol, 1, 1);
    
    Matrix m (cRow, cCol);
    pvm_upkdouble (&m.rgf[0], cRow * cCol, 1);
    
    return m;
}

void Matrix::print (std::ostream &stream)
{
    for (unsigned int i = 0; i < cRow; ++i)
    {
        for (unsigned int j = 0; j < cCol; ++j)
            stream << (*this)(i, j) << '\t';
        
        stream << std::endl;
    }
}

Matrix operator* (const Matrix &mA, const Matrix &mB)
{
    if (mB.rows () != mA.cols ())
        throw std::domain_error ("Shape mismatch");
    
    Matrix mRes (mA.rows (), mB.cols ());
    for (unsigned int i = 0; i != mRes.rows (); ++i)
        for (unsigned int j = 0; j != mRes.cols (); ++j)
        {
            double f = 0;
            for (unsigned int k = 0; k != mA.cols (); ++k)
                f += mA(i, k) * mB(k, j);
            mRes(i, j) = f;
        }

    return mRes;
}

Matrix & Matrix::operator= (const Matrix &m)
{
    if (&m == this) return *this;

    cCol = m.cCol;
    cRow = m.cRow;
    rgf = m.rgf;

    return *this;
}