
/****
    Copyright (C) 1996 McGill University.
    Copyright (C) 1996 McCAT System Group.
    Copyright (C) 1996 ACAPS Benchmark Administrator
                       benadmin@acaps.cs.mcgill.ca

    This program is free software; you can redistribute it and/or modify
    it provided this copyright notice is maintained.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
****/

/************************************************************************/
/*  author :   Mikkel Damsgaard                                         */
/*             Kirsebaerhaven 85b                                        */
/*                                                                      */
/*             DK-8520 Lystrup                                          */
/*             email mikdam@daimi.aau.dk                                */
/*                                                                      */
/*  files :                                                             */
/*  Divsol.c           QRfact.h           Divsol.h           Jacobi.c   */     
/*  Jacobi.h           Triang.c           print.c            MM.c       */   
/*  Triang.h           print.h            MM.h               QRfact.c   */   
/*  main.c             main.h                                           */
/*                                                                      */
/*  It calculates the eigenvalues for 4 different matrixes. It does not */
/*  take any input; those 4 matrixes are calculated by MakeMatrix       */
/*  function. Output is given as 4 files: val2, val3, val4, val5, that  */
/*  contains the eigenvalues for each of the matrixes.                  */
/*                                                                      */
/************************************************************************/
#include "main.h"
#include "QRfact.h"
#include "Triang.h"
#include "Jacobi.h"
#include <stdio.h>
#include <string.h>
Matrix A,Q,U;

int comp(const double *a,const double *b)
{
  if (fabs(*a)<fabs(*b)) return 1;
  else if (fabs(*a)>fabs(*b)) return -1;
  else return 0;
}

int main ()
{
  double a,b,c,d;
  int i,j,k,l,m;
  Vector v,u,z,w;
  Matrix V,T,X,Z;
  FILE *vec;
  char filename[20],num[3];

  for (l=2;l<=5;l++) 
    {
      strcpy(filename,"val");
      sprintf(num,"%i\0",l);
      strcat(filename,num);
      /* printf("filename = %s\n",filename); */

      A = MakeMatrix(l);
      /*      if (l==5) printf("%e\n",A[0][5]); */


      /* Bauer-Fike */
      /*      printf("%i : Norm af E = %e\n",l,2*A[0][l]); */

      /*      U = Trianglelise(A,l); */
      U = Jacobi(A,l); 
      QRiterate(A,U); 

      v = newVector();
      for (i=0;i<n;i++)
	v[i] = A[i][i];
      
      qsort(v,n,sizeof(double),(int (*)(const void*,const void*))comp);
      
      for (i=0;i<n;i++)
	fprintf(stdout,"%i %e\n",i,v[i]);

      /*      printf("sigma1 = %e, sigman = %e, k2 = %e\n",
	      v[0],v[n-1],v[0]/v[n-1]); */

      /*      printf("Kondition af U: %e\n",NormOne(U)*NormInf(U)); */

      
      /*      Check(A,U,l);      */
      
      if (l==6)
	{ /* Get the egenvector for the 2 largest and 
	     the 2 smallest egenvalues */
	  
	  for (i=0;i<n;i++)
	    {
/*	      j = (i>=2)?n-i+1:i; */

	      k=0;
	      while (v[i]!=A[k][k]) k++;

	      strcpy(filename,"vec");
	      sprintf(num,"%i\0",i);
	      strcat(filename,num);
	      printf("filename = %s\n",filename);
	      
	      for (m=0;m<n;m++)
		fprintf(stdout,"%i %e\n",m,U[m][k]);
	    }

	}

      freeMatrix(U);
      freeMatrix(A);
    }
  return 0;
}

void Check(Matrix A, Matrix U, int l)
{  /* TEST SUITE */

  Matrix X,T;
  double a;
  int i,j;

  X = newMatrix();
  T = MakeMatrix(l);
  

  /* Beregn U^TXU */
  matrixMult(X,T,U); 
  
  matrixTranspose(U);
  
  matrixMult(T,U,X); 

  matrixTranspose(U);
  a = 0.0;
  for (i=0;i<n;i++)
    for (j=0;j<n;j++)
      a += (A[i][j]-T[i][j])*(A[i][j]-T[i][j]);
  
  printf("Step: %i !! The frobenius norm of X-T is %e\n",l,sqrt(a));
  
  a = 0.0;
  for (i=0;i<n;i++)
    for (j=i+1;j<n;j++)
      a += fabs(A[i][j]-A[j][i]);
  
  printf("Is A symmetric? %e\n",a);

  /*
     printMatrix(A);
     
     printMatrix(T);
     
     printMatrix(U); */


  printf("\n\n");
  freeMatrix(X);
  freeMatrix(T);
}