/* -------------------------- MODULE cg.cpp -------------------------- */

/***********************************************************************
*                                                                      *
* Conjugate Gradient Method (CG Method)                                *
* -------------------------------------                                *
*                                                                      *
* Programming language: ANSI C                                         *
* Compiler:             Turbo C 2.0                                    *
* Computer:             IBM PS/2 70 with 80387                         *
* Sources:              [BUNS85], [SCHW], [MAES84], [BIBLI 11].        *
* Author:               Juergen Dietel, Computer Center, RWTH Aachen   *
* Date:                 7.31.1992                                      *
*                                                                      *
*             C++ Version By J-P Moreau (without dynamic allocations). *
*                           (www.jpmoreau.fr)                          *
***********************************************************************/
#include "cg.h"              /* for REAL, MAT, VEC ...................*/

#define  ZERO      0
#define  MACH_EPS  2e-16     /* convenient for double type on a PC    */

int cg_method                /* Conjugate Gradient Method ............*/
                (
                 int  n,     /* Size of the linear system ............*/
                 MAT  a,     /* System matrix ........................*/
                 VEC  y,     /* right hand side ......................*/
                 VEC  x      /* solution vector ......................*/
                )            /* Error code ...........................*/
/* original name: cg_verfahren()                                      */
/***********************************************************************
* CG solves the linear system                                          *
*                         A * X = Y                                    *
* for a symmetric, positive definite matrix A via the conjugate        *
* gradient method.                                                     *
*                                                                      *
* Input parameters:                                                    *
* =================                                                    *
* n  Size of the linear system                                         *
* a  [0..n-1,0..n-1] system matrix A. Only the upper triangle of A is  *
*    used.                                                             *
* y  [0..n-1] vector of the right hand side                            *
*                                                                      *
* Output parameters:                                                   *
* ==================                                                   *
* x  [0..n-1] vector giving the solution                               *
*                                                                      *
* Return value:                                                        *
* =============                                                        *
* = 0: all is ok                                                       *
* = 1: n < 2 or other disallowed input parameters                      *
* = 2: memory exceeded                                                 *
*                                                                      *
* global names used:                                                   *
* ==================                                                   *
* ZERO, MACH_EPS, MAT, REAL, VEC.                                      *
***********************************************************************/
{
  VEC  d,            /* [0..n-1] auxiliary vectors d and g            */
       g,            /*                                               */
       AmalD;        /* [0..n-1] auxiliary vector A * d               */
  REAL alpha,        /* coefficient                                   */
       beta,         /* coefficient                                   */
       dividend,     /* numerator and denominator of a fraction,      */
       divisor,      /* respectively, used to compute alpha, beta     */
       hilf,         /* auxiliary variables                           */
       hilf2,        /*                                               */
       abstand,      /* distance of two successive approximations     */
                     /* for the solution vector x (taken in the       */
                     /* euclidean norm)                               */
       xnorm;        /* euklidean norm of x                           */
  int  k, i, j;      /* loop variables                                */


  if (n < 2)                              /* invalid parameter?       */
    return 1;

  for (i = n - 1; i >= 0; i--)            /* Matrix A empty ?         */
    if (a[i] == NULL)                     /* disallowed  !            */
      return 1;

  /*------------------------------------------------------------------*/
  /* start with x at the origin                                       */
  /*------------------------------------------------------------------*/

  for (i = n - 1; i >= 0; i--)
    x[i] = ZERO;


  /*------------------------------------------------------------------*/
  /* initialize  d and g :                                            */
  /* d = -g = -(a*x - y) = y (since x = 0)                            */
  /*------------------------------------------------------------------*/

  for (i = n - 1; i >= 0; i--)
    hilf = y[i],
    d[i] = hilf,
    g[i] = -hilf;



  /*------------------------------------------------------------------*/
  /* perform at most n steps of the CG Method                         */
  /*------------------------------------------------------------------*/

  for (k = n; k > 0; k--)
  {


    /*----------------------------------------------------------------*/
    /* compute new alpha:                                             */
    /* alpha = -(d(transp) * g) / (d(transp) * (a * d))               */
    /*----------------------------------------------------------------*/


    dividend = ZERO;
    divisor  = ZERO;

    for (i = n - 1; i >= 0; i--)
    {
      dividend += d[i] * g[i];
      for (j = 0, hilf = ZERO; j < i; j++)
        hilf += a[j][i] * d[j];
      for (j = i; j < n; j++)
        hilf += a[i][j] * d[j];
      AmalD[i] = hilf;
      divisor += d[i] * hilf;
    }

    if (divisor == ZERO)
    {
      return 0;
    }

    alpha = -dividend / divisor;


    /*----------------------------------------------------------------*/
    /* compute the norm of x und  alpha * d  and find a new x:        */
    /* x  =  x + alpha * d, then check whether x is close enough,     */
    /* in order to stop the process before n complete steps           */
    /*----------------------------------------------------------------*/

    xnorm   = ZERO;
    abstand = ZERO;

    for (i = n - 1; i >= 0; i--)
      hilf    =  x[i],
      xnorm   += hilf*hilf,
      hilf2   =  alpha * d[i],
      abstand += hilf2*hilf2,
      x[i]    =  hilf + hilf2;

    if (abstand < MACH_EPS * xnorm)
    {
      return 0;
    }


    /*----------------------------------------------------------------*/
    /* compute new g:   g  =  g + alpha * (a * d)                     */
    /*----------------------------------------------------------------*/

    for (i = n - 1; i >= 0; i--)
      g[i] += alpha * AmalD[i];


    /*----------------------------------------------------------------*/
    /* compute new beta :                                             */
    /* beta = (g(transp) * (a * d)) / (d(transp) * (a * d))           */
    /*----------------------------------------------------------------*/

    dividend = ZERO;

    for (i = n - 1; i >= 0; i--)
      dividend += g[i] * AmalD[i];

    beta = dividend / divisor;


    /*----------------------------------------------------------------*/
    /* compute new d :   d  =  - g + beta * d                         */
    /*----------------------------------------------------------------*/

    for (i = n - 1; i >= 0; i--)
      d[i] = -g[i] + beta * d[i];


  }

  return 0;
}

/* ---------------------------- END cg.cpp -------------------------- */