#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <lapacke.h>

/* Auxiliary routine prototypes */
extern void print_matrix( char* desc, int m, int n, double* a, int lda );

#define N 4

/* Main program */
int main() {
        /* Locals */
        double ap[10] = {3.00, 6.00, 2.00, 7.00, 1.00, 5.00, 0.00, 4.00, 8.00, 9.00};
        double bp[10] = {3.00, 6.00, 2.00, 7.00, 1.00, 5.00, 0.00, 4.00, 8.00, 9.00};
	double wp[N], wbp[N];
        int lda = N, ldzp = N, ldzbp = N;
	double zp[N*N], zbp[N*N];
	double workopt[1];
        int n = N;
	int iworkopt[1], lwork, liwork;
        int info;
        int infob;

        /* Executable statements */
        /* LAPACKE call */
	info = LAPACKE_dspevd(LAPACK_ROW_MAJOR,'N','U', n, ap,wp,zp,ldzp);
        printf( " LAPACKE_dspevd Example Program Results\n" );
        print_matrix( "Eigenvalues", 1, N, wp, 1 );
	printf("\n");
        /* LAPACK call */
        lwork = -1;
        liwork = -1;
	dspevd_("N","L",&n,bp,wbp,zbp,&ldzp,workopt,&lwork,iworkopt,&liwork,&infob);
        lwork = (int)workopt[0];
        liwork = iworkopt[0];
        double *work = malloc(lwork*sizeof(double));
        int *iwork = malloc(liwork*sizeof(int));
	dspevd_("N","L",&n,bp,wbp,zbp,&ldzp,workopt,&lwork,iworkopt,&liwork,&infob);
        printf( " dspevd Example Program Results\n" );
        print_matrix( "Eigenvalues", 1, N, wbp, 1 );
	printf("\n");
        exit( 0 );
} /* End of DSPEVD Example */

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, double* a, int lda ) {
        int i, j;
        printf( "\n %s\n", desc );
        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.6f", a[i+j*lda] );
                printf( "\n" );
        }
}