denselutest.java

另一个功能更强大的矩阵运算软件开源代码

/*
 * Copyright (C) 2003-2006 Bjørn-Ove Heimsund
 * 
 * This file is part of MTJ.
 * 
 * This library is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation; either version 2.1 of the License, or (at your
 * option) any later version.
 * 
 * This library 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 Lesser General Public License
 * for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation,
 * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 */

package no.uib.cipr.matrix;

import no.uib.cipr.matrix.DenseLU;
import no.uib.cipr.matrix.DenseMatrix;
import no.uib.cipr.matrix.Matrices;
import no.uib.cipr.matrix.Matrix;
import no.uib.cipr.matrix.MatrixEntry;
import junit.framework.TestCase;

/**
 * Tests the dense LU decomposition
 */
public class DenseLUTest extends TestCase {

	/**
	 * Matrix to decompose
	 */
	private DenseMatrix A;

	private DenseMatrix I;

	private final int max = 100;

	public DenseLUTest(String arg0) {
		super(arg0);
	}

	@Override
	protected void setUp() throws Exception {
		int n = Utilities.getInt(1, max);
		A = new DenseMatrix(n, n);
		Utilities.populate(A);
		Utilities.addDiagonal(A, 1);
		while (Utilities.singular(A))
			Utilities.addDiagonal(A, 1);

		I = Matrices.identity(n);
	}

	@Override
	protected void tearDown() throws Exception {
		A = null;
		I = null;
	}

	public void testDenseLU() {
		int n = A.numRows();
		DenseLU lu = new DenseLU(n, n);
		lu.factor(A.copy());

		lu.solve(I);

		Matrix J = I.mult(A, new DenseMatrix(n, n));
		for (int i = 0; i < n; ++i)
			for (int j = 0; j < n; ++j)
				if (i != j)
					assertEquals(J.get(i, j), 0, 1e-10);
				else
					assertEquals(J.get(i, j), 1, 1e-10);
	}

	public void testDenseLUtranspose() {
		int n = A.numRows();
		DenseLU lu = new DenseLU(n, n);
		lu.factor(A.copy());

		lu.transSolve(I);

		Matrix J = I.transAmult(A, new DenseMatrix(n, n));
		for (int i = 0; i < n; ++i)
			for (int j = 0; j < n; ++j)
				if (i != j)
					assertEquals(J.get(i, j), 0, 1e-10);
				else
					assertEquals(J.get(i, j), 1, 1e-10);
	}

	public void testDenseLUrcond() {
		int n = A.numRows();
		DenseLU lu = new DenseLU(n, n);
		lu.factor(A.copy());

		lu.rcond(A, Matrix.Norm.One);
		lu.rcond(A, Matrix.Norm.Infinity);
	}

	public void testDenseLUToInput() {
		// MTJ bug in DenseLU code

		Matrix m = new DenseMatrix(3, 3);

		// -2 2 -3
		// -1 1 3
		// 2 0 -1
		m.set(0, 0, -2);
		m.set(0, 1, 2);
		m.set(0, 2, -3);
		m.set(1, 0, -1);
		m.set(1, 1, 1);
		m.set(1, 2, 3);
		m.set(2, 0, 2);
		m.set(2, 1, 0);
		m.set(2, 2, -1);

		// SHOULD BE:
		// L:
		// 1.000 0.000 0.000
		// -1.000 1.000 0.000
		// 0.500 0.000 1.000
		//
		// U:
		// -2.000 2.000 -3.000
		// 0.000 2.000 -4.000
		// 0.000 0.000 4.500
		//
		// Permutation matrix:
		// 1.000 0.000 0.000
		// 0.000 0.000 1.000
		// 0.000 1.000 0.000

		DenseLU dlu = DenseLU.factorize(m);

		// check that m = L . U
		Matrix lTimesU = new DenseMatrix(3, 3);
		dlu.getL().mult(dlu.getU(), lTimesU);
		int[] pivots = dlu.getPivots();
		for (MatrixEntry entry : m) {
			int row = entry.row();
			int col = entry.column();
			double val = entry.get();
			double valLU = pivots[row] * lTimesU.get(row, col);
			assert val == valLU : "Row " + row + ", Col " + col
					+ " wasn't equal! " + val + " " + valLU;
		}

		Matrix lu = dlu.getLU();
		// m == lu
		for (MatrixEntry entry : m) {
			int row = entry.row();
			int col = entry.column();
			double val = entry.get();
			double valLU = lu.get(row, col);
			assert val == valLU : "Row " + row + ", Col " + col
					+ " wasn't equal! " + val + " " + valLU;
		}
	}

}