/*********************************************************************
*        Eigenvalues and eigenvectors of a real symmetric            *
*                square matrix by JACOBI'S METHOD                    *
* ------------------------------------------------------------------ *
*   Uses function Jacobi of module ujacobi.cpp.                      *
*                                                                    *
*                                       C++ version by J-P Moreau    *
*                                           (www.jpmoreau.fr)        *
* ------------------------------------------------------------------ *
* SAMPLE RUN:                                                        *
* Input file (elpro4.dat):                                           *
*                                                                    *
*  4                                                                 *
*   4 -2 -1  0                                                       *
*      4  0 -1                                                       *
*         4 -2                                                       *
*            4                                                       *
*                                                                    *
* Output file (elpro.lst):                                           *
*                                                                    *
* ---------------------------------------------------                *
*  Eigenvalues and eigenvectors of a real, symmetric                 *
*     square matrix by Jacobi's iterative method                     *
*             (with dynamic allocations)                             *
* ---------------------------------------------------                *
*                                                                    *
*  Dimension of matrix:     4                                        *
*                                                                    *
*  Symmetric matrix is:                                              *
*                                                                    *
*     4.000000    -2.000000    -1.000000     0.000000                *
*    -2.000000     4.000000     0.000000    -1.000000                *
*    -1.000000     0.000000     4.000000    -2.000000                *
*     0.000000    -1.000000    -2.000000     4.000000                *
*                                                                    *
*  Eigenvalues:                                                      *
*     7.00000000                                                     *
*     1.00000000                                                     *
*     3.00000000                                                     *
*     5.00000000                                                     *
*                                                                    *
*  Eigenvectors (in lines):                                          *
*     1.000000    -1.000000    -1.000000     1.000000                *
*     1.000000     1.000000     1.000000     1.000000                *
*     1.000000     1.000000    -1.000000    -1.000000                *
*    -1.000000     1.000000    -1.000000     1.000000                *
*                                                                    *
*  Number of rotations:     25                                       *
* ---------------------------------------------------                *
*                                                                    *
**********************************************************************
Version with dynamic allocations                                    */

#include "ujacobi.cpp"             //for function Jacobi()


void main()  {

  int dimen;                       // Dimension of square matrix
  int nrot;                        // Number of matrix rotations
  int  i,j;                        // loop variables
  double vmax;                     // normalization value

  REAL   **Mat;                    // given square symmetric matrix
  REAL   *Eigenvalues;             // eigenvalues solution vector
  REAL   **Eigenvectors;           // eigenvectors solution matrix

  FILE *fp1, *fp2;

  void *vmblock = NULL;

  // open input/output files
  fp1 =fopen("elpro4.dat","r");    // read mode
  fp2 =fopen("elpro.lst","w");     // write mode

  //read size of matrix from input file
  fscanf(fp1,"%d", &dimen);

  //allocate vector and matrices
  vmblock = vminit();
  Mat          = (REAL **) vmalloc(vmblock, MATRIX,  dimen+1, dimen+1);  
  Eigenvalues  = (REAL *)  vmalloc(vmblock, VEKTOR,  dimen+1, 0);
  Eigenvectors = (REAL **) vmalloc(vmblock, MATRIX,  dimen+1, dimen+1);

  //read data from input file
  for (i=1; i<=dimen; i++)         // read only upper half triangle
	for (j=i; j<=dimen; j++) {
      fscanf(fp1,"%lf", &Mat[i][j]);
	  Mat[j][i]=Mat[i][j];
    }
  fclose(fp1);

  //print title and input matrix
  fprintf(fp2,"---------------------------------------------------\n");
  fprintf(fp2," Eigenvalues and eigenvectors of a real, symmetric\n");
  fprintf(fp2,"    square matrix by Jacobi's iterative method\n");
  fprintf(fp2,"             (with dynamic allocations)\n");
  fprintf(fp2,"---------------------------------------------------\n\n");
  fprintf(fp2," Dimension of matrix: %d\n\n", dimen);
  fprintf(fp2," Symmetric matrix is:\n");
  for (i=1; i<=dimen; i++) {
	for (j=1; j<=dimen; j++)  
      fprintf(fp2,"%10.6f", Mat[i][j]);
	fprintf(fp2,"\n");
  }

  Jacobi(Mat,dimen,Eigenvalues,Eigenvectors,&nrot);

  //print results to output file
  fprintf(fp2,"\n\n Eigenvalues:\n");
  for (i=1; i<=dimen; i++)
    fprintf(fp2,"%14.8f\n", Eigenvalues[i]);
  fprintf(fp2,"\n\n Eigenvectors (in lines):\n");
  for (j=1; j<=dimen; j++) {
    vmax=Eigenvectors[1][j];
	for (i=2; i<=dimen; i++)
	  if(fabs(Eigenvectors[i][j])>fabs(vmax))  vmax=Eigenvectors[i][j];
    for (i=1; i<=dimen; i++)
      fprintf(fp2,"%10.6f", Eigenvectors[i][j]/vmax);
	fprintf(fp2,"\n");
  }
  fprintf(fp2,"\n\n Number of rotations: %d\n", nrot);
  fprintf(fp2,"---------------------------------------------------\n\n");
  fclose(fp2);
  printf("\n Program done.\n");
  printf(" Results in file elpro.lst.\n\n");

  vmfree(vmblock);

}

// end of file tujacobi.cpp