#! /usr/bin/env python
#
def r83t_cg ( n, a, b, x_init ):

#*****************************************************************************80
#
## R83T_CG uses the conjugate gradient method on an R83T system.
#
#  Discussion:
#
#    The R83T storage format is used for a tridiagonal matrix.
#    The superdiagonal is stored in entries (1:N-1,3), the diagonal in
#    entries (1:N,2), and the subdiagonal in (2:N,1).  Thus, the
#    original matrix is "collapsed" horizontally into the array.
#
#    The matrix A must be a positive definite symmetric band matrix.
#
#    The method is designed to reach the solution after N computational
#    steps.  However, roundoff may introduce unacceptably large errors for
#    some problems.  In such a case, calling the routine again, using
#    the computed solution as the new starting estimate, should improve
#    the results.
#
#  Example:
#
#    Here is how an R83T matrix of order 5 would be stored:
#
#       *  A11 A12
#      A21 A22 A23
#      A32 A33 A34
#      A43 A44 A45
#      A54 A55  *
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    08 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Frank Beckman,
#    The Solution of Linear Equations by the Conjugate Gradient Method,
#    in Mathematical Methods for Digital Computers,
#    edited by John Ralston, Herbert Wilf,
#    Wiley, 1967,
#    ISBN: 0471706892,
#    LC: QA76.5.R3.
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#    N must be positive.
#
#    Input, real A(N,3), the matrix.
#
#    Input, real B(N), the right hand side vector.
#
#    Input, real X_INIT(N), an estimate for the solution, which may be 0.
#
#    Output, real X(N), the approximate solution vector.
#
  import numpy as np

  x = np.zeros ( n )
  for i in range ( 0, n ):
    x[i] = x_init[i]
#
#  Initialize
#    AP = A * x,
#    R  = b - A * x,
#    P  = b - A * x.
#
  ap = r83t_mv ( n, n, a, x )

  r = np.zeros ( n )
  for i in range ( 0, n ):
    r[i] = b[i] - ap[i]

  p = np.zeros ( n )
  for i in range ( 0, n ):
    p[i] = b[i] - ap[i]
#
#  Do the N steps of the conjugate gradient method.
#
  for it in range ( 0, n ):
#
#  Compute the matrix*vector product AP=A*P.
#
    ap = r83t_mv ( n, n, a, p )
#
#  Compute the dot products
#    PAP = P*AP,
#    PR  = P*R
#  Set
#    ALPHA = PR / PAP.
#
    pap = np.dot ( p, ap )
    pr = np.dot ( p, r )

    if ( pap == 0.0 ):
      return x

    alpha = pr / pap
#
#  Set
#    X = X + ALPHA * P
#    R = R - ALPHA * AP.
#
    for i in range ( 0, n ):
      x[i] = x[i] + alpha * p[i]

    for i in range ( 0, n ):
      r[i] = r[i] - alpha * ap[i]
#
#  Compute the vector dot product
#    RAP = R*AP
#  Set
#    BETA = - RAP / PAP.
#
    rap = np.dot ( r, ap )

    beta = - rap / pap
#
#  Update the perturbation vector
#    P = R + BETA * P.
#
    for i in range ( 0, n ):
      p[i] = r[i] + beta * p[i]

  return x

def r83t_cg_test ( ):

#*****************************************************************************80
#
## R83T_CG_TEST tests R83T_CG.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    18 June 2014
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from r8vec_norm import r8vec_norm
  from r8vec_norm_affine import r8vec_norm_affine
  from r8vec_uniform_01 import r8vec_uniform_01

  print ( '' )
  print ( 'R83T_CG_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R83T_CG applies CG to an R83T matrix.' )

  n = 10
#
#  Let A be the -1 2 -1 matrix.
#
  a = r83t_dif2 ( n, n )
#
#  Choose a random solution.
#
  seed = 123456789
  x1, seed = r8vec_uniform_01 ( n, seed )
#
#  Compute the corresponding right hand side.
#
  b = r83t_mv ( n, n, a, x1 )
#
#  Call the CG routine.
#
  x2 = np.ones ( n )
  x3 = r83t_cg ( n, a, b, x2 )
#
#  Compute the residual.
#
  r = r83t_res ( n, n, a, x3, b )
  r_norm = r8vec_norm ( n, r )
#
#  Compute the error.
#
  e_norm = r8vec_norm_affine ( n, x1, x3 )
#
#  Report.
#
  print ( '' )
  print ( '  Number of variables N = %d' % ( n ) )
  print ( '  Norm of residual ||Ax-b|| = %g' % ( r_norm ) )
  print ( '  Norm of error ||x1-x2|| = %g' % ( e_norm ) )

  return

def r83t_dif2 ( m, n ):

