eu.mihosoft.ext.j3d.javax.vecmath.GMatrix Maven / Gradle / Ivy
The newest version!
/*
* 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 eu.mihosoft.ext.j3d.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 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;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy