/*
  CLAPACK demo to illustrate spectral analysis.
  Here are solve for the eigenvectors and eigenvalues of a symmetric matrix
  in CLAPACK's single-precision real version (ssyev).
  You may explore many other variants at www.netlib.org/lapack
  J.K. Johnstone, 2012
*/

#include <iostream>
using namespace std;

// to avoid C++ incompatibilities (after all, it is expecting C)
extern "C" {               
  // if you put f2c.h and clapack.h in /usr/include (see 'installing CLAPACK')
#include "f2c.h"           
#include "clapack.h"       
  // if you didn't put them in /usr/include, use the full pathname in the 
  // above two lines, such as:
  /*
#include "/Users/jj/Downloads/CLAPACK-3.2.1/INCLUDE/f2c.h"
#include "/Users/jj/Downloads/CLAPACK-3.2.1/INCLUDE/clapack.h"
  */
}

int main (int argc, char **argv)
{
  int i,j;
  // example from p. 411 of Golub and van Loan (2nd ed)
  float *A;  // matrix as a 1d array in column-major order (a la Fortran)
  A = new float[16];
  A[0] = 1; A[1] = 1; A[2] = 1;  A[3] = 1;
  A[4] = 1; A[5] = 2; A[6] = 3;  A[7] = 4;
  A[8] = 1; A[9] = 3; A[10]= 6;  A[11]= 10;
  A[12]= 1; A[13]= 4; A[14]= 10; A[15]= 20;
  cout << "A = (";
  for (i=0; i<15; i++) cout << A[i] << ",";
  cout << A[15] << ")" << endl;

  // use call-by-reference parameters, in another adaptation to Fortran
  int n=4;
  char jobz = 'V'; // compute eigenvectors too
  char uplo = 'U'; // upper triangle of matrix is stored
  integer N=n;     // dimension of matrix
  integer LDA=n;
  float *eigvalue; 
  eigvalue=new float[n]; 
  int worksize;    // must have lwork>=3n-1, so n^2 is not big enough for n=2!
  if (n==2) worksize = 5; else worksize = n*n;
  float *work; 
  work=new float[worksize];
  integer lwork = worksize;
  integer info;    // exit status

  ssyev_ (&jobz, &uplo, &N, A, &LDA, eigvalue, work, &lwork, &info);

  cout << "info = " << info << endl;
  if (info>=0)
   {
    cout << "Eigenvalues = (" 
	 << eigvalue[0] << "," << eigvalue[1] << "," 
	 << eigvalue[2] << "," << eigvalue[3] << ")" << endl;
    cout << "Correct answer in GvL: (.0380, .4538, 2.2034 (typo?), 26.3047)" 
	 << endl;
    // eigenvectors are stored in A (overwriting the original matrix):
    // ith eigenvector = ith column of A
    for (i=0; i<4; i++)
      {
	cout << "Eigenvector " << i << " = (";
	for (j=0; j<3; j++)
	  cout << A[j+i*4] << ","; // jth elt of ith col=jth elt of ith evector
	cout << A[3+i*4] << ")" << endl;
      }
    }
  delete [] A;
  delete [] eigvalue;  // be a good boy scout
  delete [] work;
  return 0;
}