# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>

# include "cg_rc.h"

/******************************************************************************/

int cg_rc ( int n, double b[], double x[], double r[], double z[], 
  double p[], double q[], int job )

/******************************************************************************/
/*
  Purpose:

    CG_RC is a reverse communication conjugate gradient routine.

  Discussion:

    This routine seeks a solution of the linear system A*x=b
    where b is a given right hand side vector, A is an n by n
    symmetric positive definite matrix, and x is an unknown vector
    to be determined.

    Under the assumptions that the matrix A is large and sparse,
    the conjugate gradient method may provide a solution when
    a direct approach would be impractical because of excessive
    requirements of storage or even of time.

    The conjugate gradient method presented here does not require the 
    user to store the matrix A in a particular way.  Instead, it only 
    supposes that the user has a way of calculating
      y = alpha * A * x + b * y
    and of solving the preconditioned linear system
      M * x = b
    where M is some preconditioning matrix, which might be merely
    the identity matrix, or a diagonal matrix containing the
    diagonal entries of A.

    This routine was extracted from the "templates" package.
    There, it was not intended for direct access by a user;
    instead, a higher routine called "cg()" was called once by
    the user.  The cg() routine then made repeated calls to 
    cgrevcom() before returning the result to the user.

    The reverse communication feature of cgrevcom() makes it, by itself,
    a very powerful function.  It allows the user to handle issues of
    storage and implementation that would otherwise have to be
    mediated in a fixed way by the function argument list.  Therefore,
    this version of cgrecom() has been extracted from the templates
    library and documented as a stand-alone procedure.

    The user sets the value of JOB to 1 before the first call,
    indicating the beginning of the computation, and to the value of
    2 thereafter, indicating a continuation call.  
    The output value of JOB is set by cgrevcom(), which
    will return with an output value of JOB that requests a particular
    new action from the user.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    12 January 2013

  Author:

    John Burkardt

  Reference:

    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
    June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
    Charles Romine, Henk van der Vorst,
    Templates for the Solution of Linear Systems:
    Building Blocks for Iterative Methods,
    SIAM, 1994,
    ISBN: 0898714710,
    LC: QA297.8.T45.

  Parameters:

    Input, int N, the dimension of the matrix.

    Input, double B[N], the right hand side vector.

    Input/output, double X[N].  On first call, the user 
    should store an initial guess for the solution in X.  On return with
    JOB = 4, X contains the latest solution estimate.

    Input/output, double R[N], Z[N], P[N], Q[N],
    information used by the program during the calculation.  The user
    does not need to initialize these vectors.  However, specific
    return values of JOB may require the user to carry out some computation
    using data in some of these vectors.

    Input/output, int JOB, communicates the task to be done.
    The user needs to set the input value of JOB to 1, before the first call,
    and then to 2 for every subsequent call for the given problem.
    The output value of JOB indicates the requested user action.  
    * JOB = 1, compute Q = A * P;
    * JOB = 2: solve M*Z=R, where M is the preconditioning matrix;
    * JOB = 3: compute R = R - A * X;
    * JOB = 4: check the residual R for convergence.  
               If satisfactory, terminate the iteration.
               If too many iterations were taken, terminate the iteration.
*/
{
  double alpha;
  double beta;
  int i;
  static int iter;
  int job_next;
  double pdotq;
  static double rho;
  static double rho_old;
  static int rlbl;
/*
  Initialization.
  Ask the user to compute the initial residual.
*/
  if ( job == 1 )
  {
    for ( i = 0; i < n; i++ )
    {
      r[i] = b[i];
    }
    job_next = 3;
    rlbl = 2;
  }
/*
  Begin first conjugate gradient loop.
  Ask the user for a preconditioner solve.
*/
  else if ( rlbl == 2 )
  {
    iter = 1;

    job_next = 2;
    rlbl = 3;
  }
/*
  Compute the direction.
  Ask the user to compute ALPHA.
  Save A*P to Q.
*/
  else if ( rlbl == 3 )
  {
    rho = 0.0;
    for ( i = 0; i < n; i++ )
    {
      rho = rho + r[i] * z[i];
    }

    if ( 1 < iter )
    {
      beta = rho / rho_old;
      for ( i = 0; i < n; i++ )
      {
        z[i] = z[i] + beta * p[i];
      }
    }

    for ( i = 0; i < n; i++ )
    {
      p[i] = z[i];
    }

    job_next = 1;
    rlbl = 4;
  }
/*
  Compute current solution vector.
  Ask the user to check the stopping criterion.
*/
  else if ( rlbl == 4 )
  {
    pdotq = 0.0;
    for ( i = 0; i < n; i++ )
    {
      pdotq = pdotq + p[i] * q[i];
    }
    alpha = rho / pdotq;
    for ( i = 0; i < n; i++ )
    {
      x[i] = x[i] + alpha * p[i];
    }
    for ( i = 0; i < n; i++ )
    {
      r[i] = r[i] - alpha * q[i];
    }
    job_next = 4;
    rlbl = 5;
  }
/*
  Begin the next step.
  Ask for a preconditioner solve.
*/
  else if ( rlbl == 5 )
  {
    rho_old = rho;
    iter = iter + 1;

    job_next = 2;
    rlbl = 3;
  }

  return job_next;
}
/******************************************************************************/

