Any ideas on where this stuff can go???
-
+
+ /** TODO: Find a home for the following methods in the linear package */
+
+ /**
+ * <p>Uses back substitution to solve the system</p>
+ *
+ * <p>coefficients X = constants</p>
+ *
+ * <p>coefficients must upper-triangular and constants must be a column
+ * matrix. The solution is returned as a column matrix.</p>
+ *
+ * <p>The number of columns in coefficients determines the length
+ * of the returned solution vector (column matrix). If constants
+ * has more rows than coefficients has columns, excess rows are ignored.
+ * Similarly, extra (zero) rows in coefficients are ignored</p>
+ *
+ * @param coefficients upper-triangular coefficients matrix
+ * @param constants column RHS constants matrix
+ * @return solution matrix as a column matrix
+ *
+ */
+ private static RealMatrix solveUpperTriangular(RealMatrixImpl coefficients,
+ RealMatrixImpl constants) {
+ if (!isUpperTriangular(coefficients, 1E-12)) {
+ throw new IllegalArgumentException(
+ "Coefficients is not upper-triangular");
+ }
+ if (constants.getColumnDimension() != 1) {
+ throw new IllegalArgumentException(
+ "Constants not a column matrix.");
+ }
+ int length = coefficients.getColumnDimension();
+ double[][] cons = constants.getDataRef();
+ double[][] coef = coefficients.getDataRef();
+ double x[] = new double[length];
+ for (int i = 0; i < length; i++) {
+ int index = length - 1 - i;
+ double sum = 0;
+ for (int j = index + 1; j < length; j++) {
+ sum += coef[index][j] * x[j];
+ }
+ x[index] = (cons[index][0] - sum) / coef[index][index];
+ }
+ return new RealMatrixImpl(x);
+ }
+
+ /**
+ * <p>Returns true iff m is an upper-triangular matrix.</p>
+ *
+ * <p>Makes sure all below-diagonal elements are within epsilon of 0.</p>
+ *
+ * @param m matrix to check
+ * @param epsilon maximum allowable absolute value for elements below
+ * the main diagonal
+ *
+ * @return true if m is upper-triangular; false otherwise
+ * @throws NullPointerException if m is null
+ */
+ private static boolean isUpperTriangular(RealMatrixImpl m, double epsilon)
{
+ double[][] data = m.getDataRef();
+ int nCols = m.getColumnDimension();
+ int nRows = m.getRowDimension();
+ for (int r = 0; r < nRows; r++) {
+ int bound = Math.min(r, nCols);
+ for (int c = 0; c < bound; c++) {
+ if (Math.abs(data[r][c]) > epsilon) {
+ return false;
+ }
+ }
+ }
+ return true;
+ }
}
---------------------------------------------------------------------
To unsubscribe, e-mail: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
