lcsim/src/org/lcsim/recon/emid/hmatrix
diff -N HMatrix.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HMatrix.java 24 Jun 2005 18:11:01 -0000 1.1
@@ -0,0 +1,503 @@
+package org.lcsim.recon.emid.hmatrix;
+import java.io.*;
+/**
+ *This class allows one to calculate how well a measurement agrees with expectations.
+ *The vector of average values and the associated covariance matrix is used to construct
+ *the HMatrix. The degree of agreement is calculated as a chi-squared of a measurement
+ *vector against the expectations.
+ *
+ *
+ * The covariance matrix is defined as:
+ * <br clear="all" /><table border="0" width="100%"><tr><td>
+ * <table align="center"><tr><td nowrap="nowrap" align="center">
+ * <i>M</i><sub><i>ij</i></sub> = </td><td nowrap="nowrap" align="center">
+ * 1
+ * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>N</i><br /></td><td nowrap="nowrap" align="center">
+ * </td><td nowrap="nowrap" align="center">
+ * <font size="-1"><i>N</i></font><!--sup
+ * --><br /><font size="+3"><font face="symbol">å<br />
+ * </font></font><font size="-1"><i>n</i> = 1</font> <br /></td><td nowrap="nowrap" align="center">
+ * (<i>E</i><sub><i>i</i></sub><sup>(<i>n</i>)</sup> <font face="symbol">-</font
+ * ></td><td nowrap="nowrap" align="center">
+ *
+ * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>
+ * <div class="norm"><i>E</i><br /></div>
+ * <div class="comb"> </div>
+ * </td><td nowrap="nowrap" align="center">
+ * <font size="-1"></font><!--sup
+ * --><br /><font size="-1"><i>i</i></font> <br /></td><td nowrap="nowrap" align="center">
+ * )(<i>E</i><sub><i>j</i></sub><sup>(<i>n</i>)</sup> <font face="symbol">-</font
+ * ></td><td nowrap="nowrap" align="center">
+ *
+ * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>
+ * <div class="norm"><i>E</i><br /></div>
+ * <div class="comb"> </div>
+ * </td><td nowrap="nowrap" align="center">
+ * <font size="-1"></font><!--sup
+ * --><br /><font size="-1"><i>j</i></font> <br /></td><td nowrap="nowrap" align="center">
+ * )</td></tr></table>
+ * </td></tr></table>
+ *
+ * The H(essian)Matrix is the inverse of the covariance matrix.
+ * The chi-squared is then calculated <it>via</it>
+ * <br clear="all" /><table border="0" width="100%"><tr><td>
+ * <table align="center"><tr><td nowrap="nowrap" align="center">
+ * <font face="symbol">z</font
+ * ><sub><i>m</i></sub> <font face="symbol">º</font
+ * > </td><td nowrap="nowrap" align="center">
+ * <font size="-1"><i>N</i></font><!--sup
+ * --><br /><font size="+3"><font face="symbol">å<br />
+ * </font></font><font size="-1"><i>i</i>,<i>j</i> = 1</font> <br /></td><td nowrap="nowrap" align="center">
+ * (<i>E</i><sub><i>i</i></sub><sup>(<i>m</i>)</sup> <font face="symbol">-</font
+ * ></td><td nowrap="nowrap" align="center">
+ *
+ * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>
+ * <div class="norm"><i>E</i><br /></div>
+ * <div class="comb"> </div>
+ * </td><td nowrap="nowrap" align="center">
+ * <font size="-1"></font><!--sup
+ * --><br /><font size="-1"><i>i</i></font> <br /></td><td nowrap="nowrap" align="center">
+ * )<i>H</i><sub><i>ij</i></sub> (<i>E</i><sub><i>j</i></sub><sup>(<i>m</i>)</sup> <font face="symbol">-</font
+ * ></td><td nowrap="nowrap" align="center">
+ *
+ * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>
+ * <div class="norm"><i>E</i><br /></div>
+ * <div class="comb"> </div>
+ * </td><td nowrap="nowrap" align="center">
+ * <font size="-1"></font><!--sup
+ * --><br /><font size="-1"><i>j</i></font> <br /></td><td nowrap="nowrap" align="center">
+ * )</td></tr></table>
+ * </td></tr></table>
+ *
+ *Methods to persist the HMatrix are provided in both
+ *plain ASCII format and Java serialized format.
+ *
+ *@author Norman A. Graf
+ *@version 1.0
+ */
+public class HMatrix implements Serializable
+{
+ private double[][] _invcov; // the inverse covariance matrix
+ private double[] _vec; // the vector of measurements
+ private int _dim; // the dimensionality
+ private int _key; // the key for the HMatrix
+ private double[] _tmp; // temporary array
+ private double[] _tmp2; // temporary array
+
+
+ /**
+ * Constructor
+ *
+ * @param n the dimensionality of the matrix
+ * @param key the key for indexing this HMatrix
+ * @param vec The array of average values for the measurement space.
+ * @param cov The inverse of the covariance matrix for the averages.
+ */
+ public HMatrix(int dim, int key, double [] vec, double[][] cov)
+ {
+ _dim = dim;
+ _key = key;
+ _vec = new double[_dim];
+ System.arraycopy(vec, 0, _vec, 0, _dim);
+ _invcov = new double[_dim][_dim];
+ for(int i =0; i<_dim; ++i)
+ {
+ System.arraycopy(cov[i], 0, _invcov[i], 0, _dim);
+ }
+ _tmp = new double[_dim];
+ _tmp2 = new double[_dim];
+
+ }
+
+
+ public HMatrix(String resourceFileName)
+ {
+ InputStream in = this.getClass().getResourceAsStream(resourceFileName);
+ if(in!=null)
+ {
+ try
+ {
+ BufferedReader bufferedreader
+ = new BufferedReader(new InputStreamReader(in));
+ create(bufferedreader);
+ bufferedreader.close();
+
+ }
+ catch(IOException _ex)
+ {
+ System.err.println("HMatrixBuilder::read -> Error reading HMatrix from input reader.");
+ System.exit(0);
+ }
+ }
+ }
+
+ /**
+ * Calculates the chi-squared for the measurement compared to
+ * the expected values and covariances represented by the HMatrix.
+ *
+ * @param dat The array of measured values
+ * @return The value of chi-squared.
+ */
+ public double chisquared(double[] dat)
+ {
+ // vector of measured-predicted values
+ for(int i=0; i<_dim; ++i)
+ {
+ _tmp[i] = dat[i]-_vec[i];
+ }
+
+ double chisqr = 0.;
+ // covariance matrix times difference vector
+ for(int i=0; i<_dim; ++i)
+ {
+ _tmp2[i] = 0.;
+ for(int j=0; j<_dim; ++j)
+ {
+ _tmp2[i]+=_invcov[j][i]*_tmp[j];
+ }
+ chisqr += _tmp[i]*_tmp2[i];
+ }
+ return chisqr;
+ }
+
+
+ /**
+ *Output stream
+ *
+ * @return A String representation of the object.
+ */
+ public String toString()
+ {
+ StringBuffer sb = new StringBuffer("HMatrix: dimension "+_dim+" key "+_key+ "\n");
+ sb.append("vector: \n");
+ for(int i = 0; i<_dim; ++i)
+ {
+ sb.append(_vec[i]+" ");
+ }
+ sb.append("\n\nInverse covariance matrix: \n");
+ for(int i = 0; i<_dim; ++i)
+ {
+ for(int j = 0; j<_dim; ++j)
+ {
+ sb.append(_invcov[i][j]+" ");
+ }
+ sb.append("\n");
+ }
+ return sb.toString();
+ }
+
+ /**
+ * Writes the HMatrix to a Serialized file
+ */
+ public void writeSerialized(String filename)
+ {
+
+ try
+ {
+ FileOutputStream fos = new FileOutputStream(filename);
+ ObjectOutputStream objOs = new ObjectOutputStream(fos);
+ objOs.writeObject(this);
+ objOs.flush();
+ fos.close();
+ } catch (NotSerializableException se)
+ {
+ System.err.println(se);
+ } catch (FileNotFoundException fe)
+ {
+ System.err.println(fe);
+ } catch (IOException se)
+ {
+ System.err.println(se);
+ }
+
+ }
+
+ /**
+ * Reads the HMatrix from a Serialized file
+ */
+ public static HMatrix readSerialized(String filename)
+ {
+ HMatrix hm = null;
+
+ try
+ {
+ FileInputStream fis = new FileInputStream(filename);
+ ObjectInputStream objIs = new ObjectInputStream(fis);
+ hm = (HMatrix)objIs.readObject();
+ fis.close();
+ } catch (NotSerializableException se)
+ {
+ System.err.println(se);
+ } catch (FileNotFoundException fe)
+ {
+ System.err.println(fe);
+ } catch (IOException se)
+ {
+ System.err.println(se);
+ } catch (ClassNotFoundException ce)
+ {
+ System.err.println(ce);
+ }
+ return hm;
+ }
+
+ /**
+ * Writes out the HMatrix to an ASCII file
+ */
+ public void write(String filename, String comment)
+ {
+ System.out.println("Writing matrix to '" + filename + "'.");
+ try
+ {
+ PrintWriter printwriter = new PrintWriter(new BufferedWriter(new FileWriter(filename)));
+ // Print the comment
+ String s2 = "# " + comment;
+ printwriter.println(s2);
+ // Now the dimension and index of the matrix
+ String s3 = _dim+" "+_key;
+ printwriter.println(s3);
+ String s4 = "";
+ //Now the vector of average values
+ {
+ String s5 = "";
+ for(int i = 0; i < _dim; i++) s5 = s5 + _vec[i] + " ";
+ printwriter.println(s5.substring(0, s5.length() - 1));
+ }
+ //Now the covariance matrix
+ for(int i = 0; i < _dim; i++)
+ {
+ String s6 = "";
+ for(int j = 0; j < _dim; j++)
+ s6 = s6 + _invcov[i][j] + " ";
+
+ printwriter.println(s6.substring(0, s6.length() - 1));
+ }
+
+ printwriter.close();
+ }
+ catch(IOException _ex)
+ {
+ System.err.println("Matrix::write -> Error writing to '" + filename + "'.");
+ System.exit(0);
+ }
+ System.out.println("Matrix written.");
+ }
+
+ /**
+ * Reads the HMatrix from an ASCII file
+ */
+ public static HMatrix read(String filename)
+ {
+ /*
+ boolean flag = false;
+ boolean readVec = false;
+ int i = 0;
+ int ii = 0;
+ int j = 0;
+ int dim=0;
+ int key = 0;
+ double[] vec=null;
+ double[][] cov=null;
+ */
+ try
+ {
+ BufferedReader bufferedreader = new BufferedReader(new FileReader(filename));
+ HMatrix mat = create(bufferedreader);
+ /*
+ for(String s1 = new String(); (s1 = bufferedreader.readLine()) != null;)
+ if(s1.startsWith("#"))
+ System.out.println(s1);
+ else
+ if(!flag)
+ {
+ dim = Integer.parseInt(s1.substring(0, s1.indexOf(" ")));
+ key = Integer.parseInt(s1.substring(s1.indexOf(" ") + 1, s1.length()));
+ vec = new double[dim];
+ cov = new double[dim][dim];
+ flag = true;
+ }
+ else
+ {
+ if (!readVec)
+ {
+
+ // Read the vector of averages
+ while(s1.length() > 0 && s1.indexOf(" ") != -1)
+ {
+ double d = Double.valueOf(s1.substring(0, s1.indexOf(" "))).doubleValue();
+ s1 = s1.substring(s1.indexOf(" ") + 1, s1.length());
+ vec[ii] = d;
+ ii++;
+ }
+ double d1 = Double.valueOf(s1).doubleValue();
+ vec[ii] = d1;
+ readVec = true;
+ }
+ else
+ {
+ // Read the covariance matrix
+ while(s1.length() > 0 && s1.indexOf(" ") != -1)
+ {
+ double d = Double.valueOf(s1.substring(0, s1.indexOf(" "))).doubleValue();
+ s1 = s1.substring(s1.indexOf(" ") + 1, s1.length());
+ cov[i][j] = d;
+ j++;
+ }
+ double d1 = Double.valueOf(s1).doubleValue();
+ cov[i][j] = d1;
+ i++;
+ j = 0;
+ }
+ }
+ */
+ bufferedreader.close();
+ return mat;
+ }
+ catch(IOException _ex)
+ {
+ System.err.println("HMatrixBuilder::read -> Error reading '" + filename + "#.");
+ System.exit(0);
+ }
+ // return new HMatrix(dim, key, vec, cov);
+ return null;
+ }
+
+ /**
+ * Reads the HMatrix from a Reader
+ */
+ public static HMatrix read(Reader reader)
+ {
+ /*
+ boolean flag = false;
+ boolean readVec = false;
+ int i = 0;
+ int ii = 0;
+ int j = 0;
+ int dim=0;
+ int key = 0;
+ double[] vec=null;
+ double[][] cov=null;
+ */
+ try
+ {
+ BufferedReader bufferedreader = new BufferedReader(reader);
+ HMatrix mat = create(bufferedreader);
+ /*
+ for(String s1 = new String(); (s1 = bufferedreader.readLine()) != null;)
+ if(s1.startsWith("#"))
+ System.out.println(s1);
+ else
+ if(!flag)
+ {
+ dim = Integer.parseInt(s1.substring(0, s1.indexOf(" ")));
+ key = Integer.parseInt(s1.substring(s1.indexOf(" ") + 1, s1.length()));
+ vec = new double[dim];
+ cov = new double[dim][dim];
+ flag = true;
+ }
+ else
+ {
+ if (!readVec)
+ {
+
+ // Read the vector of averages
+ while(s1.length() > 0 && s1.indexOf(" ") != -1)
+ {
+ double d = Double.valueOf(s1.substring(0, s1.indexOf(" "))).doubleValue();
+ s1 = s1.substring(s1.indexOf(" ") + 1, s1.length());
+ vec[ii] = d;
+ ii++;
+ }
+ double d1 = Double.valueOf(s1).doubleValue();
+ vec[ii] = d1;
+ readVec = true;
+ }
+ else
+ {
+ // Read the covariance matrix
+ while(s1.length() > 0 && s1.indexOf(" ") != -1)
+ {
+ double d = Double.valueOf(s1.substring(0, s1.indexOf(" "))).doubleValue();
+ s1 = s1.substring(s1.indexOf(" ") + 1, s1.length());
+ cov[i][j] = d;
+ j++;
+ }
+ double d1 = Double.valueOf(s1).doubleValue();
+ cov[i][j] = d1;
+ i++;
+ j = 0;
+ }
+ }
+ */
+ bufferedreader.close();
+ return mat;
+ }
+ catch(IOException _ex)
+ {
+ System.err.println("HMatrixBuilder::read -> Error reading HMatrix from input reader.");
+ System.exit(0);
+ }
+ // return new HMatrix(dim, key, vec, cov);
+ return null;
+ }
+
+ private static HMatrix create(BufferedReader bufferedreader) throws IOException
+ {
+ boolean flag = false;
+ boolean readVec = false;
+ int i = 0;
+ int ii = 0;
+ int j = 0;
+ int dim=0;
+ int key = 0;
+ double[] vec=null;
+ double[][] cov=null;
+ for(String s1 = new String(); (s1 = bufferedreader.readLine()) != null;)
+ if(s1.startsWith("#"))
+ System.out.println(s1);
+ else
+ if(!flag)
+ {
+ dim = Integer.parseInt(s1.substring(0, s1.indexOf(" ")));
+ key = Integer.parseInt(s1.substring(s1.indexOf(" ") + 1, s1.length()));
+ vec = new double[dim];
+ cov = new double[dim][dim];
+ flag = true;
+ }
+ else
+ {
+ if (!readVec)
+ {
+
+ // Read the vector of averages
+ while(s1.length() > 0 && s1.indexOf(" ") != -1)
+ {
+ double d = Double.valueOf(s1.substring(0, s1.indexOf(" "))).doubleValue();
+ s1 = s1.substring(s1.indexOf(" ") + 1, s1.length());
+ vec[ii] = d;
+ ii++;
+ }
+ double d1 = Double.valueOf(s1).doubleValue();
+ vec[ii] = d1;
+ readVec = true;
+ }
+ else
+ {
+ // Read the covariance matrix
+ while(s1.length() > 0 && s1.indexOf(" ") != -1)
+ {
+ double d = Double.valueOf(s1.substring(0, s1.indexOf(" "))).doubleValue();
+ s1 = s1.substring(s1.indexOf(" ") + 1, s1.length());
+ cov[i][j] = d;
+ j++;
+ }
+ double d1 = Double.valueOf(s1).doubleValue();
+ cov[i][j] = d1;
+ i++;
+ j = 0;
+ }
+ }
+ return new HMatrix(dim, key, vec, cov);
+ }
+
+}
\ No newline at end of file
lcsim/src/org/lcsim/recon/emid/hmatrix
diff -N HMatrixBuilder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HMatrixBuilder.java 24 Jun 2005 18:11:01 -0000 1.1
@@ -0,0 +1,304 @@
+package org.lcsim.recon.emid.hmatrix;
+import java.io.*;
+/** A class to construct an HMatrix.
+ *
+ *@author Norman A. Graf
+ *@version 1.0
+ *@see HMatrix
+ */
+public class HMatrixBuilder
+{
+ private double[][] _cov;
+ private double[] _vec;
+ private double[] _vals;
+ private double[][] _valsq;
+ private int _dim;
+ private int _key;
+ private int _ndata;
+ private boolean _isValid;
+ private HMatrix _hmx;
+
+ /** Construct an <b>n</b> x <b>n</b> HMatrixBuilder.
+ * Leaves HMatrix in an invalid state.
+ * One must either accumulate entries to build up the HMatrix
+ * or read in the HMatrix values from a file.
+ *
+ * @param n The dimension of the measurement space.
+ * @param key The key by which to index this HMatrix.
+ */
+ public HMatrixBuilder(int n, int key)
+ {
+ _dim = n;
+ _key = key;
+ _cov = new double[_dim][_dim];
+ _vec = new double[_dim];
+ _vals = new double[_dim];
+ _valsq = new double[_dim][_dim];
+ _ndata = 0;
+ _isValid = false;
+ }
+
+ /**
+ * Add the measurement vector to the accumulated HMatrix
+ *
+ * @param dat The array of measurements.
+ */
+ public void accumulate(double[] dat)
+ {
+ //if(dat.length !=_dim) throw new InvalidArgumentException();
+ for(int i=0; i<_dim; ++i)
+ {
+ for(int j=0; j<_dim; ++j)
+ {
+ _valsq[i][j]+=dat[i]*dat[j];
+ if(i==j) _vals[i]+=dat[i];
+ //System.out.println("i= "+i+", j= "+j+"val= "+_vals[i]+", valsq= "+_valsq[i][j]);
+ }
+ }
+ _ndata++;
+ }
+
+ /**
+ * Validates the HMatrix by performing the averages
+ */
+ public void validate()
+ {
+ for(int i=0; i<_dim; ++i)
+ {
+ _vec[i] = _vals[i]/_ndata;
+ for(int j=0; j<_dim; ++j)
+ {
+ double data = (double)_ndata;
+ _cov[i][j] = _valsq[i][j]/(data) - (_vals[i]/(data)*(_vals[j]/(data)));
+ }
+ }
+
+ _hmx = new HMatrix(_dim, _key, _vec, invert(_cov,_dim));
+ _isValid = true;
+ }
+
+ /**
+ * Output Stream
+ *
+ * @return String representation of the object
+ */
+ public String toString()
+ {
+ StringBuffer sb = new StringBuffer("HMatrixBuilder: dimension "+_dim+" ");
+ if(!_isValid) return (sb.append("INVALID").toString());
+ sb.append(_hmx);
+/* sb.append("\n\nvector: \n");
+ for(int i = 0; i<_dim; ++i)
+ {
+ sb.append(_vec[i]+" ");
+ }
+ sb.append("\n\ncovariance matrix: \n");
+ for(int i = 0; i<_dim; ++i)
+ {
+ for(int j = 0; j<_dim; ++j)
+ {
+ sb.append(_cov[i][j]+" ");
+ }
+ sb.append("\n");
+ }
+ */
+ return sb.toString();
+ }
+
+ /**
+ * Writes the HMatrix to a Serialized file
+ *
+ * @param filename The file to which to write.
+ */
+ public void writeSerialized(String filename)
+ {
+ if (_isValid)
+ {
+ _hmx.writeSerialized(filename);
+ }
+ else
+ {
+ // throw new Exception();
+ System.out.println("HMatrix not valid! Cannot write!");
+ }
+ }
+
+ /**
+ * Reads the HMatrix from a Serialized file
+ *
+ * @param filename The Serialized file from which to read.
+ * @return The HMatrix
+ */
+ public HMatrix readSerialized(String filename)
+ {
+ return HMatrix.readSerialized(filename);
+ }
+
+ /**
+ * Reads the HMatrix from an ASCII file
+ *
+ * @param filename The ASCII file from which to read.
+ * @return The HMatrix
+ */
+ public HMatrix read(String filename)
+ {
+ return HMatrix.read(filename);
+ }
+
+ /**
+ * Writes out the HMatrix to an ASCII file
+ *
+ * @param filename The file to which to write.
+ * @param comment A comment for the header.
+ */
+ public void write(String filename, String comment)
+ {
+ _hmx.write(filename, comment);
+ }
+
+ /**
+ * Invert the matrix
+ */
+ private double[][] invert(double[][] mat, int dim)
+ {
+
+ int i = dim;
+ double[][] matrix = new double[i][2 * i];
+ double[][] matrix1 = new double[i][1];
+ double[][] matrix2 = new double[i][1];
+ for(int j = 0; j < i; j++)
+ {
+ for(int i2 = 0; i2 < i; i2++)
+ matrix[j][i2]= mat[j][i2];
+
+ }
+
+ for(int k = 0; k < i; k++)
+ {
+ for(int j2 = i; j2 < 2 * i; j2++)
+ matrix[k][j2]= 0.0D;
+
+ matrix[k][i + k]=1.0D;
+ }
+
+ for(int i4 = 0; i4 < i; i4++)
+ {
+ for(int l = i4; l < i; l++)
+ {
+ matrix2[l][0]=0.0D;
+ for(int k2 = i4; k2 < i; k2++)
+ matrix2[l][0]=matrix2[l][0] + matrix[l][k2];
+
+ matrix1[l][0]= Math.abs(matrix[l][i4]) / matrix2[l][0];
+ }
+
+ int k4 = i4;
+ for(int i1 = i4 + 1; i1 < i; i1++)
+ if(matrix1[i1][0] > matrix1[i1 - 1][0])
+ k4 = i1;
+
+ if(matrix1[k4][0] == 0.0D)
+ {
+ System.err.println("Matrix::inverse -> Matrix is singular.");
+ System.exit(0);
+ }
+ if(k4 != i4)
+ {
+ for(int l2 = 0; l2 < 2 * i; l2++)
+ {
+ double d = matrix[i4][l2];
+ matrix[i4][l2]=matrix[k4][l2];
+ matrix[k4][l2]=d;
+ }
+
+ }
+ for(int i3 = 2 * i - 1; i3 >= i4; i3--)
+ matrix[i4][i3]=matrix[i4][i3] / matrix[i4][i4];
+
+ for(int j1 = i4 + 1; j1 < i; j1++)
+ {
+ for(int j3 = 2 * i - 1; j3 >= i4 + 1; j3--)
+ matrix[j1][j3] =matrix[j1][j3] - matrix[i4][j3] * matrix[j1][i4];
+
+ }
+
+ }
+
+ for(int j4 = i - 1; j4 >= 0; j4--)
+ {
+ for(int k1 = j4 - 1; k1 >= 0; k1--)
+ {
+ for(int k3 = i; k3 < 2 * i; k3++)
+ matrix[k1][k3] =matrix[k1][k3] - matrix[j4][k3] * matrix[k1][j4];
+
+ }
+
+ }
+
+ double[][] matrix3 = new double[i][i];
+ for(int l1 = 0; l1 < i; l1++)
+ {
+ for(int l3 = 0; l3 < i; l3++)
+ matrix3[l1][l3] =matrix[l1][l3 + i];
+
+ }
+
+ return matrix3;
+ }
+
+/*
+// Inversion in place requires larger matrix.
+ private double[][] invert( double[][] D, int n)
+ {
+ double alpha;
+ double beta;
+ int i;
+ int j;
+ int k;
+ int error;
+
+ error = 0;
+ int n2 = 2*n;
+
+ // initialize the reduction matrix
+ for( i = 1; i <= n; i++ )
+ {
+ for( j = 1; j <= n; j++ )
+ {
+ D[i][j+n] = 0.;
+ }
+ D[i][i+n] = 1.0;
+ }
+
+ // perform the reductions
+ for( i = 1; i <= n; i++ )
+ {
+ alpha = D[i][i];
+ if( alpha == 0.0 ) // error - singular matrix
+ {
+ error = 1;
+ break;
+ }
+ else
+ {
+ for( j = 1; j <= n2; j++ )
+ {
+ D[i][j] = D[i][j]/alpha;
+ }
+ for( k = 1; k <= n; k++ )
+ {
+ if( (k-i) != 0 )
+ {
+ beta = D[k][i];
+ for( j = 1; j <= n2; j++ )
+ {
+ D[k][j] = D[k][j] - beta*D[i][j];
+ }
+ }
+ }
+ }
+ }
+ return D;
+ }
+ */
+}
