/**
 * This code was placed in the public domain by its author, Niek Sanders, on 
 * August 30, 2007.
 *
 * This code implements in inplace multplication algorithm for two 4x4 matrices.
 * Given matrices A and B, it updates B with the value AB:
 *
 *   (B) <== (A)(B)
 *
 * The algorithm is from Golub and Loan's "Matrix Computations, 3rd ed" pg 23.  
 */
#include <iostream>


void inplace( const double* A, double* B ) {

    double workspace[4];

    for ( int j = 0; j < 4; j++ ) {

        workspace[0] = 0.0; 
        workspace[1] = 0.0; 
        workspace[2] = 0.0; 
        workspace[3] = 0.0; 

        for ( int k = 0; k < 4; k++ ) {
            for ( int i = 0; i < 4; i++ ) {
                workspace[i] += A[i*4+k] * B[k*4+j];
            }
        }

        for ( int i = 0; i < 4; i++ ) {
            B[i*4+j] = workspace[i];
        }
    }
}


void prettyPrint( double* M ) {

    for ( int m = 0; m < 4; m++ ) {
        for ( int n = 0; n < 4; n++ ) {
            std::cout << M[m*4+n] << " ";
        }
        std::cout << std::endl;
    }
    std::cout << std::endl;

}


int main( int argc, char* argv[] ) {

    // Define two 4x4 matrices in row-major order
    double A[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 };
    double B[] = { 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0 };

    // Compute B=AB
    inplace( A, B );

    // Display results
    std::cout << "Updated B" << std::endl;
    prettyPrint( B ); 

    return 0;
}