#include "RMat.h"
#include "R.h"     // R functions
#include "Rmath.h" // Rmath

//-----------------------------------------------------------------------------
//- This file defines TWO functions.
//
//  The first one is a C++ function named 'CppMultiply' that directly accesses
//  the C++ RMat class.
//
//  The second one is a C function named 'CWrapper', that INdirectly accesses
//  the C++ RMat class, via a call to 'CppMultiply'.

//-----------------------------------------------------------------------------
//- CppMultiply
//
//  C++ function which uses the RMat class defined in RMat.h and RMat.cpp.
//  Since it is a C++ function, it cannot be called directly from R,
//  but instead will be called through an intermediary C function,
//  CppMultiply, defined below.
//
//  'result' is a raw pointer to the result of matrix multiplication.
//  'numRows' and 'numCols' is the number of rows and columns of the first
//  matrix; 'numCols2' is the number of columns of the second matrix.
//  The number of rows of the second matrix must be the number of
//  columns of the first matrix; this will be checked at the level
//  of the R code.
//  'mat1' and 'mat2' are raw pointers to the matrix elements of the two
//   matrices.
//
void CppMultiply(double *result,
                   int *numRows, int *numCols, int *numCols2,
                   double *mat1, double *mat2)
{
    RMat  vMat1(mat1,*numRows,*numCols);
    RMat  vMat2(mat2,*numCols,*numCols2);
    
    //- Compute matrix product using overloaded operator,
    //  just to demonstrate that we get the same answer.
    //
    RMat  vMat3 = vMat1 * vMat2;
    vMat3.RPrint("CppMultiply(): way down in the C++ code,\nvMat3 =\n");

    //- Use Constructor #3 to point vMat4's mValues pointer to the memory
    //  location provided by R.  Note that this constructor assumes that
    //  memory allocation/deallocation will be managed up in the R level,
    //  not down here in the C++ code.
    //
    //  Then, "manually" compute the matrix product.  We couldn't use the
    //  overloaded multiplication operator as we did with vMat3 above,
    //  because the first thing that the overloaded assignment operator '='
    //  does is deallocate memory.  This means that the memory location that
    //  R provided would probably become corrupted or inaccessible.
    //  (In any case, you're not supposed to use the regular C++ 'delete'
    //  statement with memory that was allocated by R.)  This is still using
    //  C++ stuff -- we're using the overloaded '()' operator.
    //
    RMat  vMat4(result,*numRows,*numCols2);

    vMat4 = 0.0; // Initialize to zero.
    for (int i=1; i<=*numRows; i++)
    {
        for (int j=1; j<=*numCols2; j++)
        {
            for (int k=1;k<=*numCols;k++)
            { 
                vMat4(i,j) += vMat1(i,k) * vMat2(k,j); 
            } 
        }
    }
    vMat4.RPrint("CppMultiply(): way down in the C++ code,\nvMat4 =\n");

    //- When program execution finishes the above three nested FOR loops,
    //  vMat4 will contain the result of multiplying the two matrices.
    //  When program execution exits this function vMat4 will go out of scope,
    //  but the result of multiplication will remain in memory.
    //  Note that the RMat destructor does NOT deallocate memory; if it did,
    //  result of multiplication would be lost when vMat4 goes out of scope!
    //  Not to mention the problem of deallocating memory in vMat1 and vMat2.
}

//-----------------------------------------------------------------------------
//- CWrapper
//
//  C function that in turn invokes the above C++ function 'CppMultiply'.
//  R can access C code but can't access C++ code directly.
//  This C function provides a C interface to the C++ code that R can access.
//  See: http://www.parashift.com/c++-faq-lite/mixing-c-and-cpp.html
//  In this C function, you must NOT include class declarations,
//  instantiate any C++ objects, or do any oher obvious C++
//  things.
//  The EXTERN statement tells the C++ compiler that the
//  enclosed function 'CWrapper' is a C function.
//  Although apparently we can insert C++ style comments and
//  we can even declare variables in the middle of the function,
//  which I thought you can't do in regular C.
//
//  'result' is a raw pointer to the result of matrix multiplication.
//  'numRows' and 'numCols' is the number of rows and columns of the first
//  matrix; 'numCols2' is the number of columns of the second matrix.
//  The number of rows of the second matrix must be the number of
//  columns of the first matrix; this will be checked at the level
//  of the R code.
//  'mat1' and 'mat2' are raw pointers to the matrix elements of the two
//   matrices.
//
extern "C" {
    void CWrapper(double *result,
                 int *numRows, int *numCols, int *numCols2,
                 double *mat1, double *mat2)
    {
        //- Invoke second function which internally can do C++ things.
        //
        CppMultiply(result,numRows,numCols,numCols2,mat1,mat2);
    }
}