void r8mat_mv ( int m, int n, double a[], double x[], double ax[] )

/******************************************************************************/
/*
  Purpose:

    R8MAT_MV multiplies a matrix times a vector.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

    For this routine, the result is returned as an argument.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 August 2011

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns of the matrix.

    Input, double A[M,N], the M by N matrix.

    Input, double X[N], the vector to be multiplied by A.

    Output, double AX[M], the product A*X.
*/
{
  int i;
  int j;

  for ( i = 0; i < m; i++ )
  {
    ax[i] = 0.0;
    for ( j = 0; j < n; j++ )
    {
      ax[i] = ax[i] + a[i+j*m] * x[j];
    }
  }

  return;
}
/******************************************************************************/

double *r8vec_uniform_01_new ( int n, int *seed )

/******************************************************************************/
/*
  Purpose:

    R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC.

  Discussion:

    This routine implements the recursion

      seed = 16807 * seed mod ( 2^31 - 1 )
      unif = seed / ( 2^31 - 1 )

    The integer arithmetic never requires more than 32 bits,
    including a sign bit.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    19 August 2004

  Author:

    John Burkardt

  Reference:

    Paul Bratley, Bennett Fox, Linus Schrage,
    A Guide to Simulation,
    Second Edition,
    Springer, 1987,
    ISBN: 0387964673,
    LC: QA76.9.C65.B73.

    Bennett Fox,
    Algorithm 647:
    Implementation and Relative Efficiency of Quasirandom
    Sequence Generators,
    ACM Transactions on Mathematical Software,
    Volume 12, Number 4, December 1986, pages 362-376.

    Pierre L'Ecuyer,
    Random Number Generation,
    in Handbook of Simulation,
    edited by Jerry Banks,
    Wiley, 1998,
    ISBN: 0471134031,
    LC: T57.62.H37.

    Peter Lewis, Allen Goodman, James Miller,
    A Pseudo-Random Number Generator for the System/360,
    IBM Systems Journal,
    Volume 8, Number 2, 1969, pages 136-143.

  Parameters:

    Input, int N, the number of entries in the vector.

    Input/output, int *SEED, a seed for the random number generator.

    Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values.
*/
{
  int i;
  int i4_huge = 2147483647;
  int k;
  double *r;

  if ( *seed == 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8VEC_UNIFORM_01_NEW - Fatal error\n" );
    fprintf ( stderr, "  Input value of SEED = 0.\n" );
    exit ( 1 );
  }

  r = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    k = *seed / 127773;

    *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

    if ( *seed < 0 )
    {
      *seed = *seed + i4_huge;
    }

    r[i] = ( double ) ( *seed ) * 4.656612875E-10;
  }

  return r;
}
/******************************************************************************/

void timestamp ( void )

/******************************************************************************/
/*
  Purpose:

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}
/******************************************************************************/

double *wathen ( int nx, int ny, int n )

/******************************************************************************/
/*
  Purpose:

    WATHEN returns the WATHEN matrix.

  Discussion:

    The Wathen matrix is a finite element matrix which is sparse.

    The entries of the matrix depend in part on a physical quantity
    related to density.  That density is here assigned random values between
    0 and 100.

    The matrix order N is determined by the input quantities NX and NY,
    which would usually be the number of elements in the X and Y directions.

    The value of N is
      N = 3*NX*NY + 2*NX + 2*NY + 1,

    and sufficient storage in A must have been set aside to hold
    the matrix.

    A is the consistent mass matrix for a regular NX by NY grid
    of 8 node serendipity elements.  

    Here is an illustration for NX = 3, NY = 2:

     23-24-25-26-27-28-29
      |     |     |     |
     19    20    21    22
      |     |     |     |
     12-13-14-15-16-17-18
      |     |     |     |
      8     9    10    11
      |     |     |     |
      1--2--3--4--5--6--7

    For this example, the total number of nodes is, as expected,

      N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29

  Properties:

    A is symmetric positive definite for any positive values of the
    density RHO(NX,NY), which is here given the value 1.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    13 January 2013

  Author:

    John Burkardt

  Reference:

    Nicholas Higham,
    Algorithm 694: A Collection of Test Matrices in MATLAB,
    ACM Transactions on Mathematical Software,
    Volume 17, Number 3, September 1991, pages 289-305.

    Andrew Wathen,
    Realistic eigenvalue bounds for the Galerkin mass matrix,
    IMA Journal of Numerical Analysis,
    Volume 7, 1987, pages 449-457.

  Parameters:

    Input, int NX, NY, values which determine the size of A.

    Input, int N, the order of the matrix.

    Output, double WATHEN[N*N], the matrix.
*/
{
  double *a;
  static double em[8*8] = {
     6.0, -6.0,  2.0, -8.0,  3.0, -8.0,  2.0, -6.0, 
    -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, 
     2.0, -6.0,  6.0, -6.0,  2.0, -8.0,  3.0, -8.0, 
    -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, 
     3.0, -8.0,  2.0, -6.0,  6.0, -6.0,  2.0, -8.0, 
    -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, 
     2.0, -8.0,  3.0, -8.0,  2.0, -6.0,  6.0, -6.0, 
    -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0 };
  int i;
  int j;
  int kcol;
  int krow;
  int node[8];
  double rho;
  
  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      a[i+j*n] = 0.0;
    }
  }

  for ( j = 1; j <= ny; j++ )
  {
    for ( i = 1; i <= nx; i++ )
    {
/*
  For the element (I,J), determine the indices of the 8 nodes.
*/
      node[0] = 3 * j * nx + 2 * j + 2 * i;
      node[1] = node[0] - 1;
      node[2] = node[0] - 2;
      node[3] = ( 3 * j - 1 ) * nx + 2 * j + i - 2;
      node[4] = ( 3 * j - 3 ) * nx + 2 * j + 2 * i - 4;
      node[5] = node[4] + 1;
      node[6] = node[4] + 2;
      node[7] = node[3] + 1;
/*
  The density RHO can also be set to a random positive value.
*/
      for ( krow = 0; krow < 8; krow++ )
      {
        for ( kcol = 0; kcol < 8; kcol++ )
        {
          rho = 1.0;

          if (      node[krow] < 0 ||
               n <= node[krow] ||
                    node[kcol] < 0 ||
               n <= node[kcol] )
          {
            fprintf ( stderr, "\n" );
            fprintf ( stderr, "WATHEN - Fatal error\n" );
            fprintf ( stderr, "  I = %d  J = %d\n", i, j );
            fprintf ( stderr, "  KROW = %d\n", krow );
            fprintf ( stderr, "  KCOL = %d\n", kcol );
            fprintf ( stderr, "  NODE[KROW] = %d\n", node[krow] );
            fprintf ( stderr, "  NODE[KCOL] = %d\n", node[kcol] );
            exit ( 1 );
          }
          a[node[krow]+node[kcol]*n] = a[node[krow]+node[kcol]*n]
            + 20.0 * rho * em[krow+kcol*8] / 9.0;
        }
      }
    }
  }
  return a;
}
/******************************************************************************/

int wathen_order ( int nx, int ny )

/******************************************************************************/
/*
  Purpose:

    WATHEN_ORDER returns the order of the WATHEN matrix.

  Discussion:

    N = 3*NX*NY + 2*NX + 2*NY + 1,

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    10 January 2013

  Author:

    John Burkardt

  Reference:

    Nicholas Higham,
    Algorithm 694: A Collection of Test Matrices in MATLAB,
    ACM Transactions on Mathematical Software,
    Volume 17, Number 3, September 1991, pages 289-305.

    Andrew Wathen,
    Realistic eigenvalue bounds for the Galerkin mass matrix,
    IMA Journal of Numerical Analysis,
    Volume 7, 1987, pages 449-457.

  Parameters:

    Input, int NX, NY, values which determine the size of A.

    Output, int WATHEN_ORDER, the order of the matrix.
*/
{
  int n;

  n = 3 * nx * ny + 2 * nx + 2 * ny + 1;

  return n;
}