# include <cmath>
# include <complex>
# include <cstdlib>
# include <string>
# include <ctime>
# include <iomanip>
# include <iostream>

using namespace std;

# include "eigs.hpp"

int main ( );
double *clement1_matrix ( int n );
double *helmert_matrix ( int n );
void c8vec_print ( int n, complex <double> a[], string title );
void r8mat_print ( int m, int n, double a[], string title );
void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, string title );
void timestamp ( );

//****************************************************************************80

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    eigs_test() tests eigs().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    06 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  double *A;
  complex <double> *lambda;
  int n;

  timestamp ( );

  cout << "\n";
  cout << "eigs_test():\n";
  cout << "  C++ version\n";
  cout << "  Test eigs():\n";

  cout << "\n";
  cout << "  The Clement1 matrix is symmetric, and has\n";
  cout << "  integer eigenvalues.\n";
  n = 5;
  A = clement1_matrix ( n );
  r8mat_print ( n, n, A, "  The matrix A:" );
  lambda = eigs ( n, A );
  c8vec_print ( n, lambda, "  The eigenvalues LAMBDA" );
  delete [] A;
  delete [] lambda;

  cout << "\n";
  cout << "  The Helmert matrix is nonsymmetric, orthogonal, and has\n";
  cout << "  complex eigenvalues of norm 1.\n";
  n = 5;
  A = helmert_matrix ( n );
  r8mat_print ( n, n, A, "  The matrix A:" );
  lambda = eigs ( n, A );
  c8vec_print ( n, lambda, "  The eigenvalues LAMBDA" );
  delete [] A;
  delete [] lambda;
//
//  Terminate.
//
  cout << "\n";
  cout << "eigs_test():\n";
  cout << "  Normal end of execution.\n";

  timestamp ( );

  return 0;
}
//****************************************************************************80

double *clement1_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    clement1() returns the CLEMENT1 matrix.
//
//  Formula:
//
//    if ( J = I + 1 )
//      A(I,J) = sqrt(I*(N-I))
//    else if ( I = J + 1 )
//      A(I,J) = sqrt(J*(N-J))
//    else
//      A(I,J) = 0
//
//  Example:
//
//    N = 5
//
//       .    sqrt(4)    .       .       .
//    sqrt(4)    .    sqrt(6)    .       .
//       .    sqrt(6)    .    sqrt(6)    .
//       .       .    sqrt(6)    .    sqrt(4)
//       .       .       .    sqrt(4)    .
//
//  Properties:
//
//    A is tridiagonal.
//
//    A is banded, with bandwidth 3.
//
//    Because A is tridiagonal, it has property A (bipartite).
//
//    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).
//
//    The diagonal of A is zero.
//
//    A is singular if N is odd.
//
//    About 64 percent of the entries of the inverse of A are zero.
//
//    The eigenvalues are plus and minus the numbers
//
//      N-1, N-3, N-5, ..., (1 or 0).
//
//    If N is even,
//
//      det ( A ) = (-1)^(N/2) * (N-1) * (N+1)^(N/2)
//
//    and if N is odd,
//
//      det ( A ) = 0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    06 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Paul Clement,
//    A class of triple-diagonal matrices for test purposes,
//    SIAM Review,
//    Volume 1, 1959, pages 50-52.
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double CLEMENT1[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  A = new double[n*n];

  for ( i = 0; i < n; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      if ( j == i + 1 )
      {
        A[i+j*n] = sqrt ( ( double ) ( ( i + 1 ) * ( n - i - 1 ) ) );
      }
      else if ( i == j + 1 )
      {
        A[i+j*n] = sqrt ( ( double ) ( ( j + 1 ) * ( n - j - 1 ) ) );
      }
      else
      {
        A[i+j*n] = 0.0;
      }
    }
  }
  return A;
}
//****************************************************************************80

