Initial creation of vecmath sources in CVS repository

git-svn-id: https://svn.java.net/svn/vecmath~svn/trunk@5 dd45e54d-f42e-c781-df72-dca083a658b1

1 files changed, 2989 insertions, 0 deletions

diff --git a/src/javax/vecmath/GMatrix.java b/src/javax/vecmath/GMatrix.java
new file mode 100644
index 0000000..5ccdbb4
--- /dev/null
+++ b/src/javax/vecmath/GMatrix.java
@@ -0,0 +1,2989 @@
+/*
+ * $RCSfile$
+ *
+ * Copyright (c) 2004 Sun Microsystems, Inc. All rights reserved.
+ *
+ * Use is subject to license terms.
+ *
+ * $Revision$
+ * $Date$
+ * $State$
+ */
+
+package javax.vecmath;
+
+import java.lang.Math;
+
+/**
+ * A double precision, general, dynamically-resizable,
+ * two-dimensional matrix class.  Row and column numbering begins with
+ * zero.  The representation is row major.
+ */
+
+public class GMatrix implements java.io.Serializable, Cloneable {
+
+    // Compatible with 1.1
+    static final long serialVersionUID = 2777097312029690941L;
+    private static final boolean debug = false;
+
+    int nRow;
+    int nCol;
+
+    // double dereference is slow 
+    double[][] values;
+
+    private static final double EPS = 1.0E-10;
+
+    /**
+     * Constructs an nRow by NCol identity matrix. 
+     * Note that because row and column numbering begins with
+     * zero, nRow and nCol will be one larger than the maximum
+     * possible matrix index values.
+     * @param nRow  number of rows in this matrix.
+     * @param nCol  number of columns in this matrix.
+     */
+    public GMatrix(int nRow, int nCol)
+    {
+        values = new double[nRow][nCol];
+	this.nRow = nRow;
+	this.nCol = nCol;
+
+	int i, j;
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = 0.0;
+	    }
+	}
+
+	int l;
+	if (nRow < nCol)
+	    l = nRow;
+        else
+	    l = nCol;
+
+	for (i = 0; i < l; i++) {
+	    values[i][i] = 1.0;
+	}
+    }
+
+    /** 
+     * Constructs an nRow by nCol matrix initialized to the values 
+     * in the matrix array.  The array values are copied in one row at
+     * a time in row major fashion.  The array should be at least 
+     * nRow*nCol in length.
+     * Note that because row and column numbering begins with 
+     * zero, nRow and nCol will be one larger than the maximum
+     * possible matrix index values.
+     * @param nRow  number of rows in this matrix. 
+     * @param nCol  number of columns in this matrix. 
+     * @param matrix  a 1D array that specifies a matrix in row major fashion
+     */ 
+    public GMatrix(int nRow, int nCol, double[] matrix)
+    {
+        values = new double[nRow][nCol];
+	this.nRow = nRow;
+	this.nCol = nCol;
+
+	int i, j;
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = matrix[i*nCol+j];
+	    }
+	}
+    }
+
+    /** 
+     * Constructs a new GMatrix and copies the initial values
+     * from the parameter matrix.
+     * @param matrix  the source of the initial values of the new GMatrix
+     */ 
+    public GMatrix(GMatrix matrix)
+    {
+        nRow = matrix.nRow;
+	nCol = matrix.nCol;
+        values = new double[nRow][nCol];
+
+	int i, j;
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = matrix.values[i][j];
+	    }
+	}
+    }
+
+    /**
+     * Sets the value of this matrix to the result of multiplying itself
+     * with matrix m1 (this = this * m1).  
+     * @param m1 the other matrix
+     */  
+    public final void mul(GMatrix m1)
+    {
+	int i, j, k;
+
+	if (nCol != m1.nRow ||  nCol != m1.nCol)
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix0"));
+
+	double [][] tmp = new double[nRow][nCol];
+
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		tmp[i][j] = 0.0;
+		for (k = 0; k < nCol; k++) {
+		    tmp[i][j] += values[i][k]*m1.values[k][j];	
+		}
+	    }
+	}
+       
+	values = tmp;
+    }
+
+    /**
+     * Sets the value of this matrix to the result of multiplying
+     * the two argument matrices together (this = m1 * m2).
+     * @param m1 the first matrix
+     * @param m2 the second matrix
+     */
+    public final void mul(GMatrix m1, GMatrix m2)
+    {
+	int i, j, k;
+
+	if (m1.nCol != m2.nRow || nRow != m1.nRow || nCol != m2.nCol)
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix1"));
+
+	double[][] tmp = new double[nRow][nCol];
+
+	for (i = 0; i < m1.nRow; i++) {
+	    for (j = 0; j < m2.nCol; j++) {
+		tmp[i][j] = 0.0;
+		for (k = 0; k < m1.nCol; k++) {
+		    tmp[i][j] += m1.values[i][k]*m2.values[k][j];	
+		}
+	    }
+	}
+       
+	values = tmp;
+    }
+
+    /**
+     * Computes the outer product of the two vectors; multiplies the
+     * the first vector by the transpose of the second vector and places
+     * the matrix result into this matrix.  This matrix must be
+     * be as big or bigger than getSize(v1)xgetSize(v2).
+     * @param v1 the first vector, treated as a row vector
+     * @param v2 the second vector, treated as a column vector
+     */
+    public final void mul(GVector v1, GVector v2)
+    {
+	int i, j;
+
+	if (nRow < v1.getSize())
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix2"));
+
+	if (nCol < v2.getSize())
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix3"));
+
+	for (i = 0; i < v1.getSize(); i++ ) {
+	    for (j = 0; j < v2.getSize(); j++ ) {
+		values[i][j] = v1.values[i]*v2.values[j];
+	    }
+	}
+    }
+
+    /**
+     * Sets the value of this matrix to sum of itself and matrix m1.
+     * @param m1 the other matrix
+     */  
+    public final void add(GMatrix m1)
+    {
+	int i, j;
+
+	if (nRow != m1.nRow) 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix4"));
+
+	if (nCol != m1.nCol) 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix5"));
+
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = values[i][j] + m1.values[i][j];
+	    }
+	}
+    }
+
+    /**
+     * Sets the value of this matrix to the matrix sum of matrices m1 and m2.
+     * @param m1 the first matrix
+     * @param m2 the second matrix
+     */  
+    public final void add(GMatrix m1, GMatrix m2)
+    {
+	int i, j;
+
+	if (m2.nRow != m1.nRow) 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix6"));
+
+	if (m2.nCol != m1.nCol) 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix7"));
+
+	if (nCol != m1.nCol  || nRow != m1.nRow) 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix8"));
+
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = m1.values[i][j] + m2.values[i][j];
+	    }
+	}
+    }
+
+    /**
+     * Sets the value of this matrix to the matrix difference of itself
+     * and matrix m1 (this = this - m1).
+     * @param m1 the other matrix
+     */  
+    public final void sub(GMatrix m1) 
+    { 
+	int i, j; 
+	if (nRow != m1.nRow)  
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix9")); 
+ 
+	if (nCol != m1.nCol) 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix28")); 
+ 
+	for (i = 0; i < nRow; i++) { 
+	    for (j = 0; j < nCol; j++) { 
+		values[i][j] = values[i][j] - m1.values[i][j];
+	    } 
+	} 
+    } 
+
+    /**
+     * Sets the value of this matrix to the matrix difference
+     * of matrices m1 and m2 (this = m1 - m2).
+     * @param m1 the first matrix
+     * @param m2 the second matrix
+     */  
+    public final void sub(GMatrix m1, GMatrix m2) 
+    {
+	int i, j;
+	if (m2.nRow != m1.nRow) 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix10"));
+
+	if (m2.nCol != m1.nCol) 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix11"));
+
+	if (nRow !=  m1.nRow || nCol != m1.nCol) 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix12"));
+	  
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = m1.values[i][j] - m2.values[i][j];
+	    }
+	}
+    }
+
+    /** 
+     * Negates the value of this matrix: this = -this.
+     */
+    public final void negate()
+    {
+	int i, j;
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0;j < nCol; j++) {
+		values[i][j] = -values[i][j];
+	    }
+	}
+    }
+
+    /**
+     *  Sets the value of this matrix equal to the negation of
+     *  of the GMatrix parameter.
+     *  @param m1  The source matrix
+     */  
+    public final void negate(GMatrix m1)
+    {
+	int i, j;
+	if (nRow != m1.nRow  || nCol != m1.nCol)
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix13"));
+
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] =  -m1.values[i][j];
+	    }
+	}
+    }
+
+    /**
+     * Sets this GMatrix to the identity matrix.
+     */
+    public final void setIdentity()
+    {
+        int i, j;
+        for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = 0.0;
+	    }
+        }
+
+        int l;
+        if (nRow < nCol)
+	    l = nRow;
+        else
+	    l = nCol;
+
+        for (i = 0; i < l; i++) {
+	    values[i][i] = 1.0;
+        }
+    }
+
+    /**
+     * Sets all the values in this matrix to zero.
+     */
+    public final void setZero()
+    {
+	int i, j;
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = 0.0; 
+	    }
+	}   
+    }
+
+    /**
+     * Subtracts this matrix from the identity matrix and puts the values
+     * back into this (this = I - this).
+     */
+    public final void identityMinus()
+    {
+	int i, j;
+ 
+	for(i = 0; i < nRow; i++) {
+	    for(j = 0; j < nCol; j++) { 
+		values[i][j] = -values[i][j]; 
+	    } 
+	}    
+
+        int l; 
+        if( nRow < nCol) 
+	    l = nRow; 
+        else 
+	    l = nCol;
+ 
+        for(i = 0; i < l; i++) { 
+	    values[i][i] += 1.0; 
+        } 
+    }
+
+   
+    /**
+     * Inverts this matrix in place.
+     */
+    public final void invert()
+    {
+	invertGeneral(this);
+    }
+
+    /**
+     * Inverts matrix m1 and places the new values into this matrix.  Matrix
+     * m1 is not modified.
+     * @param m1   the matrix to be inverted 
+     */ 
+    public final void invert(GMatrix m1)
+    {
+	invertGeneral(m1);
+    }
+
+    /**
+     * Copies a sub-matrix derived from this matrix into the target matrix.
+     * The upper left of the sub-matrix is located at (rowSource, colSource);
+     * the lower right of the sub-matrix is located at 
+     * (lastRowSource,lastColSource).  The sub-matrix is copied into the
+     * the target matrix starting at (rowDest, colDest).
+     * @param rowSource   the top-most row of the sub-matrix 
+     * @param colSource   the left-most column of the sub-matrix 
+     * @param numRow   the number of rows in the sub-matrix
+     * @param numCol  the number of columns in the sub-matrix
+     * @param rowDest  the top-most row of the position of the copied
+     *                 sub-matrix within the target matrix
+     * @param colDest  the left-most column of the position of the copied
+     *                 sub-matrix within the target matrix
+     * @param target  the matrix into which the sub-matrix will be copied
+     */
+    public final void copySubMatrix(int rowSource, int colSource, 
+				    int numRow, int numCol, int rowDest,
+				    int colDest, GMatrix target) 
+    {
+        int i, j;
+
+	if (this != target) {
+	    for (i = 0; i < numRow; i++) {
+		for (j = 0; j < numCol; j++) {
+		    target.values[rowDest+i][colDest+j] =
+			values[rowSource+i][colSource+j];
+		}
+	    }
+	} else {
+	    double[][] tmp = new double[numRow][numCol];
+	    for (i = 0; i < numRow; i++) {
+		for (j = 0; j < numCol; j++) {
+		    tmp[i][j] = values[rowSource+i][colSource+j];
+		}
+	    }
+	    for (i = 0; i < numRow; i++) {
+		for (j = 0; j < numCol; j++) {
+		    target.values[rowDest+i][colDest+j] = tmp[i][j];
+		}
+	    }
+	}
+    }
+
+    /**
+     * Changes the size of this matrix dynamically.  If the size is increased
+     * no data values will be lost.  If the size is decreased, only those data
+     * values whose matrix positions were eliminated will be lost.
+     * @param nRow  number of desired rows in this matrix
+     * @param nCol  number of desired columns in this matrix
+     */
+    public final void setSize(int nRow, int nCol)
+    {
+	double[][] tmp = new double[nRow][nCol];
+	int i, j, maxRow, maxCol;
+
+	if (this.nRow < nRow)
+	    maxRow = this.nRow;
+	else
+	    maxRow = nRow;
+
+	if (this.nCol < nCol)
+	    maxCol = this.nCol;
+	else
+	    maxCol = nCol;
+
+	for (i = 0; i < maxRow; i++) {
+	    for (j = 0; j < maxCol; j++) {
+		tmp[i][j] = values[i][j];
+	    }
+	}
+
+	this.nRow = nRow;
+	this.nCol = nCol;
+
+	values = tmp;
+    }
+
+    /**
+     * Sets the value of this matrix to the values found in the array parameter.
+     * The values are copied in one row at a time, in row major 
+     * fashion.  The array should be at least equal in length to
+     * the number of matrix rows times the number of matrix columns
+     * in this matrix.
+     * @param matrix  the row major source array
+     */  
+    public final void set(double[] matrix)
+    {
+	int i, j;
+        
+	for (i = 0; i < nRow; i++) { 
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = matrix[nCol*i+j];
+	    } 
+	}  
+    }
+
+    /**
+     * Sets the value of this matrix to that of the Matrix3f provided.
+     * @param m1 the matrix
+     */  
+    public final void set(Matrix3f m1)
+    {
+	int i, j;
+
+	if (nCol < 3 || nRow < 3) { // expand matrix if too small 
+	    nCol = 3;
+	    nRow = 3;
+	    values = new double[nRow][nCol];
+        }
+
+        values[0][0] = m1.m00;
+        values[0][1] = m1.m01;
+        values[0][2] = m1.m02;
+
+        values[1][0] = m1.m10;
+        values[1][1] = m1.m11;
+        values[1][2] = m1.m12;
+
+        values[2][0] = m1.m20;
+        values[2][1] = m1.m21;
+        values[2][2] = m1.m22;
+
+        for (i = 3; i < nRow; i++) {   // pad rest or matrix with zeros 
+	    for (j = 3; j < nCol; j++) {
+		values[i][j] = 0.0;
+	    }
+        }
+    }
+
+    /**
+     * Sets the value of this matrix to that of the Matrix3d provided.
+     * @param m1 the matrix
+     */  
+    public final void set(Matrix3d m1)
+    {
+	if (nRow < 3 || nCol < 3) {
+	    values = new double[3][3];
+	    nRow = 3;
+	    nCol = 3;
+	}
+
+        values[0][0] = m1.m00;
+        values[0][1] = m1.m01;
+        values[0][2] = m1.m02;
+
+        values[1][0] = m1.m10;
+        values[1][1] = m1.m11;
+        values[1][2] = m1.m12;
+
+        values[2][0] = m1.m20;
+        values[2][1] = m1.m21;
+        values[2][2] = m1.m22;
+
+        for (int i = 3; i < nRow; i++) {   // pad rest or matrix with zeros 
+	    for(int j = 3; j < nCol; j++) {
+		values[i][j] = 0.0;
+	    }
+        }
+
+    }
+
+    /**  
+     * Sets the value of this matrix to that of the Matrix4f provided.
+     * @param m1 the matrix
+     */  
+    public final void set(Matrix4f m1)
+    {
+	if (nRow < 4 || nCol < 4) {
+	    values = new double[4][4];
+	    nRow = 4;
+	    nCol = 4;
+	}
+
+        values[0][0] = m1.m00;
+        values[0][1] = m1.m01;
+        values[0][2] = m1.m02;
+        values[0][3] = m1.m03;
+
+        values[1][0] = m1.m10;
+        values[1][1] = m1.m11;
+        values[1][2] = m1.m12;
+        values[1][3] = m1.m13;
+
+        values[2][0] = m1.m20;
+        values[2][1] = m1.m21;
+        values[2][2] = m1.m22;
+        values[2][3] = m1.m23;
+
+        values[3][0] = m1.m30;
+        values[3][1] = m1.m31;
+        values[3][2] = m1.m32;
+        values[3][3] = m1.m33;
+
+        for (int i = 4 ; i < nRow; i++) {   // pad rest or matrix with zeros 
+	    for (int j = 4; j < nCol; j++) {
+		values[i][j] = 0.0;
+	    }
+        }
+    }
+
+    /**   
+     * Sets the value of this matrix to that of the Matrix4d provided. 
+     * @param m1 the matrix 
+     */   
+    public final void set(Matrix4d m1) 
+    { 
+	if (nRow < 4 || nCol < 4) {
+	    values = new double[4][4];
+	    nRow = 4;
+	    nCol = 4;
+	}
+
+        values[0][0] = m1.m00;
+        values[0][1] = m1.m01;
+        values[0][2] = m1.m02;
+        values[0][3] = m1.m03;
+
+        values[1][0] = m1.m10;
+        values[1][1] = m1.m11;
+        values[1][2] = m1.m12;
+        values[1][3] = m1.m13;
+
+        values[2][0] = m1.m20;
+        values[2][1] = m1.m21;
+        values[2][2] = m1.m22;
+        values[2][3] = m1.m23;
+
+        values[3][0] = m1.m30;
+        values[3][1] = m1.m31;
+        values[3][2] = m1.m32;
+        values[3][3] = m1.m33;
+
+        for (int i = 4; i < nRow; i++) {   // pad rest or matrix with zeros 
+	    for (int j = 4; j < nCol; j++) {
+		values[i][j] = 0.0;
+	    }
+        }
+    } 
+
+    /**
+     * Sets the value of this matrix to the values found in matrix m1.
+     * @param m1  the source matrix
+     */
+    public final void set(GMatrix m1)
+    {
+	int i, j;
+
+	if (nRow < m1.nRow || nCol < m1.nCol) {
+	    nRow = m1.nRow;
+	    nCol = m1.nCol;
+	    values = new double[nRow][nCol];
+	}
+
+	for (i = 0; i < Math.min(nRow, m1.nRow); i++) {
+	    for (j = 0; j < Math.min(nCol, m1.nCol); j++) {
+		values[i][j] = m1.values[i][j];
+	    }
+	}
+
+        for (i = m1.nRow; i < nRow; i++) {   // pad rest or matrix with zeros 
+	    for (j = m1.nCol; j < nCol; j++) {
+		values[i][j] = 0.0;
+	    }
+        }
+    }
+
+    /**
+     * Returns the number of rows in this matrix.
+     * @return  number of rows in this matrix
+     */
+    public final int getNumRow()
+    {
+        return(nRow);
+    }
+
+    /**
+     * Returns the number of colmuns in this matrix.
+     * @return  number of columns in this matrix
+     */
+    public final int getNumCol()
+    {
+	return(nCol);
+    }
+
+    /**
+     * Retrieves the value at the specified row and column of this matrix.
+     * @param row the row number to be retrieved (zero indexed)
+     * @param column the column number to be retrieved (zero indexed)
+     * @return the value at the indexed element
+     */  
+    public final double getElement(int row, int column)
+    {
+        return(values[row][column]);
+    }
+
+ 
+    /**  
+     * Modifies the value at the specified row and column of this matrix.
+     * @param row  the row number to be modified (zero indexed) 
+     * @param column  the column number to be modified (zero indexed) 
+     * @param value  the new matrix element value
+     */   
+    public final void setElement(int row, int column, double value)
+    {
+	values[row][column] = value;
+    }
+
+    /**
+     * Places the values of the specified row into the array parameter.
+     * @param row  the target row number
+     * @param array  the array into which the row values will be placed
+     */
+    public final void getRow(int row, double[] array)
+    {
+	for (int i = 0; i < nCol; i++) {
+            array[i] = values[row][i];
+	}
+    }
+
+    /**
+     * Places the values of the specified row into the vector parameter.
+     * @param row  the target row number
+     * @param vector  the vector into which the row values will be placed
+     */  
+    public final void getRow(int row, GVector vector)
+    {
+	if (vector.getSize() < nCol)
+	    vector.setSize(nCol);
+        
+	for (int i = 0; i < nCol; i++) {
+            vector.values[i] = values[row][i];
+	}
+    }
+
+    /**
+     * Places the values of the specified column into the array parameter.
+     * @param col  the target column number
+     * @param array  the array into which the column values will be placed
+     */
+    public final void getColumn(int col, double[] array)
+    {
+	for (int i = 0; i < nRow; i++) {
+            array[i] = values[i][col];
+	}
+
+    }
+
+    /**
+     * Places the values of the specified column into the vector parameter.
+     * @param col  the target column number
+     * @param vector  the vector into which the column values will be placed
+     */  
+    public final void getColumn(int col, GVector vector)
+    {
+	if (vector.getSize() < nRow)
+	    vector.setSize(nRow);
+                  
+	for (int i = 0; i < nRow; i++) { 
+            vector.values[i] = values[i][col];   
+	} 
+    }
+
+    /**
+     * Places the values in the upper 3x3 of this GMatrix into
+     * the matrix m1.
+     * @param m1  The matrix that will hold the new values
+     */  
+    public final void get(Matrix3d m1)
+    {
+	if (nRow < 3 || nCol < 3) {
+	    m1.setZero();
+	    if (nCol > 0) {
+		if (nRow > 0){ 
+		    m1.m00 = values[0][0];
+		    if (nRow > 1){ 
+			m1.m10 = values[1][0];
+			if( nRow > 2 ){ 
+			    m1.m20= values[2][0];
+			}
+		    }
+		}
+		if (nCol > 1) {
+		    if (nRow > 0) { 
+			m1.m01 = values[0][1];
+			if (nRow > 1){ 
+			    m1.m11 = values[1][1];
+			    if (nRow >  2){ 
+				m1.m21 = values[2][1];
+			    }
+			}
+		    }
+		    if (nCol > 2) {
+			if (nRow > 0) { 
+			    m1.m02 = values[0][2];
+			    if (nRow > 1) {
+				m1.m12 = values[1][2];
+				if (nRow > 2) { 
+				    m1.m22 = values[2][2];
+				}
+			    }
+			}
+		    }
+		}
+	    }
+	} else {
+	    m1.m00 = values[0][0];
+	    m1.m01 = values[0][1];
+	    m1.m02 = values[0][2];
+
+	    m1.m10 = values[1][0];
+	    m1.m11 = values[1][1];
+	    m1.m12 = values[1][2];
+
+	    m1.m20 = values[2][0];
+	    m1.m21 = values[2][1];
+	    m1.m22 = values[2][2];
+	}
+    }
+
+    /**
+     * Places the values in the upper 3x3 of this GMatrix into 
+     * the matrix m1. 
+     * @param m1  The matrix that will hold the new values
+     */  
+    public final void get(Matrix3f m1) 
+    { 
+
+	if (nRow < 3 || nCol < 3) {
+	    m1.setZero();
+	    if (nCol > 0) {
+		if (nRow > 0) { 
+		    m1.m00 = (float)values[0][0];
+		    if (nRow > 1) { 
+			m1.m10 = (float)values[1][0];
+			if (nRow > 2) { 
+			    m1.m20 = (float)values[2][0];
+			}
+		    }
+		}
+		if (nCol > 1) {
+		    if (nRow > 0) { 
+			m1.m01 = (float)values[0][1];
+			if (nRow >  1){ 
+			    m1.m11 = (float)values[1][1];
+			    if (nRow >  2){ 
+				m1.m21 = (float)values[2][1];
+			    }
+			}
+		    }
+		    if (nCol > 2) {
+			if (nRow > 0) { 
+			    m1.m02 = (float)values[0][2];
+			    if (nRow > 1) {
+				m1.m12 = (float)values[1][2];
+				if (nRow > 2) { 
+				    m1.m22 = (float)values[2][2];
+				}
+			    }
+			}
+		    }
+		}
+	    }
+        } else {
+	    m1.m00 = (float)values[0][0];
+	    m1.m01 = (float)values[0][1];
+	    m1.m02 = (float)values[0][2];
+
+	    m1.m10 = (float)values[1][0];
+	    m1.m11 = (float)values[1][1];
+	    m1.m12 = (float)values[1][2];
+
+	    m1.m20 = (float)values[2][0];
+	    m1.m21 = (float)values[2][1];
+	    m1.m22 = (float)values[2][2];
+	}
+    } 
+
+    /**
+     * Places the values in the upper 4x4 of this GMatrix into 
+     * the matrix m1. 
+     * @param m1  The matrix that will hold the new values
+     */  
+    public final void get(Matrix4d m1) 
+    { 
+	if (nRow < 4 || nCol < 4) {
+	    m1.setZero();
+	    if (nCol > 0) {
+		if (nRow > 0) { 
+		    m1.m00 = values[0][0];
+		    if (nRow > 1) { 
+			m1.m10 = values[1][0];
+			if (nRow > 2) { 
+			    m1.m20 = values[2][0];
+			    if (nRow > 3) { 
+				m1.m30 = values[3][0];
+			    }
+			}
+		    }
+		}
+		if (nCol > 1) {
+		    if (nRow > 0) { 
+			m1.m01 = values[0][1];
+			if (nRow > 1) { 
+			    m1.m11 = values[1][1];
+			    if (nRow > 2) { 
+				m1.m21 = values[2][1];
+				if (nRow > 3) { 
+				    m1.m31 = values[3][1];
+				}
+			    }
+			}
+		    }
+		    if (nCol > 2) {
+			if (nRow > 0) { 
+			    m1.m02 = values[0][2];
+			    if (nRow > 1) {
+				m1.m12 = values[1][2];
+				if (nRow > 2) { 
+				    m1.m22 = values[2][2];
+				    if (nRow > 3) { 
+					m1.m32 = values[3][2];
+				    }
+				}
+			    }
+			}
+			if (nCol > 3) {
+			    if (nRow > 0) { 
+				m1.m03 = values[0][3];
+				if (nRow > 1) {
+				    m1.m13 = values[1][3];
+				    if (nRow > 2) { 
+					m1.m23 = values[2][3];
+					if (nRow > 3) { 
+					    m1.m33 = values[3][3];
+					}
+				    }
+				}
+			    }
+			}
+		    }
+		}
+	    }
+        } else {
+	    m1.m00 = values[0][0];
+	    m1.m01 = values[0][1];
+	    m1.m02 = values[0][2];
+	    m1.m03 = values[0][3];
+
+	    m1.m10 = values[1][0];
+	    m1.m11 = values[1][1];
+	    m1.m12 = values[1][2];
+	    m1.m13 = values[1][3];
+
+	    m1.m20 = values[2][0];
+	    m1.m21 = values[2][1];
+	    m1.m22 = values[2][2];
+	    m1.m23 = values[2][3];
+
+	    m1.m30 = values[3][0];
+	    m1.m31 = values[3][1];
+	    m1.m32 = values[3][2];
+	    m1.m33 = values[3][3];
+	}
+
+    } 
+
+    /**
+     * Places the values in the upper 4x4 of this GMatrix into 
+     * the matrix m1. 
+     * @param m1  The matrix that will hold the new values
+     */  
+    public final void get(Matrix4f m1) 
+    { 
+
+	if (nRow < 4 || nCol < 4) {
+	    m1.setZero();
+	    if (nCol > 0) {
+		if (nRow > 0) {
+		    m1.m00 = (float)values[0][0];
+		    if (nRow > 1) {
+			m1.m10 = (float)values[1][0];
+			if (nRow > 2) {
+			    m1.m20 = (float)values[2][0];
+			    if (nRow > 3) {
+				m1.m30 = (float)values[3][0];
+			    }  
+			}
+		    }
+		}   
+		if (nCol > 1) {
+		    if (nRow > 0) {
+			m1.m01 = (float)values[0][1];
+			if (nRow > 1) {
+			    m1.m11 = (float)values[1][1];
+			    if (nRow > 2) {
+				m1.m21 = (float)values[2][1];
+				if (nRow > 3) {
+				    m1.m31 = (float)values[3][1];
+				}
+			    }   
+			}   
+		    }
+		    if (nCol > 2) {
+			if (nRow > 0) {
+			    m1.m02 = (float)values[0][2];
+			    if (nRow > 1) {
+				m1.m12 = (float)values[1][2];
+				if (nRow > 2) {
+				    m1.m22 = (float)values[2][2];
+				    if (nRow > 3) {
+					m1.m32 = (float)values[3][2];
+				    }
+				}   
+			    }   
+			}   
+			if (nCol > 3) {
+			    if (nRow > 0) {
+				m1.m03 = (float)values[0][3];
+				if (nRow > 1) {
+				    m1.m13 = (float)values[1][3];
+				    if (nRow > 2) {
+					m1.m23 = (float)values[2][3];
+					if (nRow > 3) {
+					    m1.m33 = (float)values[3][3];
+					}
+				    }   
+				}   
+			    }   
+			}  
+		    }  
+		}   
+	    }   
+        } else {
+	    m1.m00 = (float)values[0][0];
+	    m1.m01 = (float)values[0][1];
+	    m1.m02 = (float)values[0][2];
+	    m1.m03 = (float)values[0][3];
+
+	    m1.m10 = (float)values[1][0];
+	    m1.m11 = (float)values[1][1];
+	    m1.m12 = (float)values[1][2];
+	    m1.m13 = (float)values[1][3];
+
+	    m1.m20 = (float)values[2][0];
+	    m1.m21 = (float)values[2][1];
+	    m1.m22 = (float)values[2][2];
+	    m1.m23 = (float)values[2][3];
+
+	    m1.m30 = (float)values[3][0];
+	    m1.m31 = (float)values[3][1];
+	    m1.m32 = (float)values[3][2];
+	    m1.m33 = (float)values[3][3];
+	}
+    } 
+
+    /**
+     * Places the values in the this GMatrix into the matrix m1;
+     * m1 should be at least as large as this GMatrix.
+     * @param m1  The matrix that will hold the new values
+     */  
+    public final void get(GMatrix m1) 
+    { 
+	int i, j, nc, nr;
+
+	if (nCol < m1.nCol)
+	    nc = nCol;
+	else
+	    nc = m1.nCol;
+
+	if (nRow < m1.nRow)
+	    nr = nRow;
+	else
+	    nr = m1.nRow;
+
+	for (i = 0; i < nr; i++) {
+	    for (j = 0; j < nc; j++) {
+		m1.values[i][j] = values[i][j];
+	    }
+	}
+	for (i = nr; i < m1.nRow; i++) {
+	    for (j = 0; j < m1.nCol; j++) {
+		m1.values[i][j] = 0.0;
+	    }
+	}
+	for (j = nc; j < m1.nCol; j++) {
+	    for (i = 0; i < nr; i++) {
+		m1.values[i][j] = 0.0;
+	    }
+	}
+    } 
+
+    /**
+     * Copy the values from the array into the specified row of this
+     * matrix.  
+     * @param row  the row of this matrix into which the array values 
+     *             will be copied.
+     * @param array  the source array
+     */
+    public final void setRow(int row, double[] array)
+    {
+	for (int i = 0; i < nCol; i++) {
+            values[row][i] = array[i];
+	}
+    }
+
+    /**
+     * Copy the values from the vector into the specified row of this
+     * matrix.
+     * @param row  the row of this matrix into which the array values
+     *             will be copied
+     * @param vector  the source vector
+     */  
+    public final void setRow(int row, GVector vector)
+    {
+	for(int i = 0; i < nCol; i++) {
+            values[row][i] = vector.values[i];
+	}
+    }
+
+    /**
+     * Copy the values from the array into the specified column of this
+     * matrix.   
+     * @param col  the column of this matrix into which the array values 
+     *             will be copied
+     * @param array  the source array
+     */ 
+    public final void setColumn(int col, double[] array)
+    {
+	for(int i = 0; i < nRow; i++) {
+            values[i][col] = array[i];
+	}
+    }
+
+    /**
+     * Copy the values from the vector into the specified column of this
+     * matrix.
+     * @param col  the column of this matrix into which the array values
+     *             will be copied
+     * @param vector  the source vector
+     */  
+    public final void setColumn(int col, GVector vector)
+    {
+	for(int i = 0; i < nRow; i++) {
+            values[i][col] = vector.values[i];
+	}
+
+    }
+
+    /**
+     *  Multiplies the transpose of matrix m1 times the transpose of matrix
+     *  m2, and places the result into this.
+     *  @param m1  The matrix on the left hand side of the multiplication
+     *  @param m2  The matrix on the right hand side of the multiplication
+     */  
+    public final void mulTransposeBoth(GMatrix m1, GMatrix m2)
+    {    
+	int i, j, k;
+
+	if (m1.nRow != m2.nCol || nRow != m1.nCol || nCol != m2.nRow)
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix14"));
+
+	if (m1 == this || m2 == this) {
+	    double[][] tmp = new double[nRow][nCol];
+	    for (i = 0; i < nRow; i++) { 
+		for (j = 0; j < nCol; j++) {
+		    tmp[i][j] = 0.0; 
+		    for (k = 0; k < m1.nRow; k++) {
+			tmp[i][j] += m1.values[k][i]*m2.values[j][k]; 
+		    }
+		}  
+	    }
+	    values = tmp;
+	} else {
+	    for (i = 0; i < nRow; i++) { 
+		for (j = 0; j < nCol; j++) {
+		    values[i][j] = 0.0; 
+		    for (k = 0; k < m1.nRow; k++) {
+			values[i][j] += m1.values[k][i]*m2.values[j][k]; 
+		    }
+		}  
+	    }  
+	}
+    }    
+
+    /**   
+     *  Multiplies matrix m1 times the transpose of matrix m2, and
+     *  places the result into this.
+     *  @param m1  The matrix on the left hand side of the multiplication
+     *  @param m2  The matrix on the right hand side of the multiplication
+     */  
+    public final void mulTransposeRight(GMatrix m1, GMatrix m2)
+    {    
+	int i, j, k;
+
+	if (m1.nCol != m2.nCol || nCol != m2.nRow || nRow != m1.nRow)
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix15"));
+
+	if (m1 == this || m2 == this) {
+	    double[][] tmp = new double[nRow][nCol];
+	    for (i = 0; i < nRow; i++) { 
+		for (j = 0; j < nCol; j++) {
+		    tmp[i][j] = 0.0; 
+		    for (k = 0; k < m1.nCol; k++) {
+			tmp[i][j] += m1.values[i][k]*m2.values[j][k]; 
+		    }
+		}  
+	    }
+	    values = tmp;
+	} else {
+	    for (i = 0; i < nRow; i++) { 
+		for (j = 0;j < nCol; j++) {
+		    values[i][j] = 0.0; 
+		    for (k = 0; k < m1.nCol; k++) {
+			values[i][j] += m1.values[i][k]*m2.values[j][k]; 
+		    }
+		}  
+	    }  
+	}
+
+    }    
+ 
+
+    /**   
+     *  Multiplies the transpose of matrix m1 times matrix m2, and
+     *  places the result into this.
+     *  @param m1  The matrix on the left hand side of the multiplication
+     *  @param m2  The matrix on the right hand side of the multiplication
+     */  
+    public final void mulTransposeLeft(GMatrix m1, GMatrix m2)
+    {
+	int i, j, k;
+
+	if (m1.nRow != m2.nRow || nCol != m2.nCol || nRow != m1.nCol)
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix16"));
+
+	if (m1 == this || m2 == this) {
+	    double[][] tmp = new double[nRow][nCol];
+	    for (i = 0; i < nRow; i++) { 
+		for (j = 0; j < nCol; j++) {
+		    tmp[i][j] = 0.0; 
+		    for (k = 0; k < m1.nRow; k++) {
+			tmp[i][j] += m1.values[k][i]*m2.values[k][j]; 
+		    }
+		}  
+	    }
+	    values = tmp;
+	} else {
+	    for (i = 0; i < nRow; i++) { 
+		for (j = 0; j < nCol; j++) {
+		    values[i][j] = 0.0; 
+		    for (k = 0; k < m1.nRow; k++) {
+			values[i][j] += m1.values[k][i]*m2.values[k][j]; 
+		    }
+		}  
+	    }  
+	}
+    }
+   
+
+    /**
+     * Transposes this matrix in place.
+     */
+    public final void transpose()
+    {
+        int i, j;
+
+        if (nRow != nCol) {
+	    double[][] tmp;
+	    i=nRow;
+	    nRow = nCol;
+	    nCol = i;
+	    tmp = new double[nRow][nCol];
+	    for (i = 0; i < nRow; i++) { 
+		for (j = 0; j < nCol; j++) {
+		    tmp[i][j] = values[j][i]; 
+		}  
+	    }  
+	    values = tmp;
+        } else {
+	    double swap;
+	    for (i = 0; i < nRow; i++) { 
+		for (j = 0; j < i; j++) {
+		    swap = values[i][j];
+		    values[i][j] = values[j][i];
+		    values[j][i] = swap;
+		}  
+	    }  
+	}
+    }
+
+    /**
+     * Places the matrix values of the transpose of matrix m1 into this matrix.
+     * @param m1  the matrix to be transposed (but not modified)
+     */
+    public final void transpose(GMatrix m1)
+    {
+        int i, j;
+
+        if (nRow != m1.nCol || nCol != m1.nRow)
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix17"));
+
+        if (m1 != this) {
+	    for (i = 0; i < nRow; i++) {
+		for (j = 0;j < nCol; j++) {
+		    values[i][j] = m1.values[j][i];
+		}
+	    }
+	} else {
+	    transpose();
+        }
+    }
+
+    /**
+     * Returns a string that contains the values of this GMatrix.
+     * @return the String representation
+     */  
+    public String toString() 
+    {
+	StringBuffer buffer = new StringBuffer(nRow*nCol*8);
+
+	int i, j;
+
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		buffer.append(values[i][j]).append(" ");
+	    }
+	    buffer.append("\n");
+	}
+
+	return buffer.toString();
+    }
+
+    private static void checkMatrix( GMatrix m) 
+    {
+	int i, j;
+
+	for (i = 0; i < m.nRow; i++) {
+	    for (j = 0; j < m.nCol; j++) {
+		if (Math.abs(m.values[i][j]) < 0.0000000001) {
+		    System.out.print(" 0.0     ");
+		} else {
+		    System.out.print(" " + m.values[i][j]);
+		}
+	    }
+	    System.out.print("\n");
+	}
+    }
+
+
+    /**
+     * Returns a hash code value based on the data values in this
+     * object.  Two different GMatrix objects with identical data
+     * values (i.e., GMatrix.equals returns true) will return the
+     * same hash number.  Two GMatrix objects with different data
+     * members may return the same hash value, although this is not
+     * likely.
+     * @return the integer hash code value
+     */
+    public int hashCode() {
+	long bits = 1L;
+
+	bits = 31L * bits + (long)nRow;
+	bits = 31L * bits + (long)nCol;
+
+	for (int i = 0; i < nRow; i++) {
+	    for (int j = 0; j < nCol; j++) {
+		bits = 31L * bits + Double.doubleToLongBits(values[i][j]);
+	    }
+	}
+
+	return (int) (bits ^ (bits >> 32));
+    } 
+
+
+    /**
+     * Returns true if all of the data members of GMatrix m1 are
+     * equal to the corresponding data members in this GMatrix.
+     * @param m1  The matrix with which the comparison is made.
+     * @return  true or false
+     */  
+    public boolean equals(GMatrix m1)
+    {
+	try { 
+	    int i, j;
+
+	    if (nRow != m1.nRow || nCol != m1.nCol)
+		return false;
+
+	    for (i = 0;i < nRow; i++) {
+		for (j = 0; j < nCol; j++) {
+		    if (values[i][j] != m1.values[i][j])
+			return false;
+		}
+	    }
+	    return true;
+	}  
+	catch (NullPointerException e2) {
+	    return false;
+	}
+    }
+
+    /**
+     * Returns true if the Object o1 is of type GMatrix and all of the
+     * data members of o1 are equal to the corresponding data members in
+     * this GMatrix.
+     * @param o1  The object with which the comparison is made.
+     * @return  true or false
+     */  
+    public boolean equals(Object o1)
+    {
+        try { 
+	    GMatrix m2 = (GMatrix) o1;
+	    int i, j;
+	    if (nRow != m2.nRow || nCol != m2.nCol)
+		return false;
+
+	    for (i = 0; i < nRow; i++) {
+                for (j = 0; j < nCol; j++) {
+		    if (values[i][j] != m2.values[i][j])
+			return false;
+                }
+	    }
+	    return true;
+        }
+        catch (ClassCastException e1) {
+	    return false;
+	}
+        catch (NullPointerException e2) {
+	    return false;
+	}
+    }
+
+    /**
+     * @deprecated Use epsilonEquals(GMatrix, double) instead
+     */
+    public boolean epsilonEquals(GMatrix m1, float epsilon) {
+	return epsilonEquals(m1, (double)epsilon);
+    }
+
+    /**
+     * Returns true if the L-infinite distance between this matrix
+     * and matrix m1 is less than or equal to the epsilon parameter,
+     * otherwise returns false.  The L-infinite
+     * distance is equal to
+     * MAX[i=0,1,2, . . .n ; j=0,1,2, . . .n ; abs(this.m(i,j) - m1.m(i,j)]
+     * @param m1  The matrix to be compared to this matrix
+     * @param epsilon  the threshold value
+     */  
+    public boolean epsilonEquals(GMatrix m1, double epsilon)
+    {
+	int i, j;
+	double diff;
+	if (nRow != m1.nRow || nCol != m1.nCol)
+	    return false;
+
+        for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		diff = values[i][j] - m1.values[i][j];
+		if ((diff < 0 ? -diff : diff) > epsilon)
+		    return false;
+	    }
+        }
+        return true;
+    }
+
+    /**
+     * Returns the trace of this matrix.
+     * @return  the trace of this matrix
+     */
+    public final double trace()
+    {
+	int i, l;
+	double t;
+
+	if (nRow < nCol) 
+	    l = nRow;
+	else 
+	    l = nCol;
+
+	t = 0.0;
+	for (i = 0; i < l; i++) {
+	    t += values[i][i];
+	}
+	return t;
+    }
+
+    /**
+     *  Finds the singular value decomposition (SVD) of this matrix 
+     *  such that this = U*W*transpose(V); and returns the rank of
+     *  this matrix; the values of U,W,V are all overwritten.  Note 
+     *  that the matrix V is output as V, and
+     *  not transpose(V).  If this matrix is mxn, then U is mxm, W 
+     *  is a diagonal matrix that is mxn, and V is nxn.  Using the
+     *  notation W = diag(w), then the inverse of this matrix is:
+     *  inverse(this) = V*diag(1/w)*tranpose(U), where diag(1/w)
+     *  is the same matrix as W except that the reciprocal of each 
+     *  of the diagonal components is used.
+     *  @param U  The computed U matrix in the equation this = U*W*transpose(V)
+     *  @param W  The computed W matrix in the equation this = U*W*transpose(V)
+     *  @param V  The computed V matrix in the equation this = U*W*transpose(V)
+     *  @return  The rank of this matrix.
+     */
+    public final int SVD(GMatrix U, GMatrix W, GMatrix V)
+    {
+	// check for consistancy in dimensions 
+	if (nCol != V.nCol || nCol != V.nRow) {
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix18"));
+	}
+
+	if (nRow != U.nRow || nRow != U.nCol) {
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix25"));
+	}
+
+	if (nRow != W.nRow || nCol != W.nCol) {
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix26"));
+	}
+
+	// Fix ArrayIndexOutOfBounds for 2x2 matrices, which partially
+	// addresses bug 4348562 for J3D 1.2.1.
+	//
+	// Does *not* fix the following problems reported in 4348562,
+	// which will wait for J3D 1.3:
+	//
+	//   1) no output of W
+	//   2) wrong transposition of U
+	//   3) wrong results for 4x4 matrices
+	//   4) slow performance
+	if (nRow == 2 && nCol == 2) {
+	    if (values[1][0] == 0.0) {
+		U.setIdentity();
+		V.setIdentity();
+
+		if (values[0][1] == 0.0) {
+		    return 2;
+		}
+
+		double[] sinl = new double[1];
+		double[] sinr = new double[1];
+		double[] cosl = new double[1];
+		double[] cosr = new double[1];
+		double[] single_values = new double[2];
+
+		single_values[0] = values[0][0];
+		single_values[1] = values[1][1];
+
+		compute_2X2(values[0][0], values[0][1], values[1][1],
+			    single_values, sinl, cosl, sinr, cosr, 0);
+
+		update_u(0, U, cosl, sinl);
+		update_v(0, V, cosr, sinr);
+
+		return 2;
+	    }
+	    // else call computeSVD() and check for 2x2 there
+	}
+
+	return computeSVD(this, U, W, V);
+    }
+
+    /**
+     * LU Decomposition: this matrix must be a square matrix and the
+     * LU GMatrix parameter must be the same size as this matrix. 
+     * The matrix LU will be overwritten as the combination of a 
+     * lower diagonal and upper diagonal matrix decompostion of this
+     * matrix; the diagonal 
+     * elements of L (unity) are not stored.  The GVector parameter
+     * records the row permutation effected by the partial pivoting,
+     * and is used as a parameter to the GVector method LUDBackSolve
+     * to solve sets of linear equations.
+     * This method returns +/- 1 depending on whether the number
+     * of row interchanges was even or odd, respectively.
+     * @param LU  The matrix into which the lower and upper decompositions
+     * will be placed.
+     * @param permutation  The row permutation effected by the partial
+     * pivoting
+     * @return  +-1 depending on whether the number of row interchanges
+     * was even or odd respectively
+     */
+    public final int LUD(GMatrix LU, GVector permutation)
+    {
+        int size = LU.nRow*LU.nCol;
+        double[] temp = new double[size];
+	int[] even_row_exchange = new int[1];
+	int[] row_perm = new int[LU.nRow];
+	int i, j;
+
+        if (nRow != nCol) {
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix19"));
+        }
+
+        if (nRow != LU.nRow) {
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix27"));
+        }
+
+        if (nCol != LU.nCol) {
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix27"));
+        }
+
+        if (LU.nRow != permutation.getSize()) {
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix20"));
+        }
+
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		temp[i*nCol+j] = values[i][j];
+	    }
+        }
+
+        // Calculate LU decomposition: Is the matrix singular? 
+        if (!luDecomposition(LU.nRow, temp, row_perm, even_row_exchange)) {
+            // Matrix has no inverse 
+            throw new SingularMatrixException
+		(VecMathI18N.getString("GMatrix21"));
+        }
+
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		LU.values[i][j] = temp[i*nCol+j];
+	    }
+        }
+
+	for (i = 0; i < LU.nRow; i++){
+	    permutation.values[i] = (double)row_perm[i];
+        }
+
+        return even_row_exchange[0];
+    }
+
+    /** 
+     *  Sets this matrix to a uniform scale matrix; all of the  
+     *  values are reset. 
+     *  @param scale  The new scale value 
+     */   
+    public final void setScale(double scale) 
+    { 
+	int i, j, l;
+
+	if (nRow < nCol) 
+	    l = nRow;
+	else 
+	    l = nCol;
+
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] = 0.0;
+	    }
+	}
+
+	for (i = 0; i < l; i++) {
+	    values[i][i] = scale;
+	}
+    } 
+
+    /**
+     * General invert routine.  Inverts m1 and places the result in "this".
+     * Note that this routine handles both the "this" version and the
+     * non-"this" version.
+     *
+     * Also note that since this routine is slow anyway, we won't worry
+     * about allocating a little bit of garbage.
+     */
+    final void invertGeneral(GMatrix  m1) {
+        int size = m1.nRow*m1.nCol;
+	double temp[] = new double[size];
+	double result[] = new double[size];
+	int row_perm[] = new int[m1.nRow];
+	int[] even_row_exchange = new int[1];
+	int i, j;
+
+	// Use LU decomposition and backsubstitution code specifically
+	// for floating-point nxn matrices.
+	if (m1.nRow != m1.nCol) {
+	    // Matrix is either under or over determined 
+	    throw new MismatchedSizeException
+		(VecMathI18N.getString("GMatrix22"));
+	}
+
+	// Copy source matrix to temp 
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		temp[i*nCol+j] = m1.values[i][j];
+	    }
+	}
+
+	// Calculate LU decomposition: Is the matrix singular? 
+	if (!luDecomposition(m1.nRow, temp, row_perm, even_row_exchange)) {
+	    // Matrix has no inverse 
+	    throw new SingularMatrixException
+		(VecMathI18N.getString("GMatrix21"));
+	}
+
+	// Perform back substitution on the identity matrix 
+        for (i = 0; i < size; i++)
+	    result[i] = 0.0;
+
+        for (i = 0; i < nCol; i++)
+	    result[i+i*nCol] = 1.0;
+
+	luBacksubstitution(m1.nRow, temp, row_perm, result);
+
+	for (i = 0; i < nRow; i++) {
+	    for (j = 0; j < nCol; j++) {
+		values[i][j] =  result[i*nCol+j];
+	    }
+        }
+    }
+
+    /**
+     * Given a nxn array "matrix0", this function replaces it with the 
+     * LU decomposition of a row-wise permutation of itself.  The input 
+     * parameters are "matrix0" and "dim".  The array "matrix0" is also 
+     * an output parameter.  The vector "row_perm[]" is an output 
+     * parameter that contains the row permutations resulting from partial 
+     * pivoting.  The output parameter "even_row_xchg" is 1 when the 
+     * number of row exchanges is even, or -1 otherwise.  Assumes data 
+     * type is always double.
+     *
+     * @return true if the matrix is nonsingular, or false otherwise.
+     */
+    //
+    // Reference: Press, Flannery, Teukolsky, Vetterling, 
+    //	      _Numerical_Recipes_in_C_, Cambridge University Press, 
+    //	      1988, pp 40-45.
+    //
+    static boolean luDecomposition(int dim, double[] matrix0,
+				   int[] row_perm, int[] even_row_xchg) {
+
+	double row_scale[] = new double[dim];
+
+	// Determine implicit scaling information by looping over rows 
+	int i, j;
+	int ptr, rs, mtx;
+	double big, temp;
+
+	ptr = 0;
+	rs = 0;
+	even_row_xchg[0] = 1;
+
+	// For each row ... 
+	i = dim;
+	while (i-- != 0) {
+	    big = 0.0;
+
+	    // For each column, find the largest element in the row 
+	    j = dim;
+	    while (j-- != 0) {
+		temp = matrix0[ptr++];
+		temp = Math.abs(temp);
+		if (temp > big) {
+		    big = temp;
+		}
+	    }
+
+	    // Is the matrix singular? 
+	    if (big == 0.0) {
+		return false;
+	    }
+	    row_scale[rs++] = 1.0 / big;
+	}
+
+	// For all columns, execute Crout's method 
+	mtx = 0;
+	for (j = 0; j < dim; j++) {
+	    int imax, k;
+	    int target, p1, p2;
+	    double sum;
+
+	    // Determine elements of upper diagonal matrix U 
+	    for (i = 0; i < j; i++) {
+		target = mtx + (dim*i) + j;
+		sum = matrix0[target];
+		k = i;
+		p1 = mtx + (dim*i);
+		p2 = mtx + j;
+		while (k-- != 0) {
+		    sum -= matrix0[p1] * matrix0[p2];
+		    p1++;
+		    p2 += dim;
+		}
+		matrix0[target] = sum;
+	    }
+
+	    // Search for largest pivot element and calculate
+	    // intermediate elements of lower diagonal matrix L.
+	    big = 0.0;
+	    imax = -1;
+	    for (i = j; i < dim; i++) {
+		target = mtx + (dim*i) + j;
+		sum = matrix0[target];
+		k = j;
+		p1 = mtx + (dim*i);
+		p2 = mtx + j;
+		while (k-- != 0) {
+		    sum -= matrix0[p1] * matrix0[p2];
+		    p1++;
+		    p2 += dim;
+		}
+		matrix0[target] = sum;
+
+		// Is this the best pivot so far? 
+		if ((temp = row_scale[i] * Math.abs(sum)) >= big) {
+		    big = temp;
+		    imax = i;
+		}
+	    }
+
+	    if (imax < 0) {
+		throw new RuntimeException(VecMathI18N.getString("GMatrix24"));
+	    }
+
+	    // Is a row exchange necessary? 
+	    if (j != imax) {
+		// Yes: exchange rows 
+		k = dim;
+		p1 = mtx + (dim*imax);
+		p2 = mtx + (dim*j);
+		while (k-- != 0) {
+		    temp = matrix0[p1];
+		    matrix0[p1++] = matrix0[p2];
+		    matrix0[p2++] = temp;
+		}
+
+		// Record change in scale factor 
+		row_scale[imax] = row_scale[j];
+		even_row_xchg[0] = -even_row_xchg[0]; // change exchange parity
+	    }
+
+	    // Record row permutation 
+	    row_perm[j] = imax;
+
+	    // Is the matrix singular 
+	    if (matrix0[(mtx + (dim*j) + j)] == 0.0) {
+		return false;
+	    }
+
+	    // Divide elements of lower diagonal matrix L by pivot 
+	    if (j != (dim-1)) {
+		temp = 1.0 / (matrix0[(mtx + (dim*j) + j)]);
+		target = mtx + (dim*(j+1)) + j;
+		i = (dim-1) - j;
+		while (i-- != 0) {
+		    matrix0[target] *= temp;
+		    target += dim;
+		}
+	    }
+
+	}
+
+	return true;
+    }
+
+    /**
+     * Solves a set of linear equations.  The input parameters "matrix1",
+     * and "row_perm" come from luDecompostion and do not change
+     * here.  The parameter "matrix2" is a set of column vectors assembled
+     * into a nxn matrix of floating-point values.  The procedure takes each
+     * column of "matrix2" in turn and treats it as the right-hand side of the
+     * matrix equation Ax = LUx = b.  The solution vector replaces the
+     * original column of the matrix.
+     *
+     * If "matrix2" is the identity matrix, the procedure replaces its contents
+     * with the inverse of the matrix from which "matrix1" was originally
+     * derived.
+     */
+    //
+    // Reference: Press, Flannery, Teukolsky, Vetterling, 
+    //	      _Numerical_Recipes_in_C_, Cambridge University Press, 
+    //	      1988, pp 44-45.
+    //
+    static void luBacksubstitution(int dim, double[] matrix1,
+				   int[] row_perm,
+				   double[] matrix2) {
+
+	int i, ii, ip, j, k;
+	int rp;
+	int cv, rv, ri;
+	double tt;
+	
+	// rp = row_perm;
+	rp = 0;
+
+	// For each column vector of matrix2 ... 
+	for (k = 0; k < dim; k++) {
+	    // cv = &(matrix2[0][k]);
+	    cv = k;
+	    ii = -1;
+
+	    // Forward substitution 
+	    for (i = 0; i < dim; i++) {
+		double sum;
+
+		ip = row_perm[rp+i];
+		sum = matrix2[cv+dim*ip];
+		matrix2[cv+dim*ip] = matrix2[cv+dim*i];
+		if (ii >= 0) {
+		    // rv = &(matrix1[i][0]);
+		    rv = i*dim;
+		    for (j = ii; j <= i-1; j++) {
+			sum -= matrix1[rv+j] * matrix2[cv+dim*j];
+		    }
+		}
+		else if (sum != 0.0) {
+		    ii = i;
+		}
+		matrix2[cv+dim*i] = sum;
+	    }
+
+	    // Backsubstitution 
+	    for (i = 0; i < dim; i++) {
+		ri = (dim-1-i);
+		rv = dim*(ri);
+		tt = 0.0;
+		for(j=1;j<=i;j++) {
+		    tt += matrix1[rv+dim-j] * matrix2[cv+dim*(dim-j)]; 	  
+		}
+		matrix2[cv+dim*ri]= (matrix2[cv+dim*ri] - tt) / matrix1[rv+ri];
+            }
+	}
+    }
+
+    static int computeSVD(GMatrix mat, GMatrix U, GMatrix W, GMatrix V) {
+	int i, j, k;
+	int nr, nc, si;
+
+	int converged, rank;
+	double cs, sn, r, mag,scale, t;
+	int eLength, sLength, vecLength;
+
+	GMatrix tmp = new GMatrix(mat.nRow, mat.nCol);
+	GMatrix u = new GMatrix(mat.nRow, mat.nCol);
+	GMatrix v = new GMatrix(mat.nRow, mat.nCol);
+	GMatrix m = new GMatrix(mat);
+
+	// compute the number of singular values
+	if (m.nRow >= m.nCol) {
+	    sLength = m.nCol;
+	    eLength = m.nCol-1;
+	}else {
+	    sLength = m.nRow;
+	    eLength = m.nRow;
+	}
+
+	if (m.nRow > m.nCol) 
+	    vecLength = m.nRow;
+	else
+	    vecLength = m.nCol;
+
+	double[] vec = new double[vecLength];
+	double[] single_values = new double[sLength];
+	double[] e = new double[eLength];
+
+	if(debug) {
+	    System.out.println("input to compute_svd = \n"+m.toString());
+	}
+
+	rank = 0;
+
+	U.setIdentity();
+	V.setIdentity();
+
+	nr = m.nRow;
+	nc = m.nCol;
+
+	// householder reduction 
+	for (si = 0; si < sLength; si++) {
+	    // for each singular value
+
+	    if (nr > 1) {
+		// zero out column
+		if (debug)
+		    System.out.println
+			("*********************** U ***********************\n");
+
+		// compute reflector
+		mag = 0.0;
+		for (i = 0; i < nr; i++) {
+		    mag += m.values[i+si][si] * m.values[i+si][si];
+		    if (debug)
+			System.out.println
+			    ("mag = " + mag + " matrix.dot = " +
+			     m.values[i+si][si] * m.values[i+si][si]);
+		}
+
+		mag = Math.sqrt(mag);
+		if (m.values[si][si] == 0.0) {
+		    vec[0] = mag;
+		} else {
+		    vec[0] = m.values[si][si] + d_sign(mag, m.values[si][si]);
+		}
+
+		for (i = 1; i < nr; i++) {
+		    vec[i] =  m.values[si+i][si];
+		}
+
+		scale = 0.0;
+		for (i = 0; i < nr; i++) {
+		    if (debug)
+			System.out.println("vec["+i+"]="+vec[i]);
+
+		    scale += vec[i]*vec[i];
+		}
+
+		scale = 2.0/scale;
+		if (debug)
+		    System.out.println("scale = "+scale);
+
+		for (j = si; j < m.nRow; j++) {
+		    for (k = si; k < m.nRow; k++) {
+			u.values[j][k] = -scale * vec[j-si]*vec[k-si];
+		    }
+		}
+
+		for (i = si; i < m.nRow; i++){
+		    u.values[i][i] +=  1.0;
+		}
+
+		// compute s
+		t = 0.0;
+		for (i = si; i < m.nRow; i++){
+		    t += u.values[si][i] * m.values[i][si];
+		}
+		m.values[si][si] = t;
+
+		// apply reflector
+		for (j = si; j < m.nRow; j++) {
+		    for (k = si+1; k < m.nCol; k++) {
+			tmp.values[j][k] = 0.0;
+			for (i = si; i < m.nCol; i++) {
+			    tmp.values[j][k] += u.values[j][i] * m.values[i][k];
+			}
+		    }
+		}
+
+		for (j = si; j < m.nRow; j++) {
+		    for (k = si+1; k < m.nCol; k++) {
+			m.values[j][k] = tmp.values[j][k];
+		    }
+		}
+
+		if (debug) {
+		    System.out.println("U =\n" + U.toString());
+		    System.out.println("u =\n" + u.toString());
+		}
+
+		// update U matrix
+		for (j = si; j < m.nRow; j++) {
+		    for (k = 0; k < m.nCol; k++) {
+			tmp.values[j][k] = 0.0; 
+			for (i = si; i < m.nCol; i++) {
+			    tmp.values[j][k] += u.values[j][i] * U.values[i][k];
+			}
+		    }
+		}
+
+		for (j = si; j < m.nRow; j++) {
+		    for (k = 0; k < m.nCol; k++) {
+			U.values[j][k] = tmp.values[j][k];
+		    }
+		}
+
+		if (debug) {
+		    System.out.println("single_values["+si+"] =\n" +
+				       single_values[si]);
+		    System.out.println("m =\n" + m.toString());
+		    System.out.println("U =\n" + U.toString());
+		}
+
+		nr--;
+	    }
+
+	    if( nc > 2 ) {
+		// zero out row
+		if (debug)
+		    System.out.println
+			("*********************** V ***********************\n");
+
+		mag = 0.0;
+		for (i = 1; i < nc; i++){
+		    mag += m.values[si][si+i] * m.values[si][si+i];
+		}
+
+		if (debug)
+		    System.out.println("mag = " + mag);
+
+		// generate the reflection vector, compute the first entry and
+		// copy the rest from the row to be zeroed
+		mag = Math.sqrt(mag);
+		if (m.values[si][si+1] == 0.0) {
+		    vec[0] = mag;
+		} else {
+		    vec[0] = m.values[si][si+1] +
+			d_sign(mag, m.values[si][si+1]);
+		}
+
+		for (i = 1; i < nc - 1; i++){
+		    vec[i] =  m.values[si][si+i+1];
+		}
+
+		// use reflection vector to compute v matrix
+		scale = 0.0;
+		for (i = 0; i < nc - 1; i++){
+		    if( debug )System.out.println("vec["+i+"]="+vec[i]);
+		    scale += vec[i]*vec[i];
+		}
+
+		scale = 2.0/scale;
+		if (debug)
+		    System.out.println("scale = "+scale);
+
+		for (j = si + 1; j < nc; j++) {
+		    for (k = si+1; k < m.nCol; k++) {
+			v.values[j][k] = -scale * vec[j-si-1]*vec[k-si-1];
+		    }
+		}
+
+		for (i = si + 1; i < m.nCol; i++){
+		    v.values[i][i] +=  1.0;
+		}
+
+		t=0.0;
+		for (i = si; i < m.nCol; i++){
+		    t += v.values[i][si+1] * m.values[si][i];
+		}
+		m.values[si][si+1]=t;
+
+		// apply reflector
+		for (j = si + 1; j < m.nRow; j++) {
+		    for (k = si + 1; k < m.nCol; k++) {
+			tmp.values[j][k] = 0.0;
+			for (i = si + 1; i < m.nCol; i++) {
+			    tmp.values[j][k] += v.values[i][k] * m.values[j][i];
+			}
+		    }
+		}
+
+		for (j = si + 1; j < m.nRow; j++) {
+		    for (k = si + 1; k < m.nCol; k++) {
+			m.values[j][k] = tmp.values[j][k];
+		    }
+		}
+
+		if (debug) {
+		    System.out.println("V =\n" + V.toString());
+		    System.out.println("v =\n" + v.toString());
+		    System.out.println("tmp =\n" + tmp.toString());
+		}
+
+		// update V matrix
+		for (j = 0; j < m.nRow; j++) {
+		    for (k = si + 1; k < m.nCol; k++) {
+			tmp.values[j][k] = 0.0; 
+			for (i = si + 1; i < m.nCol; i++) {
+			    tmp.values[j][k] += v.values[i][k] * V.values[j][i];
+			}
+		    }
+		}
+
+		if (debug)
+		    System.out.println("tmp =\n" + tmp.toString());
+
+		for (j = 0;j < m.nRow; j++) {
+		    for (k = si + 1; k < m.nCol; k++) {
+			V.values[j][k] = tmp.values[j][k];
+		    }
+		}
+
+		if (debug) {
+		    System.out.println("m =\n" + m.toString());
+		    System.out.println("V =\n" + V.toString());
+		}
+
+		nc--;
+	    }
+	}
+
+	for (i = 0; i < sLength; i++){
+	    single_values[i] = m.values[i][i];
+	}
+
+	for (i = 0; i < eLength; i++){
+	    e[i] = m.values[i][i+1];
+	}
+
+	// Fix ArrayIndexOutOfBounds for 2x2 matrices, which partially
+	// addresses bug 4348562 for J3D 1.2.1.
+	//
+	// Does *not* fix the following problems reported in 4348562,
+	// which will wait for J3D 1.3:
+	//
+	//   1) no output of W
+	//   2) wrong transposition of U
+	//   3) wrong results for 4x4 matrices
+	//   4) slow performance
+	if (m.nRow == 2 && m.nCol == 2) {
+	    double[] cosl = new double[1];
+	    double[] cosr = new double[1];
+	    double[] sinl = new double[1];
+	    double[] sinr = new double[1];
+
+	    compute_2X2(single_values[0], e[0], single_values[1],
+			single_values, sinl, cosl, sinr, cosr, 0);
+
+	    update_u(0, U, cosl, sinl);
+	    update_v(0, V, cosr, sinr);
+
+	    return 2;
+	}
+
+	// compute_qr causes ArrayIndexOutOfBounds for 2x2 matrices
+	compute_qr (0, e.length-1, single_values, e, U, V); 
+
+	// compute rank = number of non zero singular values
+	rank = single_values.length;
+  
+	// sort by order of size of single values
+	// and check for zero's
+	return rank;
+    }
+
+    static void compute_qr(int start, int end, double[] s, double[] e,
+			   GMatrix u, GMatrix v) {
+
+	int i, j, k, n, sl;
+	boolean converged;
+	double shift, r, utemp, vtemp, f, g;
+	double[] cosl = new double[1];
+	double[] cosr = new double[1];
+	double[] sinl = new double[1];
+	double[] sinr = new double[1];
+	GMatrix m = new GMatrix(u.nCol, v.nRow);
+ 
+	final int MAX_INTERATIONS = 2;
+	final double CONVERGE_TOL = 4.89E-15;
+
+	if (debug) {
+	    System.out.println("start =" + start);
+	    System.out.println("s =\n");
+	    for(i=0;i<s.length;i++) {
+		System.out.println(s[i]);
+	    }
+
+	    System.out.println("\nes =\n");
+	    for (i = 0; i < e.length; i++) {
+		System.out.println(e[i]);
+	    }
+
+	    for (i = 0; i < s.length; i++) {
+		m.values[i][i] = s[i];
+	    }
+
+	    for (i = 0; i < e.length; i++) {
+		m.values[i][i+1] = e[i];
+	    }
+	    System.out.println("\nm =\n" + m.toString());
+	}
+
+	double c_b48 =  1.0;
+	double c_b71 = -1.0;
+	converged = false;
+
+	if (debug)
+	    print_svd(s, e, u, v);
+
+	f = 0.0;
+	g = 0.0;
+   
+	for (k = 0; k < MAX_INTERATIONS && !converged;k++) {
+	    for (i = start; i <= end; i++) {
+
+		// if at start of iterfaction compute shift 
+		if (i == start) { 
+		    if (e.length == s.length)
+			sl = end;
+		    else
+			sl = end + 1;
+
+		    shift = compute_shift(s[sl-1], e[end], s[sl]);
+
+		    f = (Math.abs(s[i]) - shift) *
+			(d_sign(c_b48, s[i]) + shift/s[i]);
+		    g = e[i];
+		} 
+
+		r = compute_rot(f, g, sinr, cosr);
+		if (i != start)
+		    e[i-1] = r;
+
+		f = cosr[0] * s[i] + sinr[0] * e[i];
+		e[i] = cosr[0] * e[i] - sinr[0] * s[i];
+		g = sinr[0] * s[i+1];
+		s[i+1] = cosr[0] * s[i+1];
+
+		// if (debug) print_se(s,e);
+		update_v (i, v, cosr, sinr);
+		if (debug)
+		    print_m(m,u,v);
+
+		r = compute_rot(f, g, sinl, cosl);
+		s[i] = r;
+		f = cosl[0] * e[i] + sinl[0] * s[i+1];
+		s[i+1] = cosl[0] * s[i+1] - sinl[0] * e[i];
+
+		if( i < end) {
+		    // if not last
+		    g = sinl[0] * e[i+1];
+		    e[i+1] =  cosl[0] * e[i+1];
+		}
+		//if (debug) print_se(s,e);
+
+		update_u(i, u, cosl, sinl);
+		if (debug)
+		    print_m(m,u,v);
+	    }
+
+	    // if extra off diagonal perform one more right side rotation
+	    if (s.length == e.length) {
+		r = compute_rot(f, g, sinr, cosr);
+		f = cosr[0] * s[i] + sinr[0] * e[i];
+		e[i] = cosr[0] * e[i] - sinr[0] * s[i];
+		s[i+1] = cosr[0] * s[i+1];
+
+		update_v(i, v, cosr, sinr);
+		if (debug)
+		    print_m(m,u,v);
+	    } 
+
+	    if (debug) {
+		System.out.println
+		    ("\n*********************** iteration #" + k +
+		     " ***********************\n");
+		print_svd(s, e, u, v);
+	    }
+
+	    // check for convergence on off diagonals and reduce 
+	    while ((end-start > 1) && (Math.abs(e[end]) < CONVERGE_TOL)) { 
+		end--;
+	    }
+
+	    // check if need to split
+	    for (n = end - 2; n > start; n--) {
+		if (Math.abs(e[n]) < CONVERGE_TOL) {     // split
+		    compute_qr(n + 1, end, s, e, u, v);  // do lower matrix
+		    end = n - 1;                         // do upper matrix
+
+		    // check for convergence on off diagonals and reduce 
+		    while ((end - start > 1) &&
+			   (Math.abs(e[end]) < CONVERGE_TOL)) { 
+			end--;
+		    }
+		}
+	    }
+
+	    if (debug)
+		System.out.println("start = " + start);
+
+	    if ((end - start <= 1) && (Math.abs(e[start+1]) < CONVERGE_TOL)) {
+		converged = true;
+	    } else {
+		// check if zero on the diagonal
+	    }
+
+	}
+
+	if (debug)
+	    System.out.println("\n****call compute_2X2 ********************\n");
+
+	if (Math.abs(e[1]) < CONVERGE_TOL) {
+	    compute_2X2(s[start], e[start], s[start+1], s,
+			sinl, cosl, sinr, cosr, 0);
+	    e[start] = 0.0;
+	    e[start+1] = 0.0;
+	} else {
+	}
+
+	i = start;
+	update_u(i, u, cosl, sinl);
+	update_v(i, v, cosr, sinr);
+
+	if(debug) {
+	    System.out.println
+		("\n*******after call compute_2X2 **********************\n");
+	    print_svd(s, e, u, v);
+	}
+
+	return;
+    }
+
+    private static void print_se(double[] s, double[] e) {
+	System.out.println("\ns =" + s[0] + " " + s[1] + " " + s[2]);
+	System.out.println("e =" + e[0] + " " + e[1]);
+    }
+
+    private static void update_v(int index, GMatrix v,
+				 double[] cosr, double[] sinr) {
+	int j;
+	double vtemp;
+
+	for (j = 0; j < v.nRow; j++) {
+	    vtemp = v.values[j][index];
+	    v.values[j][index] =
+		cosr[0]*vtemp + sinr[0]*v.values[j][index+1];
+	    v.values[j][index+1] =
+	       -sinr[0]*vtemp + cosr[0]*v.values[j][index+1];
+	}
+    }
+
+    private static void chase_up(double[] s, double[] e, int k, GMatrix v) {
+	double f, g, r;
+	double[] cosr = new double[1];
+	double[] sinr = new double[1];
+	int i;
+	GMatrix t = new GMatrix(v.nRow, v.nCol);
+	GMatrix m = new GMatrix(v.nRow, v.nCol);
+
+	if (debug) {
+	    m.setIdentity();
+	    for (i = 0; i < s.length; i++) {
+		m.values[i][i] = s[i];
+	    }
+	    for (i = 0; i < e.length; i++) {
+		m.values[i][i+1] = e[i];
+	    }
+	}
+
+	f = e[k];
+	g = s[k];
+
+	for (i = k; i > 0; i--) {
+	    r = compute_rot(f, g, sinr, cosr);
+	    f = -e[i-1] * sinr[0];
+	    g = s[i-1];
+	    s[i] = r;
+	    e[i-1] = e[i-1] * cosr[0];
+	    update_v_split(i, k+1, v, cosr, sinr, t, m);
+	}
+
+	s[i+1] = compute_rot(f, g, sinr, cosr);
+	update_v_split(i, k+1, v, cosr, sinr, t, m);
+    }
+    
+    private static void chase_across(double[] s, double[] e, int k, GMatrix u) {
+	double f, g, r;
+	double[] cosl = new double[1];
+	double[] sinl = new double[1];
+	int i;
+	GMatrix t = new GMatrix(u.nRow, u.nCol);
+	GMatrix m = new GMatrix(u.nRow, u.nCol);
+
+	if (debug) {
+	    m.setIdentity();
+	    for (i = 0; i < s.length; i++) {
+		m.values[i][i] = s[i];
+	    }
+	    for (i = 0; i < e.length; i++) {
+		m.values[i][i+1] = e[i];
+	    }
+	}
+
+	g = e[k];
+	f = s[k+1];
+
+	for (i = k; i < u.nCol-2; i++){
+	    r = compute_rot(f, g, sinl, cosl);
+	    g = -e[i+1] * sinl[0];
+	    f = s[i+2];
+	    s[i+1] = r;
+	    e[i+1] = e[i+1] * cosl[0];
+	    update_u_split(k, i + 1, u, cosl, sinl, t, m);
+	}
+
+	s[i+1] = compute_rot(f, g, sinl, cosl);
+	update_u_split(k, i + 1, u, cosl, sinl, t, m);
+    }
+
+    private static void update_v_split(int topr, int bottomr, GMatrix v,
+				       double[] cosr, double[] sinr,
+				       GMatrix t, GMatrix m) {
+	int j;
+	double vtemp;
+
+	for (j = 0; j < v.nRow; j++) {
+	    vtemp = v.values[j][topr];
+	    v.values[j][topr] = cosr[0]*vtemp - sinr[0]*v.values[j][bottomr];
+	    v.values[j][bottomr] = sinr[0]*vtemp + cosr[0]*v.values[j][bottomr];
+	}
+
+	if (debug) {
+	    t.setIdentity();
+	    for (j = 0; j < v.nRow; j++) {
+		vtemp = t.values[j][topr];
+		t.values[j][topr] =
+		    cosr[0]*vtemp - sinr[0]*t.values[j][bottomr];
+		t.values[j][bottomr] =
+		    sinr[0]*vtemp + cosr[0]*t.values[j][bottomr];
+	    }
+	}
+
+	System.out.println("topr    =" + topr);
+	System.out.println("bottomr =" + bottomr);
+	System.out.println("cosr =" + cosr[0]);
+	System.out.println("sinr =" + sinr[0]);
+	System.out.println("\nm =");
+	checkMatrix(m);
+	System.out.println("\nv =");
+	checkMatrix(t);
+	m.mul(m,t);
+	System.out.println("\nt*m =");
+	checkMatrix(m);
+    }
+
+    private static void update_u_split(int topr, int bottomr, GMatrix u,
+				       double[] cosl, double[] sinl,
+				       GMatrix t, GMatrix m) {
+	int j;
+	double utemp;
+
+	for (j = 0; j < u.nCol; j++) {
+	    utemp = u.values[topr][j];
+	    u.values[topr][j]    = cosl[0]*utemp - sinl[0]*u.values[bottomr][j];
+	    u.values[bottomr][j] = sinl[0]*utemp + cosl[0]*u.values[bottomr][j];
+	}
+
+	if(debug) {
+	    t.setIdentity();
+	    for (j = 0;j < u.nCol; j++) {
+		utemp = t.values[topr][j];
+		t.values[topr][j] =
+		    cosl[0]*utemp - sinl[0]*t.values[bottomr][j];
+		t.values[bottomr][j] =
+		    sinl[0]*utemp + cosl[0]*t.values[bottomr][j];
+	    }
+	}
+	System.out.println("\nm=");
+	checkMatrix(m);
+	System.out.println("\nu=");
+	checkMatrix(t);
+	m.mul(t,m);
+	System.out.println("\nt*m=");
+	checkMatrix(m);
+    }
+
+    private static void update_u(int index, GMatrix u,
+				 double[] cosl, double[] sinl) {
+	int j;
+	double utemp;
+
+	for (j = 0; j < u.nCol; j++) {
+	    utemp = u.values[index][j];
+	    u.values[index][j] =
+		cosl[0]*utemp + sinl[0]*u.values[index+1][j];
+	    u.values[index+1][j] =
+	       -sinl[0]*utemp + cosl[0]*u.values[index+1][j];
+	}
+    }
+
+    private static void print_m(GMatrix m, GMatrix u, GMatrix v) {
+	GMatrix mtmp = new GMatrix(m.nCol, m.nRow);
+
+	mtmp.mul(u, mtmp);
+	mtmp.mul(mtmp, v);
+	System.out.println("\n m = \n" + mtmp.toString(mtmp));
+            
+    }
+
+    private static String toString(GMatrix m) 
+    {
+	StringBuffer buffer = new StringBuffer(m.nRow * m.nCol * 8);
+	int i, j;
+
+	for (i = 0; i < m.nRow; i++) {
+	    for(j = 0; j < m.nCol; j++) {
+		if (Math.abs(m.values[i][j]) < .000000001) {
+		    buffer.append("0.0000 ");
+		} else {
+		    buffer.append(m.values[i][j]).append(" ");
+		}
+	    }
+	    buffer.append("\n");
+	}
+	return buffer.toString();
+    }
+
+    private static void print_svd(double[] s, double[] e,
+				  GMatrix u, GMatrix v) {
+	int i;
+	GMatrix mtmp = new GMatrix(u.nCol, v.nRow);
+
+	System.out.println(" \ns = ");
+	for (i = 0; i < s.length; i++) {
+	    System.out.println(" " + s[i]);	
+	}
+
+	System.out.println(" \ne = ");
+	for (i = 0; i < e.length; i++) {
+	    System.out.println(" " + e[i]);	
+	}
+
+	System.out.println(" \nu  = \n" + u.toString());
+	System.out.println(" \nv  = \n" + v.toString());
+
+	mtmp.setIdentity();
+	for (i = 0; i < s.length; i++) {
+	    mtmp.values[i][i] = s[i];
+	}
+	for (i = 0; i < e.length; i++) {
+	    mtmp.values[i][i+1] = e[i];
+	}
+	System.out.println(" \nm  = \n"+mtmp.toString());
+
+	mtmp.mulTransposeLeft(u, mtmp);
+	mtmp.mulTransposeRight(mtmp, v);
+
+	System.out.println(" \n u.transpose*m*v.transpose  = \n" +
+			   mtmp.toString());
+    }
+
+    static double max(double a, double b) {
+	if (a > b) 
+	    return a;
+	else
+	    return b;
+    }
+
+    static double min(double a, double b) {
+	if (a < b) 
+	    return a;
+	else
+	    return b;
+    }
+
+    static double compute_shift(double f, double g, double h) {
+	double d__1, d__2;
+	double fhmn, fhmx, c, fa, ga, ha, as, at, au;
+	double ssmin;
+ 
+	fa = Math.abs(f);
+	ga = Math.abs(g);
+	ha = Math.abs(h);
+	fhmn = min(fa,ha);
+	fhmx = max(fa,ha);
+
+	if (fhmn == 0.0) {
+	    ssmin = 0.0;
+	    if (fhmx == 0.0) {
+	    } else {
+		d__1 = min(fhmx,ga) / max(fhmx,ga);
+	    }
+	} else {
+	    if (ga < fhmx) {
+		as = fhmn / fhmx + 1.0;
+		at = (fhmx - fhmn) / fhmx;
+		d__1 = ga / fhmx;
+		au = d__1 * d__1;
+		c = 2.0 / (Math.sqrt(as * as + au) + Math.sqrt(at * at + au));
+		ssmin = fhmn * c;
+	    } else {
+		au = fhmx / ga;
+		if (au == 0.0) {
+		    ssmin = fhmn * fhmx / ga;
+		} else {
+		    as = fhmn / fhmx + 1.0;
+		    at = (fhmx - fhmn) / fhmx;
+		    d__1 = as * au;
+		    d__2 = at * au;
+		    c = 1.0 / (Math.sqrt(d__1 * d__1 + 1.0) +
+			       Math.sqrt(d__2 * d__2 + 1.0));
+		    ssmin = fhmn * c * au;
+		    ssmin += ssmin;
+		}
+	    }
+	}    
+
+	return ssmin;
+    }
+
+    static int compute_2X2(double f, double g, double h,
+			   double[] single_values, double[] snl, double[] csl,
+			   double[] snr, double[] csr, int index) {
+
+	double c_b3 = 2.0;
+	double c_b4 = 1.0;
+ 
+	double d__1;
+	int pmax;
+	double temp;
+	boolean swap;
+	double a, d, l, m, r, s, t, tsign, fa, ga, ha;
+	double ft, gt, ht, mm;
+	boolean gasmal;
+	double tt, clt, crt, slt, srt;
+	double ssmin,ssmax;
+ 
+	ssmax = single_values[0];
+	ssmin = single_values[1];
+	clt = 0.0;
+	crt = 0.0;
+	slt = 0.0;
+	srt = 0.0;
+	tsign = 0.0;
+
+	ft = f;
+	fa = Math.abs(ft);
+	ht = h;
+	ha = Math.abs(h);
+ 
+	pmax = 1;
+	if (ha > fa) 
+	    swap = true;
+	else 
+	    swap = false;
+
+	if (swap) {
+	    pmax = 3;
+	    temp = ft;
+	    ft = ht;
+	    ht = temp;
+	    temp = fa;
+	    fa = ha;
+	    ha = temp;
+ 
+	}
+
+	gt = g;
+	ga = Math.abs(gt);
+	if (ga == 0.0) {
+	    single_values[1] = ha;
+	    single_values[0] = fa;
+	    clt = 1.0;
+	    crt = 1.0;
+	    slt = 0.0;
+	    srt = 0.0;
+	} else {
+	    gasmal = true;
+	    if (ga > fa) {
+		pmax = 2;
+		if (fa / ga < EPS) {
+		    gasmal = false;
+		    ssmax = ga;
+
+		    if (ha > 1.0) {
+			ssmin = fa / (ga / ha);
+		    } else {
+			ssmin = fa / ga * ha;
+		    }
+		    clt = 1.0;
+		    slt = ht / gt;
+		    srt = 1.0;
+		    crt = ft / gt;
+		}
+	    }
+	    if (gasmal) {
+		d = fa - ha;
+		if (d == fa) {
+ 
+		    l = 1.0;
+		} else {
+		    l = d / fa;
+		}
+ 
+		m = gt / ft;
+ 		t = 2.0 - l;
+ 		mm = m * m;
+		tt = t * t;
+		s = Math.sqrt(tt + mm);
+ 
+		if (l == 0.0) {
+		    r = Math.abs(m);
+		} else {
+		    r = Math.sqrt(l * l + mm);
+		}
+ 
+		a = (s + r) * 0.5;
+		if (ga > fa) {
+		    pmax = 2;
+		    if (fa / ga < EPS) {
+			gasmal = false;
+			ssmax = ga;
+			if (ha > 1.0) {
+			    ssmin = fa / (ga / ha);
+			} else {
+			    ssmin = fa / ga * ha;
+			}
+			clt = 1.0;
+			slt = ht / gt;
+			srt = 1.0;
+			crt = ft / gt;
+		    }
+		}
+		if (gasmal) {
+		    d = fa - ha;
+		    if (d == fa) {
+			l = 1.0;
+		    } else {
+			l = d / fa;
+		    }
+ 
+		    m = gt / ft;
+		    t = 2.0 - l;
+ 
+		    mm = m * m;
+		    tt = t * t;
+		    s = Math.sqrt(tt + mm);
+ 
+		    if (l == 0.) {
+			r = Math.abs(m);
+		    } else {
+			r = Math.sqrt(l * l + mm);
+		    }
+ 
+		    a = (s + r) * 0.5;
+		    ssmin = ha / a;
+		    ssmax = fa * a;
+
+		    if (mm == 0.0) {
+			if (l == 0.0) {
+			    t = d_sign(c_b3, ft) * d_sign(c_b4, gt);
+			} else {
+			    t = gt / d_sign(d, ft) + m / t;
+			}
+		    } else {
+			t = (m / (s + t) + m / (r + l)) * (a + 1.0);
+		    }
+
+		    l = Math.sqrt(t * t + 4.0);
+		    crt = 2.0 / l;
+		    srt = t / l;
+		    clt = (crt + srt * m) / a;
+		    slt = ht / ft * srt / a;
+		}
+	    }
+	    if (swap) {
+		csl[0] = srt;
+		snl[0] = crt;
+		csr[0] = slt;
+		snr[0] = clt;
+	    } else {
+		csl[0] = clt;
+		snl[0] = slt;
+		csr[0] = crt;
+		snr[0] = srt;
+	    }
+       
+	    if (pmax == 1) {
+		tsign = d_sign(c_b4, csr[0]) *
+		    d_sign(c_b4, csl[0]) * d_sign(c_b4, f);
+	    }
+	    if (pmax == 2) {
+		tsign = d_sign(c_b4, snr[0]) *
+		    d_sign(c_b4, csl[0]) * d_sign(c_b4, g);
+	    }
+	    if (pmax == 3) {
+		tsign = d_sign(c_b4, snr[0]) *
+		    d_sign(c_b4, snl[0]) * d_sign(c_b4, h);
+	    }
+
+	    single_values[index] = d_sign(ssmax, tsign);
+	    d__1 = tsign * d_sign(c_b4, f) * d_sign(c_b4, h);
+	    single_values[index+1] = d_sign(ssmin, d__1);
+	}  
+
+	return 0;
+    }
+
+    static double compute_rot(double f, double g, double[] sin, double[] cos) {
+	int i__1;
+	double d__1, d__2;
+	double cs, sn;
+	int i;
+	double scale;
+	int count;
+	double f1, g1;
+	double r;
+	final double safmn2 = 2.002083095183101E-146;
+	final double safmx2 = 4.994797680505588E+145;
+ 
+	if (g == 0.0) {
+	    cs = 1.0;
+	    sn = 0.0;
+	    r = f;
+	} else if (f == 0.0) {
+	    cs = 0.0;
+	    sn = 1.0;
+	    r = g;
+	} else {
+	    f1 = f;
+	    g1 = g;
+	    scale = max(Math.abs(f1),Math.abs(g1));
+	    if (scale >= safmx2) {
+		count = 0;
+		while(scale >= safmx2) {
+		    ++count;
+		    f1 *= safmn2;
+		    g1 *= safmn2;
+		    scale = max(Math.abs(f1), Math.abs(g1));
+		}
+		r = Math.sqrt(f1*f1 + g1*g1);
+		cs = f1 / r;
+		sn = g1 / r;
+		i__1 = count;
+		for (i = 1; i <= count; ++i) {
+		    r *= safmx2;
+		}
+	    } else if (scale <= safmn2) {
+		count = 0;
+		while(scale <= safmn2) {
+		    ++count;
+		    f1 *= safmx2;
+		    g1 *= safmx2;
+		    scale = max(Math.abs(f1), Math.abs(g1));
+		}
+		r = Math.sqrt(f1*f1 + g1*g1);
+		cs = f1 / r;
+		sn = g1 / r;
+		i__1 = count;
+		for (i = 1; i <= count; ++i) {
+		    r *= safmn2;
+		}
+	    } else {
+		r = Math.sqrt(f1*f1 + g1*g1);
+		cs = f1 / r;
+		sn = g1 / r;
+	    }
+	    if (Math.abs(f) > Math.abs(g) && cs < 0.0) {
+		cs = -cs;
+		sn = -sn;
+		r = -r;
+	    }
+	}    
+	sin[0] = sn;
+	cos[0] = cs;
+	return r;
+    }
+
+    static double d_sign(double a, double b) {
+	double x;
+	x = (a >= 0 ? a : - a);
+	return (b >= 0 ? x : -x);
+    }
+
+    /**
+     * Creates a new object of the same class as this object.
+     *
+     * @return a clone of this instance.
+     * @exception OutOfMemoryError if there is not enough memory.
+     * @see java.lang.Cloneable
+     * @since Java 3D 1.3
+     */
+    public Object clone() {
+	GMatrix m1 = null;
+	try {
+	    m1 = (GMatrix)super.clone();
+	} catch (CloneNotSupportedException e) {
+	    // this shouldn't happen, since we are Cloneable
+	    throw new InternalError();
+	}
+
+	// Also need to clone array of values
+        m1.values = new double[nRow][nCol];
+	for (int i = 0; i < nRow; i++) {
+	   for(int j = 0; j < nCol; j++) {
+	       m1.values[i][j] = values[i][j];
+	   }
+	}
+
+	return m1;
+    }
+
+}