/* Test program to read in matrix and pass to subroutine HICHOL
   for Cholesky decomposition */

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

#define FITS 20         /* Maximum number of nights */
#define STARS 60        /* Maximum number of stars */
#define NSTARS 40       /* Max. number of stars per night */
#define FITPARMS 8      /* Max. number of arc parameters per night */ 
#define STARPARMS 2     /* Max. number of parameters per star */

#define NV ((long)FITS*(FITPARMS*(FITPARMS+1)/2+FITPARMS*NSTARS*STARPARMS))
#define NT ((long)STARS*STARPARMS*(STARS*STARPARMS+1)/2)
#define NP ((long)FITS*FITPARMS+STARS*STARPARMS)

#define ST1 2000000.0
#define ST2 2000000.0

/*
#define NVD FITS*(FITPARMS*(FITPARMS+1)/2+FITPARMS*NSTARS*STARPARMS) \
          + STARS*STARPARMS*(STARS*STARPARMS+1)/2
*/
#define NVD NV+NT

int hichol(float v[],unsigned long iv[],unsigned short jv[],
           unsigned short ni[],unsigned short ns,unsigned short nm, 
           float t[]);
int hi_prep(double vd[],unsigned long iv[],unsigned short jv[],double rd[],
            unsigned short ni[],unsigned short ns,unsigned short nm, 
            float v[],float t[],float r[],float s[]);
int hi_chol(float v[],unsigned long iv[],unsigned short jv[],
            unsigned short ir,unsigned short nr,unsigned short mb,
            float t[],unsigned short nm,unsigned short nb);
int cholesky(float t[],unsigned short n);
int hi_solve(float v[],unsigned long iv[],unsigned short jv[],
             unsigned short ni[],unsigned short ns,unsigned short nm,
	     float t[],float r[],float s[]);
int hi_res(double vd[],unsigned long iv[],unsigned short jv[],double rd[],
	   double rn[],float s[],float r[],unsigned short nm);
int hi_pow(double vd[],unsigned long iv[],unsigned short jv[],double u[],
	   double w[],double rn[],unsigned short nm);

main()
        {
        char filename[80];
        FILE *fp;

        unsigned long i,j,ii,jj;
        unsigned short nm,ns;

        /* Input matrix in row-wise format */
        double vd[NVD+1];          /* data */
        unsigned short jv[NVD+1];  /* column numbers */
        unsigned long iv[NP+2];    /* row positions */

 	/* Further input information */
        unsigned short ni[FITS+1]; /* dimensions of sub-matrices */
        double rd[NP+1];           /* right hand side */

	/* Work space */
        double rn[NP+1],u[NP+1],w[NP+1];
        float v[NV+1],t[NT+1],r[NP+1],s[NP+1];

	/* Miscellaneous */
	float sum,ev,tiny;

printf("%f %f %lf\n",ST1,ST2,ST1*ST2);
        puts("\nThis program reads a matrix and ");
        puts("calculates the Cholesky decomposition.");

/*--------------------------------------------------------------------*/
        printf("%s","\nPlease enter filename of sparse matrix: ");
        scanf("%s",filename);
        if((fp=fopen(filename,"r"))==NULL){
                fprintf(stderr,"Error opening file.");
                exit(1);       
                }
        i=1;
        while(i<=NVD){
                if(fscanf(fp,"%lf %hu",&vd[i],&jv[i])==EOF)break;
                i+=1;
                }
        fclose(fp);

        printf("%s","Please enter filename for array iv: ");
        scanf("%s",filename);
        if((fp=fopen(filename,"r"))==NULL){
                fprintf(stderr, "Error opening file.");
                exit(1);       
                }
        i=1;
        while(i<=NP+2){
                if(fscanf(fp,"%lu",&iv[i])==EOF)break;
                i+=1;
                }
        nm=i-2;
        fclose(fp);

        printf("%s","Please enter filename for right hand side: ");
        scanf("%s",filename);
        if((fp=fopen(filename,"r"))==NULL){
                fprintf(stderr, "Error opening file.");
                exit(1);       
                }
        i=1;
        while(i<=NP+1){
                if(fscanf(fp,"%lf",&rd[i])==EOF)break;
                i+=1;
                }
        fclose(fp);

        printf("Dimension of matrix is: %d",nm);
        printf("%s","\nEnter number of fits: ");
        scanf("%hu",&ns);
/*--------------------------------------------------------------------*/

        /* Set array of sub-matrix dimensions */
        for(i=1;i<=ns;i++){ni[i]=4;}

	/* Prepare for hichol */
	puts("Setup...");
        hi_prep(vd,iv,jv,rd,ni,ns,nm,v,t,r,s);

	/* Calculate Cholesky decomposition */
	puts("Call hichol now...");
        hichol(v,iv,jv,ni,ns,nm,t);

	/* Solve for right hand side */
	puts("Now solve for right hand side...");
        hi_solve(v,iv,jv,ni,ns,nm,t,r,s);

	/* Improve solution by repeatedly calling these functions */
	puts("Improving solution...");
	for(i=1;i<=2;i++){
		hi_res(vd,iv,jv,rd,rn,s,r,nm);
        	hi_solve(v,iv,jv,ni,ns,nm,t,r,s); 
		}

	puts("Write output...");
        if((fp=fopen("s.out","w"))==NULL){
                fprintf(stderr, "Error opening file.");
                exit(1);
                }
        for(i=1;i<=nm;i++){
                fprintf(fp,"%4.13lf \n ",s[i]);
                }
        fclose(fp);

	/* Calculate smallest eigenvalue */
  	for(i=1;i<=nm;i++){r[i]=1.0;}
	sum=1.0;
	tiny=1.0e-6;
	do{
		ev=sum;
		for(i=1;i<=nm;i++){rd[i]=r[i];}
		hi_solve(v,iv,jv,ni,ns,nm,t,r,s);
	/*	hi_res(vd,iv,jv,rd,rn,s,r,nm);
        	hi_solve(v,iv,jv,ni,ns,nm,t,r,s); */
		hi_res(vd,iv,jv,rd,rn,s,r,nm);
        	hi_solve(v,iv,jv,ni,ns,nm,t,r,s); 
		sum=0.0;
		for(j=1;j<=nm;j++){sum+=s[j]*s[j];}
		sum=sqrt(sum);
		for(j=1;j<=nm;j++){r[j]=s[j]/sum;}
		} while(fabs(ev-sum)>tiny);
	sum=1.0/sum;
	printf("Smallest eigen value = %f",sum);   

/*      puts("\nNow start calculating inverse...");
        for(i=1;i<=nm;i++){
                for(j=1;j<=nm;j++){r[j]=0.0;}
                r[i]=1.0;
                hi_solve(v,iv,jv,ni,ns,nm,t,r,s);
                }
        puts("Finished calculating inverse");
*/
    	/* Calculate largest eigenvalue */
	for(i=1;i<=nm;i++){u[i]=1.0;}
	for(i=1;i<=1;i++){hi_pow(vd,iv,jv,u,w,rn,nm);}
	puts("\nHITEST finished.");
  
        return 0;
        }