# include <cmath>
# include <cstdlib>
# include <ctime>
# include <iostream>

extern "C" void dgeev_ ( char *jobvl, char *jobvr, int *n, double *A, int *lda, 
  double *wr, double *wi, double *vl, int *ldvl, double *vr,  int *ldvr, 
  double *work, int *lwork, int *info );

extern "C" void dgesv_ ( int *n, int *nrhs, double  *a, int *lda, int *ipiv, 
  double *b, int *lbd, int *info  );

using namespace std;

int main ( );
void dgeev_test ( );
void dgesv_test ( );
double *clement1_matrix ( int n );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    lapack_test() tests lapack().
//
//  Discussion:
//
//    This example shows how a C++ calling program can interact with the
//    lapack() library, which is written in Fortran.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    06 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );

  cout << "\n";
  cout << "lapack_test():\n";
  cout << "  C++ version\n";
  cout << "  lapack() is a Fortran library of linear algebra utilities.\n";
  cout << "  A C++ program can call it, if careful.\n";

  dgeev_test ( );
  dgesv_test ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "lapack_test():\n";
  cout << "  Normal end of execution,\n";
  timestamp ( );

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

void dgeev_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    dgeev_test() tests dgeev().
//
//  Discussion:
//
//    dgeev() computes the eigenvalues and eigenvectors of a matrix A,
//    where the matrix A is stored in general format.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    06 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  int n = 5;

  double *A;
  int i;
  int info;
  char jobvl;
  char jobvr;
  int lda = n;
  int ldvl = 1;
  int ldvr = 1;
  int lwork;
  double *vl;
  double *vr;
  double *wi;
  double *work;
  double *wr;

  cout << "\n";
  cout << "dgeev_test():\n";
  cout << "  dgeev() computes eigenvalues and eigenvectors of a matrix A,\n";
  cout << "  with matrix A in general storage.\n";

  A = clement1_matrix ( n );

  jobvl = 'N';
  jobvr = 'N';
  lda = n;
  wr = new double[n];
  wi = new double[n];
  vl = NULL;
  ldvl = 1;
  vr = NULL;
  ldvr = 1;
  work = new double[3*n];
  lwork = 3 * n;

  dgeev_ ( &jobvl, &jobvr, &n, A, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
    work, &lwork, &info );

  cout << "\n";
  cout << "  Computed eigenvalues lambda:\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << wr[i] << " + " << wi[i] << " * i\n";
  }

  delete [] A;
  delete [] wi;
  delete [] work;
  delete [] wr;

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

void dgesv_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    dgesv_test() tests dgesv().
//
//  Discussion:
//
//    dgesv() computes the solution x of a linear system A*x=b,
//    where the matrix A is stored in general format.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    06 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  int n = 3;

  double *A;
  double *b;
  int i;
  int info;
  int *ipvt;
  int lda = n;
  int ldb = n;
  int nrhs = 1;
//
//  Matrix is given by columns.
//
  A = new double[n*n];
  A[0+0*3] = 1.0;
  A[1+0*3] = 4.0;
  A[2+0*3] = 7.0;
  A[0+1*3] = 2.0;
  A[1+1*3] = 5.0;
  A[2+1*3] = 8.0;
  A[0+2*3] = 3.0;
  A[1+2*3] = 6.0;
  A[2+2*3] = 0.0;

  b = new double[n];
  b[0] = 14.0;
  b[1] = 32.0;
  b[2] = 23.0;

  ipvt = new int[n];

  cout << "\n";
  cout << "dgesv_test():\n";
  cout << "  dgesv() solves a linear system A*x=b,\n";
  cout << "  with matrix A in general storage.\n";

  dgesv_ ( &n, &nrhs, A, &lda, ipvt, b, &ldb, &info );

  cout << "\n";
  cout << "  Computed solution x:\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << b[i];
  }
  cout << "\n";

  delete [] A;
  delete [] b;
  delete [] ipvt;

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

double *clement1_matrix ( int n )

//****************************************************************************80
//
//  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 = 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

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
}