double *helmert_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    helmert() returns the HELMERT matrix.
//
//  Formula:
//
//    If I = 1 then
//      A(I,J) = 1 / sqrt ( N )
//    else if J < I then
//      A(I,J) = 1 / sqrt ( I * ( I - 1 ) )
//    else if J = I then
//      A(I,J) = (1-I) / sqrt ( I * ( I - 1 ) )
//    else
//      A(I,J) = 0
//
//  Discussion:
//
//    The matrix given above by Helmert is the classic example of
//    a family of matrices which are now called Helmertian or
//    Helmert matrices.
//
//    A matrix is a (standard) Helmert matrix if it is orthogonal,
//    and the elements which are above the diagonal and below the
//    first row are zero.
//
//    If the elements of the first row are all strictly positive,
//    then the matrix is a strictly Helmertian matrix.
//
//    It is possible to require in addition that all elements below
//    the diagonal be strictly positive, but in the reference, this
//    condition is discarded as being cumbersome and not useful.
//
//    A Helmert matrix can be regarded as a change of basis matrix
//    between a pair of orthonormal coordinate bases.  The first row
//    gives the coordinates of the first new basis vector in the old
//    basis.  Each later row describes combinations of (an increasingly
//    extensive set of) old basis vectors that constitute the new
//    basis vectors.
//
//    Helmert matrices have important applications in statistics.
//
//  Example:
//
//    N = 5
//
//    0.4472    0.4472    0.4472    0.4472    0.4472
//    0.7071   -0.7071         0         0         0
//    0.4082    0.4082   -0.8165         0         0
//    0.2887    0.2887    0.2887   -0.8660         0
//    0.2236    0.2236    0.2236    0.2236   -0.8944
//
//  Properties:
//
//    A is generally not symmetric: A' /= A.
//
//    A is orthogonal: A' * A = A * A' = I.
//
//    Because A is orthogonal, it is normal: A' * A = A * A'.
//
//    A is not symmetric: A' /= A.
//
//    det ( A ) = (-1)^(N+1)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 June 2008
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    HO Lancaster,
//    The Helmert Matrices,
//    American Mathematical Monthly,
//    Volume 72, 1965, pages 4-12.
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double HELMERT[N*N], the matrix.
//
{
  double *a;
  int i;
  int j;

  a = new double[n*n];
//
//  A begins with the first row, diagonal, and lower triangle set to 1.
//
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i == 0 )
      {
        a[i+j*n] = 1.0 / sqrt ( ( double ) ( n ) );
      }
      else if ( j < i )
      {
        a[i+j*n] = 1.0 / sqrt ( ( double ) ( i * ( i + 1 ) ) );
      }
      else if ( i == j )
      {
        a[i+j*n] = - sqrt ( ( double ) ( i     ) )
                   / sqrt ( ( double ) ( i + 1 ) );
      }
      else
      {
        a[i+j*n] = 0.0;
      }
    }
  }
  return a;
}
//****************************************************************************80

void c8vec_print ( int n, complex <double> a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    c8vec_print() prints a C8VEC.
//
//  Discussion:
//
//    A C8VEC is a vector of complex <double> values.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    10 September 2009
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of components of the vector.
//
//    complex <double> A[N], the vector to be printed.
//
//    string TITLE, a title.
//
{
  int i;

  cout << "\n";
  cout << title << "\n";

  cout << "\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << setw(8) << i
         << ": " << real ( a[i] )
         << "  " << imag ( a[i] ) << "\n";
  }

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

void r8mat_print ( int m, int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_print() prints an R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//    Entry A(I,J) is stored as A[I+J*M]
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 September 2009
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, the number of rows in A.
//
//    int N, the number of columns in A.
//
//    double A[M*N], the M by N matrix.
//
//    string TITLE, a title.
//
{
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

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

void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_print_some() prints some of an R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 June 2013
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, the number of rows of the matrix.
//    M must be positive.
//
//    int N, the number of columns of the matrix.
//    N must be positive.
//
//    double A[M*N], the matrix.
//
//    int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    string TITLE, a title.
//
{
# define INCX 5

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  cout << "\n";
  cout << title << "\n";

  if ( m <= 0 || n <= 0 )
  {
    cout << "\n";
    cout << "  (None)\n";
    return;
  }
//
//  Print the columns of the matrix, in strips of 5.
//
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    if ( n < j2hi )
    {
      j2hi = n;
    }
    if ( jhi < j2hi )
    {
      j2hi = jhi;
    }
    cout << "\n";
//
//  For each column J in the current range...
//
//  Write the header.
//
    cout << "  Col:    ";
    for ( j = j2lo; j <= j2hi; j++ )
    {
      cout << setw(7) << j - 1 << "       ";
    }
    cout << "\n";
    cout << "  Row\n";
    cout << "\n";
//
//  Determine the range of the rows in this strip.
//
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( ihi < m )
    {
      i2hi = ihi;
    }
    else
    {
      i2hi = m;
    }

    for ( i = i2lo; i <= i2hi; i++ )
    {
//
//  Print out (up to) 5 entries in row I, that lie in the current strip.
//
      cout << setw(5) << i - 1 << ": ";
      for ( j = j2lo; j <= j2hi; j++ )
      {
        cout << setw(12) << a[i-1+(j-1)*m] << "  ";
      }
      cout << "\n";
    }
  }

  return;
# undef INCX
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  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:
//
//    19 March 2018
//
//  Author:
//
//    John Burkardt
//
{
# define TIME_SIZE 40

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

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

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

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}