/*************************************************************
*       Program to test Cholesky Method for symmetric        *
*               and positive definite matrices               *
*                                                            *
*                        C++ version by J-P Moreau, Paris    *
*                               (www.jpmoreau.fr)            *
* ---------------------------------------------------------- *
* Reference:                                                 *
*                                                            *
*   "Numerical Algorithms with C, By Gisela Engeln-Muellges  *
*    and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].     *
* ---------------------------------------------------------- *
* SAMPLE RUN:                                                *
*                                                            *
* Input file: matsym.dat                                     *
*                                                            *
* 4                                                          *
* 5.0 -1.0 -1.0 -1.0  1.0                                    *
*      5.0 -1.0 -1.0  1.0                                    *
*           5.0 -1.0  1.0                                    *
*                5.0  1.0                                    *
*                                                            *
* Output file: matsym.lst                                    *
*                                                            *
* ---------------------------------------------------------- *
*      Solving a linear system by Cholesky Method            *
* ---------------------------------------------------------- *
*  The matrix of left hand coefficients must be a symmetric  *
*  positive definite matrix (else error code = 2)            *
*                                                            *
*  Size of system: 4                                         *
*                                                            *
*  Input Linear System:                                      *
*                                                            *
*   5.000000   -1.000000   -1.000000   -1.000000    1.000000 *
*                                                            *
*  -1.000000    5.000000   -1.000000   -1.000000    1.000000 *
*                                                            *
*  -1.000000   -1.000000    5.000000   -1.000000    1.000000 *
*                                                            *
*  -1.000000   -1.000000   -1.000000    5.000000    1.000000 *
*                                                            *
*  Solution:                                                 *
*                                                            *
*   0.500000    0.500000    0.500000    0.500000             *
* ---------------------------------------------------------- *
*                                                            *
*************************************************************/
// Program to link with basis_r.cpp,vmblock.cpp, fcholy.cpp
#include <basis.h>
#include <vmblock.h>  

#include "Fcholy.h"
                      

// Only diagonal and upper-diagonal elements are read

int main ()  {

  REAL **a, *b, *x, **org;
  REAL tmp;
  int n, i, j, cas, rc;
  void *vmblock;

  FILE *fp1, *fp2;
  char input[20], output[20], s[20];
           
  printf("\n Input data file name (without *.dat): "); 
  scanf("%s",s);
  
  strcpy(input,s);
  strcat(input,".dat");
  strcpy(output,s);
  strcat(output,".lst");

  fp1=fopen(input,"r");
  fp2=fopen(output,"w");

  WriteHead(fp2,"     Solving a linear system by Cholesky Method");
  fprintf(fp2," The matrix of left hand coefficients must be a symmetric\n");
  fprintf(fp2," positive definite matrix (else error code = 2)\n\n");

  if (fscanf (fp1,"%d", &n) <= 0)
  {
    LogError ("Input Stream", 0, __FILE__, __LINE__);
    return 1;
  }

  if (n < 1)
  {
    LogError ("n must be > 0", 0, __FILE__, __LINE__);
    return 1;
  }

  //dynamic allocations
  vmblock = vminit();
  a   = (REAL **)vmalloc(vmblock, MATRIX, n, n);
  org = (REAL **)vmalloc(vmblock, MATRIX, n, n);
  x   = (REAL *) vmalloc(vmblock, VEKTOR,  n, 0);
  b   = (REAL *) vmalloc(vmblock, VEKTOR,  n, 0);

  if (! vmcomplete(vmblock))
  {
    LogError ("No Memory", 0, __FILE__, __LINE__);
    return 1;
  }

  for (i = 0; i < n; i++)
  {
	//read left hand coefficients of line i
    for (j = i; j < n; j++)
    {
      if (fscanf (fp1,FORMAT_IN, &tmp) <= 0)
      {
        LogError ("Input Stream", 0, __FILE__, __LINE__);
        return 1;
      }
      org[i][j] = tmp;    
      org[j][i] = tmp;
    }
    //read right hand coefficient of line i
    fscanf(fp1,FORMAT_IN,&tmp); b[i]=tmp;
  }
  fclose(fp1);                     
     
  fprintf(fp2," Size of system: %d\n\n",n);
                     
  fprintf(fp2," Input Linear System:\n\n");
  for (i = 0; i < n; i++)
  {
    for (j = 0; j < n; j++)
    {
      a[i][j] = org[i][j];
      fprintf(fp2,FORMAT_126LF, a[i][j]);
    }
    fprintf(fp2,FORMAT_126LF, b[i]);
    fprintf(fp2,"\n");
  }

  cas = 0;

  //call Cholesky subroutine
  rc = choly ( cas, n, a, b, x ); 
    
  if (rc == 0)
  {
    fprintf(fp2,"\n Solution:\n\n");
    WriteVec(fp2,n,x);
  }
  else {
    fprintf(fp2,"\n ERROR choly: code=%d\n",rc);
  }

  WriteEnd(fp2);
  fclose(fp2);

  printf("\n Results in %s.\n\n",output); 
  return (0);
}

// end of file tcholy.cpp