#include "math.h"
#include "stdio.h"
#include "/home/2008/nlamee/numerics/include/clapack.h"
#define M 13	/*No of blocks*/
#define N 10	/*No of iterations*/

void print(double a[M][N], char * blocks[M])
{
  int i, j;
  for (i = 0; i < M; ++i)
  {
    printf("%s [", blocks[i]);
    for (j = 0; j < N; ++j)
      printf("%10.2f ", a[i][j]);
    printf("]\n");
  }
}


void printToFile(double a[M][N], char * blocks[M], FILE *fptr)
{
  int i, j;
  for (i = 0; i < M; ++i)
  {
    fprintf(fptr, "%s [", blocks[i]);
    for (j = 0; j < N; ++j)
      fprintf(fptr, "%10.2f ", a[i][j]);
    fprintf(fptr, "]\n");
  }
}

int main(int argc, char *args[])
{
  double result[M][N];	/*result*/
  char * labels[M] = {"five", "neg", "six", "adder", "adder2", "neg2", "nine", "seven", "inv", "prod", "adder3", "adder4", "neg3"  };
/*  double 0:five = 5.0;
  double 1:neg = -5.0;
  double 2:six = 6.0;
  double 3:adder = _b1[0];
  double 4:adder2 = _b1[1];
  double 5:neg2 = _b1[2];
  double 6:nine = 9.0;
  double 7:seven = 7.0;
  double 8:inv = 1/9.0;
  double 9:prod = result[4][0]*7.0;
  double 10:adder3 = _b2[0];
  double 11:adder4 = _b2[1];
  double 12:neg3 = _b2[2];
*/
  result[0][0] = 5.0;
  result[1][0] = -5.0;
  result[2][0] = 6.0;

 /* solve the system of linear equations */
  double _a1[9] = {1,-1,0,0,1,-1,1,0,1};
  double _b1[3] = {0+6.0,0+result[1][0],0};
  int ipiv_b1[3];
  int _b1_res = clapack_dgesv(CblasRowMajor, 3, 1, _a1, 3, ipiv_b1, _b1, 3);
  if (_b1_res != 0 )
  {
    printf("Matrix:_b1 Error code:: %d; cannot solve the system of linear equations",_b1_res);
    exit(1);
  }
  result[3][0] = _b1[0];
  result[4][0] = _b1[1];
  result[5][0] = _b1[2];
  result[6][0] = 9.0;
  result[7][0] = 7.0;
  result[8][0] = 1/9.0;
  result[9][0] = result[4][0]*7.0;

 /* solve the system of linear equations */
  double _a2[9] = {1,-1,0,0,1,-1,1,0,1};
  double _b2[3] = {0+result[9][0],0+result[8][0],0};
  int ipiv_b2[3];
  int _b2_res = clapack_dgesv(CblasRowMajor, 3, 1, _a2, 3, ipiv_b2, _b2, 3);
  if (_b2_res != 0 )
  {
    printf("Matrix:_b2 Error code:: %d; cannot solve the system of linear equations",_b2_res);
    exit(1);
  }
  result[10][0] = _b2[0];
  result[11][0] = _b2[1];
  result[12][0] = _b2[2];

  int i;
  for (i = 1; i < N; ++i)
  {
/*    double 0:five = 5.0;
    double 1:neg = result[1][0];
    double 2:six = 6.0;
    double 3:adder = result[3][0];
    double 4:adder2 = result[4][0];
    double 5:neg2 = result[5][0];
    double 6:nine = 9.0;
    double 7:seven = 7.0;
    double 8:inv = result[8][0];
    double 9:prod = result[9][0];
    double 10:adder3 = result[10][0];
    double 11:adder4 = result[11][0];
    double 12:neg3 = result[12][0];
*/
    result[0][i] = 5.0;
    result[1][i] = result[1][0];
    result[2][i] = 6.0;
    result[3][i] = result[3][0];
    result[4][i] = result[4][0];
    result[5][i] = result[5][0];
    result[6][i] = 9.0;
    result[7][i] = 7.0;
    result[8][i] = result[8][0];
    result[9][i] = result[9][0];
    result[10][i] = result[10][0];
    result[11][i] = result[11][0];
    result[12][i] = result[12][0];
  }
  if (argc > 1)
  {  
    FILE *fptr = fopen(args[1],"w");
    printToFile(result, labels, fptr);
    fclose(fptr);
  }
  else
    print(result, labels);	 /*print the results*/

  return 0;
}