Relocate package prefix to org.jogamp.vecmath

Updating the package prefix avoids clashes with old versions of Java 3D.
This is especially important on OS X, where Java 3D 1.3 is sometimes
present on the java.ext.path, taking precedence over the classpath.

1 files changed, 3006 insertions, 0 deletions

diff --git a/src/org/jogamp/vecmath/GMatrix.java b/src/org/jogamp/vecmath/GMatrix.java
new file mode 100644
index 0000000..2c7a6ae
--- /dev/null
+++ b/src/org/jogamp/vecmath/GMatrix.java
@@ -0,0 +1,3006 @@
+/*
+ * Copyright 1997-2008 Sun Microsystems, Inc.  All Rights Reserved.
+ * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
+ *
+ * This code is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License version 2 only, as
+ * published by the Free Software Foundation.  Sun designates this
+ * particular file as subject to the "Classpath" exception as provided
+ * by Sun in the LICENSE file that accompanied this code.
+ *
+ * This code is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+ * version 2 for more details (a copy is included in the LICENSE file that
+ * accompanied this code).
+ *
+ * You should have received a copy of the GNU General Public License version
+ * 2 along with this work; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * Please contact Sun Microsystems, Inc., 4150 Network Circle, Santa Clara,
+ * CA 95054 USA or visit www.sun.com if you need additional information or
+ * have any questions.
+ *
+ */
+
+package org.jogamp.vecmath;
+
+
+/**
+ * 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
+     */
+    @Override
+    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
+     */
+    @Override
+    public int hashCode() {
+	long bits = 1L;
+
+	bits = VecMathUtil.hashLongBits(bits, nRow);
+	bits = VecMathUtil.hashLongBits(bits, nCol);
+
+	for (int i = 0; i < nRow; i++) {
+		for (int j = 0; j < nCol; j++) {
+			bits = VecMathUtil.hashDoubleBits(bits, values[i][j]);
+		}
+	}
+
+	return VecMathUtil.hashFinish(bits);
+    }
+
+
+    /**
+     * 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
+     */
+    @Override
+    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 vecmath 1.3
+     */
+    @Override
+    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;
+    }
+
+}