#include <stdio.h>
#include <math.h>
#include "opt.h"

test1()
{
/* 	dbl	a[3][3] = {{2, 1, 0}, {1, 3, 0}, {0, 0, 5}}; */
 	dbl	a[3][3] = {{9, 2, -3}, {2, 6, 1}, {-3, 1, 6}};
	dbl	aorg[3][3], b[3][3], c[3][3];
	
	printf("matrix A\n");
	matrixfprint(stdout, 3, 3, a);	printf("\n");
	matrixcopy(3, 3, aorg, a);
	
	matrixinverse(3, b, a, 1e-6);
/*	matrixfprint(stdout, 3, 3, a);	printf("\n"); */
	printf("B = inverse of A\n");
	matrixfprint(stdout, 3, 3, b);	printf("\n");
	printf("check if AB = I or not\n");
	matrixmatrix(3, 3, 3, c, aorg, b);
	matrixfprint(stdout, 3, 3, c);	printf("\n");
}

test11()
{
 	dbl	a[6][6] = 
	{{ 9,  2, -3,  1,  0,  0},
	 { 2,  6,  1,  0,  1,  0},
	 {-3,  1,  6,  0,  0,  1},
	 { 9,  2, -3,  2,  0,  0},
	 { 0,  1,  0,  2,  6,  1},
	 { 0,  0,  1, -3,  1,  6}};
	dbl	aorg[6][6], b[6][6], c[6][6];
	
	printf("matrix A\n");
	matrixfprint(stdout, 6, 6, a);	printf("\n");
	matrixcopy(6, 6, aorg, a);
	
	matrixinverse(6, b, a, 1e-6);
/*	matrixfprint(stdout, 3, 3, a);	printf("\n"); */
	printf("B = inverse of A\n");
	matrixfprint(stdout, 6, 6, b);	printf("\n");
	printf("check if AB = I or not\n");
	matrixmatrix(6, 6, 6, c, aorg, b);
	matrixfprint(stdout, 6, 6, c);	printf("\n");
}

test2()
{
  dbl a[12][12] =
    {
      {2, 0,  -1, 0,  0, 0,  0,  0,  -0.5, -0.5,  0, 0},
      {0, 3,   0, 0,  0, 0,  0, -1,  -0.5, -0.5,  0, 0},
      
      {-1, 0,  4, 0, -1, 0, -0.5,  0.5,  0,  0,  -0.5, -0.5},
      { 0, 0,  0, 4,  0, 0,  0.5, -0.5,  0, -1,  -0.5, -0.5},
      
      {0, 0,  -1, 0,  2, 0,  0, 0, -0.5,  0.5,  0,  0},
      {0, 0,   0, 0,  0, 3,  0, 0,  0.5, -0.5,  0, -1},
      
      {0,  0, -0.5,  0.5,  0, 0,  1.5, -0.5,  -1, 0,  0, 0},
      {0, -1,  0.5, -0.5,  0, 0, -0.5,  1.5,   0, 0,  0, 0},
      
      {-0.5, -0.5,  0,  0, -0.5,  0.5, -1, 0,  3, 0, -1, 0},
      {-0.5, -0.5,  0, -1,  0.5, -0.5,  0, 0,  0, 2,  0, 0},
      
      {0, 0,  -0.5, -0.5,  0,  0,  0, 0, -1, 0,  1.5, 0.5},
      {0, 0,  -0.5, -0.5,  0, -1,  0, 0,  0, 0,  0.5, 1.5}
    };
  dbl x[12] = {
    -1./4, -8./31,    0, -10./31,      1./4, -8./31,
     1./62, -17./62,  0, -1,	    -1./62, -17./62
  };
  dbl y[12];
  
  matrixvector(12, 12, y, a, x);
  vectorfprint(stdout, 12, x);
  vectorfprint(stdout, 12, y);
}

test3()
{
  double a[2][3] = {{1,2,3},
		    {5,6,7}};
  double b[3][2];
  
  matrixfprint(stdout,2,3,a);
  matrixtranspose(3,2,b,a);
  matrixfprint(stdout,3,2,b);
}

main()
{
  test11();
}