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

# include "eigs.h"

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

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

int main ( )

/******************************************************************************/
/*
  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 ( );

  printf ( "\n" );
  printf ( "eigs_test():\n" );
  printf ( "  C version\n" );
  printf ( "  Test eigs():\n" );

  printf ( "\n" );
  printf ( "  The Clement1 matrix is symmetric, and has\n" );
  printf ( "  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" );
  free ( A );
  free ( lambda );

  printf ( "\n" );
  printf ( "  The Helmert matrix is nonsymmetric, orthogonal, and has\n" );
  printf ( "  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" );
  free ( A );
  free ( lambda );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "eigs_test():\n" );
  printf ( "  Normal end of execution.\n" );

  timestamp ( );

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

double *clement1_matrix ( int n )

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

    clement1_matrix() 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 = ( double * ) malloc ( n * n * sizeof ( double ) );

  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;
}
/******************************************************************************/

double *helmert_matrix ( int n )

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

    helmert_matrix() 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 = ( double * ) malloc ( n * n * sizeof ( double ) );
/*
  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;
}
/******************************************************************************/

void c8vec_print ( int n, double complex a[], char *title )

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

    c8vec_print() prints a C8VEC.

  Discussion:

    A C8VEC is a vector of double complex values.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 January 2013

  Author:

    John Burkardt

  Input:

    int N, the number of components of the vector.

    double complex A[N], the vector to be printed.

    char *TITLE, a title.
*/
{
  int i;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );
  fprintf ( stdout, "\n" );

  for ( i = 0; i < n; i++ )
  {
    fprintf ( stdout, "  %8d: %14f  %14f\n", i, creal ( a[i] ), cimag ( a[i] ) );
  }

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

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

/******************************************************************************/
/*
  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:

    28 May 2008

  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.

    char *TITLE, a title.
*/
{
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

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

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

/******************************************************************************/
/*
  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.

    char *TITLE, a title.
*/
{
# define INCX 5

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

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );

  if ( m <= 0 || n <= 0 )
  {
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  (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;
    }

    fprintf ( stdout, "\n" );
/*
  For each column J in the current range...

  Write the header.
*/
    fprintf ( stdout, "  Col:  ");
    for ( j = j2lo; j <= j2hi; j++ )
    {
      fprintf ( stdout, "  %7d     ", j - 1 );
    }
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  Row\n" );
    fprintf ( stdout, "\n" );
/*
  Determine the range of the rows in this strip.
*/
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( m < ihi )
    {
      i2hi = m;
    }
    else
    {
      i2hi = ihi;
    }

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

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

void timestamp ( )

/******************************************************************************/
/*
  Purpose:
 
    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# 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
}