MatrixMult <-
function(A,B) {
    # Make sure A and B are numeric.
    if ( ( ! is.numeric(A) ) || ( ! is.numeric(B) ) ) {
        print("MatrixMult(): input must be numeric!")
        return(NaN);
    }

    # Make sure A and B are matrices.
    if ( ( ! is.matrix(A) ) || ( ! is.matrix(B) ) ) {
        print("MatrixMult(): input must be matrices!")
        return(NaN);
    }

    # Check inner dimensions of matrices.
    ncolA = ncol(A)
    nrowB = nrow(B)
    if ( ncolA != nrowB ) {
        print("MatrixMult(): improperly dimensioned matrices!")
        return(NaN);
    }

    # Create a matrix of zeros which will hold the result
    # of multiplication.  This step allocates memory that
    # the C++ code will use.
    nrowA      = nrow(A)
    ncolB      = ncol(B)
    multResult = matrix(rep(0,nrowA*ncolB),nrow=nrowA)

    # Invoke C function 'CWrapper'.
    # The first argument of .C() is the name of the C wrapper function.
    # The rest of the arguments provide pointers to six R variables, in order:
    # (1) The matrix 'multResult', which will contain the result of multiplication
    # (2) The number of rows in the first matrix.
    # (3) The number of columns in the first matrix.
    # (4) The number of columns in the second matrix.
    # (5) The actual values of the first matrix.
    # (6) The actual values of the second matrix.
    output =.C("CWrapper",
        product = as.double(multResult),
        nRows   = as.integer(nrow(A)),
        nCols   = as.integer(ncolA),
        nCols2  = as.integer(ncol(B)),
        matrix1 = as.double(A),
        matrix2 = as.double(B)
    )

    # Reshape 'product' back into an nrowA by ncolB matrix.
    # Then return the matrix.
    multResult = matrix(output$product,nrow=nrowA)	
    return(multResult)
}