/********************************************************************
*               Test program for Gauss Seidel method                *
* ----------------------------------------------------------------- *
* SAMPLE RUN:                                                       *
* (Solve a linear system by Gauss Seidel iterative method).         *
*                                                                   *
* Input file tseidel.dat contains:                                  *
*                                                                   *
* 16                                                                *
*  4 -1  0  0  -1  0  0  0   0  0  0  0   0  0  0  0                *
* -1  4 -1  0   0 -1  0  0   0  0  0  0   0  0  0  0                *
*  0 -1  4 -1   0  0 -1  0   0  0  0  0   0  0  0  0                *
*  0  0 -1  4   0  0  0 -1   0  0  0  0   0  0  0  0                *
* -1  0  0  0   4 -1  0  0  -1  0  0  0   0  0  0  0                *
*  0 -1  0  0  -1  4 -1  0   0 -1  0  0   0  0  0  0                *
*  0  0 -1  0   0 -1  4 -1   0  0 -1  0   0  0  0  0                *
*  0  0  0 -1   0  0 -1  4   0  0  0 -1   0  0  0  0                *
*  0  0  0  0  -1  0  0  0   4 -1  0  0  -1  0  0  0                *
*  0  0  0  0   0 -1  0  0  -1  4 -1  0   0 -1  0  0                *
*  0  0  0  0   0  0 -1  0   0 -1  4 -1   0  0 -1  0                *
*  0  0  0  0   0  0  0 -1   0  0 -1  4   0  0  0 -1                *
*  0  0  0  0   0  0  0  0  -1  0  0  0   4 -1  0  0                *
*  0  0  0  0   0  0  0  0   0 -1  0  0  -1  4 -1  0                *
*  0  0  0  0   0  0  0  0   0  0 -1  0   0 -1  4 -1                *
*  0  0  0  0   0  0  0  0   0  0  0 -1   0  0 -1  4                *
*  2                                                                *
*  1                                                                *
*  1                                                                *
*  2                                                                *
*  1                                                                *
*  0                                                                *
*  0                                                                *
*  1                                                                *
*  1                                                                *
*  0                                                                *
*  0                                                                *
*  1                                                                *
*  2                                                                *
*  1                                                                *
*  1                                                                *
*  2                                                                *
*                                                                   *
* See full results in output file tseidel2.lst.                     *
*                                                                   *
*  Last number of iterations: 49                                    *
*                                                                   *
*  Solution Vector:                                                 *
*   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                                                        *
*                                                                   *
* ----------------------------------------------------------------- *
* Ref.: "Numerical algorithms with C, By Gisela Engeln-Muellges and *
*        Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].            *
*                                                                   *
*                          Visual C++ Release By J-P Moreau, Paris. *
*                                     (www.jpmoreau.fr)             *
* ----------------------------------------------------------------- *
* Use files: Basis_r.cpp, Fseidel.cpp, Vmblock.cpp                  *
*            basis.h,                  vmblock.h                    *
*-------------------------------------------------------------------*/
#include "basis.h"        //for REAL, etc.
#include "vmblock.h"      //for vmalloc, etc.

  FILE *fp;

//header file for seidel method (see Fseidel.cpp)
int seidel              /* Gauss Seidel Method with relaxation .......*/
           (
            int     crit,         /* crit = 0, 1, 2, 3 ...............*/
            int     n,            /* size of matrix ..................*/
            REAL *  mat[],        /* matrix ..........................*/
            REAL    b[],          /* Right hand side .................*/
            REAL    omega,        /* Relaxaktion coefficient .........*/
            REAL    x[],          /* solution ........................*/
            REAL    residu[],     /* Residuum vector .................*/
            int *   iter          /* # of iterations .................*/
           );


//test program
int main ()  {

  REAL **a, *b, *x, *residu;
  REAL omega = (REAL) 1.5;
  int n, rc, krit = 0, iter;
  void *vmblock;
  
  fp=fopen("tseidel2.dat","r");         
  
  if (fscanf (fp,"%d ", &n) <= 0)
  {
    LogError ("Input stream", 0,  __FILE__, __LINE__);
    return 1;
  }

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

  vmblock = vminit();
  a      = (REAL **)vmalloc(vmblock, MATRIX, n, n);
  b      = (REAL *) vmalloc(vmblock, VEKTOR, n, 0);
  x      = (REAL *) vmalloc(vmblock, VEKTOR, n, 0);
  residu = (REAL *) vmalloc(vmblock, VEKTOR, n, 0);

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

  if (ReadMat (fp,n, n, a) != 0)
  {
    LogError ("Input stream", 0,  __FILE__, __LINE__);
    return 1;
  }
  
  if (ReadVec (fp,n, b) != 0)
  {
    LogError ("Input stream", 0,  __FILE__, __LINE__);
    return 1;
  }

  fclose(fp);
                           
  // open output file                         
  fp=fopen("tseidel2.lst","w");  
  
  WriteHead (fp," Gauss Seidel method");  
  fprintf (fp," Dimension of the input matrix = %d\n\n", n);
  fprintf (fp," Input matrix:\n");
  WriteMat (fp,n, n, a);
  fprintf (fp," Second member:\n");
  WriteVec (fp,n, b);

  fprintf (fp,"\n Solution vector:\n");

  if ( (rc = seidel (krit, n, a, b, omega, x, residu, &iter)) != 0)
  {
    LogError ("seidel", rc,  __FILE__, __LINE__);
    return 1;
  }
  WriteVec (fp,n, x);

  fprintf(fp,"\n Residual vector (must be near zero):\n");
  WriteVec (fp,n, residu);
  
  fprintf(fp,"\n Number of iterations: %d\n",iter);
  WriteEnd(fp);  
  
  printf("\n\n Results in tseidel2.lst...\n\n");  
  fclose(fp);
  return (0);
}            

// J-P Moreau April, 2005