static int CountErrors = 0; static int DontDump =0; /* GaussFit - A System for Least Squares and Robust Estimation Source Code Copyright (C) 1987 by William H. Jefferys, Michael J. Fitzpatrick and Barbara E. McArthur All Rights Reserved. */ /******************************************/ /*Householder Reduction System Functions*/ /******************************************/ #include #include #include #include #include "alloc.h" #include "house.h" #include "array.h" #include "datum.h" #include "simpledefs.h" #include "files.h" #include "prototypes.p" #include "symboltable.h" #include "defines.h" #ifdef THINK_C #include #endif extern double envhuber, envtukey, envfair, envtrim; extern int envminsum, envirls, envorm, envprec; char *PtoCstr(); unsigned char *CtoPstr(); double MAD(); /* Median Absolute Deviation Function */ FITS *getenvtabpntr(); #define Matrix(i,j)matrix[RowPtr[i]][ColPtr[j]] /* All arrays are alloc- */ #define RowType(i)rowtype[RowPtr[i]]/* ated dynamically and*/ #define ColType(i)coltype[ColPtr[i]]/* are accessed indirect-*/ #define Deltap(i) deltap[ColPtr[i]]/* ly through arrays*/ #define CondModulus(i)condmodulus[ColPtr[i]] /* ColPtr and RowPtr.*/ #define ConstModulus(i) constmodulus[ColPtr[i]] #define ConstTest(i)consttest[ColPtr[i]] #define CondTest(i)condtest[ColPtr[i]] #define Number(i,j) number[ColPtr[i]][j] #define NumberX(i) number[ColPtr[i]] #define VarName(i)varname[ColPtr[i]] extern double lambda; extern int MARQ; /* Allocate pointers for each matrix & vector */ static double **matrix = NULL; /* Matrix for Eqns of Condition */ static int *RowPtr = NULL; /* RowPtr matrix */ static int *ColPtr = NULL; /* ColPtr matrix */ static int *rowtype = NULL; /* Type of each row (condition, constraint) */ static int *coltype = NULL; /* Type of each col (parameter,indexed parameter, right hand side) */ static short (*number)[5] = NULL; /* index of a columns parameter */ static double *condmodulus = NULL;/* modulus of condition eqns in a column */ static double *constmodulus = NULL;/* modulus of condition eqns in a row */ static double *condtest = NULL; static double *consttest = NULL; static double *deltap = NULL; /* vector of parameter corrections */ static double *vector = NULL; static double *colvec = NULL; static double *sigvec = NULL; static double *sigval = NULL; static int matcol=MINCOLS; static int matrow=MINROWS; /* extern int totsize; extern int totrealloc; extern int totfree; */ double *resdum = NULL;/* For minsum routine */ int numresiduals; /* For minsum routine */ static char **varname = NULL; #ifdef THINK_C #define WIDTH 60 extern unsigned char *CtoPstr(); #else #define WIDTH 120 #endif typedef int INDEX[3]; /* arrays for "sort" of columns by pivot ability*/ static INDEX *indexarray=NULL; double *matstart,*matend; static MEMSIZE matsize; /* pointer to current matrix size */ static int LastNonzero; /* Last possible pivot row */ static int CurrentRow; /* The row being eliminated */ static int Recalculate; /* Flag to recalculate moduli */ static int LastCol; /* Last column active in matrix */ static double Tol; /* Acceptable Tolerance */ static int PrintFlag; /* Print intermediate results */ static int LastRow; /* Last row active in matrix */ static int FirstPivots; /* 0 = pivot plates first; 1 = pivot stars first; -1 = let program choose pivots */ static int Constraints; /* Number of constraints entered */ static int Touched = 0; /* Whether more info has been added */ static int noindex = 0; /* there are no indexed parameters */ static int lbl = 0; /* the indexes are large; put out line by line*/ static double DOF; /* Count of Degrees of Freedom */ static double NEQN; /* Count of Eqs of Condition */ double Sigma; /* Current value of Sigma */ double Sigma1; /* (extern) Current value of Sigma1 */ double Sigma1Sq; /* (extern) Current value of Sigma1Sq */ double SumRho = 0.0; /* (extern) Running Sum for Rho */ double SumPsi = 0.0; /* (extern) Running Sum for Psi */ double SumPsiSq = 0.0; /* (extern) Running Sum for PsiSq */ double SumPsiP = 0.0; /* (extern) Running Sum for PsiP */ double SumPsiPSq = 0.0; /* (extern) Running Sum for PsiPSq */ double ScaleFac = 1.0; /* (extern) ScaleFac */ double DeltaV = 0.0; /* (extern) Largest Residual Change */ double ChiSqr; /* Chi Square for real variances*/ static int StartAt = 0; /* Column to start pivoting */ int firstprint = 1; /* for printing out doubles */ FILE *fp,*fpc,*fpcv; double integral(); #define BOUNDS #ifdef BOUNDS int rowptr(i) int i; { int k; if((i>=matsize.nrows)||(i<0)) printf("Row Index Error %d, Rows = %d\n",i,matsize.nrows); k = RowPtr[i]; if((k>=matsize.nrows)||(k<0)) printf("Row Pointer Error %d, Rows = %d\n",k,matsize.nrows); return k; } int colptr(i) int i; { int k; if((i>=matsize.ncols)||(i<0)) printf("Column Index Error %d, Cols = %d\n",i,matsize.ncols); k = ColPtr[i]; if((k>=matsize.ncols)||(k<0)) printf("Column Pointer Error %d, Cols = %d\n",k,matsize.ncols); return k; } #else #define rowptr(i) RowPtr[i] #define colptr(i) ColPtr[i] #endif #define TOOSMALL 30 long totalalloc=0, totalfree=0, freememry0; char *MemAlloc(s,size)/* Allocate 'size' bytes */ char *s; long size; { char *x, **z; long freememry; unsigned usize; usize = (unsigned)size; /* error if not enough memory */ #ifdef DEBUGMEM if (totalalloc == 0) { freememry0 = FreeMem(); } totalalloc ++; freememry = FreeMem(); if (size>TOOSMALL) printf("memalloc : %s : size %ld %ld - %ld = %ld : used %ld free %ld", s, size, totalalloc, totalfree, totalalloc-totalfree, freememry0-freememry, freememry); #endif if((x = (char *)malloc(usize)) == NULL) /* Allocate the memory */ fatalerror("Memory Manager Error in Allocation--%s\n",s); /* error if not enough memory */ #ifdef DEBUGMEM if (size>TOOSMALL) printf(".\n"); #endif return x; } char *Reallocate(s,size,x)/* use for both initial and reallocation */ char *s; char *x; long size; { char *z; unsigned usize; void *vx; usize = (unsigned)size; #ifdef DEBUGMEM if (size>TOOSMALL) printf("realloc : %s : size %ld #[%ld]", s, size, totalalloc); #endif if (x == NULL) { fatalerror("Attempt to reallocate a NULL ptr in %s", s); } #ifdef USEANSI vx = (void*)x; if ((z = (char*)realloc(vx,usize)) == NULL)/* reallocate rowarray */ fatalerror("Memory Manager Error in Reallocation. Allocated pointer is NULL. --%s\n",s); #else if ((z = realloc(x,usize)) == NULL)/* reallocate rowarray */ fatalerror("Memory Manager Error in Reallocation. Allocated pointer is NULL. --%s\n",s); #endif #ifdef DEBUGMEM if (size>TOOSMALL) printf(".\n"); #endif return z; } freemem(s, ptr) char *s; char *ptr; { #ifdef DEBUGMEM totalfree ++; #endif #ifdef DEBUGFREE printf("freemem : %s : %ld - %ld = %ld : %ld", s, totalalloc, totalfree, totalalloc-totalfree, freememry); #endif /*totfree = totfree + sizeof(ptr);*/ if (ptr != NULL) { free((void *)ptr); } } double LSB()/*Computes the size of the LSB*/ /*of a word assuming the MSB is 1.0*/ { double x, y; x = 1.0; do { x = 0.5 * x; /* Divide by 2 repeatedly until no change is */ y = 1.0 + x; /* seen when result is added to 1 */ } while( y != 1.0); return 2.0 * x; } CountVars(type) /* Count the number of variables of type 'type' */ int type; { int k,count; count=0; /* initialize */ for(k=0;k<=LastCol;k++) /* do for each columsn that has a variable */ if(type==ColType(k)) /* if the desired type */ count++; /* increment count */ return count; } DumpColNames() { int i; fprintf(stderr, "\n"); fprintf(fp, "\n"); for (i=0; i<=LastCol; i++) { fprintf(fp, "%s %d %d[%d,%d,%d,%d]\n", VarName(i), ColType(i), Number(i,4), Number(i,0), Number(i,1), Number(i,2), Number(i,3)); fprintf(stderr, "%s %d %d[%d,%d,%d,%d]\n", VarName(i), ColType(i), Number(i,4), Number(i,0), Number(i,1), Number(i,2), Number(i,3)); } fprintf(stderr, "\n"); fprintf(fp, "\n"); } GetDeltaValue(k,name,index,type,value) /* Get data which was computed for */ int k; /* k'th parameter */ char **name; short *index; int *type; double *value; { if(k= matsize.ncols) { /* column doesn't exist; add to table */ matsize.lastcol = matsize.ncols; matsize.ncols = LastCol + 10;/*maxval(1.0,0.2*LastCol);*/ colspace(); matspace(); } ColType(LastCol) = type; /* enter data into tables- coltype filled*/ NumPtr = (short *)(number + ColPtr[LastCol]); NumPtr[0] = indices[0]; /* number filled */ NumPtr[1] = indices[1]; /* number filled */ NumPtr[2] = indices[2]; /* number filled */ NumPtr[3] = indices[3]; /* number filled */ NumPtr[4] = indices[4]; /* number filled */ VarName(LastCol) = name; /* varname filled */ /* Correct degrees of freedom if not RAS */ if(type == ParameterType) { DOF--; } /* DOF = n - k + r where n = # of eqs of cond k = # of parameters r = # of constraints */ return LastCol; /* return pointer to new column */ } double getdeltas(type,name,i) /* gets parameter correction from array */ char *name; int type; short *i; { int k; k = getcolumn(type,name,i); /* get the column number */ return Deltap(k); /* fetch parameter correction */ } PrintMatrix() { /* Print current matrix of eqns of condition */ int i, begin, width,inc; int colwid; colwid = getmaxprmname(); if (colwid <6) colwid = 6; width = (WIDTH-colwid)/(colwid+2); if (PrintFlag /*&& Touched*/) { /* Only print if matrix has been altered */ /* and user wants it */ begin = 0; inc = width; i = 1; while (inc < LastCol) { PrinttheMatrix(begin,inc,colwid); i++; begin = inc + 1; inc = i * width; } PrinttheMatrix(begin,LastCol,colwid); } } PrinttheMatrix(begin,end,colwid) int begin, end,colwid; { int i, j,ndim; double matint; fprintf(fp,"%*s ",colwid," "); for (j = begin;j<=end;j++) /* print column number */ fprintf(fp,"%*d ",colwid,colptr(j)); fprintf(fp,"\n"); fprintf(fp,"%*s ",colwid,"Type"); /* print name of row */ for (j = begin;j<=end;j++) /* print each item of vector */ fprintf(fp,"%*d ",colwid,coltype[colptr(j)]); fprintf(fp,"\n"); /* print carriage return */ fprintf(fp,"%*s ",colwid,"Name"); /* print name of each variable*/ for (j = begin;j<=end;j++) fprintf(fp,"%*s ",colwid,VarName(j)); fprintf(fp,"\n"); ndim = getdimnum(); for (j = 0; j= matsize.nrows) /* Check that row exists*/ { matsize.lastrow = matsize.nrows; matsize.nrows = LastRow + 10;/*LastRow*maxval(1.0,0.2*LastRow);*/ rowspace(); matspace(); } if(eqtype == Constraint)/* increment constraints if it's a constraint */ Constraints += 1; RowType(LastRow) = eqtype;/* fix row type */ for (i = 0;i= ThisCol */ /*This is an order(n) sort for the special case we have here.*/ { int n[3]; /* number of rows of each type */ int i, j, k; for (i = 0;i<= 2;i++) /* zero out number vector */ n[i] = 0; /* collect the possible pivots of each type */ for (i = ThisCol;i<=LastRow;i++) { j = Modulus(ThisCol, i); indexarray[n[j]][j] = rowptr(i);/*put row pointers into array*/ n[j]++; /* increment number of pivots of this type */ } k = ThisCol; /* Order rows so nonzero pivot constraints come before nonzero pivot conditions which come before zero pivot rows */ for (j = 2;j>=0;j--) for (i = 0;i 0) LastNonzero = k; k++; } } ApplyConstraintTransform (column) int column; /*This is actually just an elementary Gaussian elimination */ /*on each row that has a leading nonzero element.This*/ /*is the limit of a Givens transformation on that row as */ /*the weight of the pivot row tends to infinity. */ { int j, k, colj,rowclm; double divisor, s, Old, New, Delta,f1; if (notzero(Matrix(column,column))) divisor = 1.0 / Matrix(column,column); /* compute reciprocal of pivot*/ else fatalerror("Apply Constraint Transformation Division by zero--\n"," "); /* error if attempt to divide by zero */ for (j = column + 1;j<=LastCol;j++) /* Do for all columns */ { s = divisor * Matrix(column,j); /* compute normalized element from pivot row */ if (s != 0.0) /* only need to do it for nonzero element */ { for (k = column + 1;k<=LastNonzero;k++)/*do for each col */ { Old = Matrix(k,j);/*subtract element from row */ Delta = s * Matrix(k,column); New = Old - Delta; Matrix(k,j) = New; if (RowType(k) == Condition) CondModulus(j) -= /*CondModulus(j)update condition Modulus -*/ Delta * (Old + New); else ConstModulus(j) -= /*ConstModulus(j)update constraint modulus -*/ Delta * (Old + New); } /* see Lawson and Hanson for this step */ ConstModulus(j) = RoundToZero(ConstModulus(j) - sqr(Matrix(column,j)),ConstTest(j)); } } } ApplyHouseholderTransform (column) int column; /*This is straight out of Lawson and Hanson, except that */ /*we check to see if s==0 and don"t transform a column */ /*if that is the case. */ { int j, k; int columnCP, columnRP, jCP; double divisor, s, s2, h, b, Result,f1, f2; double *matRow; columnCP = ColPtr[column]; columnRP = RowPtr[column]; s = 0.0; for (j = column;j<=LastNonzero;j++) { /* Form modulus of the pivot column */ s2 = matrix[RowPtr[j]][columnCP]; s += s2 * s2; } s = sqrt(s); /* square root of the sum of the squares */ h = matrix[columnRP][columnCP]; /* compute b (reflection coefficient) */ if (h > 0.0) s = -s; h -= s; matrix[columnRP][columnCP] = s; b = s * h; if (notzero(b)) { /* ONly do it if b != 0 */ divisor = 1.0 / b; /* Zero out the column and replace it with the orthogonal transformation vector */ for (j = column + 1;j<=LastCol;j++) { jCP = ColPtr[j]; s = h * matrix[columnRP][jCP]; for (k = column + 1;k<=LastNonzero;k++) { matRow = matrix[RowPtr[k]]; s += matRow[columnCP] * matRow[jCP]; } s = s * divisor; if (s != 0.0) { /*no need to transform the column if s==0*/ Result = matrix[columnRP][jCP] += s * h; f1 = condmodulus[jCP] - sqr(Result); f2 = condtest[jCP]; condmodulus[jCP] = RoundToZero(f1,f2); for (k = column + 1;k<=LastNonzero;k++) { matRow = matrix[RowPtr[k]]; matRow[jCP] += s * matRow[columnCP]; } } } } } ApplyTransformation (ThisCol) /* Apply transformation to each column */ int ThisCol; { SortPivotColumn(ThisCol);/* First sort the column to get the pivot */ if (RowType(ThisCol) == Constraint) /* if pivot is constrained */ ApplyConstraintTransform(ThisCol);/* apply constraint transform */ else /* otherwise appy HouseHolder transform */ ApplyHouseholderTransform(ThisCol); } CalculateModuli() /* Modulus of each column */ { double ConstSum, CondSum; int i, j; for (j = CurrentRow;j-for pivoting. */ { double x; if(x = TypeModulus(i) - TypeModulus(j)) return x; /* Constraint >- Condition >- _RHS */ else if(x = NameModulus(i) - NameModulus(j)) return x;/* Subscripted >- Globals */ else return CCModulus(i) - CCModulus(j); /* largest >- smallest */ } FindPivotColumn (ThisCol) int ThisCol; /*Pivoting strategy is:*/ /*Never pivot on RHS */ /*Constraints before condition equations*/ /*Low-numbered subscripts after high-numbered*/ /*Globals after subscripted variables*/ /*Large column (measured by Sum-Of-Squares) before small */ { int i, maximum; maximum = ThisCol;/* limits of search */ CurrentRow = ThisCol; if (Recalculate)/* recalculate moduli if significance lost */ CalculateModuli(); for (i = ThisCol + 1;i<=LastCol;i++) /* get maximum pivot */ if (CompareCol(maximum,i) < 0.0) maximum = i; SwitchCol(ThisCol, maximum);/* switch it into first place */ } double sign(dvalu) double dvalu; { if (dvalu < 0.0) return (-1.0); else return (1.0); } SolveLinearSystem() { int i, k; double s, smallest, largest,whichsign; if(LastRow=0;i--) /* Go from last to first column */ { s = Matrix(i,LastCol); /* get RHS */ for (k = LastCol-1;k>i;k--) /* Subtract terms already solved for */ s -= Matrix(i,k) * Deltap(k); if (Matrix(i,i) == 0.0) /* Check if matrix is singular */ { fatalerror("Matrix Singular!\n",""); Deltap(i) = 0.0; /* This line can never be activated. */ } else { /*if nonsingular, compute parameter correction */ Deltap(i) = s / Matrix(i,i);/* deltap filled */ } /* update smallest and largest matrix element*/ /*smallest = minval(smallest, fabs(Matrix(i,i))); largest = maxval(largest, fabs(Matrix(i,i)));*/ } for (i = LastCol-1;i>=0;i--) /* Go from last to first column */ { if (colvec[i] != 0.0) Deltap(i) = Deltap(i)/colvec[i]; } /* estimate condition as ratio */ /*fprintf(fp,"Condition = %7.1le\n\n", largest/maxval(smallest,Tol));*/ } SumIt() { double KFAC; if(DOF <= 0.0) /* can't solve if DOF <= 0/0 */ fatalerror("The observation set is underdetermined!\n",""); SumRho *= 2.0/DOF; /* Sums for this iteration */ SumPsi /= (NEQN - Constraints); SumPsiSq /= DOF; SumPsiP /= (NEQN - Constraints); SumPsiPSq /= (NEQN - Constraints); /* compute new sigma and scale factor */ KFAC = 1.0 + (NEQN - Constraints-DOF)/(NEQN - Constraints)*(SumPsiPSq - sqr(SumPsiP)); Sigma1Sq = KFAC*SumPsiSq/SumPsiP; Sigma1 = sqrt(Sigma1Sq); Sigma = Sigma1*(ScaleFac==0.0?1.0:ScaleFac); putenvval("sigma",Sigma); ScaleFac = (ScaleFac==0.0?1.0:ScaleFac)*sqrt(SumRho); putenvval("scale",ScaleFac); /* Print results of iteration */ printiter_results(stdout,KFAC); printiter_results(fp,KFAC); } printiter_results(fp,KFAC) FILE *fp; double KFAC; { fprintf(fp,"scale =\t"); printdouble(fp,ScaleFac,0); fprintf(fp,"\n"); fprintf(fp,"SumRho =\t"); printdouble(fp,SumRho,0); fprintf(fp,"\n"); fprintf(fp,"SumPsi =\t"); printdouble(fp,SumPsi,0); fprintf(fp,"\n"); fprintf(fp,"SumPsiSq =\t"); printdouble(fp,SumPsiSq,0); fprintf(fp,"\n"); fprintf(fp,"SumPsiP =\t"); printdouble(fp,SumPsiP,0); fprintf(fp,"\n"); fprintf(fp,"SumPsiPSq =\t"); printdouble(fp,SumPsiPSq,0); fprintf(fp,"\n\n"); fprintf(fp,"NEQN = %5d, DOF = %5d\n",(int)NEQN,(int)DOF); fprintf(fp,"KFAC =\t "); printdouble(fp,KFAC,0); fprintf(fp,"\n"); fprintf(fp,"Sigma =\t "); printdouble(fp,Sigma,0); fprintf(fp,", DOF = %d\n\n",(int)DOF); } CovarianceMtx() { int i,j,k,width,begin,inc; int indsize, mxprnam, colwid; double s; if(envminsum) { /* see if covariance matrix can be provided */ printf("Sorry, can't get covariance matrix when doing minsum\n"); return 0; /*Can't get cov mtx when doing minsum */ } if(envirls) { if ((envfair) || (envtrim) || (envtukey) || (envhuber) || (envminsum)) { printf("Sorry, covariance matrix not implemented for robust estimation using IRLS\n"); return 0; /* Can't get covariance mtx when doing IRLS */ } } printf("Calculating Covariance Matrix\n"); for(k=LastCol-1; k>=0; k--) { /* Do for each column */ for (i=0; i=0;i--) { /* compute inverse of upper triang. matrix */ s = deltafn(i,k)*deltafn(RowType(i),Condition); /* Handle covariance of condition eqn */ /* rhs vector is zero for a constraint, 1 on diagonal for a condition */ for (j=k; j>i; j--) s -= Matrix(i,j) * vector[j]; if (Matrix(i,i) == 0.0) { /* Back substitute the column */ fprintf(fp,"Matrix Singular!\n"); vector[i] = 0.0; /* vector filled */ } else /* recompute column of inverse matrix */ vector[i] = s / Matrix(i,i); /* vector filled */ } for(i=0;i=0; k--) /* multiply matrix by it's transpose */ for(i=0; i<=k; i++) { s = 0.0; for(j=k; j 8) lbl = 1; if (indsize == 0) { indsize = mxprnam; noindex = 1; } if (lbl) colwid = 8; else colwid = indsize; if (lbl) width = WIDTH-(mxprnam+8)/colwid+2; else width = WIDTH-(mxprnam+indsize)/colwid+2; begin = 0; inc = width; i = 1; prSigma(indsize,mxprnam); while (inc j) { fprintf(fptr,"%*s\t", colwid,""); } else { if (fabs(matint) < 1.0E9) { if (mult >10) fprintf(fptr,"%*d\t",colwid,(int)matint); else fprintf(fptr,"%le\t",matint); } } } else { /* it's a covariance matrix */ if (i >j) { matint = (double)mult * (Matrix(j,i)/div); fprintf(fptr,"%le\t",matint*sigval[i]*sigval[j]); } else { fprintf(fptr,"%le\t",matint*sigval[i]*sigval[j]); } } } fprintf(fptr,"\n"); } fprintf(fptr,"\n"); } prSigma(indsize,mxprnam) int indsize, mxprnam; { int i,j; short snum[5]; fprintf(fp,"\n"); fprintf(fp, "Sigma Values \n"); fprintf(fp,"\n"); for(i=0;i 4) { fprintf(stderr, "\nprintIndex: Number of dimensions = %d\n",indx[4]); exit(-1); } */ strcpy(string,"["); for (i=0; i 4) { fprintf(stderr, "\nprintIndex: Number of dimensions = %d\n",indx[4]); exit(-1); } sprintf(tiny, "[%d]",indx[indx[4]-1-level]); fprintf(fd,"%*s",indsize,tiny); } printdouble(fp,dvalue,preon) double dvalue; FILE *fp; int preon; { double floor() ; double truncnum,diffround,tnum; char dnum[30],dnum2[30]; int i, j; static int prec = 18; static int file_prec = 18; if ((dvalue < UUND) && (dvalue > LUND)) { fprintf(fp,"#N/A\t"); return; } if (firstprint) { firstprint = 0; prec = envprec; file_prec = machine_prec(); if ((file_prec < 14) || (file_prec > 30)) file_prec = IEEE_DBL_DIG; if ((prec < 4) || (prec > 30)) prec = file_prec; } if ((prec != file_prec) && (preon == 0)) fprintf(fp,"% .*g\t",prec,dvalue); else { if (((fabs(dvalue) < 1.0E-10) || (fabs(dvalue) > 1.0E+9)) && (dvalue != 0.0)) fprintf(fp,"% .*e\t",file_prec,dvalue); else { sprintf(dnum,"% .*f",file_prec,dvalue); i = strlen(dnum)-1; if (IEEE_DBL_DIG < i) i = IEEE_DBL_DIG; if ((dnum[i] != '0') && (dnum[i-1] == '0')) i--; while (dnum[i] == '0') i--; i++; if (dnum[i-1] == '.') i++; if (dvalue < 0.0) i++; dnum[i] = '\0'; fprintf(fp,"%s\t",dnum); } } } getindexsz() { int indxsiz,i,j,k; long pow10; int max[4]; static int sumsize = -1; if (sumsize == -1) { sumsize = 0; for (i=0;i<4; i++) max[i] = -1; for (i=0;i max[j]) max[j] = Number(i,j); /* find size of indices */ for (i=0; i<4;i++) { if (max[i] > -1) { k = 1; pow10 = 10; while (max[i] >= pow10) { k++; pow10 *= 10; } sumsize = sumsize + k; } } /* add space for commas*/ for (i=0;i<4;i++) if (max[i] > -1) sumsize++; if (sumsize > 0) sumsize ++; } return (sumsize); } getmaxprmname() { int length,i; static int maxp = -1; if (maxp == -1) { maxp = 0; for (i = 0; i< LastCol; i++) if (maxp < (length = strlen(varname[i]))) maxp =length; } return maxp; } getdimnum() { int i; static int maxd = -1; if (maxd == -1) { maxd = 0; for (i = 0; i< LastCol; i++) if (maxd < (Number(i,4))) maxd = Number(i,4); } return maxd; } Solve1() /* Solve matrix as it exists so far */ { int ThisCol; int i,j; int EndAt; if(!Touched) return; /* Matrix has already been transformed */ FirstPivots = 0;/* we now always start at zero */ Recalculate = 1; /* recalculate moduli after solution */ /*fprintf(fp,"\nTransformation Number ");*/ EndAt = minval(LastRow,LastCol); /* last row to transform */ for (ThisCol = StartAt;ThisCol<=EndAt;ThisCol++) /* transform each column */ { /*fprintf(fp,", %3d", ThisCol);*/ if ((ThisCol % 20) == 0) fprintf(fp,"\n"); FindPivotColumn(ThisCol); /* Find pivot columns and */ ApplyTransformation(ThisCol);/* apply appropriate transformation */ } for(i=1;i<=EndAt;i++) /* zero out the lower triangle */ for(j=0;ji && RowType(j) >= pivtype;) j = j-1; if(i= matsize.nrows) { matsize.lastrow = matsize.nrows; matsize.nrows = numresiduals + 10;/*numresiduals*maxval(1.0,0.2*numresiduals);*/ rowspace(); matspace(); } resdum[numresiduals]= x; numresiduals++; } dump(s) char *s; { int i, j,k,l; int collimit, rowlimit; if(DontDump) return; if(CountErrors>100) return; CountErrors++; fprintf(fp,s); fprintf(fp,", Entry # %d\n",CountErrors); collimit = matsize.ncols; rowlimit = minval(LastRow,30); fprintf(fp,"Matrix\n"); for (j = 0; j