/*****************************************************************
*       Solving a symmetric linear system by Gauss method        *
*                                                                *
*                               C++ version by J-P Moreau        *
*                                   (www.jpmoreau.fr)            *
* -------------------------------------------------------------- *
* Reference:                                                     *
*                                                                *
*   "ALGEBRE Algorithmes et programmes en Pascal                 *
*    de Jean-Louis Jardrin - Dunod BO-PRE 1988" [BIBLI 10].      *
* -------------------------------------------------------------- *
* SAMPLE RUN:                                                    *
*                                                                *
* Input file (matsym.dat):                                       *
*                                                                *
* 4                                                              *
* 8  2  3  12     25                                             *
*    4  7   0.25  13.25                                          *
*       3   5     18                                             *
*           2     19.25                                          *
*                                                                *
* Output file (matsym.lst):                                      *
*                                                                *
* -------------------------------------------------------------- *
*  SYMMETRIC SYSTEM TO BE SOLVED:                                *
* -------------------------------------------------------------- *
*                                                                *
*  N=4                                                           *
*                                                                *
*   8.000000     2.000000     3.000000    12.000000    25.000000 *
*   2.000000     4.000000     7.000000     0.250000    13.250000 * 
*   3.000000     7.000000     3.000000     5.000000    18.000000 * 
*  12.000000     0.250000     5.000000     2.000000    19.250000 * 
*                                                                *
* -------------------------------------------------------------- *
* System solution:                                               *
* -------------------------------------------------------------- *
*                                                                *
*  X1 =   1.000000                                               *
*  X2 =   1.000000                                               *
*  X3 =   1.000000                                               *
*  X4 =   1.000000                                               *
*                                                                *
* -------------------------------------------------------------- *
* End of file matsym.lst.                                        *
*                                                                *
*****************************************************************/
#include <basis.h>
#include <vmblock.h>  


/*******************************************************************
*  This function solves a symmetric linear system by Gauss method  *
* ---------------------------------------------------------------- *
* The input matrix includes the right hand column, only the upper  *
* triangle terms are given.                                        *
*                                                                  *
* INPUTS:                                                          *
*           eps:    desired precision (real)                       *
*             n:    size of system (integer)                       *
*             A:    extended matrix of system of type MAT          *
* OUTPUTS:                                                         *
*            it:    error indicator (0=singular system, 1=OK)      *
*             X:    system solution (x1,...,xn) of type VEC        *
* ---------------------------------------------------------------- *
* LIMITATION:       No term in main diagonal must equal zero.      *
*******************************************************************/
  void RSLSG (REAL eps, int n, REAL *A[], int *it, REAL X[]) {
    REAL q0, s;
    int i,j,k;
    *it=1; k=1;
    while ((*it != 0) && (k < n))  {
      if (fabs(A[k][k]) < eps)  *it=0;
      else {
	for (i=k+1; i<n+1; i++) {
          q0=A[k][i]/A[k][k];
	  for (j=i; j<n+2; j++) { 
	    A[i][j] = A[i][j]-q0*A[k][j];
	  }
        }
        k++;
      }
    } //while

    if ((*it==1) && (fabs(A[n][n])>=eps)) {
      X[n] = A[n][n+1]/A[n][n];
      for (i=n-1; i>0; i--) { 
        s=0.0;
        for (j=i+1; j<n+1; j++) {  
	  s = s + A[i][j]*X[j];
	}
        X[i] = (A[i][n+1]-s) / A[i][i];
      }  
    }
    else 
      *it=0;   
  } //RSLSG 


  int main()  {

    int  i,it,j,n;
    REAL **A, *X, eps, tmp;
    void *vmblock = NULL; 

    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");       

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

    eps = 1e-12;

    vmblock = vminit();  
    A       = (REAL **)vmalloc(vmblock, MATRIX,  n+1, n+2);  
    X       = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);

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

    WriteHead(fp2,"  SYMMETRIC SYSTEM TO BE SOLVED:");

    fprintf(fp2,"  N=%d\n",n);

    for (i=1; i<n+1; i++)
      for (j=i; j<n+2; j++)  {
	fscanf(fp1,"%lf",&tmp);
	A[i][j] = (REAL) tmp;
	if (j<n+1) A[j][i] = (REAL) tmp;
      }

    fclose(fp1);	

    WriteMat1 (fp2, n, n+1, A); 

    RSLSG(eps,n,A,&it,X);

    if (it==0) 
      fprintf(fp2,"  At least one diagonal term equals zero !");
    else {
      WriteHead(fp2,"  System solution:");
      for (i=1; i<n+1; i++) 
        fprintf(fp2,"  X%d = %10.6f\n",i,X[i]); 
    }
    WriteEnd(fp2);
    fprintf(fp2," End of file %s.\n",output);
    fclose(fp2);
    printf("\n Results in file %s.\n\n",output);
    return 0;

 } //main
 
// end of file syslin.cpp