/***********************************************************************
* This is a test program for the function cg_method() from the module  *
* cg.cpp to solve a linear system                                      *
*                       A * X  =  Y                                    *
* for a symmetric positive definite matrix A using the conjugate       *
* gradient method.                                                     *
* -------------------------------------------------------------------- *
* Ref.: "Numerical algorithms with C, By Gisela Engeln-Muellges and    *
*        Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].               *
* -------------------------------------------------------------------- *                                                                       *
* SAMPLE RUN:                                                          *
*                                                                      *
* Solve the following linear system:                                   *
*                                                                      *
* 4                     size N of matrix                               *
* 5.0 -1.0 -1.0 -1.0    the upper triangle (with diagonal)             *
*      5.0 -1.0 -1.0    of the positive definite matrix A              *
*           5.0 -1.0                                                   *
*                5.0                                                   *
* 1.0                   right hand side of the linear system           * 
* 1.0                                                                  *
* 1.0                                                                  * 
* 1.0                                                                  *
*                                                                      *
* -----------------------------------------------------                *
*  Symmetric Linear System (Conjugate Gradient Method)                 *
* -----------------------------------------------------                *
*  Test data (Matrix and right hand side):                             *
*    5.00000  -1.00000  -1.00000  -1.00000     1.00000                 *
*   -1.00000   5.00000  -1.00000  -1.00000     1.00000                 *
*   -1.00000  -1.00000   5.00000  -1.00000     1.00000                 *
*   -1.00000  -1.00000  -1.00000   5.00000     1.00000                 *
*                                                                      *
*  Solution for the linear system:                                     *
*    0.50000   0.50000   0.50000   0.50000                             *
* -----------------------------------------------------                *
*                                                                      *
*        C++ version without dynamic allocations by J-P Moreau, Paris. *
*                            (www.jpmoreau.fr)                         *
***********************************************************************/
#include <stdio.h>
#include <math.h>   

#include "cg.cpp"       /*  for  cg_method() (with header)            */


int main()  {

  MAT  a;         /* the upper triangle of a positive definite real   */
                  /* matrix.                                          */
  VEC  y,         /* right hand side of the linear system             */
       x;         /* solution vector of the system                    */
  int  n,         /* size of matrix A                                 */
       i,         /* row index                                        */
       j,         /* column index                                     */
       fehler;    /* return value of cg_method()                      */

/* ----------------------------- print header --------------------- */
  
  printf("-----------------------------------------------------\n");
  printf(" Symmetric Linear System (Conjugate Gradient Method) \n");
  printf("-----------------------------------------------------\n");
  
/* ----------------------------- Set inputs ----------------------- */
  n=4;                                 /* size of system to solve   */

  /* define a matrix */
  a[0][0]=5.0; a[0][1]=-1.0; a[0][2]=-1.0; a[0][3]=-1.0;
               a[1][1]= 5.0; a[1][2]=-1.0; a[1][3]=-1.0;
                             a[2][2]= 5.0; a[2][3]=-1.0;
			                   a[3][3]= 5.0;

  /*define b vector */
  y[0]=1.0; y[1]=1.0; y[2]=1.0; y[3]=1.0; 

  /* -------------- Print input data as a safeguard ---------------- */

  printf(" Test data (Matrix and right hand side):\n");
  for (i = 0; i < n; i++)
  {
    for (j = 0; j < i; j++)
      printf(" %9.5f", a[j][i]);
    for (j = i; j < n; j++)
      printf(" %9.5f", a[i][j]);
    printf("   %9.5f\n", y[i]);
  }  

  /* ---------------- solve linear system --------------------------- */

  fehler = cg_method(n, a, y, x);              /* perform  CG method  */


  /* -------------------- Print output results ---------------------- */

  if (fehler)                                  /* Error in CG method? */
  {
    printf("\n Error in executing CG method.\n\n");
    return 4 + fehler;
  }

  printf("\n Solution for the linear system:\n");
  for (i = 0; i < n; i++)
    printf(" %9.5f", x[i]); 
  printf("\n");      
    
  printf("-----------------------------------------------------\n\n");
  return 0;
}

// end of file cgtst1.cpp