// * Row-wise Cholesky Factorization in C++ with Timing Routines
// 
// * Note: This code will overwrite the input matrix! Create a copy if needed.
//
// * Implemented in ANSI C++. 
//   - To compile with g++             : g++ -o chol_row chol_row.cpp -lm
//   - Intriguing performance results  : g++ -O3 -o chol_row chol_row.cpp -lm
//
// * Author: Michael Mazack, <mazack @ yahoo . com>
//
// * Date: April 27th, 2010
//
// * License: Public Domain. Redistribution and modification without 
// restriction is granted. If you find this code helpful, please let me know.

#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;

// For gettimeofday().
#include <sys/time.h>

// Timing routine.
double tictoc();

int main(int argc, char *argv[])
{
  int i, j, k, l;

  int n;

  double t;

  double **a;
  
  // Get the dimensions from the command-line if available.
  if(argc > 1)
    n = atoi(argv[1]);
  else
    n = 5;

  // Allocate the memory for the matrix.
  a = new double *[n];
  for(i = 0; i < n; i++)
    a[i] = new double [n];

  // Initialize the matrix to a tridiagonal system.
  for(i = 0; i < n; i++)
    for(j = 0; j < n; j++)
      {
	if(i == j)
	  a[i][j] = 4;
	else if(abs(i - j) == 1)
	  a[i][j] = -1;
	else
	  a[i][j] = 0;
      }

  // Start the timer.
  t = tictoc();

  for(i = 0; i < n; i++)
    {
      // Square root of the diagonal.
      a[i][i] = sqrt(a[i][i]);

      // Scale the row below. 
      for(j = (i+1); j < n; j++)
	a[i][j] = a[i][j]/a[i][i];
	
      // Subtract the matrix formed by row'*row from the
      // minor matrix.
      for(l = (i+1); l < n; l++)
	for(k = (i+1); k < n; k++)
	  a[l][k] = a[l][k] - a[l][i]*a[k][i];
    }

  // Stop the timer.
  t = tictoc() - t;

  // Display the result. Comment out to supress.
  for(i = 0; i < n; i++)
    {
      for(j = 0; j < n; j++)
	{
	  if(j >= i)
	    cout << a[i][j] << "\t\t";
	  else
	    cout << "0\t\t";
	}
      cout << endl;
    }

  // Statistical results.
  cout << "Flop: " << n*n*n/3.0 << " operations." << endl;
  cout << "Time: " << t << " seconds. " << endl;
  cout << "Flop Rate: " << 1.0e-6*(n*n*n/3.0)/t << " MFlops" << endl;

  // Collect garbage.
  for(i = 0; i < n; i++)
    delete [] a[i];
  delete [] a;

  return 0;
}

// Timing routine. Returns the number of seconds elapsed as a double.
double tictoc() 
{
  struct timeval tv;

  gettimeofday(&tv, NULL);

  // Add the number of seconds to the number of (re-scaled) microseconds.
  return (tv.tv_sec + 1.0e-6*tv.tv_usec);
}