/*
  Sposob podlaczenie LAPACKA w ANSI C - wg Rubin H. Landau, zob:    
  http://physics.oregonstate.edu/~landaur/nacphy/lapack/codes/eigen-c.html 
  ORAZ: 
  http://nicolas.limare.net/pro/notes/2014/10/31_cblas_clapack_lapacke/ 
*/

/* 
   Jesli "liblapackck.a" oraz "libblas.a" sa >>na swoich miejscach<<  
   ==> KOMPILUJEMY KOMENDA: 
   $ gcc -lblas -llapack lapack-example.c
*/


#include <stdio.h>

#define size 3	     /* dimension of matrix in THIS example */
struct complex {double re; double im;};    /* a complex number */



/* Deklaracja f. zewnetrznej - mozna pominac w ANSI C (!!) */
extern void zgeev_(char * jobvl, char * jobvr, int * n,
		   struct complex * A, int * lda,
		   struct complex * w, struct complex * vl, int * ldvl,
		   struct complex * vr, int * ldvr,
		   struct complex * work, int * lwork,
		   double * rwork, int * info);



/* Eigenvalues of a generix complex matrix with ZGEEV */
int main()
{
  /*                        ( 1  -i   0 )   */
  /*   Input matrix:   A =  ( i  -1   0 )   */
  /*                        ( 0   0   7 )   */  
  struct complex A[3][3] = {
    { 1.0, 0.0,   0.0,-1.0,   0.0, 0.0 },
    { 0.0, 1.0,  -1.0, 0.0,   0.0, 0.0 },
    { 0.0, 0.0,   0.0, 0.0,   7.0, 0.0 },
  };

  struct complex b[3], DUMMY[1], WORK[6];  
  struct complex AT[size*size];   /* ==> input for the LAPACK routine */
  int i, j, ok, c1, c2, c3;
  char c4;


 
  for (i=0; i<size; i++)       /*  to call a Fortran routine from C   */
    {			       /*  we have to transform the matrix ...  */
      for(j=0; j<size; j++) {
	AT[j + size*i].re = A[j][i].re;
	AT[j + size*i].im = A[j][i].im;
      }		
    }

  c1 = size;		       /*  ... and put all numbers and characters */ 
  c2 = 2*size;    	       /*  we want to pass  */
  c3 = 1;		       /*  to the routine in variables  */
  c4 = 'N';

  /*  find solution using LAPACK routine ZGEEV, all the arguments have to */
  /*  be pointers and you have to add an underscore to the routine name  */
  zgeev_(&c4, &c4, &c1, AT, &c1, b, DUMMY, &c3, DUMMY, &c3,
	 WORK, &c2, (double *)WORK, &ok);      

  /*
    parameters in the order as they appear in the function call
    no left eigenvectors, no right eigenvectors, order of input matrix A,
    input matrix A, leading dimension of A, array for eigenvalues, 
    array for left eigenvalue, leading dimension of DUMMY, 
    array for right eigenvalues, leading dimension of DUMMY,
    workspace array dim>=2*order of A, dimension of WORK
    workspace array dim=2*order of A, return value 
  */

  if (ok==0) 		      	/* print eigenvalues */
    for (i=0; i<size; i++) {
      printf("%lf\t%lf\n", b[i].re, b[i].im);
    }
  else
    printf("An error occured\n");
  
  /* The correct output is:  
     1.414214      0.000000
     -1.414214	   0.000000
     7.000000	   0.000000
  */
  return 0; 
}