You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mahout.apache.org by gs...@apache.org on 2009/11/23 16:14:38 UTC
svn commit: r883365 [45/47] - in /lucene/mahout/trunk: ./ examples/ matrix/
matrix/src/ matrix/src/main/ matrix/src/main/java/
matrix/src/main/java/org/ matrix/src/main/java/org/apache/
matrix/src/main/java/org/apache/mahout/ matrix/src/main/java/org/a...
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Property.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Property.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Property.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Property.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,1151 @@
+/*
+Copyright 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.colt.matrix.linalg;
+
+import org.apache.mahout.colt.matrix.DoubleMatrix1D;
+import org.apache.mahout.colt.matrix.DoubleMatrix2D;
+import org.apache.mahout.colt.matrix.DoubleMatrix3D;
+/**
+Tests matrices for linear algebraic properties (equality, tridiagonality, symmetry, singularity, etc).
+<p>
+Except where explicitly indicated, all methods involving equality tests (<tt>==</tt>) allow for numerical instability, to a degree specified upon instance construction and returned by method {@link #tolerance()}.
+The public static final variable <tt>DEFAULT</tt> represents a default Property object with a tolerance of <tt>1.0E-9</tt>.
+The public static final variable <tt>ZERO</tt> represents a Property object with a tolerance of <tt>0.0</tt>.
+The public static final variable <tt>TWELVE</tt> represents a Property object with a tolerance of <tt>1.0E-12</tt>.
+As long as you are happy with these tolerances, there is no need to construct Property objects.
+Simply use idioms like <tt>Property.DEFAULT.equals(A,B)</tt>, <tt>Property.ZERO.equals(A,B)</tt>, <tt>Property.TWELVE.equals(A,B)</tt>.
+<p>
+To work with a different tolerance (e.g. <tt>1.0E-15</tt> or <tt>1.0E-5</tt>) use the constructor and/or method {@link #setTolerance(double)}.
+Note that the public static final Property objects are immutable: Is is not possible to alter their tolerance.
+Any attempt to do so will throw an Exception.
+<p>
+Note that this implementation is not synchronized.
+<p>
+Example: <tt>equals(DoubleMatrix2D A, DoubleMatrix2D B)</tt> is defined as follows
+<table>
+<td class="PRE">
+<pre>
+{ some other tests not related to tolerance go here }
+double epsilon = tolerance();
+for (int row=rows; --row >= 0;) {
+ for (int column=columns; --column >= 0;) {
+ //if (!(A.getQuick(row,column) == B.getQuick(row,column))) return false;
+ if (Math.abs(A.getQuick(row,column) - B.getQuick(row,column)) > epsilon) return false;
+ }
+}
+return true;
+</pre>
+</td>
+</table>
+Here are some example properties
+<table border="1" cellspacing="0">
+ <tr align="left" valign="top">
+ <td valign="middle" align="left"><tt>matrix</tt></td>
+ <td> <tt>4 x 4 <br>
+ 0 0 0 0<br>
+ 0 0 0 0<br>
+ 0 0 0 0<br>
+ 0 0 0 0 </tt></td>
+ <td><tt>4 x 4<br>
+ 1 0 0 0<br>
+ 0 0 0 0<br>
+ 0 0 0 0<br>
+ 0 0 0 1 </tt></td>
+ <td><tt>4 x 4<br>
+ 1 1 0 0<br>
+ 1 1 1 0<br>
+ 0 1 1 1<br>
+ 0 0 1 1 </tt></td>
+ <td><tt> 4 x 4<br>
+ 0 1 1 1<br>
+ 0 1 1 1<br>
+ 0 0 0 1<br>
+ 0 0 0 1 </tt></td>
+ <td><tt> 4 x 4<br>
+ 0 0 0 0<br>
+ 1 1 0 0<br>
+ 1 1 0 0<br>
+ 1 1 1 1 </tt></td>
+ <td><tt>4 x 4<br>
+ 1 1 0 0<br>
+ 0 1 1 0<br>
+ 0 1 0 1<br>
+ 1 0 1 1 </tt><tt> </tt> </td>
+ <td><tt>4 x 4<br>
+ 1 1 1 0<br>
+ 0 1 0 0<br>
+ 1 1 0 1<br>
+ 0 0 1 1 </tt> </td>
+ </tr>
+ <tr align="center" valign="middle">
+ <td><tt>upperBandwidth</tt></td>
+ <td>
+ <div align="center"><tt>0</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>0</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td><tt>3</tt></td>
+ <td align="center" valign="middle"><tt>0</tt></td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>2</tt></div>
+ </td>
+ </tr>
+ <tr align="center" valign="middle">
+ <td><tt>lowerBandwidth</tt></td>
+ <td>
+ <div align="center"><tt>0</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>0</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td><tt>0</tt></td>
+ <td align="center" valign="middle"><tt>3</tt></td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>3</tt></div>
+ </td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>2</tt></div>
+ </td>
+ </tr>
+ <tr align="center" valign="middle">
+ <td><tt>semiBandwidth</tt></td>
+ <td>
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>2</tt></div>
+ </td>
+ <td><tt>4</tt></td>
+ <td align="center" valign="middle"><tt>4</tt></td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>4</tt></div>
+ </td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>3</tt></div>
+ </td>
+ </tr>
+ <tr align="center" valign="middle">
+ <td><tt>description</tt></td>
+ <td>
+ <div align="center"><tt>zero</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>diagonal</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>tridiagonal</tt></div>
+ </td>
+ <td><tt>upper triangular</tt></td>
+ <td align="center" valign="middle"><tt>lower triangular</tt></td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>unstructured</tt></div>
+ </td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>unstructured</tt></div>
+ </td>
+ </tr>
+</table>
+
+@author wolfgang.hoschek@cern.ch
+@version 1.1, 28/May/2000 (fixed strange bugs involving NaN, -inf, inf)
+*/
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class Property extends org.apache.mahout.colt.PersistentObject {
+ /**
+ * The default Property object; currently has <tt>tolerance()==1.0E-9</tt>.
+ */
+ public static final Property DEFAULT = new Property(1.0E-9);
+
+ /**
+ * A Property object with <tt>tolerance()==0.0</tt>.
+ */
+ public static final Property ZERO = new Property(0.0);
+
+ /**
+ * A Property object with <tt>tolerance()==1.0E-12</tt>.
+ */
+ public static final Property TWELVE = new Property(1.0E-12);
+
+ protected double tolerance;
+/**
+ * Not instantiable by no-arg constructor.
+ */
+private Property() {
+ this(1.0E-9); // just to be on the safe side
+}
+/**
+ * Constructs an instance with a tolerance of <tt>Math.abs(newTolerance)</tt>.
+ */
+public Property(double newTolerance) {
+ tolerance = Math.abs(newTolerance);
+}
+/**
+ * Returns a String with <tt>length</tt> blanks.
+ */
+protected static String blanks(int length) {
+ if (length <0 ) length = 0;
+ StringBuffer buf = new StringBuffer(length);
+ for (int k = 0; k < length; k++) {
+ buf.append(' ');
+ }
+ return buf.toString();
+}
+/**
+ * Checks whether the given matrix <tt>A</tt> is <i>rectangular</i>.
+ * @throws IllegalArgumentException if <tt>A.rows() < A.columns()</tt>.
+ */
+public void checkRectangular(DoubleMatrix2D A) {
+ if (A.rows() < A.columns()) {
+ throw new IllegalArgumentException("Matrix must be rectangular: "+ org.apache.mahout.colt.matrix.doublealgo.Formatter.shape(A));
+ }
+}
+/**
+ * Checks whether the given matrix <tt>A</tt> is <i>square</i>.
+ * @throws IllegalArgumentException if <tt>A.rows() != A.columns()</tt>.
+ */
+public void checkSquare(DoubleMatrix2D A) {
+ if (A.rows() != A.columns()) throw new IllegalArgumentException("Matrix must be square: "+ org.apache.mahout.colt.matrix.doublealgo.Formatter.shape(A));
+}
+/**
+ * Returns the matrix's fraction of non-zero cells; <tt>A.cardinality() / A.size()</tt>.
+ */
+public double density(DoubleMatrix2D A) {
+ return A.cardinality() / (double) A.size();
+}
+/**
+ * Returns whether all cells of the given matrix <tt>A</tt> are equal to the given value.
+ * The result is <tt>true</tt> if and only if <tt>A != null</tt> and
+ * <tt>! (Math.abs(value - A[i]) > tolerance())</tt> holds for all coordinates.
+ * @param A the first matrix to compare.
+ * @param value the value to compare against.
+ * @return <tt>true</tt> if the matrix is equal to the value;
+ * <tt>false</tt> otherwise.
+ */
+public boolean equals(DoubleMatrix1D A, double value) {
+ if (A==null) return false;
+ double epsilon = tolerance();
+ for (int i = A.size(); --i >= 0;) {
+ //if (!(A.getQuick(i) == value)) return false;
+ //if (Math.abs(value - A.getQuick(i)) > epsilon) return false;
+ double x = A.getQuick(i);
+ double diff = Math.abs(value - x);
+ if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
+ if (!(diff <= epsilon)) return false;
+ }
+ return true;
+}
+/**
+ * Returns whether both given matrices <tt>A</tt> and <tt>B</tt> are equal.
+ * The result is <tt>true</tt> if <tt>A==B</tt>.
+ * Otherwise, the result is <tt>true</tt> if and only if both arguments are <tt>!= null</tt>,
+ * have the same size and
+ * <tt>! (Math.abs(A[i] - B[i]) > tolerance())</tt> holds for all indexes.
+ * @param A the first matrix to compare.
+ * @param B the second matrix to compare.
+ * @return <tt>true</tt> if both matrices are equal;
+ * <tt>false</tt> otherwise.
+ */
+public boolean equals(DoubleMatrix1D A, DoubleMatrix1D B) {
+ if (A==B) return true;
+ if (! (A != null && B != null)) return false;
+ int size = A.size();
+ if (size != B.size()) return false;
+
+ double epsilon = tolerance();
+ for (int i=size; --i >= 0;) {
+ //if (!(getQuick(i) == B.getQuick(i))) return false;
+ //if (Math.abs(A.getQuick(i) - B.getQuick(i)) > epsilon) return false;
+ double x = A.getQuick(i);
+ double value = B.getQuick(i);
+ double diff = Math.abs(value - x);
+ if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
+ if (!(diff <= epsilon)) return false;
+ }
+ return true;
+}
+/**
+ * Returns whether all cells of the given matrix <tt>A</tt> are equal to the given value.
+ * The result is <tt>true</tt> if and only if <tt>A != null</tt> and
+ * <tt>! (Math.abs(value - A[row,col]) > tolerance())</tt> holds for all coordinates.
+ * @param A the first matrix to compare.
+ * @param value the value to compare against.
+ * @return <tt>true</tt> if the matrix is equal to the value;
+ * <tt>false</tt> otherwise.
+ */
+public boolean equals(DoubleMatrix2D A, double value) {
+ if (A==null) return false;
+ int rows = A.rows();
+ int columns = A.columns();
+
+ double epsilon = tolerance();
+ for (int row=rows; --row >= 0;) {
+ for (int column=columns; --column >= 0;) {
+ //if (!(A.getQuick(row,column) == value)) return false;
+ //if (Math.abs(value - A.getQuick(row,column)) > epsilon) return false;
+ double x = A.getQuick(row,column);
+ double diff = Math.abs(value - x);
+ if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
+ if (!(diff <= epsilon)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * Returns whether both given matrices <tt>A</tt> and <tt>B</tt> are equal.
+ * The result is <tt>true</tt> if <tt>A==B</tt>.
+ * Otherwise, the result is <tt>true</tt> if and only if both arguments are <tt>!= null</tt>,
+ * have the same number of columns and rows and
+ * <tt>! (Math.abs(A[row,col] - B[row,col]) > tolerance())</tt> holds for all coordinates.
+ * @param A the first matrix to compare.
+ * @param B the second matrix to compare.
+ * @return <tt>true</tt> if both matrices are equal;
+ * <tt>false</tt> otherwise.
+ */
+public boolean equals(DoubleMatrix2D A, DoubleMatrix2D B) {
+ if (A==B) return true;
+ if (! (A != null && B != null)) return false;
+ int rows = A.rows();
+ int columns = A.columns();
+ if (columns != B.columns() || rows != B.rows()) return false;
+
+ double epsilon = tolerance();
+ for (int row=rows; --row >= 0;) {
+ for (int column=columns; --column >= 0;) {
+ //if (!(A.getQuick(row,column) == B.getQuick(row,column))) return false;
+ //if (Math.abs((A.getQuick(row,column) - B.getQuick(row,column)) > epsilon) return false;
+ double x = A.getQuick(row,column);
+ double value = B.getQuick(row,column);
+ double diff = Math.abs(value - x);
+ if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
+ if (!(diff <= epsilon)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * Returns whether all cells of the given matrix <tt>A</tt> are equal to the given value.
+ * The result is <tt>true</tt> if and only if <tt>A != null</tt> and
+ * <tt>! (Math.abs(value - A[slice,row,col]) > tolerance())</tt> holds for all coordinates.
+ * @param A the first matrix to compare.
+ * @param value the value to compare against.
+ * @return <tt>true</tt> if the matrix is equal to the value;
+ * <tt>false</tt> otherwise.
+ */
+public boolean equals(DoubleMatrix3D A, double value) {
+ if (A==null) return false;
+ int rows = A.rows();
+ int columns = A.columns();
+
+ double epsilon = tolerance();
+ for (int slice=A.slices(); --slice >= 0;) {
+ for (int row=rows; --row >= 0;) {
+ for (int column=columns; --column >= 0;) {
+ //if (!(A.getQuick(slice,row,column) == value)) return false;
+ //if (Math.abs(value - A.getQuick(slice,row,column)) > epsilon) return false;
+ double x = A.getQuick(slice,row,column);
+ double diff = Math.abs(value - x);
+ if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
+ if (!(diff <= epsilon)) return false;
+ }
+ }
+ }
+ return true;
+}
+/**
+ * Returns whether both given matrices <tt>A</tt> and <tt>B</tt> are equal.
+ * The result is <tt>true</tt> if <tt>A==B</tt>.
+ * Otherwise, the result is <tt>true</tt> if and only if both arguments are <tt>!= null</tt>,
+ * have the same number of columns, rows and slices, and
+ * <tt>! (Math.abs(A[slice,row,col] - B[slice,row,col]) > tolerance())</tt> holds for all coordinates.
+ * @param A the first matrix to compare.
+ * @param B the second matrix to compare.
+ * @return <tt>true</tt> if both matrices are equal;
+ * <tt>false</tt> otherwise.
+ */
+public boolean equals(DoubleMatrix3D A, DoubleMatrix3D B) {
+ if (A==B) return true;
+ if (! (A != null && B != null)) return false;
+ int slices = A.slices();
+ int rows = A.rows();
+ int columns = A.columns();
+ if (columns != B.columns() || rows != B.rows() || slices != B.slices()) return false;
+
+ double epsilon = tolerance();
+ for (int slice=slices; --slice >= 0;) {
+ for (int row=rows; --row >= 0;) {
+ for (int column=columns; --column >= 0;) {
+ //if (!(A.getQuick(slice,row,column) == B.getQuick(slice,row,column))) return false;
+ //if (Math.abs(A.getQuick(slice,row,column) - B.getQuick(slice,row,column)) > epsilon) return false;
+ double x = A.getQuick(slice,row,column);
+ double value = B.getQuick(slice,row,column);
+ double diff = Math.abs(value - x);
+ if ((diff!=diff) && ((value!=value && x!=x) || value==x)) diff = 0;
+ if (!(diff <= epsilon)) return false;
+ }
+ }
+ }
+ return true;
+}
+/**
+Modifies the given matrix square matrix <tt>A</tt> such that it is diagonally dominant by row and column, hence non-singular, hence invertible.
+For testing purposes only.
+@param A the square matrix to modify.
+@throws IllegalArgumentException if <tt>!isSquare(A)</tt>.
+*/
+public void generateNonSingular(DoubleMatrix2D A) {
+ checkSquare(A);
+ org.apache.mahout.jet.math.Functions F = org.apache.mahout.jet.math.Functions.functions;
+ int min = Math.min(A.rows(), A.columns());
+ for (int i = min; --i >= 0; ) {
+ A.setQuick(i,i, 0);
+ }
+ for (int i = min; --i >= 0; ) {
+ double rowSum = A.viewRow(i).aggregate(F.plus,F.abs);
+ double colSum = A.viewColumn(i).aggregate(F.plus,F.abs);
+ A.setQuick(i,i, Math.max(rowSum,colSum) + i+1);
+ }
+}
+/**
+ */
+protected static String get(org.apache.mahout.colt.list.ObjectArrayList list, int index) {
+ return ((String) list.get(index));
+}
+/**
+ * A matrix <tt>A</tt> is <i>diagonal</i> if <tt>A[i,j] == 0</tt> whenever <tt>i != j</tt>.
+ * Matrix may but need not be square.
+ */
+public boolean isDiagonal(DoubleMatrix2D A) {
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int row = rows; --row >=0; ) {
+ for (int column = columns; --column >= 0; ) {
+ //if (row!=column && A.getQuick(row,column) != 0) return false;
+ if (row!=column && !(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>diagonally dominant by column</i> if the absolute value of each diagonal element is larger than the sum of the absolute values of the off-diagonal elements in the corresponding column.
+ * <tt>returns true if for all i: abs(A[i,i]) > Sum(abs(A[j,i])); j != i.</tt>
+ * Matrix may but need not be square.
+ * <p>
+ * Note: Ignores tolerance.
+ */
+public boolean isDiagonallyDominantByColumn(DoubleMatrix2D A) {
+ org.apache.mahout.jet.math.Functions F = org.apache.mahout.jet.math.Functions.functions;
+ double epsilon = tolerance();
+ int min = Math.min(A.rows(), A.columns());
+ for (int i = min; --i >= 0; ) {
+ double diag = Math.abs(A.getQuick(i,i));
+ diag += diag;
+ if (diag <= A.viewColumn(i).aggregate(F.plus,F.abs)) return false;
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>diagonally dominant by row</i> if the absolute value of each diagonal element is larger than the sum of the absolute values of the off-diagonal elements in the corresponding row.
+ * <tt>returns true if for all i: abs(A[i,i]) > Sum(abs(A[i,j])); j != i.</tt>
+ * Matrix may but need not be square.
+ * <p>
+ * Note: Ignores tolerance.
+ */
+public boolean isDiagonallyDominantByRow(DoubleMatrix2D A) {
+ org.apache.mahout.jet.math.Functions F = org.apache.mahout.jet.math.Functions.functions;
+ double epsilon = tolerance();
+ int min = Math.min(A.rows(), A.columns());
+ for (int i = min; --i >= 0; ) {
+ double diag = Math.abs(A.getQuick(i,i));
+ diag += diag;
+ if (diag <= A.viewRow(i).aggregate(F.plus,F.abs)) return false;
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is an <i>identity</i> matrix if <tt>A[i,i] == 1</tt> and all other cells are zero.
+ * Matrix may but need not be square.
+ */
+public boolean isIdentity(DoubleMatrix2D A) {
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int row = rows; --row >=0; ) {
+ for (int column = columns; --column >= 0; ) {
+ double v = A.getQuick(row,column);
+ if (row==column) {
+ if (!(Math.abs(1-v) < epsilon)) return false;
+ }
+ else if (!(Math.abs(v) <= epsilon)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>lower bidiagonal</i> if <tt>A[i,j]==0</tt> unless <tt>i==j || i==j+1</tt>.
+ * Matrix may but need not be square.
+ */
+public boolean isLowerBidiagonal(DoubleMatrix2D A) {
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int row = rows; --row >=0; ) {
+ for (int column = columns; --column >= 0; ) {
+ if (!(row==column || row==column+1)) {
+ //if (A.getQuick(row,column) != 0) return false;
+ if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
+ }
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>lower triangular</i> if <tt>A[i,j]==0</tt> whenever <tt>i < j</tt>.
+ * Matrix may but need not be square.
+ */
+public boolean isLowerTriangular(DoubleMatrix2D A) {
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int column = columns; --column >= 0; ) {
+ for (int row = Math.min(column,rows); --row >= 0; ) {
+ //if (A.getQuick(row,column) != 0) return false;
+ if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>non-negative</i> if <tt>A[i,j] >= 0</tt> holds for all cells.
+ * <p>
+ * Note: Ignores tolerance.
+ */
+public boolean isNonNegative(DoubleMatrix2D A) {
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int row = rows; --row >=0; ) {
+ for (int column = columns; --column >= 0; ) {
+ if (! (A.getQuick(row,column) >= 0)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * A square matrix <tt>A</tt> is <i>orthogonal</i> if <tt>A*transpose(A) = I</tt>.
+ * @throws IllegalArgumentException if <tt>!isSquare(A)</tt>.
+ */
+public boolean isOrthogonal(DoubleMatrix2D A) {
+ checkSquare(A);
+ return equals(A.zMult(A,null,1,0,false,true), org.apache.mahout.colt.matrix.DoubleFactory2D.dense.identity(A.rows()));
+}
+/**
+ * A matrix <tt>A</tt> is <i>positive</i> if <tt>A[i,j] > 0</tt> holds for all cells.
+ * <p>
+ * Note: Ignores tolerance.
+ */
+public boolean isPositive(DoubleMatrix2D A) {
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int row = rows; --row >=0; ) {
+ for (int column = columns; --column >= 0; ) {
+ if (!(A.getQuick(row,column) > 0)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>singular</i> if it has no inverse, that is, iff <tt>det(A)==0</tt>.
+ */
+public boolean isSingular(DoubleMatrix2D A) {
+ return !(Math.abs(Algebra.DEFAULT.det(A)) >= tolerance());
+}
+/**
+ * A square matrix <tt>A</tt> is <i>skew-symmetric</i> if <tt>A = -transpose(A)</tt>, that is <tt>A[i,j] == -A[j,i]</tt>.
+ * @throws IllegalArgumentException if <tt>!isSquare(A)</tt>.
+ */
+public boolean isSkewSymmetric(DoubleMatrix2D A) {
+ checkSquare(A);
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int row = rows; --row >=0; ) {
+ for (int column = rows; --column >= 0; ) {
+ //if (A.getQuick(row,column) != -A.getQuick(column,row)) return false;
+ if (!(Math.abs(A.getQuick(row,column) + A.getQuick(column,row)) <= epsilon)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>square</i> if it has the same number of rows and columns.
+ */
+public boolean isSquare(DoubleMatrix2D A) {
+ return A.rows() == A.columns();
+}
+/**
+ * A matrix <tt>A</tt> is <i>strictly lower triangular</i> if <tt>A[i,j]==0</tt> whenever <tt>i <= j</tt>.
+ * Matrix may but need not be square.
+ */
+public boolean isStrictlyLowerTriangular(DoubleMatrix2D A) {
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int column = columns; --column >= 0; ) {
+ for (int row = Math.min(rows,column+1); --row >= 0; ) {
+ //if (A.getQuick(row,column) != 0) return false;
+ if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>strictly triangular</i> if it is triangular and its diagonal elements all equal 0.
+ * Matrix may but need not be square.
+ */
+public boolean isStrictlyTriangular(DoubleMatrix2D A) {
+ if (!isTriangular(A)) return false;
+
+ double epsilon = tolerance();
+ for (int i = Math.min(A.rows(), A.columns()); --i >= 0; ) {
+ //if (A.getQuick(i,i) != 0) return false;
+ if (!(Math.abs(A.getQuick(i,i)) <= epsilon)) return false;
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>strictly upper triangular</i> if <tt>A[i,j]==0</tt> whenever <tt>i >= j</tt>.
+ * Matrix may but need not be square.
+ */
+public boolean isStrictlyUpperTriangular(DoubleMatrix2D A) {
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int column = columns; --column >= 0; ) {
+ for (int row = rows; --row >= column; ) {
+ //if (A.getQuick(row,column) != 0) return false;
+ if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>symmetric</i> if <tt>A = tranpose(A)</tt>, that is <tt>A[i,j] == A[j,i]</tt>.
+ * @throws IllegalArgumentException if <tt>!isSquare(A)</tt>.
+ */
+public boolean isSymmetric(DoubleMatrix2D A) {
+ checkSquare(A);
+ return equals(A,A.viewDice());
+}
+/**
+ * A matrix <tt>A</tt> is <i>triangular</i> iff it is either upper or lower triangular.
+ * Matrix may but need not be square.
+ */
+public boolean isTriangular(DoubleMatrix2D A) {
+ return isLowerTriangular(A) || isUpperTriangular(A);
+}
+/**
+ * A matrix <tt>A</tt> is <i>tridiagonal</i> if <tt>A[i,j]==0</tt> whenever <tt>Math.abs(i-j) > 1</tt>.
+ * Matrix may but need not be square.
+ */
+public boolean isTridiagonal(DoubleMatrix2D A) {
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int row = rows; --row >=0; ) {
+ for (int column = columns; --column >= 0; ) {
+ if (Math.abs(row-column) > 1) {
+ //if (A.getQuick(row,column) != 0) return false;
+ if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
+ }
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>unit triangular</i> if it is triangular and its diagonal elements all equal 1.
+ * Matrix may but need not be square.
+ */
+public boolean isUnitTriangular(DoubleMatrix2D A) {
+ if (!isTriangular(A)) return false;
+
+ double epsilon = tolerance();
+ for (int i = Math.min(A.rows(), A.columns()); --i >= 0; ) {
+ //if (A.getQuick(i,i) != 1) return false;
+ if (!(Math.abs(1 - A.getQuick(i,i)) <= epsilon)) return false;
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>upper bidiagonal</i> if <tt>A[i,j]==0</tt> unless <tt>i==j || i==j-1</tt>.
+ * Matrix may but need not be square.
+ */
+public boolean isUpperBidiagonal(DoubleMatrix2D A) {
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int row = rows; --row >=0; ) {
+ for (int column = columns; --column >= 0; ) {
+ if (!(row==column || row==column-1)) {
+ //if (A.getQuick(row,column) != 0) return false;
+ if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
+ }
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>upper triangular</i> if <tt>A[i,j]==0</tt> whenever <tt>i > j</tt>.
+ * Matrix may but need not be square.
+ */
+public boolean isUpperTriangular(DoubleMatrix2D A) {
+ double epsilon = tolerance();
+ int rows = A.rows();
+ int columns = A.columns();
+ for (int column = columns; --column >= 0; ) {
+ for (int row = rows; --row > column; ) {
+ //if (A.getQuick(row,column) != 0) return false;
+ if (!(Math.abs(A.getQuick(row,column)) <= epsilon)) return false;
+ }
+ }
+ return true;
+}
+/**
+ * A matrix <tt>A</tt> is <i>zero</i> if all its cells are zero.
+ */
+public boolean isZero(DoubleMatrix2D A) {
+ return equals(A,0);
+}
+/**
+The <i>lower bandwidth</i> of a square matrix <tt>A</tt> is the maximum <tt>i-j</tt> for which <tt>A[i,j]</tt> is nonzero and <tt>i > j</tt>.
+A <i>banded</i> matrix has a "band" about the diagonal.
+Diagonal, tridiagonal and triangular matrices are special cases.
+
+@param A the square matrix to analyze.
+@return the lower bandwith.
+@throws IllegalArgumentException if <tt>!isSquare(A)</tt>.
+@see #semiBandwidth(DoubleMatrix2D)
+@see #upperBandwidth(DoubleMatrix2D)
+*/
+public int lowerBandwidth(DoubleMatrix2D A) {
+ checkSquare(A);
+ double epsilon = tolerance();
+ int rows = A.rows();
+
+ for (int k=rows; --k >= 0; ) {
+ for (int i = rows-k; --i >= 0; ) {
+ int j = i+k;
+ //if (A.getQuick(j,i) != 0) return k;
+ if (!(Math.abs(A.getQuick(j,i)) <= epsilon)) return k;
+ }
+ }
+ return 0;
+}
+/**
+Returns the <i>semi-bandwidth</i> of the given square matrix <tt>A</tt>.
+A <i>banded</i> matrix has a "band" about the diagonal.
+It is a matrix with all cells equal to zero,
+with the possible exception of the cells along the diagonal line,
+the <tt>k</tt> diagonal lines above the diagonal, and the <tt>k</tt> diagonal lines below the diagonal.
+The <i>semi-bandwith l</i> is the number <tt>k+1</tt>.
+The <i>bandwidth p</i> is the number <tt>2*k + 1</tt>.
+For example, a tridiagonal matrix corresponds to <tt>k=1, l=2, p=3</tt>,
+a diagonal or zero matrix corresponds to <tt>k=0, l=1, p=1</tt>,
+<p>
+The <i>upper bandwidth</i> is the maximum <tt>j-i</tt> for which <tt>A[i,j]</tt> is nonzero and <tt>j > i</tt>.
+The <i>lower bandwidth</i> is the maximum <tt>i-j</tt> for which <tt>A[i,j]</tt> is nonzero and <tt>i > j</tt>.
+Diagonal, tridiagonal and triangular matrices are special cases.
+<p>
+Examples:
+<table border="1" cellspacing="0">
+ <tr align="left" valign="top">
+ <td valign="middle" align="left"><tt>matrix</tt></td>
+ <td> <tt>4 x 4 <br>
+ 0 0 0 0<br>
+ 0 0 0 0<br>
+ 0 0 0 0<br>
+ 0 0 0 0 </tt></td>
+ <td><tt>4 x 4<br>
+ 1 0 0 0<br>
+ 0 0 0 0<br>
+ 0 0 0 0<br>
+ 0 0 0 1 </tt></td>
+ <td><tt>4 x 4<br>
+ 1 1 0 0<br>
+ 1 1 1 0<br>
+ 0 1 1 1<br>
+ 0 0 1 1 </tt></td>
+ <td><tt> 4 x 4<br>
+ 0 1 1 1<br>
+ 0 1 1 1<br>
+ 0 0 0 1<br>
+ 0 0 0 1 </tt></td>
+ <td><tt> 4 x 4<br>
+ 0 0 0 0<br>
+ 1 1 0 0<br>
+ 1 1 0 0<br>
+ 1 1 1 1 </tt></td>
+ <td><tt>4 x 4<br>
+ 1 1 0 0<br>
+ 0 1 1 0<br>
+ 0 1 0 1<br>
+ 1 0 1 1 </tt><tt> </tt> </td>
+ <td><tt>4 x 4<br>
+ 1 1 1 0<br>
+ 0 1 0 0<br>
+ 1 1 0 1<br>
+ 0 0 1 1 </tt> </td>
+ </tr>
+ <tr align="center" valign="middle">
+ <td><tt>upperBandwidth</tt></td>
+ <td>
+ <div align="center"><tt>0</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>0</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td><tt>3</tt></td>
+ <td align="center" valign="middle"><tt>0</tt></td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>2</tt></div>
+ </td>
+ </tr>
+ <tr align="center" valign="middle">
+ <td><tt>lowerBandwidth</tt></td>
+ <td>
+ <div align="center"><tt>0</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>0</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td><tt>0</tt></td>
+ <td align="center" valign="middle"><tt>3</tt></td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>3</tt></div>
+ </td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>2</tt></div>
+ </td>
+ </tr>
+ <tr align="center" valign="middle">
+ <td><tt>semiBandwidth</tt></td>
+ <td>
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>1</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>2</tt></div>
+ </td>
+ <td><tt>4</tt></td>
+ <td align="center" valign="middle"><tt>4</tt></td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>4</tt></div>
+ </td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>3</tt></div>
+ </td>
+ </tr>
+ <tr align="center" valign="middle">
+ <td><tt>description</tt></td>
+ <td>
+ <div align="center"><tt>zero</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>diagonal</tt></div>
+ </td>
+ <td>
+ <div align="center"><tt>tridiagonal</tt></div>
+ </td>
+ <td><tt>upper triangular</tt></td>
+ <td align="center" valign="middle"><tt>lower triangular</tt></td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>unstructured</tt></div>
+ </td>
+ <td align="center" valign="middle">
+ <div align="center"><tt>unstructured</tt></div>
+ </td>
+ </tr>
+</table>
+
+@param A the square matrix to analyze.
+@return the semi-bandwith <tt>l</tt>.
+@throws IllegalArgumentException if <tt>!isSquare(A)</tt>.
+@see #lowerBandwidth(DoubleMatrix2D)
+@see #upperBandwidth(DoubleMatrix2D)
+*/
+public int semiBandwidth(DoubleMatrix2D A) {
+ checkSquare(A);
+ double epsilon = tolerance();
+ int rows = A.rows();
+
+ for (int k=rows; --k >= 0; ) {
+ for (int i = rows-k; --i >= 0; ) {
+ int j = i+k;
+ //if (A.getQuick(j,i) != 0) return k+1;
+ //if (A.getQuick(i,j) != 0) return k+1;
+ if (!(Math.abs(A.getQuick(j,i)) <= epsilon)) return k+1;
+ if (!(Math.abs(A.getQuick(i,j)) <= epsilon)) return k+1;
+ }
+ }
+ return 1;
+}
+/**
+ * Sets the tolerance to <tt>Math.abs(newTolerance)</tt>.
+ * @throws UnsupportedOperationException if <tt>this==DEFAULT || this==ZERO || this==TWELVE</tt>.
+ */
+public void setTolerance(double newTolerance) {
+ if (this==DEFAULT || this==ZERO || this==TWELVE) {
+ throw new IllegalArgumentException("Attempted to modify immutable object.");
+ //throw new UnsupportedOperationException("Attempted to modify object."); // since JDK1.2
+ }
+ tolerance = Math.abs(newTolerance);
+}
+/**
+ * Returns the current tolerance.
+ */
+public double tolerance() {
+ return tolerance;
+}
+/**
+Returns summary information about the given matrix <tt>A</tt>.
+That is a String with (propertyName, propertyValue) pairs.
+Useful for debugging or to quickly get the rough picture of a matrix.
+For example,
+<pre>
+density : 0.9
+isDiagonal : false
+isDiagonallyDominantByRow : false
+isDiagonallyDominantByColumn : false
+isIdentity : false
+isLowerBidiagonal : false
+isLowerTriangular : false
+isNonNegative : true
+isOrthogonal : Illegal operation or error: Matrix must be square.
+isPositive : true
+isSingular : Illegal operation or error: Matrix must be square.
+isSkewSymmetric : Illegal operation or error: Matrix must be square.
+isSquare : false
+isStrictlyLowerTriangular : false
+isStrictlyTriangular : false
+isStrictlyUpperTriangular : false
+isSymmetric : Illegal operation or error: Matrix must be square.
+isTriangular : false
+isTridiagonal : false
+isUnitTriangular : false
+isUpperBidiagonal : false
+isUpperTriangular : false
+isZero : false
+lowerBandwidth : Illegal operation or error: Matrix must be square.
+semiBandwidth : Illegal operation or error: Matrix must be square.
+upperBandwidth : Illegal operation or error: Matrix must be square.
+</pre>
+*/
+public String toString(DoubleMatrix2D A) {
+ final org.apache.mahout.colt.list.ObjectArrayList names = new org.apache.mahout.colt.list.ObjectArrayList();
+ final org.apache.mahout.colt.list.ObjectArrayList values = new org.apache.mahout.colt.list.ObjectArrayList();
+ String unknown = "Illegal operation or error: ";
+
+ // determine properties
+ names.add("density");
+ try { values.add(String.valueOf(density(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ // determine properties
+ names.add("isDiagonal");
+ try { values.add(String.valueOf(isDiagonal(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ // determine properties
+ names.add("isDiagonallyDominantByRow");
+ try { values.add(String.valueOf(isDiagonallyDominantByRow(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ // determine properties
+ names.add("isDiagonallyDominantByColumn");
+ try { values.add(String.valueOf(isDiagonallyDominantByColumn(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isIdentity");
+ try { values.add(String.valueOf(isIdentity(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isLowerBidiagonal");
+ try { values.add(String.valueOf(isLowerBidiagonal(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isLowerTriangular");
+ try { values.add(String.valueOf(isLowerTriangular(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isNonNegative");
+ try { values.add(String.valueOf(isNonNegative(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isOrthogonal");
+ try { values.add(String.valueOf(isOrthogonal(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isPositive");
+ try { values.add(String.valueOf(isPositive(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isSingular");
+ try { values.add(String.valueOf(isSingular(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isSkewSymmetric");
+ try { values.add(String.valueOf(isSkewSymmetric(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isSquare");
+ try { values.add(String.valueOf(isSquare(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isStrictlyLowerTriangular");
+ try { values.add(String.valueOf(isStrictlyLowerTriangular(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isStrictlyTriangular");
+ try { values.add(String.valueOf(isStrictlyTriangular(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isStrictlyUpperTriangular");
+ try { values.add(String.valueOf(isStrictlyUpperTriangular(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isSymmetric");
+ try { values.add(String.valueOf(isSymmetric(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isTriangular");
+ try { values.add(String.valueOf(isTriangular(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isTridiagonal");
+ try { values.add(String.valueOf(isTridiagonal(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isUnitTriangular");
+ try { values.add(String.valueOf(isUnitTriangular(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isUpperBidiagonal");
+ try { values.add(String.valueOf(isUpperBidiagonal(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isUpperTriangular");
+ try { values.add(String.valueOf(isUpperTriangular(A)));}
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("isZero");
+ try { values.add(String.valueOf(isZero(A))); }
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("lowerBandwidth");
+ try { values.add(String.valueOf(lowerBandwidth(A))); }
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("semiBandwidth");
+ try { values.add(String.valueOf(semiBandwidth(A))); }
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+ names.add("upperBandwidth");
+ try { values.add(String.valueOf(upperBandwidth(A))); }
+ catch (IllegalArgumentException exc) { values.add(unknown+exc.getMessage()); }
+
+
+ // sort ascending by property name
+ org.apache.mahout.colt.function.IntComparator comp = new org.apache.mahout.colt.function.IntComparator() {
+ public int compare(int a, int b) {
+ return get(names,a).compareTo(get(names,b));
+ }
+ };
+ org.apache.mahout.colt.Swapper swapper = new org.apache.mahout.colt.Swapper() {
+ public void swap(int a, int b) {
+ Object tmp;
+ tmp = names.get(a); names.set(a,names.get(b)); names.set(b,tmp);
+ tmp = values.get(a); values.set(a,values.get(b)); values.set(b,tmp);
+ }
+ };
+ org.apache.mahout.colt.GenericSorting.quickSort(0,names.size(),comp,swapper);
+
+ // determine padding for nice formatting
+ int maxLength = 0;
+ for (int i = 0; i < names.size(); i++) {
+ int length = ((String) names.get(i)).length();
+ maxLength = Math.max(length, maxLength);
+ }
+
+ // finally, format properties
+ StringBuffer buf = new StringBuffer();
+ for (int i = 0; i < names.size(); i++) {
+ String name = ((String) names.get(i));
+ buf.append(name);
+ buf.append(blanks(maxLength - name.length()));
+ buf.append(" : ");
+ buf.append(values.get(i));
+ if (i < names.size() - 1)
+ buf.append('\n');
+ }
+
+ return buf.toString();
+}
+/**
+The <i>upper bandwidth</i> of a square matrix <tt>A</tt> is the
+maximum <tt>j-i</tt> for which <tt>A[i,j]</tt> is nonzero and <tt>j > i</tt>.
+A <i>banded</i> matrix has a "band" about the diagonal.
+Diagonal, tridiagonal and triangular matrices are special cases.
+
+@param A the square matrix to analyze.
+@return the upper bandwith.
+@throws IllegalArgumentException if <tt>!isSquare(A)</tt>.
+@see #semiBandwidth(DoubleMatrix2D)
+@see #lowerBandwidth(DoubleMatrix2D)
+*/
+public int upperBandwidth(DoubleMatrix2D A) {
+ checkSquare(A);
+ double epsilon = tolerance();
+ int rows = A.rows();
+
+ for (int k=rows; --k >= 0; ) {
+ for (int i = rows-k; --i >= 0; ) {
+ int j = i+k;
+ //if (A.getQuick(i,j) != 0) return k;
+ if (!(Math.abs(A.getQuick(i,j)) <= epsilon)) return k;
+ }
+ }
+ return 0;
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Property.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/QRDecomposition.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/QRDecomposition.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/QRDecomposition.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/QRDecomposition.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,266 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.colt.matrix.linalg;
+
+import org.apache.mahout.colt.matrix.DoubleMatrix1D;
+import org.apache.mahout.colt.matrix.DoubleMatrix2D;
+
+/**
+For an <tt>m x n</tt> matrix <tt>A</tt> with <tt>m >= n</tt>, the QR decomposition is an <tt>m x n</tt>
+orthogonal matrix <tt>Q</tt> and an <tt>n x n</tt> upper triangular matrix <tt>R</tt> so that
+<tt>A = Q*R</tt>.
+<P>
+The QR decompostion always exists, even if the matrix does not have
+full rank, so the constructor will never fail. The primary use of the
+QR decomposition is in the least squares solution of nonsquare systems
+of simultaneous linear equations. This will fail if <tt>isFullRank()</tt>
+returns <tt>false</tt>.
+*/
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class QRDecomposition implements java.io.Serializable {
+ static final long serialVersionUID = 1020;
+ /** Array for internal storage of decomposition.
+ @serial internal array storage.
+ */
+ private DoubleMatrix2D QR;
+ //private double[][] QR;
+
+ /** Row and column dimensions.
+ @serial column dimension.
+ @serial row dimension.
+ */
+ private int m, n;
+
+ /** Array for internal storage of diagonal of R.
+ @serial diagonal of R.
+ */
+ private DoubleMatrix1D Rdiag;
+/**
+Constructs and returns a new QR decomposition object; computed by Householder reflections;
+The decomposed matrices can be retrieved via instance methods of the returned decomposition object.
+@param A A rectangular matrix.
+@return a decomposition object to access <tt>R</tt> and the Householder vectors <tt>H</tt>, and to compute <tt>Q</tt>.
+@throws IllegalArgumentException if <tt>A.rows() < A.columns()</tt>.
+*/
+
+public QRDecomposition (DoubleMatrix2D A) {
+ Property.DEFAULT.checkRectangular(A);
+
+ org.apache.mahout.jet.math.Functions F = org.apache.mahout.jet.math.Functions.functions;
+ // Initialize.
+ QR = A.copy();
+ m = A.rows();
+ n = A.columns();
+ Rdiag = A.like1D(n);
+ //Rdiag = new double[n];
+ org.apache.mahout.colt.function.DoubleDoubleFunction hypot = Algebra.hypotFunction();
+
+ // precompute and cache some views to avoid regenerating them time and again
+ DoubleMatrix1D[] QRcolumns = new DoubleMatrix1D[n];
+ DoubleMatrix1D[] QRcolumnsPart = new DoubleMatrix1D[n];
+ for (int k = 0; k < n; k++) {
+ QRcolumns[k] = QR.viewColumn(k);
+ QRcolumnsPart[k] = QR.viewColumn(k).viewPart(k,m-k);
+ }
+
+ // Main loop.
+ for (int k = 0; k < n; k++) {
+ //DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k);
+ // Compute 2-norm of k-th column without under/overflow.
+ double nrm = 0;
+ //if (k<m) nrm = QRcolumnsPart[k].aggregate(hypot,F.identity);
+
+ for (int i = k; i < m; i++) { // fixes bug reported by hong.44@osu.edu
+ nrm = Algebra.hypot(nrm,QR.getQuick(i,k));
+ }
+
+
+ if (nrm != 0.0) {
+ // Form k-th Householder vector.
+ if (QR.getQuick(k,k) < 0) nrm = -nrm;
+ QRcolumnsPart[k].assign(org.apache.mahout.jet.math.Functions.div(nrm));
+ /*
+ for (int i = k; i < m; i++) {
+ QR[i][k] /= nrm;
+ }
+ */
+
+ QR.setQuick(k,k, QR.getQuick(k,k) + 1);
+
+ // Apply transformation to remaining columns.
+ for (int j = k+1; j < n; j++) {
+ DoubleMatrix1D QRcolj = QR.viewColumn(j).viewPart(k,m-k);
+ double s = QRcolumnsPart[k].zDotProduct(QRcolj);
+ /*
+ // fixes bug reported by John Chambers
+ DoubleMatrix1D QRcolj = QR.viewColumn(j).viewPart(k,m-k);
+ double s = QRcolumnsPart[k].zDotProduct(QRcolumns[j]);
+ double s = 0.0;
+ for (int i = k; i < m; i++) {
+ s += QR[i][k]*QR[i][j];
+ }
+ */
+ s = -s / QR.getQuick(k,k);
+ //QRcolumnsPart[j].assign(QRcolumns[k], F.plusMult(s));
+
+ for (int i = k; i < m; i++) {
+ QR.setQuick(i,j, QR.getQuick(i,j) + s*QR.getQuick(i,k));
+ }
+
+ }
+ }
+ Rdiag.setQuick(k, -nrm);
+ }
+}
+/**
+Returns the Householder vectors <tt>H</tt>.
+@return A lower trapezoidal matrix whose columns define the householder reflections.
+*/
+public DoubleMatrix2D getH () {
+ return Algebra.DEFAULT.trapezoidalLower(QR.copy());
+}
+/**
+Generates and returns the (economy-sized) orthogonal factor <tt>Q</tt>.
+@return <tt>Q</tt>
+*/
+public DoubleMatrix2D getQ () {
+ org.apache.mahout.jet.math.Functions F = org.apache.mahout.jet.math.Functions.functions;
+ DoubleMatrix2D Q = QR.like();
+ //double[][] Q = X.getArray();
+ for (int k = n-1; k >= 0; k--) {
+ DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k);
+ Q.setQuick(k,k, 1);
+ for (int j = k; j < n; j++) {
+ if (QR.getQuick(k,k) != 0) {
+ DoubleMatrix1D Qcolj = Q.viewColumn(j).viewPart(k,m-k);
+ double s = QRcolk.zDotProduct(Qcolj);
+ s = -s / QR.getQuick(k,k);
+ Qcolj.assign(QRcolk, F.plusMult(s));
+ }
+ }
+ }
+ return Q;
+}
+/**
+Returns the upper triangular factor, <tt>R</tt>.
+@return <tt>R</tt>
+*/
+public DoubleMatrix2D getR() {
+ DoubleMatrix2D R = QR.like(n,n);
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ if (i < j)
+ R.setQuick(i,j, QR.getQuick(i,j));
+ else if (i == j)
+ R.setQuick(i,j, Rdiag.getQuick(i));
+ else
+ R.setQuick(i,j, 0);
+ }
+ }
+ return R;
+}
+/**
+Returns whether the matrix <tt>A</tt> has full rank.
+@return true if <tt>R</tt>, and hence <tt>A</tt>, has full rank.
+*/
+public boolean hasFullRank() {
+ for (int j = 0; j < n; j++) {
+ if (Rdiag.getQuick(j) == 0) return false;
+ }
+ return true;
+}
+/**
+Least squares solution of <tt>A*X = B</tt>; <tt>returns X</tt>.
+@param B A matrix with as many rows as <tt>A</tt> and any number of columns.
+@return <tt>X</tt> that minimizes the two norm of <tt>Q*R*X - B</tt>.
+@exception IllegalArgumentException if <tt>B.rows() != A.rows()</tt>.
+@exception IllegalArgumentException if <tt>!this.hasFullRank()</tt> (<tt>A</tt> is rank deficient).
+*/
+public DoubleMatrix2D solve(DoubleMatrix2D B) {
+ org.apache.mahout.jet.math.Functions F = org.apache.mahout.jet.math.Functions.functions;
+ if (B.rows() != m) {
+ throw new IllegalArgumentException("Matrix row dimensions must agree.");
+ }
+ if (!this.hasFullRank()) {
+ throw new IllegalArgumentException("Matrix is rank deficient.");
+ }
+
+ // Copy right hand side
+ int nx = B.columns();
+ DoubleMatrix2D X = B.copy();
+
+ // Compute Y = transpose(Q)*B
+ for (int k = 0; k < n; k++) {
+ for (int j = 0; j < nx; j++) {
+ double s = 0.0;
+ for (int i = k; i < m; i++) {
+ s += QR.getQuick(i,k)*X.getQuick(i,j);
+ }
+ s = -s / QR.getQuick(k,k);
+ for (int i = k; i < m; i++) {
+ X.setQuick(i,j, X.getQuick(i,j) + s*QR.getQuick(i,k));
+ }
+ }
+ }
+ // Solve R*X = Y;
+ for (int k = n-1; k >= 0; k--) {
+ for (int j = 0; j < nx; j++) {
+ X.setQuick(k,j, X.getQuick(k,j) / Rdiag.getQuick(k));
+ }
+ for (int i = 0; i < k; i++) {
+ for (int j = 0; j < nx; j++) {
+ X.setQuick(i,j, X.getQuick(i,j) - X.getQuick(k,j)*QR.getQuick(i,k));
+ }
+ }
+ }
+ return X.viewPart(0,0,n,nx);
+}
+/**
+Returns a String with (propertyName, propertyValue) pairs.
+Useful for debugging or to quickly get the rough picture.
+For example,
+<pre>
+rank : 3
+trace : 0
+</pre>
+*/
+public String toString() {
+ StringBuffer buf = new StringBuffer();
+ String unknown = "Illegal operation or error: ";
+
+ buf.append("-----------------------------------------------------------------\n");
+ buf.append("QRDecomposition(A) --> hasFullRank(A), H, Q, R, pseudo inverse(A)\n");
+ buf.append("-----------------------------------------------------------------\n");
+
+ buf.append("hasFullRank = ");
+ try { buf.append(String.valueOf(this.hasFullRank()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ buf.append("\n\nH = ");
+ try { buf.append(String.valueOf(this.getH()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ buf.append("\n\nQ = ");
+ try { buf.append(String.valueOf(this.getQ()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ buf.append("\n\nR = ");
+ try { buf.append(String.valueOf(this.getR()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ buf.append("\n\npseudo inverse(A) = ");
+ try { buf.append(String.valueOf(this.solve(org.apache.mahout.colt.matrix.DoubleFactory2D.dense.identity(QR.rows()))));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ return buf.toString();
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/QRDecomposition.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SeqBlas.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SeqBlas.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SeqBlas.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SeqBlas.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,233 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.colt.matrix.linalg;
+
+import org.apache.mahout.colt.matrix.DoubleMatrix1D;
+import org.apache.mahout.colt.matrix.DoubleMatrix2D;
+/**
+Sequential implementation of the Basic Linear Algebra System.
+
+@author wolfgang.hoschek@cern.ch
+@version 0.9, 16/04/2000
+*/
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class SeqBlas implements Blas {
+ /**
+ Little trick to allow for "aliasing", that is, renaming this class.
+ Time and again writing code like
+ <p>
+ <tt>SeqBlas.blas.dgemm(...);</tt>
+ <p>
+ is a bit awkward. Using the aliasing you can instead write
+ <p>
+ <tt>Blas B = SeqBlas.blas; <br>
+ B.dgemm(...);</tt>
+ */
+ public static final Blas seqBlas = new SeqBlas();
+
+ private static final org.apache.mahout.jet.math.Functions F = org.apache.mahout.jet.math.Functions.functions;
+/**
+Makes this class non instantiable, but still let's others inherit from it.
+*/
+protected SeqBlas() {}
+public void assign(DoubleMatrix2D A, org.apache.mahout.colt.function.DoubleFunction function) {
+ A.assign(function);
+}
+public void assign(DoubleMatrix2D A, DoubleMatrix2D B, org.apache.mahout.colt.function.DoubleDoubleFunction function) {
+ A.assign(B,function);
+}
+public double dasum(DoubleMatrix1D x) {
+ return x.aggregate(F.plus, F.abs);
+}
+public void daxpy(double alpha, DoubleMatrix1D x, DoubleMatrix1D y) {
+ y.assign(x,F.plusMult(alpha));
+}
+public void daxpy(double alpha, DoubleMatrix2D A, DoubleMatrix2D B) {
+ B.assign(A, F.plusMult(alpha));
+}
+public void dcopy(DoubleMatrix1D x, DoubleMatrix1D y) {
+ y.assign(x);
+}
+public void dcopy(DoubleMatrix2D A, DoubleMatrix2D B) {
+ B.assign(A);
+}
+public double ddot(DoubleMatrix1D x, DoubleMatrix1D y) {
+ return x.zDotProduct(y);
+}
+public void dgemm(boolean transposeA, boolean transposeB, double alpha, DoubleMatrix2D A, DoubleMatrix2D B, double beta, DoubleMatrix2D C) {
+ A.zMult(B,C,alpha,beta,transposeA,transposeB);
+}
+public void dgemv(boolean transposeA, double alpha, DoubleMatrix2D A, DoubleMatrix1D x, double beta, DoubleMatrix1D y) {
+ A.zMult(x,y,alpha,beta,transposeA);
+}
+public void dger(double alpha, DoubleMatrix1D x, DoubleMatrix1D y, DoubleMatrix2D A) {
+ org.apache.mahout.jet.math.PlusMult fun = org.apache.mahout.jet.math.PlusMult.plusMult(0);
+ for (int i=A.rows(); --i >= 0; ) {
+ fun.multiplicator = alpha * x.getQuick(i);
+ A.viewRow(i).assign(y,fun);
+
+ }
+}
+public double dnrm2(DoubleMatrix1D x) {
+ return Math.sqrt(Algebra.DEFAULT.norm2(x));
+}
+public void drot(DoubleMatrix1D x, DoubleMatrix1D y, double c, double s) {
+ x.checkSize(y);
+ DoubleMatrix1D tmp = x.copy();
+
+ x.assign(F.mult(c));
+ x.assign(y,F.plusMult(s));
+
+ y.assign(F.mult(c));
+ y.assign(tmp,F.minusMult(s));
+}
+public void drotg(double a, double b, double rotvec[]) {
+ double c,s,roe,scale,r,z,ra,rb;
+
+ roe = b;
+
+ if (Math.abs(a) > Math.abs(b)) roe = a;
+
+ scale = Math.abs(a) + Math.abs(b);
+
+ if (scale != 0.0) {
+
+ ra = a/scale;
+ rb = b/scale;
+ r = scale*Math.sqrt(ra*ra + rb*rb);
+ r = sign(1.0,roe)*r;
+ c = a/r;
+ s = b/r;
+ z = 1.0;
+ if (Math.abs(a) > Math.abs(b)) z = s;
+ if ((Math.abs(b) >= Math.abs(a)) && (c != 0.0)) z = 1.0/c;
+
+ } else {
+
+ c = 1.0;
+ s = 0.0;
+ r = 0.0;
+ z = 0.0;
+
+ }
+
+ a = r;
+ b = z;
+
+ rotvec[0] = a;
+ rotvec[1] = b;
+ rotvec[2] = c;
+ rotvec[3] = s;
+
+}
+public void dscal(double alpha, DoubleMatrix1D x) {
+ x.assign(F.mult(alpha));
+}
+
+public void dscal(double alpha, DoubleMatrix2D A) {
+ A.assign(F.mult(alpha));
+}
+
+public void dswap(DoubleMatrix1D x, DoubleMatrix1D y) {
+ y.swap(x);
+}
+public void dswap(DoubleMatrix2D A, DoubleMatrix2D B) {
+ //B.swap(A); not yet implemented
+ A.checkShape(B);
+ for(int i = A.rows(); --i >= 0;) A.viewRow(i).swap(B.viewRow(i));
+}
+public void dsymv(boolean isUpperTriangular, double alpha, DoubleMatrix2D A, DoubleMatrix1D x, double beta, DoubleMatrix1D y) {
+ if (isUpperTriangular) A = A.viewDice();
+ Property.DEFAULT.checkSquare(A);
+ int size = A.rows();
+ if (size != x.size() || size!=y.size()) {
+ throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort() + ", " + y.toStringShort());
+ }
+ DoubleMatrix1D tmp = x.like();
+ for (int i = 0; i < size; i++) {
+ double sum = 0;
+ for (int j = 0; j <= i; j++) {
+ sum += A.getQuick(i,j) * x.getQuick(j);
+ }
+ for (int j = i + 1; j < size; j++) {
+ sum += A.getQuick(j,i) * x.getQuick(j);
+ }
+ tmp.setQuick(i, alpha * sum + beta * y.getQuick(i));
+ }
+ y.assign(tmp);
+}
+public void dtrmv(boolean isUpperTriangular, boolean transposeA, boolean isUnitTriangular, DoubleMatrix2D A, DoubleMatrix1D x) {
+ if (transposeA) {
+ A = A.viewDice();
+ isUpperTriangular = !isUpperTriangular;
+ }
+
+ Property.DEFAULT.checkSquare(A);
+ int size = A.rows();
+ if (size != x.size()) {
+ throw new IllegalArgumentException(A.toStringShort() + ", " + x.toStringShort());
+ }
+
+ DoubleMatrix1D b = x.like();
+ DoubleMatrix1D y = x.like();
+ if (isUnitTriangular) {
+ y.assign(1);
+ }
+ else {
+ for (int i = 0; i < size; i++) {
+ y.setQuick(i, A.getQuick(i,i));
+ }
+ }
+
+ for (int i = 0; i < size; i++) {
+ double sum = 0;
+ if (!isUpperTriangular) {
+ for (int j = 0; j < i; j++) {
+ sum += A.getQuick(i,j) * x.getQuick(j);
+ }
+ sum += y.getQuick(i) * x.getQuick(i);
+ }
+ else {
+ sum += y.getQuick(i) * x.getQuick(i);
+ for (int j = i + 1; j < size; j++) {
+ sum += A.getQuick(i,j) * x.getQuick(j); }
+ }
+ b.setQuick(i,sum);
+ }
+ x.assign(b);
+}
+public int idamax(DoubleMatrix1D x) {
+ int maxIndex = -1;
+ double maxValue = Double.MIN_VALUE;
+ for (int i=x.size(); --i >= 0; ) {
+ double v = Math.abs(x.getQuick(i));
+ if (v > maxValue) {
+ maxValue = v;
+ maxIndex = i;
+ }
+ }
+ return maxIndex;
+}
+/**
+Implements the FORTRAN sign (not sin) function.
+See the code for details.
+@param a a
+@param b b
+*/
+private double sign(double a, double b) {
+ if (b < 0.0) {
+ return -Math.abs(a);
+ } else {
+ return Math.abs(a);
+ }
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SeqBlas.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SingularValueDecomposition.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SingularValueDecomposition.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SingularValueDecomposition.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SingularValueDecomposition.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,573 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.colt.matrix.linalg;
+
+import org.apache.mahout.colt.matrix.DoubleFactory2D;
+import org.apache.mahout.colt.matrix.DoubleMatrix2D;
+/**
+For an <tt>m x n</tt> matrix <tt>A</tt> with <tt>m >= n</tt>, the singular value decomposition is
+an <tt>m x n</tt> orthogonal matrix <tt>U</tt>, an <tt>n x n</tt> diagonal matrix <tt>S</tt>, and
+an <tt>n x n</tt> orthogonal matrix <tt>V</tt> so that <tt>A = U*S*V'</tt>.
+<P>
+The singular values, <tt>sigma[k] = S[k][k]</tt>, are ordered so that
+<tt>sigma[0] >= sigma[1] >= ... >= sigma[n-1]</tt>.
+<P>
+The singular value decomposition always exists, so the constructor will
+never fail. The matrix condition number and the effective numerical
+rank can be computed from this decomposition.
+*/
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class SingularValueDecomposition implements java.io.Serializable {
+ static final long serialVersionUID = 1020;
+ /** Arrays for internal storage of U and V.
+ @serial internal storage of U.
+ @serial internal storage of V.
+ */
+ private double[][] U, V;
+
+ /** Array for internal storage of singular values.
+ @serial internal storage of singular values.
+ */
+ private double[] s;
+
+ /** Row and column dimensions.
+ @serial row dimension.
+ @serial column dimension.
+ */
+ private int m, n;
+/**
+Constructs and returns a new singular value decomposition object;
+The decomposed matrices can be retrieved via instance methods of the returned decomposition object.
+@param Arg A rectangular matrix.
+@return A decomposition object to access <tt>U</tt>, <tt>S</tt> and <tt>V</tt>.
+@throws IllegalArgumentException if <tt>A.rows() < A.columns()</tt>.
+*/
+public SingularValueDecomposition(DoubleMatrix2D Arg) {
+ Property.DEFAULT.checkRectangular(Arg);
+
+ // Derived from LINPACK code.
+ // Initialize.
+ double[][] A = Arg.toArray();
+ m = Arg.rows();
+ n = Arg.columns();
+ int nu = Math.min(m,n);
+ s = new double [Math.min(m+1,n)];
+ U = new double [m][nu];
+ V = new double [n][n];
+ double[] e = new double [n];
+ double[] work = new double [m];
+ boolean wantu = true;
+ boolean wantv = true;
+
+ // Reduce A to bidiagonal form, storing the diagonal elements
+ // in s and the super-diagonal elements in e.
+
+ int nct = Math.min(m-1,n);
+ int nrt = Math.max(0,Math.min(n-2,m));
+ for (int k = 0; k < Math.max(nct,nrt); k++) {
+ if (k < nct) {
+
+ // Compute the transformation for the k-th column and
+ // place the k-th diagonal in s[k].
+ // Compute 2-norm of k-th column without under/overflow.
+ s[k] = 0;
+ for (int i = k; i < m; i++) {
+ s[k] = Algebra.hypot(s[k],A[i][k]);
+ }
+ if (s[k] != 0.0) {
+ if (A[k][k] < 0.0) {
+ s[k] = -s[k];
+ }
+ for (int i = k; i < m; i++) {
+ A[i][k] /= s[k];
+ }
+ A[k][k] += 1.0;
+ }
+ s[k] = -s[k];
+ }
+ for (int j = k+1; j < n; j++) {
+ if ((k < nct) & (s[k] != 0.0)) {
+
+ // Apply the transformation.
+
+ double t = 0;
+ for (int i = k; i < m; i++) {
+ t += A[i][k]*A[i][j];
+ }
+ t = -t/A[k][k];
+ for (int i = k; i < m; i++) {
+ A[i][j] += t*A[i][k];
+ }
+ }
+
+ // Place the k-th row of A into e for the
+ // subsequent calculation of the row transformation.
+
+ e[j] = A[k][j];
+ }
+ if (wantu & (k < nct)) {
+
+ // Place the transformation in U for subsequent back
+ // multiplication.
+
+ for (int i = k; i < m; i++) {
+ U[i][k] = A[i][k];
+ }
+ }
+ if (k < nrt) {
+
+ // Compute the k-th row transformation and place the
+ // k-th super-diagonal in e[k].
+ // Compute 2-norm without under/overflow.
+ e[k] = 0;
+ for (int i = k+1; i < n; i++) {
+ e[k] = Algebra.hypot(e[k],e[i]);
+ }
+ if (e[k] != 0.0) {
+ if (e[k+1] < 0.0) {
+ e[k] = -e[k];
+ }
+ for (int i = k+1; i < n; i++) {
+ e[i] /= e[k];
+ }
+ e[k+1] += 1.0;
+ }
+ e[k] = -e[k];
+ if ((k+1 < m) & (e[k] != 0.0)) {
+
+ // Apply the transformation.
+
+ for (int i = k+1; i < m; i++) {
+ work[i] = 0.0;
+ }
+ for (int j = k+1; j < n; j++) {
+ for (int i = k+1; i < m; i++) {
+ work[i] += e[j]*A[i][j];
+ }
+ }
+ for (int j = k+1; j < n; j++) {
+ double t = -e[j]/e[k+1];
+ for (int i = k+1; i < m; i++) {
+ A[i][j] += t*work[i];
+ }
+ }
+ }
+ if (wantv) {
+
+ // Place the transformation in V for subsequent
+ // back multiplication.
+
+ for (int i = k+1; i < n; i++) {
+ V[i][k] = e[i];
+ }
+ }
+ }
+ }
+
+ // Set up the final bidiagonal matrix or order p.
+
+ int p = Math.min(n,m+1);
+ if (nct < n) {
+ s[nct] = A[nct][nct];
+ }
+ if (m < p) {
+ s[p-1] = 0.0;
+ }
+ if (nrt+1 < p) {
+ e[nrt] = A[nrt][p-1];
+ }
+ e[p-1] = 0.0;
+
+ // If required, generate U.
+
+ if (wantu) {
+ for (int j = nct; j < nu; j++) {
+ for (int i = 0; i < m; i++) {
+ U[i][j] = 0.0;
+ }
+ U[j][j] = 1.0;
+ }
+ for (int k = nct-1; k >= 0; k--) {
+ if (s[k] != 0.0) {
+ for (int j = k+1; j < nu; j++) {
+ double t = 0;
+ for (int i = k; i < m; i++) {
+ t += U[i][k]*U[i][j];
+ }
+ t = -t/U[k][k];
+ for (int i = k; i < m; i++) {
+ U[i][j] += t*U[i][k];
+ }
+ }
+ for (int i = k; i < m; i++ ) {
+ U[i][k] = -U[i][k];
+ }
+ U[k][k] = 1.0 + U[k][k];
+ for (int i = 0; i < k-1; i++) {
+ U[i][k] = 0.0;
+ }
+ } else {
+ for (int i = 0; i < m; i++) {
+ U[i][k] = 0.0;
+ }
+ U[k][k] = 1.0;
+ }
+ }
+ }
+
+ // If required, generate V.
+
+ if (wantv) {
+ for (int k = n-1; k >= 0; k--) {
+ if ((k < nrt) & (e[k] != 0.0)) {
+ for (int j = k+1; j < nu; j++) {
+ double t = 0;
+ for (int i = k+1; i < n; i++) {
+ t += V[i][k]*V[i][j];
+ }
+ t = -t/V[k+1][k];
+ for (int i = k+1; i < n; i++) {
+ V[i][j] += t*V[i][k];
+ }
+ }
+ }
+ for (int i = 0; i < n; i++) {
+ V[i][k] = 0.0;
+ }
+ V[k][k] = 1.0;
+ }
+ }
+
+ // Main iteration loop for the singular values.
+
+ int pp = p-1;
+ int iter = 0;
+ double eps = Math.pow(2.0,-52.0);
+ while (p > 0) {
+ int k,kase;
+
+ // Here is where a test for too many iterations would go.
+
+ // This section of the program inspects for
+ // negligible elements in the s and e arrays. On
+ // completion the variables kase and k are set as follows.
+
+ // kase = 1 if s(p) and e[k-1] are negligible and k<p
+ // kase = 2 if s(k) is negligible and k<p
+ // kase = 3 if e[k-1] is negligible, k<p, and
+ // s(k), ..., s(p) are not negligible (qr step).
+ // kase = 4 if e(p-1) is negligible (convergence).
+
+ for (k = p-2; k >= -1; k--) {
+ if (k == -1) {
+ break;
+ }
+ if (Math.abs(e[k]) <= eps*(Math.abs(s[k]) + Math.abs(s[k+1]))) {
+ e[k] = 0.0;
+ break;
+ }
+ }
+ if (k == p-2) {
+ kase = 4;
+ } else {
+ int ks;
+ for (ks = p-1; ks >= k; ks--) {
+ if (ks == k) {
+ break;
+ }
+ double t = (ks != p ? Math.abs(e[ks]) : 0.) +
+ (ks != k+1 ? Math.abs(e[ks-1]) : 0.);
+ if (Math.abs(s[ks]) <= eps*t) {
+ s[ks] = 0.0;
+ break;
+ }
+ }
+ if (ks == k) {
+ kase = 3;
+ } else if (ks == p-1) {
+ kase = 1;
+ } else {
+ kase = 2;
+ k = ks;
+ }
+ }
+ k++;
+
+ // Perform the task indicated by kase.
+
+ switch (kase) {
+
+ // Deflate negligible s(p).
+
+ case 1: {
+ double f = e[p-2];
+ e[p-2] = 0.0;
+ for (int j = p-2; j >= k; j--) {
+ double t = Algebra.hypot(s[j],f);
+ double cs = s[j]/t;
+ double sn = f/t;
+ s[j] = t;
+ if (j != k) {
+ f = -sn*e[j-1];
+ e[j-1] = cs*e[j-1];
+ }
+ if (wantv) {
+ for (int i = 0; i < n; i++) {
+ t = cs*V[i][j] + sn*V[i][p-1];
+ V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1];
+ V[i][j] = t;
+ }
+ }
+ }
+ }
+ break;
+
+ // Split at negligible s(k).
+
+ case 2: {
+ double f = e[k-1];
+ e[k-1] = 0.0;
+ for (int j = k; j < p; j++) {
+ double t = Algebra.hypot(s[j],f);
+ double cs = s[j]/t;
+ double sn = f/t;
+ s[j] = t;
+ f = -sn*e[j];
+ e[j] = cs*e[j];
+ if (wantu) {
+ for (int i = 0; i < m; i++) {
+ t = cs*U[i][j] + sn*U[i][k-1];
+ U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1];
+ U[i][j] = t;
+ }
+ }
+ }
+ }
+ break;
+
+ // Perform one qr step.
+
+ case 3: {
+
+ // Calculate the shift.
+
+ double scale = Math.max(Math.max(Math.max(Math.max(
+ Math.abs(s[p-1]),Math.abs(s[p-2])),Math.abs(e[p-2])),
+ Math.abs(s[k])),Math.abs(e[k]));
+ double sp = s[p-1]/scale;
+ double spm1 = s[p-2]/scale;
+ double epm1 = e[p-2]/scale;
+ double sk = s[k]/scale;
+ double ek = e[k]/scale;
+ double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
+ double c = (sp*epm1)*(sp*epm1);
+ double shift = 0.0;
+ if ((b != 0.0) | (c != 0.0)) {
+ shift = Math.sqrt(b*b + c);
+ if (b < 0.0) {
+ shift = -shift;
+ }
+ shift = c/(b + shift);
+ }
+ double f = (sk + sp)*(sk - sp) + shift;
+ double g = sk*ek;
+
+ // Chase zeros.
+
+ for (int j = k; j < p-1; j++) {
+ double t = Algebra.hypot(f,g);
+ double cs = f/t;
+ double sn = g/t;
+ if (j != k) {
+ e[j-1] = t;
+ }
+ f = cs*s[j] + sn*e[j];
+ e[j] = cs*e[j] - sn*s[j];
+ g = sn*s[j+1];
+ s[j+1] = cs*s[j+1];
+ if (wantv) {
+ for (int i = 0; i < n; i++) {
+ t = cs*V[i][j] + sn*V[i][j+1];
+ V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1];
+ V[i][j] = t;
+ }
+ }
+ t = Algebra.hypot(f,g);
+ cs = f/t;
+ sn = g/t;
+ s[j] = t;
+ f = cs*e[j] + sn*s[j+1];
+ s[j+1] = -sn*e[j] + cs*s[j+1];
+ g = sn*e[j+1];
+ e[j+1] = cs*e[j+1];
+ if (wantu && (j < m-1)) {
+ for (int i = 0; i < m; i++) {
+ t = cs*U[i][j] + sn*U[i][j+1];
+ U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1];
+ U[i][j] = t;
+ }
+ }
+ }
+ e[p-2] = f;
+ iter = iter + 1;
+ }
+ break;
+
+ // Convergence.
+
+ case 4: {
+
+ // Make the singular values positive.
+
+ if (s[k] <= 0.0) {
+ s[k] = (s[k] < 0.0 ? -s[k] : 0.0);
+ if (wantv) {
+ for (int i = 0; i <= pp; i++) {
+ V[i][k] = -V[i][k];
+ }
+ }
+ }
+
+ // Order the singular values.
+
+ while (k < pp) {
+ if (s[k] >= s[k+1]) {
+ break;
+ }
+ double t = s[k];
+ s[k] = s[k+1];
+ s[k+1] = t;
+ if (wantv && (k < n-1)) {
+ for (int i = 0; i < n; i++) {
+ t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t;
+ }
+ }
+ if (wantu && (k < m-1)) {
+ for (int i = 0; i < m; i++) {
+ t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t;
+ }
+ }
+ k++;
+ }
+ iter = 0;
+ p--;
+ }
+ break;
+ }
+ }
+ }
+/**
+Returns the two norm condition number, which is <tt>max(S) / min(S)</tt>.
+*/
+public double cond() {
+ return s[0]/s[Math.min(m,n)-1];
+}
+/**
+Returns the diagonal matrix of singular values.
+@return S
+*/
+public DoubleMatrix2D getS() {
+ double[][] S = new double[n][n];
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ S[i][j] = 0.0;
+ }
+ S[i][i] = this.s[i];
+ }
+ return DoubleFactory2D.dense.make(S);
+}
+/**
+Returns the diagonal of <tt>S</tt>, which is a one-dimensional array of singular values
+@return diagonal of <tt>S</tt>.
+*/
+public double[] getSingularValues() {
+ return s;
+}
+/**
+Returns the left singular vectors <tt>U</tt>.
+@return <tt>U</tt>
+*/
+public DoubleMatrix2D getU() {
+ //return new DoubleMatrix2D(U,m,Math.min(m+1,n));
+ return DoubleFactory2D.dense.make(U).viewPart(0,0,m,Math.min(m+1,n));
+}
+/**
+Returns the right singular vectors <tt>V</tt>.
+@return <tt>V</tt>
+*/
+public DoubleMatrix2D getV() {
+ return DoubleFactory2D.dense.make(V);
+}
+/**
+Returns the two norm, which is <tt>max(S)</tt>.
+*/
+public double norm2() {
+ return s[0];
+}
+/**
+Returns the effective numerical matrix rank, which is the number of nonnegligible singular values.
+*/
+public int rank() {
+ double eps = Math.pow(2.0,-52.0);
+ double tol = Math.max(m,n)*s[0]*eps;
+ int r = 0;
+ for (int i = 0; i < s.length; i++) {
+ if (s[i] > tol) {
+ r++;
+ }
+ }
+ return r;
+}
+/**
+Returns a String with (propertyName, propertyValue) pairs.
+Useful for debugging or to quickly get the rough picture.
+For example,
+<pre>
+rank : 3
+trace : 0
+</pre>
+*/
+public String toString() {
+ StringBuffer buf = new StringBuffer();
+ String unknown = "Illegal operation or error: ";
+
+ buf.append("---------------------------------------------------------------------\n");
+ buf.append("SingularValueDecomposition(A) --> cond(A), rank(A), norm2(A), U, S, V\n");
+ buf.append("---------------------------------------------------------------------\n");
+
+ buf.append("cond = ");
+ try { buf.append(String.valueOf(this.cond()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ buf.append("\nrank = ");
+ try { buf.append(String.valueOf(this.rank()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ buf.append("\nnorm2 = ");
+ try { buf.append(String.valueOf(this.norm2()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ buf.append("\n\nU = ");
+ try { buf.append(String.valueOf(this.getU()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ buf.append("\n\nS = ");
+ try { buf.append(String.valueOf(this.getS()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ buf.append("\n\nV = ");
+ try { buf.append(String.valueOf(this.getV()));}
+ catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
+
+ return buf.toString();
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/SingularValueDecomposition.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Smp.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Smp.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Smp.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Smp.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,195 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.colt.matrix.linalg;
+
+import EDU.oswego.cs.dl.util.concurrent.FJTask;
+import EDU.oswego.cs.dl.util.concurrent.FJTaskRunnerGroup;
+import org.apache.mahout.colt.matrix.DoubleMatrix2D;
+/*
+*/
+class Smp {
+ protected FJTaskRunnerGroup taskGroup; // a very efficient and light weight thread pool
+
+ protected int maxThreads;
+/**
+Constructs a new Smp using a maximum of <tt>maxThreads<tt> threads.
+*/
+protected Smp(int maxThreads) {
+ maxThreads = Math.max(1,maxThreads);
+ this.maxThreads = maxThreads;
+ if (maxThreads>1) {
+ this.taskGroup = new FJTaskRunnerGroup(maxThreads);
+ }
+ else { // avoid parallel overhead
+ this.taskGroup = null;
+ }
+}
+/**
+ * Clean up deamon threads, if necessary.
+ */
+public void finalize() {
+ if (this.taskGroup!=null) this.taskGroup.interruptAll();
+}
+protected void run(final DoubleMatrix2D[] blocksA, final DoubleMatrix2D[] blocksB, final double[] results, final Matrix2DMatrix2DFunction function) {
+ final FJTask[] subTasks = new FJTask[blocksA.length];
+ for (int i=0; i<blocksA.length; i++) {
+ final int k = i;
+ subTasks[i] = new FJTask() {
+ public void run() {
+ double result = function.apply(blocksA[k],blocksB != null ? blocksB[k] : null);
+ if (results!=null) results[k] = result;
+ //System.out.print(".");
+ }
+ };
+ }
+
+ // run tasks and wait for completion
+ try {
+ this.taskGroup.invoke(
+ new FJTask() {
+ public void run() {
+ coInvoke(subTasks);
+ }
+ }
+ );
+ } catch (InterruptedException exc) {}
+}
+protected DoubleMatrix2D[] splitBlockedNN(DoubleMatrix2D A, int threshold, long flops) {
+ /*
+ determine how to split and parallelize best into blocks
+ if more B.columns than tasks --> split B.columns, as follows:
+
+ xx|xx|xxx B
+ xx|xx|xxx
+ xx|xx|xxx
+ A
+ xxx xx|xx|xxx C
+ xxx xx|xx|xxx
+ xxx xx|xx|xxx
+ xxx xx|xx|xxx
+ xxx xx|xx|xxx
+
+ if less B.columns than tasks --> split A.rows, as follows:
+
+ xxxxxxx B
+ xxxxxxx
+ xxxxxxx
+ A
+ xxx xxxxxxx C
+ xxx xxxxxxx
+ --- -------
+ xxx xxxxxxx
+ xxx xxxxxxx
+ --- -------
+ xxx xxxxxxx
+
+ */
+ //long flops = 2L*A.rows()*A.columns()*A.columns();
+ int noOfTasks = (int) Math.min(flops / threshold, this.maxThreads); // each thread should process at least 30000 flops
+ boolean splitHoriz = (A.columns() < noOfTasks);
+ //boolean splitHoriz = (A.columns() >= noOfTasks);
+ int p = splitHoriz ? A.rows() : A.columns();
+ noOfTasks = Math.min(p,noOfTasks);
+
+ if (noOfTasks < 2) { // parallelization doesn't pay off (too much start up overhead)
+ return null;
+ }
+
+ // set up concurrent tasks
+ int span = p/noOfTasks;
+ final DoubleMatrix2D[] blocks = new DoubleMatrix2D[noOfTasks];
+ for (int i=0; i<noOfTasks; i++) {
+ final int offset = i*span;
+ if (i==noOfTasks-1) span = p - span*i; // last span may be a bit larger
+
+ final DoubleMatrix2D AA,BB,CC;
+ if (!splitHoriz) { // split B along columns into blocks
+ blocks[i] = A.viewPart(0,offset, A.rows(), span);
+ }
+ else { // split A along rows into blocks
+ blocks[i] = A.viewPart(offset,0,span,A.columns());
+ }
+ }
+ return blocks;
+}
+protected DoubleMatrix2D[][] splitBlockedNN(DoubleMatrix2D A, DoubleMatrix2D B, int threshold, long flops) {
+ DoubleMatrix2D[] blocksA = splitBlockedNN(A,threshold, flops);
+ if (blocksA==null) return null;
+ DoubleMatrix2D[] blocksB = splitBlockedNN(B,threshold, flops);
+ if (blocksB==null) return null;
+ DoubleMatrix2D[][] blocks = {blocksA,blocksB};
+ return blocks;
+}
+protected DoubleMatrix2D[] splitStridedNN(DoubleMatrix2D A, int threshold, long flops) {
+ /*
+ determine how to split and parallelize best into blocks
+ if more B.columns than tasks --> split B.columns, as follows:
+
+ xx|xx|xxx B
+ xx|xx|xxx
+ xx|xx|xxx
+ A
+ xxx xx|xx|xxx C
+ xxx xx|xx|xxx
+ xxx xx|xx|xxx
+ xxx xx|xx|xxx
+ xxx xx|xx|xxx
+
+ if less B.columns than tasks --> split A.rows, as follows:
+
+ xxxxxxx B
+ xxxxxxx
+ xxxxxxx
+ A
+ xxx xxxxxxx C
+ xxx xxxxxxx
+ --- -------
+ xxx xxxxxxx
+ xxx xxxxxxx
+ --- -------
+ xxx xxxxxxx
+
+ */
+ //long flops = 2L*A.rows()*A.columns()*A.columns();
+ int noOfTasks = (int) Math.min(flops / threshold, this.maxThreads); // each thread should process at least 30000 flops
+ boolean splitHoriz = (A.columns() < noOfTasks);
+ //boolean splitHoriz = (A.columns() >= noOfTasks);
+ int p = splitHoriz ? A.rows() : A.columns();
+ noOfTasks = Math.min(p,noOfTasks);
+
+ if (noOfTasks < 2) { // parallelization doesn't pay off (too much start up overhead)
+ return null;
+ }
+
+ // set up concurrent tasks
+ int span = p/noOfTasks;
+ final DoubleMatrix2D[] blocks = new DoubleMatrix2D[noOfTasks];
+ for (int i=0; i<noOfTasks; i++) {
+ final int offset = i*span;
+ if (i==noOfTasks-1) span = p - span*i; // last span may be a bit larger
+
+ final DoubleMatrix2D AA,BB,CC;
+ if (!splitHoriz) {
+ // split B along columns into blocks
+ blocks[i] = A.viewPart(0,i,A.rows(),A.columns()-i).viewStrides(1,noOfTasks);
+ }
+ else {
+ // split A along rows into blocks
+ blocks[i] = A.viewPart(i,0,A.rows()-i,A.columns()).viewStrides(noOfTasks,1);
+ }
+ }
+ return blocks;
+}
+/**
+ * Prints various snapshot statistics to System.out; Simply delegates to {@link EDU.oswego.cs.dl.util.concurrent.FJTaskRunnerGroup#stats}.
+ */
+public void stats() {
+ if (this.taskGroup!=null) this.taskGroup.stats();
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/matrix/matrix/linalg/Smp.java
------------------------------------------------------------------------------
svn:eol-style = native