#*****************************************************************************80
#
## R83T_DIF2 returns the DIF2 matrix in R83T format.
#
#  Example:
#
#    N = 5
#
#    2 -1  .  .  .
#   -1  2 -1  .  .
#    . -1  2 -1  .
#    .  . -1  2 -1
#    .  .  . -1  2
#
#  Properties:
#
#    A is banded, with bandwidth 3.
#
#    A is tridiagonal.
#
#    Because A is tridiagonal, it has property A (bipartite).
#
#    A is a special case of the TRIS or tridiagonal scalar matrix.
#
#    A is integral, therefore det ( A ) is integral, and 
#    det ( A ) * inverse ( A ) is integral.
#
#    A is Toeplitz: constant along diagonals.
#
#    A is symmetric: A' = A.
#
#    Because A is symmetric, it is normal.
#
#    Because A is normal, it is diagonalizable.
#
#    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
#
#    A is positive definite.
#
#    A is an M matrix.
#
#    A is weakly diagonally dominant, but not strictly diagonally dominant.
#
#    A has an LU factorization A = L * U, without pivoting.
#
#      The matrix L is lower bidiagonal with subdiagonal elements:
#
#        L(I+1,I) = -I/(I+1)
#
#      The matrix U is upper bidiagonal, with diagonal elements
#
#        U(I,I) = (I+1)/I
#
#      and superdiagonal elements which are all -1.
#
#    A has a Cholesky factorization A = L * L', with L lower bidiagonal.
#
#      L(I,I) =    sqrt ( (I+1) / I )
#      L(I,I-1) = -sqrt ( (I-1) / I )
#
#    The eigenvalues are
#
#      LAMBDA(I) = 2 + 2 * COS(I*PI/(N+1))
#                = 4 SIN^2(I*PI/(2*N+2))
#
#    The corresponding eigenvector X(I) has entries
#
#       X(I)(J) = sqrt(2/(N+1)) * sin ( I*J*PI/(N+1) ).
#
#    Simple linear systems:
#
#      x = (1,1,1,...,1,1),   A*x=(1,0,0,...,0,1)
#
#      x = (1,2,3,...,n-1,n), A*x=(0,0,0,...,0,n+1)
#
#    det ( A ) = N + 1.
#
#    The value of the determinant can be seen by induction,
#    and expanding the determinant across the first row:
#
#      det ( A(N) ) = 2 * det ( A(N-1) ) - (-1) * (-1) * det ( A(N-2) )
#                = 2 * N - (N-1)
#                = N + 1
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Robert Gregory, David Karney,
#    A Collection of Matrices for Testing Computational Algorithms,
#    Wiley, 1969,
#    ISBN: 0882756494,
#    LC: QA263.68
#
#    Morris Newman, John Todd,
#    Example A8,
#    The evaluation of matrix inversion programs,
#    Journal of the Society for Industrial and Applied Mathematics,
#    Volume 6, Number 4, pages 466-476, 1958.
#
#    John Todd,
#    Basic Numerical Mathematics,
#    Volume 2: Numerical Algebra,
#    Birkhauser, 1980,
#    ISBN: 0817608117,
#    LC: QA297.T58.
#
#    Joan Westlake,
#    A Handbook of Numerical Matrix Inversion and Solution of 
#    Linear Equations,
#    John Wiley, 1968,
#    ISBN13: 978-0471936756,
#    LC: QA263.W47.
#
#  Parameters:
#
#    Input, integer M, N, the order of the matrix.
#
#    Output, real A(M,3), the matrix.
#
  import numpy as np

  a = np.zeros ( [ m, 3 ] )

  mn = min ( m, n )

  for i in range ( 1, mn ):
    a[i,0] = -1.0
  for i in range ( 0, mn ):
    a[i,1] = 2.0
  for i in range ( 0, mn - 1 ):
    a[i,2] = -1.0

  if ( m < n ):
    a[mn-1,2] = -1.0
  elif ( n < m ):
    a[mn,0] = -1.0
  
  return a

def r83t_mv ( m, n, a, x ):

#*****************************************************************************80
#
## R83T_MV multiplies an R83T matrix times an R8VEC.
#
#  Discussion:
#
#    The R83T storage format is used for a tridiagonal matrix.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    08 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer M, N, the number of rows and columns.
#
#    Input, real A(M,3), the matrix.
#
#    Input, real X(N), the vector to be multiplied by A.
#
#    Output, real B(M), the product A * x.
#
  import numpy as np

  b = np.zeros ( m )

  mn = min ( m, n )

  if ( n == 1 ):
    b[0] = a[0,1] * x[0]
    if ( 1 < m ):
      b[1] = a[1,0] * x[0]
    return b

  b[0]      = a[0,1]       * x[0] \
            + a[0,2]       * x[1]

  for i in range ( 1, mn - 1 ):
    b[i] = a[i,0] * x[i-1] \
         + a[i,1] * x[i] \
         + a[i,2] * x[i+1]

  b[mn-1]   = a[mn-1,0]    * x[mn-2] \
            + a[mn-1,1]    * x[mn-1]

  if ( n < m ):
    b[mn] = b[mn] + a[mn,0] * x[mn-1]
  elif ( m < n ):
    b[mn-1] = b[mn-1] + a[mn-1,2] * x[mn]

  return b

def r83t_res ( m, n, a, x, b ):

#*****************************************************************************80
#
## R83T_RES computes the residual R = B-A*X for R83T matrices.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    08 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer M, the number of rows of the matrix.
#    M must be positive.
#
#    Input, integer N, the number of columns of the matrix.
#    N must be positive.
#
#    Input, real A(M,3), the matrix.
#
#    Input, real X(N), the vector to be multiplied by A.
#
#    Input, real B(M), the desired result A * x.
#
#    Output, real R(M), the residual R = B - A * X.
#
  import numpy as np

  r = r83t_mv ( m, n, a, x )

  for i in range ( 0, m ):
    r[i] = b[i] - r[i]

  return r

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  r83t_cg_test ( )
  timestamp ( )