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

#*****************************************************************************80
#
## R8PBU_CG uses the conjugate gradient method on an R8PBU system.
#
#  Discussion:
#
#    The R8PBU storage format is for a symmetric positive definite band matrix.
#
#    To save storage, only the diagonal and upper triangle of A is stored,
#    in a compact diagonal format that preserves columns.
#
#    The diagonal is stored in row MU+1 of the array.
#    The first superdiagonal in row MU, columns 2 through N.
#    The second superdiagonal in row MU-1, columns 3 through N.
#    The MU-th superdiagonal in row 1, columns MU+1 through N.
#
#    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.
#
#  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, integer MU, the number of superdiagonals.
#    MU must be at least 0, and no more than N-1.
#
#    Input, real A(MU+1,N), the R8PBU 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 = r8pbu_mv ( n, n, mu, 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 = r8pbu_mv ( n, n, mu, 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 r8pbu_cg_test ( ):

#*****************************************************************************80
#
## R8PBU_CG_TEST tests R8PBU_CG.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    08 July 2015
#
#  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 ( 'R8PBU_CG_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8PBU_CG applies CG to an R8PBU matrix.' )

  seed = 123456789
  n = 10
  mu = 1
#
#  Let A be the -1 2 -1 matrix.
#
  a = r8pbu_dif2 ( n, n, mu )
#
#  Choose a random solution.
#
  x1, seed = r8vec_uniform_01 ( n, seed )
#
#  Compute the corresponding right hand side.
#
  b = r8pbu_mv ( n, n, mu, a, x1 )
#
#  Call the CG routine.
#
  x2 = np.ones ( n )
  x3 = r8pbu_cg ( n, mu, a, b, x2 )
#
#  Compute the residual.
#
  r = r8pbu_res ( n, n, mu, 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 r8pbu_dif2 ( m, n, mu ):

#*****************************************************************************80
#
## R8PBU_DIF2 returns the DIF2 matrix in R8PBU 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 number of rows and columns.
#
#    Input, integer MU, the number of superdiagonals.
#    MU must be at least 0, and no more than N-1.
#
#    Output, real A(MU+1,N), the matrix.
#
  import numpy as np

  a = np.zeros ( [ mu+1, n ] )

  for j in range ( 1, n ):
    a[mu-1,j] = -1.0

  for j in range ( 0, n ):
    a[mu,j] = +2.0
 
  return a

def r8pbu_mv ( m, n, mu, a, x ):

#*****************************************************************************80
#
## R8PBU_MV multiplies an R8PBU matrix by an R8VEC.
#
#  Discussion:
#
#    The R8PBU storage format is for a symmetric positive definite band matrix.
#
#    To save storage, only the diagonal and upper triangle of A is stored,
#    in a compact diagonal format that preserves columns.
#
#    The diagonal is stored in row MU+1 of the array.
#    The first superdiagonal in row MU, columns 2 through N.
#    The second superdiagonal in row MU-1, columns 3 through N.
#    The MU-th superdiagonal in row 1, columns MU+1 through N.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    08 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the order of the matrix.
#    N must be positive.
#
#    Input, integer MU, the number of superdiagonals in the matrix.
#    MU must be at least 0 and no more than N-1.
#
#    Input, real A(MU+1,N), the R8PBU matrix.
#
#    Input, real X(N), the vector to be multiplied by A.
#
#    Output, real B(N), the result vector A * x.
#
  import numpy as np

  b = np.zeros ( n )
#
#  Multiply X by the diagonal of the matrix.
#
  for i in range ( 0, n ):
    b[i] = a[mu,i] * x[i]
#
#  Multiply X by the superdiagonals of the matrix.
#
  for i in range ( mu - 1, -1, - 1 ):
    for j in range ( mu - i, n ):
      ieqn = i + j - mu
      b[ieqn] = b[ieqn] + a[i,j] * x[j]
      b[j] = b[j] + a[i,j] * x[ieqn]

  return b

def r8pbu_res ( m, n, mu, a, x, b ):

#*****************************************************************************80
#
## R8PBU_RES computes the residual R = B-A*X for R8PBU 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, integer MU, the number of superdiagonals in the matrix.
#    MU must be at least 0 and no more than N-1.
#
#    Input, real A(MU+1,N), 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.
#
  r = r8pbu_mv ( m, n, mu, a, x )
 
  for i in range ( 0, m ):
    r[i] = b[i] - r[i]

  return r

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