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

#define pi     3.14159265359
#define m      50

/*  The following commands compile and link the code.        */

/*  Sun's C compiler:                                        */
/*  cc pack.c -xlic_lib=sunperf              (creates a.out)  */
/*               or                                          */
/*  cc pack.c -xlic_lib=sunperf -o dot.exe   (creates dot.exe)*/

/*  Gnu's C compiler: (Currently does NOT work for sgesv!!!) */
/*  gcc pack.c -lsunperf -lF77 -lsunmath -lM77                */ 
/*      -L/nau/share/lang/SUNWspro/SC4.2/lib -lSUNWPro_lic   */

/* prototypes of LAPACK and BLAS function calls */
/* see httP://www.netlib.org */
float sdot(long int n, float *sx, long int incx, float *sy, long int incy);
void  sgesv(long int n1, long int n2, float a[3][3], long int n3, int *piv, 
            float *b, long int n4, int *info);
void ssyev(char jobz, char uplo, long int n5, float a[3][3], long int n6,          
           float w[3], float work[10],
			  int lwork, int info);

/* prototypes of user define functions */
void		tofile(float u[m+1], FILE *fp);
float	   f(float x);

float aa = 0, bb=pi;     
float dx;

const  char fname[10]="data.txt"; /* store vector in this file */
float  u[m+1];
FILE   *fp;

/***************************************************************/
main ()
{
   /* variable declarations, array dimensioning */
   float sx[3], sy[3];
   float result;
   float a[3][3], b[3], w[3];
   float work[10];
   float x, sum;
   int   piv[3];
   int   info;
   int   i,j;

   /* initialization */
   
   a[0][0]=1; a[0][1]=2; a[0][2]=0;
   a[1][0]=2; a[1][1]=2; a[1][2]=4;
   a[2][0]=0; a[2][1]=4; a[2][2]=3;

   b[0]=3; b[1]=8; b[2]=7;

   /* compute eigenvalue/vector and print */
	ssyev('V', 'L', 3, a, 3, w, work, 10, info);
	printf ("info from ssyev = %d\n",info);
	printf("Eigenvectors (now stored in a):\n\n");
	for(i=0;i<=2;i++) 
		{
		for(j=0;j<=2;j++)
			printf("%6.3f ",a[i][j]);
		printf("\n");
		}
   printf("\n");
	
	printf("w variable, the eigenvalues in ascending order: \n\n");
	for(i=0;i<=2;i++)
		printf("w[%d]=%6.3f \n", i, w[i]);

	printf("\nDone with eigenvalues and eigen vector.\n\n");

   /* reinitialize */
   a[0][0]=1; a[0][1]=2; a[0][2]=0;
   a[1][0]=2; a[1][1]=2; a[1][2]=4;
   a[2][0]=0; a[2][1]=4; a[2][2]=3;

   printf("after reinitializing, a is as it should be\n");
   for(i=0;i<=2;i++) 
   {
      for(j=0;j<=2;j++)
         printf("%6.3f ",a[i][j]);
      printf("\n");
   }
   printf("\n");

   printf("Here is right hand side of system to solve\n");
   for(j=0;j<=2;j++)
      printf("%6.3f ",b[j]);
   printf("\n"); 
   printf("\n");

   /* solve linear system ax=b, answer->b 
   /* Do this before ssyev changes a or reinitialize */
   sgesv(3,1,a,3,piv,b,3,&info);
   printf ("info from sgesv = %d\n",info);
   printf("solution of system overwrites right hand side\n");
   for(j=0;j<=2;j++)
      printf("%6.3f ",b[j]);
   printf("\n\n");

   /* now do dot product */
   sx[0] = 1.0;  sx[1] = 2.0; sx[2] = 3.0;
   sy[0] = 3.0;  sy[1] = 4.0; sy[2] = -1.0;

   /* compute dot product of sx and sy */
   result = sdot (3, sx, 1, sy, 1);
   printf ("Sdot result is %f\n", result);

   /* open file, populate array with values from f, write to file */

   /* initialize u with f */
   dx=(bb-aa)/m;
   for (i=0; i<=m; i++)
	{
		x = aa + dx * (float) i;
		u[i] = f( x );
	}

	/* write u to file */
   if ((fp= fopen(&fname[0],"w")) == NULL) 
      printf("Error open write file\n");
	tofile(u, fp);
	fclose(fp);

	/* Trapezoid rule */
	sum = 0;
   for (i=1; i<=m; i++)
		sum += (u[i-1] + u[i])/2.0;
	sum *= dx;
	printf("the TRAP rule gives: %g\n",sum);

} /* end of main program */
/**********************************************************************/
float f(float x)
{
	float ret;

	ret = sin(x);
	
	return(ret);
}
/***************************************/
/* send array to FILE                  */
/***************************************/
void tofile(float v[m+1], FILE *fp)
{
   int i;

   for(i=0;i<=m;i++)
	{
	   fprintf(fp, "%9.5f \n", v[i]);
	}
}