/* Test program to read in matrix and pass to subroutine HICHOL for Cholesky decomposition */ #include #include #include 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; /* data */ unsigned short *jv; /* column numbers */ unsigned long *iv; /* row positions */ /* Further input information */ unsigned short *ni; /* dimensions of sub-matrices */ double *rd; /* right hand side */ /* Work space */ double *rn,*u,*w; float *v,*t,*r,*s; /* Miscellaneous */ float sum,ev,tiny; unsigned long FITS=200; /* Maximum number of nights */ unsigned long STARS=600; /* Maximum number of stars */ unsigned long NSTARS=600; /* Max. number of stars per night */ unsigned long FITPARMS=8; /* Max. number of parameters per night */ unsigned long STARPARMS=2; /* Max. number of parameters per star */ unsigned long NV=FITS*(FITPARMS*(FITPARMS+1)/2+FITPARMS*NSTARS*STARPARMS); unsigned long NT=STARS*STARPARMS*(STARS*STARPARMS+1)/2; unsigned long NP=FITS*FITPARMS+STARS*STARPARMS; unsigned long NVD=NV+NT; vd=(double*)malloc((NVD+1)*sizeof(double)); jv=(unsigned short*)malloc((NVD+1)*sizeof(unsigned short)); iv=(unsigned long*)malloc((NP+2)*sizeof(unsigned long)); ni=(unsigned short*)malloc((FITS+1)*sizeof(unsigned short)); rd=(double*)malloc((NP+1)*sizeof(double)); rn=(double*)malloc((NP+1)*sizeof(double)); u=(double*)malloc((NP+1)*sizeof(double)); w=(double*)malloc((NP+1)*sizeof(double)); v=(float*)malloc((NV+1)*sizeof(float)); t=(float*)malloc((NT+1)*sizeof(float)); r=(float*)malloc((NP+1)*sizeof(float)); s=(float*)malloc((NP+1)*sizeof(float)); 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]=FITPARMS;} /* 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; }