# include <iostream>
# include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main ( );
void jacobiSvd_test ( );
void colPivHouseholderQr_test ( );
void normal_equations_test ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for EIGEN_TEST01.
//
//  Modified:
//
//    07 March 2017
//
{
  cout << "\n";
  cout <<"EIGEN_TEST03:\n";
  cout << "  Least squares solutions of A*x=b\n";

  jacobiSvd_test ( );
  colPivHouseholderQr_test ( );
  normal_equations_test ( );
//
//  Terminate.
//
  cout << "\n";
  cout <<"EIGEN_TEST03:\n";
  cout << "  Normal end of execution.\n";

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

void jacobiSvd_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    JACOBISVD_TEST uses jacobiSvd.
//
//  Modified:
//
//    07 March 2017
//
{
  MatrixXf A;
  VectorXf b;
  VectorXf r;
  VectorXf x;

  cout << "\n";
  cout <<"JACOBISVD_TEST:\n";
  cout << "  Least squares solution of A*x=b using jacobiSvd.\n";

  srand ( 123456789 );

  A = MatrixXf::Random ( 3, 2 );

  cout << "\n";
  cout << "  A:\n";
  cout << "\n";
  cout << A << "\n";;

  b = VectorXf::Random ( 3 );

  cout << "\n";
  cout << "  b:\n";
  cout << "\n";
  cout << b << "\n";;

  x = A.jacobiSvd ( ComputeThinU | ComputeThinV ).solve ( b );

  cout << "\n";
  cout << "  x:\n";
  cout << "\n";
  cout << x << "\n";

  r = A * x - b;

  cout << "\n";
  cout << "  r = A * x - b:\n";
  cout << "\n";
  cout << r << "\n";;

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

void colPivHouseholderQr_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    COLPIVHOUSEHOLDERQR_TEST uses colPivHouseholderQr.
//
//  Modified:
//
//    07 March 2017
//
{
  MatrixXf A;
  VectorXf b;
  VectorXf r;
  VectorXf x;

  cout << "\n";
  cout <<"COLPIVHOUSEHOLDERQR_TEST:\n";
  cout << "  Least squares solution of A*x=b using colPivHouseholderQr.\n";

  srand ( 123456789 );

  A = MatrixXf::Random ( 3, 2 );

  cout << "\n";
  cout << "  A:\n";
  cout << "\n";
  cout << A << "\n";;

  b = VectorXf::Random ( 3 );

  cout << "\n";
  cout << "  b:\n";
  cout << "\n";
  cout << b << "\n";;

  x = A.colPivHouseholderQr ( ).solve ( b );

  cout << "\n";
  cout << "  x:\n";
  cout << "\n";
  cout << x << "\n";

  r = A * x - b;

  cout << "\n";
  cout << "  r = A * x - b:\n";
  cout << "\n";
  cout << r << "\n";;

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

void normal_equations_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    NORMAL_EQUATIONS_TEST uses the normal equations.
//
//  Modified:
//
//    07 March 2017
//
{
  MatrixXf A;
  VectorXf b;
  VectorXf r;
  VectorXf x;

  cout << "\n";
  cout <<"NORMAL_EQUATIONS_TEST:\n";
  cout << "  Least squares solution of A*x=b using normal equations.\n";

  srand ( 123456789 );

  A = MatrixXf::Random ( 3, 2 );

  cout << "\n";
  cout << "  A:\n";
  cout << "\n";
  cout << A << "\n";;

  b = VectorXf::Random ( 3 );

  cout << "\n";
  cout << "  b:\n";
  cout << "\n";
  cout << b << "\n";;

  x = ( A.transpose ( ) * A ).ldlt().solve ( A.transpose ( ) * b );

  cout << "\n";
  cout << "  x:\n";
  cout << "\n";
  cout << x << "\n";

  r = A * x - b;

  cout << "\n";
  cout << "  r = A * x - b:\n";
  cout << "\n";
  cout << r << "\n";;

  return;
}