// ................................................................... diago.cc
#include <curses.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "corso.h"

extern XVideoData xvd;

main(int argc,char *argv[])
{
	int i,j,k,nn,salta;
	double *diag,*prod,*subdiag;
	double **mat,**transf;
	char answ[256],buff[80],MatFile[80];
	FILE *hnd;
	
	if (argc<2 || argc>2)
	{ // ............................................... la chiamata e' corretta?
		printf("\n Uso: diago file\n");
		if (argc<2) printf("\n Hai dimenticato il nome del file con la matrice!!!\n");
		else if (argc>2) printf("\n Hai scritto troppo!!!\n");
		printf("\n");
		exit(0);
	}
	// ................................................... copia il nome del file
	strcpy(MatFile,argv[1]);
	// ............................................................. apri il file
	if ((hnd=fopen(MatFile,"r"))==NULL)
	{
		printf("\n non riesco ad aprire %s\n",MatFile);
		exit(0);
	}
	// ........................................................... conta le righe
	nn=0;
	for (;;)
	{
		fgets(answ,256,hnd);
		if (feof(hnd)) break;
		fieldcpy(buff,answ,0);
		if (atof(buff)==0.0 && buff[0]!='0') continue;
		nn++;
	}
	rewind(hnd);
	// ........................................................ alloca la memoria
	mat=new double *[nn];
	for (i=0;i<nn;i++) mat[i]=new double [nn];
	transf=new double *[nn];
	for (i=0;i<nn;i++) transf[i]=new double [nn];
	diag=new double [nn];
	subdiag=new double [nn];
	prod=new double [nn];
	// ......................................................... copia la matrice
	i=0;
	for (;;)
	{
		j=0;
		fgets(answ,256,hnd);
		if (feof(hnd)) break;
		salta=FALSE;
		for (j=0;j<nn;j++)
		{
			fieldcpy(buff,answ,j);
			if (atof(buff)==0.0 && buff[0]!='0') {salta=TRUE;break;}
			mat[i][j]=transf[i][j]=atof(buff);
		}
		if (salta) continue;
		i++;
	}
	fclose(hnd);
	// .......................................................... tridiagonalizza
	tred2(transf,diag,subdiag,nn);
	// ............................................... apri il file dei risultati
	hnd=fopen("risultati","w");
	fprintf(hnd,"matrice originale:\n\n");
	for (i=0;i<nn;i++)
	{ // ...................................................... scrivi la matrice
		for (j=0;j<nn;j++) fprintf(hnd,"%15.5le",mat[i][j]);
		fprintf(hnd,"\n");
	}
	// ............................................................. diagonalizza
	if (!tqli(diag,subdiag,transf,nn)) printf("Problemi!!!\n");
	// .............................................. scrivi la diagonalizzazione
	fprintf(hnd,"Autovalori:\n\n");
	for (i=0;i<nn;i++) fprintf(hnd,"%4d%15.5le\n",i,diag[i]);
	fprintf(hnd,"\nColonne di autovettori:\n\n");
	for (i=0;i<nn;i++)
	{ // .......................... transf contiene gli autovettori nelle colonne
		for (j=0;j<nn;j++) fprintf(hnd,"%15.5le",transf[i][j]);
		fprintf(hnd,"\n");
	}
	// ................................ moltiplica la matrice per gli autovettori
	for (k=0;k<nn;k++)
	{ // ......................................... matrice l'autovettore numero k
		fprintf(hnd,"\n\nProdotto matrice per l'autovettore n. %d\n\n",k);
		for (i=0;i<nn;i++)
		{
			prod[i]=0.0;
			for (j=0;j<nn;j++) prod[i]+=mat[i][j]*transf[j][k];
			fprintf(hnd,"%15.5le %15.5le\n",prod[i],prod[i]/transf[i][k]);
		}
	}
	// ............................................. chiudi il file dei risultati
	fprintf(hnd,"\n");
	// ........................................................ libera la memoria
	delete [] prod;
	delete [] subdiag;
	delete [] diag;
	for (i=nn-1;i>=0;i--) delete [] transf[i];
	delete [] transf;
	for (i=nn-1;i>=0;i--) delete [] mat[i];
	delete [] mat;
}