summaryrefslogtreecommitdiffstats
path: root/src/javax/vecmath/GMatrix.java
diff options
context:
space:
mode:
Diffstat (limited to 'src/javax/vecmath/GMatrix.java')
-rw-r--r--src/javax/vecmath/GMatrix.java3006
1 files changed, 0 insertions, 3006 deletions
diff --git a/src/javax/vecmath/GMatrix.java b/src/javax/vecmath/GMatrix.java
deleted file mode 100644
index 8b6a1cc..0000000
--- a/src/javax/vecmath/GMatrix.java
+++ /dev/null
@@ -1,3006 +0,0 @@
-/*
- * 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 javax.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;
- }
-
-}
generated by cgit v1.2.3 (git 2.39.1) at 2024-12-11 12:45:50 +0000