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

#*****************************************************************************80
#
## R8SD_CG uses the conjugate gradient method on an R8SD linear system.
#
#  Discussion:
#
#    The R8SD storage format is for symmetric matrices whose only nonzero
#    entries occur along a few diagonals, but for which these diagonals are 
#    not all close enough to the main diagonal for band storage to be efficient.
#
#    In that case, we assign the main diagonal the offset value 0, and 
#    each successive superdiagonal gets an offset value 1 higher, until
#    the highest superdiagonal (the A(1,N) entry) is assigned the offset N-1.
#
#    Assuming there are NDIAG nonzero diagonals (ignoring subdiagonals#),
#    we then create an array B that has N rows and NDIAG columns, and simply
#    "collapse" the matrix A to the left:
#
#    For the conjugate gradient method to be applicable, the matrix A must 
#    be a positive definite symmetric matrix.
#
#    The method is designed to reach the solution to the linear system
#      A * x = b
#    after N computational steps.  However, roundoff may introduce
#    unacceptably large errors for some problems.  In such a case,
#    calling the routine a second time, using the current solution estimate
#    as the new starting guess, should result in improved results.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    09 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 NDIAG, the number of diagonals that are stored.
#    NDIAG must be at least 1 and no more than N.
#
#    Input, integer OFFSET(NDIAG), the offsets for the diagonal
#    storage.
#
#    Input, real A(N,NDIAG), the R8SD 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 = r8sd_mv ( n, n, ndiag, offset, 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 = r8sd_mv ( n, n, ndiag, offset, 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 r8sd_cg_test ( ):

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

  n = 10

  ndiag = 2
  offset = np.zeros ( ndiag, dtype = np.int32 )
  offset[0] = 0
  offset[1] = 1
#
#  Set A to the [-1 2 -1] matrix.
#
  a = r8sd_dif2 ( n, n, ndiag, offset )
#
#  Choose a random solution.
#
  seed = 123456789
  x1, seed = r8vec_uniform_01 ( n, seed )
#
#  Compute the corresponding right hand side.
#
  b = r8sd_mv ( n, n, ndiag, offset, a, x1 )
#
#  Call the CG routine.
#
  x2 = np.ones ( n )
  x2 = r8sd_cg ( n, ndiag, offset, a, b, x2 )
#
#  Compute the residual.
#
  r = r8sd_res ( n, n, ndiag, offset, a, x2, b )
  r_norm = r8vec_norm ( n, r )
#
#  Compute the error.
#
  e_norm = r8vec_norm_affine ( n, x1, x2 )
#
#  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 r8sd_dif2 ( m, n, ndiag, offset ):

#*****************************************************************************80
#
## R8SD_DIF2 returns the DIF2 matrix in R8SD 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:
#
#    09 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 NDIAG, the number of diagonals that are stored.
#    NDIAG must be at least 2.
#
#    Input, integer OFFSET(NDIAG), the offsets for the diagonal
#    storage.  It is simply assumed that OFFSET(1) = 0 and OFFSET(2) = 1.
#
#    Output, real A(N,NDIAG), the matrix.
#
  import numpy as np

  a = np.zeros ( [ n, ndiag ] )

  for i in range ( 0, n ):
    a[i,0] = 2.0

  for i in range ( 0, n - 1 ):
    a[i,1] = -1.0
 
  return a

def r8sd_mv ( m, n, ndiag, offset, a, x ):

#*****************************************************************************80
#
## R8SD_MV multiplies an R8SD matrix by an R8VEC.
#
#  Discussion:
#
#    The R8SD storage format is for symmetric matrices whose only nonzero 
#    entries occur along a few diagonals, but for which these diagonals are not 
#    all close enough to the main diagonal for band storage to be efficient.
#
#    In that case, we assign the main diagonal the offset value 0, and 
#    each successive superdiagonal gets an offset value 1 higher, until
#    the highest superdiagonal (the A(1,N) entry) is assigned the offset N-1.
#
#    Assuming there are NDIAG nonzero diagonals (ignoring subdiagonals#),
#    we then create an array B that has N rows and NDIAG columns, and simply
#    "collapse" the matrix A to the left:
#
#  Example:
#
#    The "offset" value is printed above each column.
#
#    Original matrix               New Matrix
#
#       0   1   2   3   4   5       0   1   3   5
#
#      11  12   0  14   0  16      11  12  14  16
#      21  22  23   0  25   0      22  23  25  --
#       0  32  33  34   0  36      33  34  36  --
#      41   0  43  44  45   0      44  45  --  --
#       0  52   0  54  55  56      55  56  --  --
#      61   0  63   0  65  66      66  --  --  --
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    09 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer M, N, the number of rows and columns.
#
#    Input, integer NDIAG, the number of diagonals that are stored.
#    NDIAG must be at least 1 and no more than N.
#
#    Input, integer OFFSET(NDIAG), the offsets for the diagonal
#    storage.
#
#    Input, real A(N,NDIAG), the R8SD matrix.
#
#    Input, real X(N), the vector to be multiplied by A.
#
#    Output, real B(N), the product A * x.
#
  import numpy as np

  b = np.zeros ( n )

  for i in range ( 0, n ):
    for jdiag in range ( 0, ndiag ):
      if ( 0 <= offset[jdiag] ):
        j = i + offset[jdiag]
        if ( 0 <= j and j < n ):
          b[i] = b[i] + a[i,jdiag] * x[j]
          if ( offset[jdiag] != 0 ):
            b[j] = b[j] + a[i,jdiag] * x[i]

  return b

def r8sd_res ( m, n, ndiag, offset, a, x, b ):

#*****************************************************************************80
#
## R8SD_RES computes the residual R = B-A*X for R8SD matrices.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    09 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 NDIAG, the number of diagonals that are stored.
#    NDIAG must be at least 1 and no more than N.
#
#    Input, integer OFFSET(NDIAG), the offsets for the diagonal
#    storage.
#
#    Input, real A(N,NDIAG), the R8SD 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 = r8sd_mv ( m, n, ndiag, offset, a, x )

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

  return r

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