From ab78365644f52a0330b430bf314accdb77e31d42 Mon Sep 17 00:00:00 2001
From: Daniel Weschke <daniel.weschke@gmail.com>
Date: Fri, 14 Mar 2014 22:10:43 +0100
Subject: [PATCH] Deploy matrix

---
 src/math/matrix/Matrix.java | 2603 +++++++++++++++++++++++++++++++++++
 1 file changed, 2603 insertions(+)
 create mode 100644 src/math/matrix/Matrix.java

diff --git a/src/math/matrix/Matrix.java b/src/math/matrix/Matrix.java
new file mode 100644
index 0000000..e45cd53
--- /dev/null
+++ b/src/math/matrix/Matrix.java
@@ -0,0 +1,2603 @@
+package math.matrix;
+
+import java.io.BufferedReader;
+import java.io.IOException;
+import java.io.PrintWriter;
+import java.io.Serializable;
+import java.io.StreamTokenizer;
+import java.text.DecimalFormat;
+import java.text.DecimalFormatSymbols;
+import java.text.NumberFormat;
+import java.util.Locale;
+import java.util.Scanner;
+
+import math.Maths;
+import string.WString;
+import thisandthat.WObject;
+import exception.ComplexException;
+import exception.IllegalDimensionException;
+import exception.NoSquareException;
+import exception.SingularException;
+
+/**
+ * Matrix object M&isin;&#x211D;
+ * with m rows and n columns
+ * @author Daniel Weschke
+ * @version 01 November 2013
+ */
+public class Matrix extends WObject implements Cloneable, Serializable{
+	
+	enum Flag{Square, Symmetric} // TODO: subclass?
+	
+	/**
+	 * UID
+	 */
+	private static final long serialVersionUID = 6314556279121894887L;
+
+	/**
+	 * Row dimension.
+	 * @serial m rows
+	 */
+	protected int m;
+	
+	/**
+	 * Column dimension.
+	 * @serial n columns
+	 */
+	protected int n;
+	
+	/**
+	 * @serial data array for internal storage of elements 
+	 */
+//	private double[][] data;
+	protected double[] data;
+	
+	/**
+	 * @serial solution of eigenvalue decomposition
+	 */
+	private EigenvalueDecomposition eig;
+
+	/**
+	 * Create null-object matrix
+	 */
+	public Matrix(){
+	}
+	
+	/**
+	 * Create an m-by-n matrix of zeros.
+	 * @param m number of rows
+	 * @param n number of columns
+	 */
+	public Matrix(int m, int n){
+		this.m = m;
+		this.n = n;
+//		data = new double[m][n];
+		data = new double[m*n];
+	}
+	
+	/**
+	 * Create an m-by-n matrix of scalar value.
+	 * @param m number of rows
+	 * @param n number of columns
+	 * @param s scalar value
+	 */
+	public Matrix(int m, int n, double s){
+		this(m, n);
+		int i,j;
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				set(i, j, s);
+	}
+	
+	/**
+	 * Create a 1-by-1 matrix based on a scalar.
+	 * @param s scalar
+	 */
+	public Matrix(double s){
+		this(1,1);
+		set(0,0,s);
+	}
+	
+	/**
+	 * Create an m-by-1 matrix based on an array.
+	 * @param a double array
+	 */
+	public Matrix(double[] a){
+		this(a.length,1);
+		int i;
+		for(i=0; i<m; i++)
+			set(i, 0, a[i]);
+	}
+	
+	/**
+	 * Create an m-by-2 matrix based two arrays.
+	 * @param a double array
+	 * @param b double array
+	 */
+	public Matrix(double[] a, double[] b){
+		this(a.length>b.length ? a.length : b.length, 2);
+		int i;
+		for(i=0; i<a.length; i++)
+			set(i, 0, a[i]);
+		for(i=0; i<b.length; i++)
+			set(i, 1, b[i]);
+	}
+
+	/** Create an m-by-n matrix from a one-dimensional packed array.
+	 * {a<sub>1</sub>, a<sub>2</sub>, ..., a<sub>mn</sub>} →
+	 * {{a<sub>1</sub>, a<sub>2</sub>, ..., a<sub>1n</sub>}, 
+	 * {a<sub>1n+1</sub>, a<sub>1n+2</sub>, ..., a<sub>2n</sub>}
+	 * @param n number of columns.
+	 * @param a one-dimensional array of doubles, packed by rows.
+	 * @exception IllegalArgumentException Array length must be a multiple of m.
+	 */
+	public Matrix(int n, double[] a){
+		this((n!=0 ? a.length/n : 0),n);
+		if(m*n != a.length)
+			throw new IllegalArgumentException("Array length must be a multiple of n.");
+		int i,j;
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				set(i, j, a[i*n+j]);
+	}
+
+	/** Create an m-by-n matrix from a one-dimensional packed array.
+	 * {a<sub>1</sub>, a<sub>2</sub>, ..., a<sub>mn</sub>} →
+	 * {{a<sub>1</sub>, a<sub>2</sub>, ..., a<sub>1n</sub>}, 
+	 * {a<sub>1n+1</sub>, a<sub>1n+2</sub>, ..., a<sub>2n</sub>}
+	 * @param a one-dimensional array of doubles, packed by columns.
+	 * @param m number of rows.
+	 * @exception IllegalArgumentException Array length must be a multiple of m.
+	 */
+	public Matrix(double[] a, int m){
+		this(m,(m!=0 ? a.length/m : 0));
+		if(m*n != a.length)
+			throw new IllegalArgumentException("Array length must be a multiple of m.");
+		int i,j;
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				set(i, j, a[i+j*m]);
+	}
+	
+	/**
+	 * Create an m-by-n matrix based on a 2d array.
+	 * @param a 2d int array
+	 * @exception IllegalArgumentException All rows must have the same length
+	 */
+	public Matrix(int[][] a){
+		this(a.length, a[0].length);
+		int i,j;
+        for(i=0; i<m; i++){
+        	if(a[i].length != n)
+        		throw new IllegalArgumentException("All rows must have the same length.");
+            for(j=0; j<n; j++)
+            	set(i, j, a[i][j]);
+        }
+	}
+	
+	/**
+	 * Create an m-by-n matrix based on a 2d array.
+	 * @param a double array
+	 */
+	public Matrix(double[][] a){
+		this(a.length, a[0].length);
+		int i,j;
+        for(i=0; i<m; i++){
+        	if(a[i].length != n)
+        		throw new IllegalArgumentException("All rows must have the same length.");
+            for(j=0; j<n; j++)
+            	set(i, j, a[i][j]);
+        }
+	}
+	
+	/**
+	 * Create an m-by-1 matrix based on a vector.
+	 * @param a vector
+	 * @see #Matrix(double[])
+	 * @see #Matrix(Vector, Vector)
+	 */
+    public Matrix(Vector a){
+    	this(a.getArray());
+    }
+	
+	/**
+	 * Create an m-by-2 matrix based on two vectors.
+	 * @param a vector
+	 * @param b vector
+	 * @see #Matrix(double[])
+	 * @see #Matrix(Vector)
+	 */
+    public Matrix(Vector a, Vector b){
+    	this(a.getArray(), b.getArray());
+    }
+	
+	/**
+	 * Create an m-by-n matrix based on a vector array.
+	 * @param a vector array
+	 * @exception IllegalArgumentException All rows must have the same length
+	 */
+    public Matrix(Vector[] a){
+    	this(a.length, a[0].n());
+        int i,j;
+        for(i=0; i<m; i++){
+        	if(a[i].n() != n)
+        		throw new IllegalArgumentException("All rows must have the same length.");
+        	for(j=0; j<n; j++)
+        		set(i, j, a[i].get(j));
+        }
+    }
+	
+	/**
+	 * Copy constructor.
+	 * @param a matrix
+	 */
+    public Matrix(Matrix a){
+    	this(a.n, a.get());
+    }
+	
+	public Matrix create(int m, int n){
+		if(this instanceof TensorII2D) return new TensorII2D();
+		if(this instanceof TensorII) return new TensorII();
+		if(this instanceof Matrix3D) return new Matrix3D(n);
+		if(this instanceof Matrix2D) return new Matrix2D(n);
+		if(this instanceof Matrix3D) return new Matrix3D(n);
+		if(this instanceof Diagonal) return new Diagonal(m,n);
+		else return new Matrix(m,n);
+	}
+	
+	public Matrix create(int m, int n, Matrix B){
+        if(m == 2 && B.n == 2) return new TensorII2D();
+        if(m == 3 && B.n == 3) return new TensorII();
+		if(this instanceof Matrix3D) return new Matrix3D(n);
+		if(this instanceof Matrix2D) return new Matrix2D(n);
+		if(this instanceof Matrix3D) return new Matrix3D(n);
+		if(this instanceof Diagonal) return new Diagonal(m,n);
+		else return new Matrix(m,n);
+	}
+	
+	public Matrix create(Matrix B){
+		if(this instanceof TensorII2D) return new TensorII2D(B);
+		if(this instanceof Matrix2D) return new Matrix2D(B);
+		if(this instanceof Matrix3D) return new Matrix3D(B);
+		if(this instanceof Diagonal) return new Diagonal(B);
+		else return new Matrix(B);
+	}
+	
+	private Vector create(int n){
+		if(this instanceof Matrix2D) return new Vector2D();
+		if(this instanceof Matrix3D) return new Vector3D();
+		else return new Vector(n);
+	}
+	
+	private Vector create(double[] v){
+		if(this instanceof Matrix2D) return new Vector2D(v);
+		if(this instanceof Matrix3D) return new Vector3D(v);
+		else return new Vector(v);
+	}
+	
+	private Vector create(Vector v){
+		if(this instanceof Matrix2D) return new Vector2D(v);
+		if(this instanceof Matrix3D) return new Vector3D(v);
+		else return new Vector(v);
+	}
+	
+	/**
+	 * Generate an n-by-n identity matrix.
+	 * An m-by-n matrix with ones on the diagonal and zeros elsewhere.
+	 * @param n rows/columns
+	 * @return n-by-n identity matrix
+	 */
+	public static Matrix identity(int n){
+        return Diagonal.identity(n);
+	}
+	
+	/**
+	 * Generate an m-by-n identity matrix.
+	 * An m-by-n matrix with ones on the diagonal and zeros elsewhere.
+	 * @param m rows
+	 * @param n columns
+	 * @return m-by-n identity matrix
+	 */
+	public static Matrix identity(int m, int n){
+        return Diagonal.identity(m, n);
+	}
+	
+	/**
+	 * Create an m-by-m diagonal matrix.
+	 * @return m-by-m diagonal matrix
+	 */
+	public static Matrix diag(double... d){
+		int i;
+		int n = d.length;
+        Matrix c = new Matrix(n, n);
+        for(i=0; i<n; i++)
+        	c.set(i, i, d[i]);
+        return c;
+	}
+	
+    /**
+     * Generate matrix with random elements.
+     * @param m rows and columns
+     * @return random m-by-m matrix with values between 0 and 1
+     */
+    public static Matrix random(int m) {
+        return random(m, m);
+    }
+	
+    /**
+     * Generate matrix with random elements
+     * @param m rows
+     * @param n columns
+     * @return random m-by-n matrix with values between 0 and 1
+     */
+    public static Matrix random(int m, int n) {
+		int i,j;
+        Matrix c = new Matrix(m, n);
+        for(i=0; i<m; i++)
+            for(j=0; j<n; j++)
+                c.set(i, j, Math.random());
+        return c;
+    }
+	
+    /**
+     * Element-wise cosine of argument in radians.
+     * @param a matrix
+     * @return the cosine for each element of the matrix
+     */
+    public static Matrix cos(Matrix a) {
+		int i,j;
+		int m = a.getM();
+		int n = a.getN();
+        Matrix c = new Matrix(m, n);
+        for(i=0; i<m; i++)
+            for(j=0; j<n; j++)
+                c.set(i, j, Math.cos(a.get(i, j)));
+        return c;
+    }
+	
+    /**
+     * Element-wise cosine of argument in radians.
+     * @param a matrix
+     * @return the cosine for each element of the matrix
+     */
+    public static Matrix sin(Matrix a) {
+		int i,j;
+		int m = a.getM();
+		int n = a.getN();
+        Matrix c = new Matrix(m, n);
+        for(i=0; i<m; i++)
+            for(j=0; j<n; j++)
+                c.set(i, j, Math.sin(a.get(i, j)));
+        return c;
+    }
+    
+    /**
+     * Clone this object.
+     * Make a deep copy of this matrix.
+     * @return copy matrix
+     */
+    public Matrix clone(){
+		int i,j;
+        Matrix c = create(m, n);
+        for(i=0; i<m; i++)
+            for(j=0; j<m; j++)
+                c.set(i, j, get(i,j));
+        return c;
+    }
+
+    /**
+     * @return the number of rows
+     */
+	public int getM(){
+		return m;
+	}
+
+	/**
+	 * @param i index
+	 * @return the row vector
+	 */
+	public Vector getM(int i){
+		Vector c = create(n);
+		int j;
+		for(j=0; j<n; j++)
+			c.set(j, get(i,j));
+		return c;
+	}
+
+	/**
+	 * @param indices
+	 * @return the row vector
+	 */
+	public Matrix getM(int[] indices){
+		int ni = indices.length;
+		Matrix c = new Matrix(indices.length,n);
+		int i,j;
+		for(i=0; i<ni; i++)
+			for(j=0; j<n; j++)
+				c.set(i, j, get(indices[i],j));
+		return c;
+	}
+
+	/**
+	 * @return the number of columns
+	 */
+	public int getN(){
+		return n;
+	}
+
+	/**
+	 * @param i index
+	 * @return the column vector
+	 * @throws ArrayIndexOutOfBoundsException
+	 */
+	public Vector getN(int i){
+		Vector c = create(m);
+		int j;
+		for(j=0; j<m; j++)
+			c.set(j, get(j,i));
+		return c;
+	}
+
+	/**
+	 * @return reference to the double array of the matrix
+	 */
+	public double[][] getD(){
+		int i, j;
+		double[][] result = new double[m][n];
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++){
+				result[i][j] = get(i,j);
+			}
+		return result;
+	}
+	
+	/**
+	 * @return reference to the double array of the matrix
+	 */
+	public double[] get(){
+		return data;
+	}
+
+	/**
+	 * @return copy of the double array of the matrix
+	 */
+	public double[][] getCopy(){
+		int i,j;
+		double[][] c = new double[m][n];
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				c[i][j] = get(i,j);
+		return c;
+	}
+
+	/**
+	 * Get the element A<sub>ij</sub> at given row i and column j
+	 * @param i row
+	 * @param j column
+	 * @return the value at the given row i and column j
+	 * @throws IllegalDimensionException 
+	 */
+	public double get(int i, int j){
+//		return data[i][j];
+		int index;
+		index = i*n + j;
+		return data[index];
+	}
+
+	/**
+	 * Get a sub matrix.
+	 * @param is first row index
+	 * @param ie last row index
+	 * @param js first column index
+	 * @param je last column index
+	 * @return sub matrix
+	   */
+	public Matrix get(int is, int ie, int js, int je) throws ArrayIndexOutOfBoundsException{
+		int i,j;
+		Matrix c = new Matrix(ie-is+1,je-js+1);
+		for(i=is; i<=ie; i++)
+			for(j=js; j<=je; j++)
+				c.set(i-is, j-js, get(i,j));
+		return c;
+	}
+    
+    /**
+     * Get a sub matrix only with chosen rows and columns.
+	 * @param i row indices
+	 * @param j column indices
+	 * @return sub matrix
+     */
+    public Matrix get(int[] i, int[] j){
+		return get(i, j, false);
+    }
+
+	/**
+	 * Get a sub matrix only with chosen rows and columns.
+	 * @param i row indices
+	 * @param j column indices
+	 * @return sub matrix
+	 */
+	public Matrix get(int[] i, int[] j, boolean positions){
+		int ii,jj;
+		int ni = i.length;
+		int mj = j.length;
+    	int d = positions?1:0;
+		Matrix a = new Matrix(ni,mj);
+		for(ii=0; ii<ni; ii++)
+			for(jj=0; jj<mj; jj++)
+				a.set(ii, jj, get(i[ii]-d, j[jj]-d));
+		return a;
+	}
+
+	/**
+	 * Get a sub matrix.
+	 * @param is first row index
+	 * @param ie last row index
+	 * @param j column indices
+	 * @return sub matrix
+	 */
+	public Matrix get(int is, int ie, int[] j){
+		int i,jj;
+		int mj = j.length;
+		Matrix a = new Matrix(ie-is+1,mj);
+		for(i=is; i<=ie; i++)
+			for(jj=0; jj<mj; jj++)
+				a.set(i-is, jj, get(i,j[jj]));
+		return a;
+	}
+    
+    /**
+     * Creates a sub matrix only with chosen rows and column
+     * @param i rows indices
+     * @param j index
+     * @return the sub matrix
+     */
+    public Vector get(int[] i, int j){
+    	int m = i.length;
+    	int ii;
+		Vector red = new Vector(m);
+		for(ii=0; ii<m; ii++)
+			red.set(ii, get(i[ii],j));
+		return red;
+    }
+
+	/**
+	 * Get a sub matrix only with chosen rows and columns.
+	 * @param i row indices
+	 * @param js starting column index
+	 * @param je ending column index
+	 * @return sub matrix
+	 */
+	public Matrix get(int[] i, int js, int je){
+		int ii,j;
+		int mi = i.length;
+		Matrix c = new Matrix(mi,je-js+1);
+		for(ii=0; ii<mi; ii++)
+			for(j=js; j<=je; j++)
+				c.set(ii, j-js, get(i[ii],j));
+		return c;
+	}
+    
+    /**
+     * Creates a sub matrix only with chosen rows and columns
+     * @param indices of the type which are defined in the boolean row 
+     * @param row boolean, true for the whole row and false for the whole column 
+     * @return the sub matrix
+     */
+    public Matrix get(int[] indices, boolean row){
+    	int m = this.m;
+    	int n = this.n;
+    	if(row) m = indices.length;
+    	else n = indices.length;
+    	int i,j;
+		Matrix c = new Matrix(m,n);
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++){
+				if(row)	c.set(i,j, get(indices[i],j));
+				else c.set(i,j, get(i,indices[j]));
+			}
+		return c;
+    }
+
+	/**
+	 * @return matrix elements packed in an 1d array by columns.
+	 */
+	public double[] getColumnPacked(){
+		int i,j;
+		double[] c = new double[m*n];
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				c[i+j*m] = get(i,j);
+		return c;
+	}
+
+	/**
+	 * @return matrix elements packed in a 1d array by rows.
+	 */
+	public double[] getRowPacked(){
+		int i,j;
+		double[] c = new double[m*n];
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				c[i*n+j] = get(i,j);
+		return c;
+	}
+
+	/**
+	 * Set value at given row i and column j
+	 * @param i row
+	 * @param j column
+	 * @param a value
+	 * @return filled matrix
+	 */
+	public Matrix set(int i, int j, double a){
+//		data[i][j] = a;
+		int index = i*n + j;
+		data[index] = a;
+//		System.out.println(i+" "+j+" = "+index+" : "+a);
+		return this;
+	}
+	
+	/**
+	 * Set the matrix with given vector at given row
+	 * @param i row
+	 * @param v vector
+	 * @return filled matrix
+	 */
+	public Matrix setM(int i, Vector v){
+		int ii = i;
+		int jj;
+		for(jj=0; jj<n; jj++)
+			set(ii, jj, v.get(jj));
+		return this;
+	}
+	
+	/**
+	 * Set the matrix with given array at given column
+	 * @param j column
+	 * @param vector values
+	 * @return filled matrix
+	 */
+	public Matrix setN(int j, double[] vector){
+//		if(m < vector.length || n<j){
+//			m = vector.length;
+//			n = j;
+//			data = new double[m][n];
+//		}
+		int i;
+		for(i=0; i<m; i++)
+			set(i, n, vector[i]);
+		return this;
+	}
+	
+	/**
+	 * Set the matrix with given vector at given column
+	 * @param j column
+	 * @param v vector
+	 * @return filled matrix
+	 */
+	public Matrix setN(int j, Vector v){
+		int jj = j;
+		int ii;
+		for(ii=0; ii<m; ii++)
+			set(ii, jj, v.get(ii));
+		return this;
+	}
+
+	/**
+	 * Set the matrix with given double array
+	 * @param B matrix
+	 * @return filled matrix
+	 */
+	public Matrix set(double[][] B){
+		return set(0, 0, B);
+	}
+
+	/**
+	 * Set the matrix with given sub matrix.
+	 * @param i first row index
+	 * @param j first column index
+	   @param B sub matrix
+	   @return filled matrix
+	 */
+	public Matrix set(int i, int j, double[][] B){
+		int ii, jj;
+		for(ii=i; ii<=B.length; ii++)
+			for(jj= j; jj<=B[0].length; jj++)
+				set(ii,jj, B[ii-i][jj-j]);
+		return this;
+	}
+
+	/**
+	 * Set the matrix with given sub matrix.
+	 * @param i first row index
+	 * @param j first column index
+	 * @param B sub matrix
+	 * @return filled matrix
+	 */
+	public Matrix set(int i, int j, Matrix B){
+		int ii, jj;
+		for(ii=i; ii<(i+B.m); ii++)
+			for(jj=j; jj<(j+B.n); jj++)
+				set(ii,jj, B.get(ii-i,jj-j));
+		return this;
+	}
+
+	/**
+	 * Set a submatrix.
+	 * @param i row indices
+	 * @param j column indices
+	 * @param B submatrix
+	 * @return filled matrix
+	 */
+	public Matrix set(int[] i, int[] j, Matrix B){
+		int ii, jj;
+		for(ii=0; ii<i.length; ii++)
+			for (jj=0; jj<j.length; jj++)
+				set(i[ii],j[jj], B.get(ii,jj));
+		return this;
+	}
+
+	/**
+	 * Set matrix with a sub matrix only with chosen rows and columns.
+	 * @param i row indices
+	 * @param j starting column index
+	 * @param B sub matrix
+	 * @return filled matrix
+	 */
+	public Matrix set(int[] i, int j, Matrix B){
+		int ii,jj;
+		for(ii=0; ii<i.length; ii++)
+			for(jj=j; jj<(j+B.n); jj++)
+				set(i[ii],jj, B.get(ii,jj-j));
+		return this;
+	}
+
+	/** Set a submatrix.
+	 * @param i starting row index
+	 * @param j column indices
+	 * @param B sub matrix
+	 * @return filled matrix
+	 */
+	public Matrix set(int i, int[] j, Matrix B){
+		int ii,jj;
+		for(ii=i; ii<(i+B.m); ii++)
+			for(jj=0; jj<j.length; jj++)
+				set(ii,j[jj], B.get(ii-i,jj));
+		return this;
+	}
+	
+    /**
+     * Set values from a sub matrix only with chosen rows and columns
+     * @param indices of the type which are defined in the boolean row
+     * @param row boolean, true for the whole row and false for the whole column 
+     * @param B sub matrix 
+     * @return filled matrix
+     */
+    public Matrix set(int[] indices, boolean row, Matrix B){
+    	int m = B.m;
+    	int n = B.n;
+    	if(row) m = indices.length;
+    	else n = indices.length;
+    	int i,j;
+		if(row)
+			for(i=0; i<m; i++)
+				for(j=0; j<n; j++)
+					set(indices[i],j, B.get(i,j));
+		else
+			for(i=0; i<m; i++)
+				for(j=0; j<n; j++)
+					set(i,indices[j], B.get(i,j));
+		return this;
+    }
+    
+    /**
+     * a<sub>ij</sub> = a<sub>ij</sub> + b
+     * @param i row
+     * @param j column
+     * @param b scalar
+     */
+    public void plus(int i, int j, double b){
+    	set(i,j, get(i,j)+b);
+    }
+
+    /**
+     * Adds two matrices of the same dimension (entrywise addition).
+     * The addition of two m-by-n matrices <b>A</b> and <b>B</b>,
+     * is again an m-by-n matrix computed by adding corresponding elements.
+     * @param B matrix
+     * @return <b>C</b> = <b>A</b> + <b>B</b>, entrywise summarized matrix
+     * @throws IllegalDimensionException
+     */
+    public Matrix plus(Matrix B) throws IllegalDimensionException{
+    	checkDimensions(B);
+    	int i,j;
+        Matrix c = create(this);
+        if(B instanceof Diagonal)
+        	for(i=0; i<n; i++)
+        		c.plus(i, i, ((Diagonal) B).get(i));
+        else
+        	for(i=0; i<m; i++)
+        		for(j=0; j<n; j++)
+        			c.plus(i,j, B.get(i,j));
+        return c;
+    }
+
+    /**
+     * Adds two matrices of the same dimension (entrywise addition).
+     * The addition of two m-by-n matrices <b>A</b> and <b>B</b>,
+     * is again an m-by-n matrix computed by adding corresponding elements.
+     * @param B matrix
+     * @return <b>A</b> = <b>A</b> + <b>B</b>, entrywise summarized matrix
+     * @throws IllegalDimensionException
+     */
+    public Matrix plusEquals(Matrix B) throws IllegalDimensionException {
+    	checkDimensions(B);
+    	int i, j;
+        if(B instanceof Diagonal)
+        	for(i=0; i<n; i++)
+        		plus(i, i, ((Diagonal) B).get(i));
+        else
+	    	for(i=0; i<m; i++)
+	    		for(j=0; j<n; j++)
+	                plus(i,j, B.get(i,j));
+    	return this;
+    }
+
+    /**
+     * Unary minus
+     * @return -A
+     */
+    public Matrix uminus(){
+    	Matrix X = create(m,n);
+    	int i, j;
+        if(this instanceof Diagonal)
+        	for(i=0; i<n; i++)
+        		X.set(i, i, -((Diagonal) this).get(i));
+        else
+	    	for(i=0; i<m; i++)
+	    		for(j=0; j<n; j++)
+	    			X.set(i,j, -get(i,j));
+    	return X;
+    }
+    
+    public void minus(int i, int j, double a){
+    	plus(i,j,-a);
+    }
+    
+    /**
+     * Subtract two matrices of the same dimension (entrywise subtraction).
+     * The subtraction of two m-by-n matrices <b>A</b> and <b>B</b>,
+     * is again an m-by-n matrix computed by subtracting corresponding elements.
+     * @param B matrix
+     * @return <b>C</b> = <b>A</b> - <b>B</b>, entrywise subtracted matrix
+     * @throws IllegalDimensionException
+     */
+    public Matrix minus(Matrix B) throws IllegalDimensionException{
+    	checkDimensions(B);
+    	int i,j;
+        Matrix c = create(B);
+        System.out.println(this);
+        System.out.println(B);
+        if(B instanceof Diagonal)
+        	for(i=0; i<n; i++)
+        		c.minus(i, i, ((Diagonal) B).get(i));
+        else
+	        for(i=0; i<m; i++)
+	            for(j=0; j<n; j++)
+	                c.minus(i,j, B.get(i,j));
+        return c;
+    }
+    
+    /**
+     * Subtract two matrices of the same dimension (entrywise subtraction).
+     * The subtraction of two m-by-n matrices <b>A</b> and <b>B</b>,
+     * is again an m-by-n matrix computed by subtracting corresponding elements.
+     * @param B matrix
+     * @throws IllegalDimensionException
+     */
+    public Matrix minusEquals(Matrix B) throws IllegalDimensionException {
+    	checkDimensions(B);
+    	int i,j;
+        if(B instanceof Diagonal)
+        	for(i=0; i<n; i++)
+        		minus(i, i, ((Diagonal) B).get(i));
+        else
+	        for(i=0; i<m; i++)
+	            for(j=0; j<n; j++)
+	                minus(i,j, B.get(i,j));
+        return this;
+    }
+    
+    /**
+     * Direct sum.
+     * The direct sum of any pair of matrices <b>A</b> of size m x n and
+     * <b>B</b> of size p x q is a matrix of size (m + p) x (n + q).
+     * @param B matrix
+     * @return <b>C</b> = <b>A</b> &oplus; <b>B</b>  matrix of size (m + p) × (n + q)
+     */
+    public Matrix sumDirect(Matrix B){
+    	int i,j = 0,k,l;
+    	Matrix c = new Matrix(m+B.m,n+B.n);
+    	for(i=0; i<m; i++) // first, old matrix, top left
+    		for(j=0; j<n; j++)
+        		c.set(i,j, get(i,j));
+        if(B instanceof Diagonal) // second, bottom right
+        	for(k=0; k<B.n; k++)
+        		c.set(i+k, j+k, ((Diagonal) B).get(i));
+        else
+	    	for(k=0; k<B.m; k++) // second, bottom right
+	    		for(l=0; l<B.n; l++)
+	    			c.set(i+k,j+l, B.get(k,l));
+    	return c;
+    }
+    
+    /**
+     * Scalar multiplication. Multiply a matrix element-wise by a scalar.
+     * Scaling by s in all directions.
+     * @param s scalar
+     * @return scaled matrix <b>C</b> = s <b>A</b>
+     */
+    public Matrix times(double s){
+    	int i,j;
+    	Matrix c = create(m,n);
+        if(this instanceof Diagonal)
+        	for(i=0; i<m; i++)
+        		c.set(i, i, s * ((Diagonal) this).get(i));
+        else
+	        for(i=0; i<m; i++)
+	            for(j=0; j<n; j++)
+	                c.set(i,j, s * get(i,j));
+        return c;
+    }
+    
+    /**
+     * c<sub>ij</sub> = sc<sub>ij</sub> 
+     * @param i
+     * @param j
+     * @param s
+     */
+    public void times(int i, int j, double s){
+    	set(i,j, s*get(i,j));
+    }
+    
+    /**
+     * Scalar multiplication. Multiply a matrix element-wise by a scalar.
+     * @param s scalar
+     * @return <b>A</b> = s <b>A</b>
+     */
+    public Matrix timesEquals(double s){
+    	int i,j;
+        if(this instanceof Diagonal)
+        	for(i=0; i<m; i++)
+        		set(i, i, s * ((Diagonal) this).get(i));
+        else
+	        for(i=0; i<m; i++)
+	            for(j=0; j<n; j++)
+	                set(i,j, s * get(i,j));
+        return this;
+    }
+	
+    /**
+     * Multiply by a vector right. matrix-vector multiplication
+     * @param b vector
+     * @return <b>C</b> = <b>A</b> • <b>b</b> = <b>A</b> <b>b</b>
+     * @throws IllegalDimensionException Illegal vector dimensions.
+     */
+	public Vector times(Vector b) throws IllegalDimensionException{
+		int i,j;
+        if(n != b.n()) throw new IllegalDimensionException("Illegal vector dimensions.");
+		Vector c = create(m);
+        if(this instanceof Diagonal)
+        	for(i=0; i<m; i++)
+        		c.set(i, ((Diagonal) this).get(i) * b.get(i));
+        else
+			for(i=0; i<m; i++)
+				for(j=0; j<n; j++)
+					c.plus(i, (get(i,j) * b.get(j)));
+		return c;
+	}
+	
+    /**
+     * Multiply a matrix right. Linear algebraic matrix multiplication.
+     * @param B matrix
+     * @return <b>C</b> = <b>A</b> <b>B</b>
+     * @throws IllegalDimensionException Inner matrix dimensions must agree.
+     */
+	public Matrix times(Matrix B) throws IllegalDimensionException{
+		int i,j,k;
+        if(n != B.m) throw new IllegalDimensionException("Inner matrix dimensions must agree.");
+        Matrix c = create(m,B.n,B);
+        if(B instanceof Diagonal)
+        	for(i=0; i<c.m; i++)
+        		for(j=0; j<c.n; j++)
+					c.set(i,j, get(i,j) * ((Diagonal) B).get(j));
+        else
+			for(i=0; i<c.m; i++)
+				for(j=0; j<c.n; j++)
+					for(k=0; k<n; k++)
+						c.plus(i,j, get(i,k) * B.get(k,j));
+		return c;
+	}
+	
+    /**
+     * Multiply a matrix right. Linear algebraic matrix multiplication.
+     * Plus scalar multiplication, multiply a matrix element-wise by a scalar.
+     * Scaling by s in all directions.
+     * @param s scalar
+     * @param B matrix
+     * @return <b>C</b> = <b>A</b> s<b>B</b>
+     * @throws IllegalDimensionException Inner matrix dimensions must agree.
+     */
+    public Matrix times(double s, Matrix B) throws IllegalDimensionException{
+    	return times(s).times(B);
+    }
+	
+	// TODO: too lame, check 0 then calculate only sub matrix
+    /**
+     * Multiply a matrix right
+     * @param B matrix
+     * @return <b>C</b> = <b>A</b> <b>B</b>
+     * @throws IllegalDimensionException Inner matrix dimensions must agree.
+     */
+	public Matrix timesV2(Matrix B) throws IllegalDimensionException{
+		int i,j,k;
+		Vector p,q;
+        if(n != B.m) throw new IllegalDimensionException("Inner matrix dimensions must agree.");
+		Matrix c = create(m,B.n);
+		for(i=0; i<c.m; i++){
+			p = getM(i);
+			if(p.isZero()) continue;
+			for(j=0; j<c.n; j++){
+				q = B.getN(j);
+				if(q.isZero()) continue;
+				for(k=0; k<n; k++)
+					c.plus(i,j, p.get(k) * q.get(k));
+			}
+		}
+		return c;
+	}
+	
+	/**
+	 * Multiply a vector left and a vector right.
+	 * @param u vector
+	 * @param v vector
+	 * @return c = <b>u</b> • <b>A</b> • <b>v</b> = <b>u</b><sup>T</sup> <b>A</b> <b>v</b>
+	 * @throws IllegalDimensionException
+	 */
+	public double times(Vector u, Vector v) throws IllegalDimensionException{
+		return u.times(this).dot(v);
+	}
+	
+	/**
+	 * Multiply a matrix left and a matrix right
+	 * @param U matrix
+	 * @param V matrix
+	 * @return <b>C</b> = <b>U</b> <b>A</b> <b>V</b>
+	 * @throws IllegalDimensionException
+	 */
+	public Matrix times(Matrix U, Matrix V) throws IllegalDimensionException{
+		return U.times(this).times(V);
+	}
+	
+    /**
+     * Hadamard / Schur Product. Element-wise multiplication with a matrix B
+     * @param B matrix
+     * @return c<sub>i,j</sub> = a<sub>i,j</sub> b<sub>i,j</sub>
+     * @throws IllegalDimensionException
+     */
+	public Matrix hadamardProduct(Matrix B) throws IllegalDimensionException {
+    	checkDimensions(B);
+		int i,j;
+        if(this instanceof Diagonal){
+        	Diagonal c = new Diagonal(m,B.n);
+        	for(i=0; i<m; i++)
+        		c.set(i, ((Diagonal) this).get(i) * B.get(i, i));
+        	return c;
+        } else if(B instanceof Diagonal){
+        	Diagonal c = new Diagonal(m,B.n);
+        	for(i=0; i<m; i++)
+        		c.set(i, get(i,i) * ((Diagonal) B).get(i));
+        	return c;
+        } else {
+        	Matrix c = create(m,B.n);
+			for(i=0; i<c.m; i++)
+				for(j=0; j<c.n; j++)
+					c.set(i,j, get(i,j) * B.get(i,j));
+			return c;
+        }
+	}
+	
+    /**
+     * Hadamard / Schur Product. Element-wise multiplication with a matrix B
+     * @param B matrix
+     * @return a<sub>i,j</sub> = a<sub>i,j</sub> b<sub>i,j</sub>
+     * @throws IllegalDimensionException
+     */
+	public Matrix hadamardProductEquals(Matrix B) throws IllegalDimensionException {
+    	checkDimensions(B);
+		int i,j;
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				set(i,j, get(i,j) * B.get(i,j));
+		return this;
+	}
+	
+	/**
+	 * Frobenius product.
+	 * @param B matrix
+	 * @return <b>A</b>:<b>B</b>=
+	 * &Sigma;<sub>i</sub>&Sigma;<sub>j</sub>A<sub>ij</sub>B<sub>ij</sub>=
+	 * vec(<b>A</b>)<sup>T</sup>vec(<b>B</b>)=
+	 * tr(<b>A</b><sup>T</sup><b>B</b>)=tr(<b>A</b><b>B</b><sup>T</sup>)
+	 * @throws IllegalDimensionException
+	 * @throws NoSquareException
+	 */
+	public double frobeniusProduct(Matrix B) throws IllegalDimensionException, NoSquareException{
+		return this.transpose().times(B).trace();
+	}
+	
+	/**
+	 * Kronecker product.
+	 * @param B matrix
+	 * @return (A&otimes;B)<sub>ij</sub>=(A)<sub>ij</sub><b>B</b>
+	 */
+	public Matrix kroneckerProduct(Matrix B){
+		int i, j;
+		Matrix C = new Matrix(m*B.m,n*B.n);
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				C.set(i*B.m, j*B.n, B.times(this.get(i, j)));
+		return C;
+	}
+
+    /**
+     * Hadamard / Schur Division. Element-wise division with a matrix B
+     * @param B matrix
+     * @return c<sub>i,j</sub> = a<sub>i,j</sub> / b<sub>i,j</sub>
+     * @throws IllegalDimensionException
+     */
+	public Matrix hadamardDivision(Matrix B) throws IllegalDimensionException {
+    	checkDimensions(B);
+		int i,j;
+		Matrix c = create(m,n);
+		for(i=0; i<c.m; i++)
+			for(j=0; j<c.n; j++)
+				c.set(i,j, get(i,j) / B.get(i,j));
+		return c;
+	}
+
+    /**
+     * Hadamard / Schur Division. Element-wise division with a matrix B
+     * @param B matrix
+     * @return a<sub>i,j</sub> = a<sub>i,j</sub> / b<sub>i,j</sub>
+     * @throws IllegalDimensionException
+     */
+	public Matrix hadamardDivisionRightEquals(Matrix B) throws IllegalDimensionException {
+    	checkDimensions(B);
+		int i,j;
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				set(i,j, get(i,j) / B.get(i,j));
+		return this;
+	}
+
+    /**
+     * Hadamard / Schur Division. Element-wise division with a matrix B
+     * @param B matrix
+     * @return a<sub>i,j</sub> = b<sub>i,j</sub> / a<sub>i,j</sub>
+     * @throws IllegalDimensionException
+     */
+	public Matrix hadamardDivisionLeftEquals(Matrix B) throws IllegalDimensionException {
+    	checkDimensions(B);
+		int i,j;
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				set(i,j, B.get(i,j) / get(i,j));
+		return this;
+	}
+    
+    /**
+     * Scalar division
+     * @param s scalar
+     * @return <b>C</b> = <b>A</b> / s
+     */
+    public Matrix over(double s){
+        return times(1/s);
+    }
+
+	/**
+	 * Divide this matrix <b>A</b> by a matrix <b>B</b>
+	 * @param B matrix
+	 * @return <b>C</b> = <b>A</b> / <b>B</b>
+	 * @throws NoSquareException
+	 * @throws SingularException
+	 * @throws IllegalDimensionException
+	 */
+	public Matrix over(Matrix B) throws NoSquareException, SingularException, IllegalDimensionException{
+    	return times(B.inverse());
+	}
+	
+	/**
+	 * @param e exponent
+	 * @return the matrix to the power of e
+	 * @throws IllegalDimensionException
+	 */
+	public Matrix pow(double e) throws IllegalDimensionException{
+		int i;
+		Matrix C = create(m,n);
+		if(this instanceof Diagonal)
+			for(i=0; i<(m<n?m:n); i++)
+				C.set(i,i, Math.pow(((Diagonal) this).get(i), e));
+		else if(isDiagonal()){
+			for(i=0; i<(m<n?m:n); i++)
+				C.set(i,i, Math.pow(get(i,i), e));
+		} else {
+			throw new IllegalArgumentException("Operation only defined for diagonal matrix.");
+		}
+		return C;
+	}
+	
+	/**
+	 * Transpose of a matrix by swapping the columns with rows
+	 * @return the transposed matrix A<sup>T</sup>
+	 */
+	public Matrix transpose(){
+		int i, j;
+		if(this instanceof Diagonal || isDiagonal())
+			return new Diagonal(this);
+		Matrix c = create(n,m);
+//		System.out.println(m + " " + n);
+//		System.out.println(c.m + " " + c.n);
+		for(i=0; i<(m<n?m:n); i++) // copy main diagonal
+			c.set(i,i, get(i,i));
+		if(isSquare()){
+	        for(i=0; i<n; i++)
+	            for(j=i+1; j<n; j++){
+	                c.set(i,j, get(j,i));
+	                c.set(j,i, get(i,j));
+	            }
+		}
+		else if(m<=n)
+			for(i=0; i<m; i++)
+				for(j=i+1; j<n; j++){
+					if(j<m&&i<n) c.set(i,j, get(j,i));
+					c.set(j,i, get(i,j));
+				}
+		else
+			for(j=0; j<n; j++)
+				for(i=j+1; i<m; i++){
+					if(j<m&&i<n) c.set(i,j, get(j,i));
+					c.set(j,i, get(i,j));
+				}
+		return c;
+	}
+
+    /**
+     * swap rows i and j
+     * @param i row
+     * @param j row
+     */
+    public Matrix swap(int i, int j){
+        Vector temp = getM(i);
+        setM(i,getM(j));
+        setM(j, temp);
+        return this;
+    }
+    
+    /**
+     * swap columns or raws i and j
+     * @param i column
+     * @param j column
+     * @param columns true (false: raws)
+     */
+    public Matrix swap(int i, int j, boolean columns){
+    	if(columns){
+            Vector temp = getN(i);
+            setN(i,getN(j));
+            setN(j, temp);
+    	} else swap(i, j);
+    	return this;
+    }
+
+    /**
+     * One norm
+     * @return maximum column sum
+     */
+    public double norm1(){
+    	int i, j;
+    	double f = 0;
+    	for(j=0; j<n; j++){
+    		double s = 0;
+    		for(i=0; i<m; i++)
+    			s += Math.abs(get(i,j));
+    		f = Math.max(f,s);
+    	}
+    	return f;
+    }
+
+    /**
+     * Two norm
+     * @return maximum singular value
+     */
+    public double norm2(){
+    	return(new SingularValueDecomposition(this).norm2());
+    }
+
+    /** Infinity norm
+     * @return maximum row sum
+     */
+    public double normInf(){
+    	int i, j;
+    	double f = 0;
+    	for(i=0; i<m; i++){
+    		double s = 0;
+    		for(j=0; j<n; j++)
+    			s += Math.abs(get(i,j));
+    		f = Math.max(f,s);
+    	}
+    	return f;
+    }
+
+    /**
+     * Frobenius norm
+     * @return sqrt of sum of squares of all elements.
+     */
+    public double normF(){
+    	int i, j;
+    	double f = 0;
+    	for(i=0; i<m; i++)
+    		for(j=0; j<n; j++)
+    			f = Maths.hypot(f,get(i,j));
+    	return f;
+    }
+
+	/**
+	 * Inverse of a square matrix <b>A</b> is the matrix <b>A</b><sup>-1</sup>
+	 * where <b>A</b><sup>-1</sup> <b>A</b> = <b>I</b> ,
+	 * with <b>I</b> as the identity matrix.
+	 * A matrix that have inverse is called non-singular or invertible
+	 * otherwise it is called singular.
+	 * @return the inverse matrix. For a singular matrix the values of the
+	 * inverted matrix are either NaN or Infinity
+     * @throws SingularException
+     * @see #inverse()
+	 */
+	public Matrix inverse() throws SingularException{
+        if(n==1 && n==m){
+            if(get(0,0)==0.0) throw new SingularException();
+            return new Matrix(1.0/det());
+        }
+		return cofactor().transpose().times(1/det());
+	}
+
+	/**
+	 * Inverse of a square matrix <b>A</b> is the matrix <b>A</b><sup>-1</sup>
+	 * where <b>A</b><sup>-1</sup> <b>A</b> = <b>I</b> ,
+	 * with <b>I</b> as the identity matrix.
+	 * A matrix that have inverse is called non-singular or invertible
+	 * otherwise it is called singular.
+	 * @return the inverse matrix. For a singular matrix the values of the
+	 * inverted matrix are either NaN or Infinity
+	 * @throws NoSquareException
+     * @throws SingularException
+     * @see #inverse()
+	 */
+	public Matrix inverse(int analytic) throws NoSquareException, SingularException{
+        if(n==1 && n==m){
+            if(get(0,0)==0.0) throw new SingularException();
+            return new Matrix(1.0/det(analytic));
+        }
+		return cofactor(1).transpose().times(1/det(analytic));
+	}
+
+	/**
+	 * Inverse or pseudoinverse of a square matrix
+	 * NOTE: No good in FachwerkLinear: K<sup>-1</sup>
+	 * @return inverse(A) if A is square, pseudoinverse otherwise.
+	 * @see #inverse()
+	 */
+	public Matrix inverse2() {
+		return solve(identity(m,m));
+	}
+	
+	/**
+	 * Cofactor of a matrix<br />
+	 * The cofactor of a matrix <b>A</b> is matrix <b>C</b>
+	 * that the value of element C<sub>ij</sub> equals the determinant of a matrix
+	 * created by removing row i and column j from matrix <b>A</b>.
+	 * @return the cofactor of the matrix
+	 */
+	public Matrix cofactor(){
+		if(m==1) return new Matrix(get());
+		int i, j;
+		Matrix mat = create(m, n);
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				mat.set(i, j, Maths.changeSign(i)*Maths.changeSign(j)*sub(i,j).det());
+//		System.out.println(mat);
+		return mat;
+	}
+	
+	/**
+	 * Cofactor of a matrix<br />
+	 * The cofactor of a matrix <b>A</b> is matrix <b>C</b>
+	 * that the value of element C<sub>ij</sub> equals the determinant of a matrix
+	 * created by removing row i and column j from matrix <b>A</b>.
+	 * @return the cofactor of the matrix
+	 * @throws NoSquareException
+	 */
+	public Matrix cofactor(int analytic) throws NoSquareException {
+		if(m==1) return new Matrix(get());
+		int i, j;
+		Matrix mat = create(m, n);
+		for(i=0; i<m; i++)
+			for(j=0; j<n; j++)
+				mat.set(i, j, Maths.changeSign(i)*Maths.changeSign(j)*sub(i,j).det(analytic));
+//		System.out.println(mat);
+		return mat;
+	}
+	
+	/**
+	 * @return the trace
+	 */
+	public double trace() throws NoSquareException{
+		int i;
+		double trace = 0;
+		for(i=0; i<Math.min(m,n); i++){
+			trace += get(i,i);
+		}
+		return trace;
+	}
+
+	/**
+	 * @return the determinant of the matrix
+	 */
+	public double det(){
+		return new LUDecomposition(this).det();
+	}
+    
+    /**
+     * @return the determinant of the matrix
+     * @throws NoSquareException
+     */
+    public double det(int analytic) throws NoSquareException{
+    	if(m==1) return get(0,0);
+        if(!isSquare()) throw new NoSquareException();
+    	else{
+    		double determinant = 0;
+    		for(int i=0; i<m; i++) {
+    			determinant += Maths.changeSign(i)*get(i,0)*sub(i,0).det(analytic); // over rows
+//or   			determinant += WMath.changeSign(i)*get(0,i)*sub(0,i).det(analytic); // over columns
+    		}
+    		return determinant;
+    	}
+    }
+    
+    /**
+     * Creates a sub matrix excluding the given row and column
+     * @param row the excluded row
+     * @param col the excluded column
+     * @return the sub matrix
+     */
+    public Matrix sub(int row, int col){
+    	double [][] subMatrix = new double[m-1][m-1];
+    	int r, c, rowOffset=0, colOffset=0;
+     
+    	for(r=0; r<subMatrix.length; r++){
+    		colOffset=0;
+    		if(row==r) rowOffset=1;
+    		for(c=0; c<subMatrix.length; c++){
+    			if(col==c) colOffset=1;
+    			subMatrix[r][c]= get(r+rowOffset,c+colOffset);
+    		}
+    	}
+    	return new Matrix(subMatrix);
+    }
+    
+    /**
+     * @return vector with all diagonal values
+     */
+    public Vector diag(){
+    	int nn = n>m?n:m;
+    	Vector d = create(nn);
+    	int i;
+    	for(i=0; i<nn; i++)
+    		d.set(i, get(i,i));
+    	return d;
+    }
+
+
+    /**
+     * @return effective numerical rank, obtained from SVD.
+     * @see #rank(int)
+     */
+ 	  public int rank() {
+ 		  return new SingularValueDecomposition(this).rank();
+ 	  }
+    
+    /**
+     * @return the rank of the matrix, assuming the matrix is square
+     * @throws NoSquareException
+     */
+    public int rank(int analytic) throws NoSquareException{
+    	if(m==1) return 1;
+    	double determinant = det(analytic);
+    	if(determinant!=0) return m;
+    	else return subRank();
+    }
+
+    /**
+     * @return the rank of the matrix, assuming the matrix is square
+     * @throws NoSquareException
+     */
+    private int subRank() throws NoSquareException{
+    	int subRank = 0;
+    	for(int i=0; i<m; i++){
+    		if(i==m-1) subRank = subRank(true);
+    		else subRank = subRank(false);
+        	if(subRank!=0)
+        		return subRank;
+    	}
+		return subRank;
+    }
+    
+    /**
+     * @param deeper false: get the under rank; true: even deeper rank
+     * @return the rank of the matrix, assuming the matrix is square
+     * @throws NoSquareException
+     */
+    private int subRank(boolean deeper) throws NoSquareException{
+    	int i;
+    	for(i=0; i<m; i++){
+    		double determinant = sub(i,0).det();
+        	if(determinant!=0)
+        		return m-1;
+    	}
+    	int subRank = 0;
+    	if(deeper)
+    	for(i=0; i<m; i++){
+    		if(i==m-1) subRank = sub(i,0).subRank(true);
+    		else subRank = sub(i,0).subRank(false);
+        	if(subRank!=0)
+        		return subRank;
+    	}
+    	return subRank;
+    }
+
+    /**
+     * Matrix condition (2 norm)
+     * @return ratio of largest to smallest singular value s<sub>max</sub>/s<sub>min</sub>.
+     */
+    public double cond(){
+    	return new SingularValueDecomposition(this).cond();
+    }
+
+    /**
+     * Solve A X = B
+     * @param B right hand side
+     * @return solution if A is square, least squares solution otherwise
+     */
+    public Matrix solve(Matrix B){
+    	return (m == n ? (new LUDecomposition(this)).solve(B) :
+                         (new QRDecomposition(this)).solve(B));
+    }
+
+    /**
+     * Solve X A = B, which is also A' X' = B'
+     * @param B right hand side
+     * @return solution if A is square, least squares solution otherwise.
+    */
+    public Matrix solveTranspose(Matrix B){
+       return transpose().solve(B.transpose());
+    }
+	
+    /**
+     * <b>A</b> <b>x</b> = <b>b</b>,
+     * assuming <b>A</b> is square and has full rank
+     * @param b matrix (as vector)
+     * @return <b>x</b> = <b>A</b><sup>-1</sup> <b>b</b>
+     */
+    public Vector solve(Vector b){
+        if (m != n || b.n() != n)
+        	 throw new IllegalArgumentException("Illegal matrix dimensions.");
+
+        // create copies of the data
+        Matrix A = create(this);
+        Vector B = create(b);
+
+        // Gaussian elimination with partial pivoting
+        for (int i = 0; i < n; i++) {
+
+            // find pivot row and swap
+            int max = i;
+            for (int j = i + 1; j < n; j++)
+                if(Math.abs(A.get(j,i)) > Math.abs(A.get(max,i)))
+                    max = j;
+            A.swap(i, max);
+            B.swap(i, max);
+
+            // singular
+            if (A.get(i,i) == 0.0){
+                System.err.println("Matrix is singular.");
+                return null;
+            }
+
+            // pivot within b
+            for (int j = i + 1; j < n; j++)
+                B.set(j, B.get(j) - B.get(i) * A.get(j,i) / A.get(i,i));
+
+            // pivot within A
+            for (int j = i + 1; j < n; j++) {
+                double m = A.get(j,i) / A.get(i,i);
+                for (int k = i+1; k < n; k++) {
+                    A.minus(j,k, A.get(i,k) * m);
+                }
+                A.set(j,i, 0.0);
+            }
+        }
+
+        // back substitution
+        Vector x = create(n);
+        for (int j = n - 1; j >= 0; j--) {
+            double t = 0.0;
+            for (int k = j + 1; k < n; k++)
+                t += A.get(j,k) * x.get(k);
+            x.set(j, (B.get(j) - t) / A.get(j,j));
+        }
+        return x;
+    }
+
+    /**
+     * @return LUDecomposition
+     * @see LUDecomposition
+     */
+    public LUDecomposition lu(){
+    	return new LUDecomposition(this);
+    }
+
+    /**
+     * @return QRDecomposition
+     * @see QRDecomposition
+     */
+    public QRDecomposition qr(){
+    	return new QRDecomposition(this);
+    }
+
+    /**
+     * @return CholeskyDecomposition
+     * @see CholeskyDecomposition
+     */
+    public CholeskyDecomposition chol(){
+    	return new CholeskyDecomposition(this);
+    }
+
+    /**
+     * @return SingularValueDecomposition
+     * @see SingularValueDecomposition
+     */
+    public SingularValueDecomposition svd(){
+    	return new SingularValueDecomposition(this);
+    }
+
+    /**
+     * @return EigenvalueDecomposition
+     * @see EigenvalueDecomposition
+     */
+    public EigenvalueDecomposition eig(){
+    	eig = new EigenvalueDecomposition(this);
+    	return eig;
+    }
+
+    /**
+     * Eigenvalues calculated with the QR algorithm.
+     * @return eigenvalues λ<sub>i</sub>
+     * @see EigenvalueDecomposition
+     */
+    public Matrix eigenvalues(){
+    	if(eig == null) eig();
+    	return create(eig.getD());
+    }
+
+    /**
+     * Real eigenvalues calculated with the QR algorithm.
+     * @return eigenvalues λ<sub>i</sub>
+     * @see EigenvalueDecomposition
+     */
+    public Vector eigenvaluesRe(){
+    	if(eig == null) eig();
+    	return create(eig.getRealEigenvalues());
+    }
+
+    /**
+     * Imaginary eigenvalues calculated with the QR algorithm.
+     * @return eigenvalues λ<sub>i</sub>
+     * @see EigenvalueDecomposition
+     */
+    public VectorComplex eigenvaluesIm(){
+    	if(eig == null) eig();
+    	return new VectorComplex(eig.getRealEigenvalues(), eig.getImagEigenvalues());
+    }
+    
+    /**
+     * @return eigenvectors x<sub>i</sub> in matrix Λ
+     */
+    public Matrix eigenvectors(){
+    	if(eig == null) eig();
+    	return create(eig.getV());
+    }
+	
+    /**
+     * Real eigenvalue calculated with the root formula of second and third degree.
+     * @return eigenvalues λ<sub>i</sub>
+     * @throws NoSquareException 
+     * @throws IllegalDimensionException 
+     * @throws ComplexException 
+     */
+	public Vector eigenvaluesRe(int analytic) throws NoSquareException, IllegalDimensionException, ComplexException{
+		// TODO also n > 3
+		// TODO: also Imaginary
+		if(m!=n && m>3) throw new IllegalDimensionException("Illegal matrix dimensions.");
+		Vector lambda = new Vector(m);
+		if(m==1) lambda.set(0, get(0, 0));
+		else if(isR2()) lambda = new Matrix2D(this).eigenvaluesRe(analytic);
+		else if(isR3())	lambda = new Matrix3D(this).eigenvaluesRe(analytic);
+		return lambda;
+	}
+
+	/**
+	 * Compute Cholesky decomposition of symmetric positive definite
+	 * matrix A = LL^T.
+	 * @return Cholesky factor L of psd matrix A = L L^T
+	 */
+    public Matrix cholesky() {
+        if(!isSquare()) throw new RuntimeException("Matrix is not square");
+        if(!isSymmetric()) throw new RuntimeException("Matrix is not symmetric");
+        int i, j, k;
+        Matrix L = create(n,n);
+        for(i=0; i<n; i++){
+            for(j=0; j<=i; j++){
+                double sum = 0.0;
+                for(k=0; k<j; k++)
+                    sum += L.get(i, k) * L.get(j, k);
+                if(i == j) L.set(i, i, Math.sqrt(get(i, i) - sum));
+                else       L.set(i, j, 1.0 / L.get(j, j) * (get(i, j) - sum));
+            }
+            if(L.get(i, i) <= Maths.ε)
+                throw new RuntimeException("Matrix not positive definite");
+        }
+        return L;
+    }
+	
+	/**
+	 * @return Mirror matrix (orthogonal). Mirrored vector due to a plane defined through a vector v. 
+	 * @throws IllegalDimensionException
+	 */
+	public static Matrix mirror(Vector v) throws IllegalDimensionException{
+		// P = I - alpha v v'; alpha = 2 / (v'v)
+		return Diagonal.identity(v.n()).minus(v.tensorProduct(v).times(2/v.dot(v)));
+	}
+	
+	/**
+	 * Find value s in row.
+	 * @param s value
+	 * @param row
+	 * @return the index of the row where the value s was found
+	 */
+	public int findInRow(double s, int row){
+		int i;
+		int index = 0;
+		for(i=0; i<m; i++)
+			if(get(i,row)==s){
+				index = i;
+				break;
+			}
+		return index;
+	}
+
+	/**
+	 * @return maximum value in matrix, -∞ if no such value.
+	 */
+	public double max(){
+		return Maths.max(get());
+	}
+
+	/**
+	 * @return minimum value in matrix, +∞ if no such value.
+	 */
+	public double min(){
+		return Maths.min(get());
+	}
+	
+	/**
+	 * @return minimum and maximum value in array,
+	 * +∞ for minimum or -∞ for maximum if no such value
+	 */
+	public double[] minmax(){
+		return Maths.minmax(get());
+	}
+	
+	/**
+	 * sum each column
+	 * @return a row vector
+	 */
+	public Vector sum(){
+		int i,j;
+		double sum;
+		Vector c = create(getN());
+		for(j=0;j<getN();j++){
+			sum = 0;
+			for(i=0;i<getM();i++)
+				sum += get(i, j);
+			c.set(j, sum);
+		}
+		return c;
+	}
+
+	/**
+	 * m == n ?
+	 * @return true or false
+	 */
+	public boolean isSquare(){
+		return m == n;
+	}
+
+	/**
+	 * m == n ?
+	 * @param A double array
+	 * @return true or false
+	 */
+    public static boolean isSquare(double[][] A) {
+    	int i;
+        int N = A.length;
+        for(i=0; i<N; i++)
+            if(A[i].length != N) return false;
+        return true;
+    }
+	
+	/**
+	 * Values in euclidean plane &#x211D;<sup>2</sup>
+	 * 2x2 matrix
+	 * @return true or false
+	 */
+	public boolean isR2(){
+		return (m==2 && n==2);
+	}
+	
+	/**
+	 * Vectors or Matrix in euclidean plane &#x211D;<sup>2</sup>
+	 * 2xn matrix
+	 * @return true or false
+	 */
+	public boolean isM2(){
+		return (m==2);
+	}
+	
+	/**
+	 * Values in euclidean space &#x211D;<sup>3</sup>
+	 * 3x3 matrix
+	 * @return true or false
+	 */
+	public boolean isR3(){
+		return (m==3 && n==3);
+	}
+	
+	/**
+	 * Check if size(A) == size(B)
+	 * @throws IllegalDimensionException 
+	 **/
+	private void checkDimensions(Matrix b) throws IllegalDimensionException {
+		if(b.m != m || b.n != n)
+			throw new IllegalDimensionException("Two matrices have to be the same dimension.");
+	}
+	
+	/**
+	 * x<sub>ij</sub> == x<sub>ji</sub>?
+	 * @return true or false
+	 */
+	public boolean isSymmetric(){
+		int i, j;
+		for(i=0; i<m; i++)
+			for(j=i; j<n; j++)
+				if(i!=j && get(i,j)!=get(j,i))
+					return false;
+		return true;
+	}
+	
+	/**
+	 * a<sub>ij</sub> == a<sub>ji</sub>?
+	 * @param A double array
+	 * @return true or false
+	 */
+    public static boolean isSymmetric(double[][] A) {
+    	int i, j;
+        int N = A.length;
+        for(i=0; i<N; i++)
+            for(j=0; j<i; j++)
+                if(A[i][j] != A[j][i]) return false;
+        return true;
+    }
+	
+	/**
+	 * x<sub>i</sub>&#x22C5;x<sub>j</sub>=δ<sub>ij</sub>. x<sub>i</sub> ⊥ x<sub>j</sub> , i ≠ j
+	 * @return true or false
+	 * @throws SingularException 
+	 * @throws NoSquareException 
+	 * @throws IllegalDimensionException 
+	 */
+	public boolean isOrthogonal() throws IllegalDimensionException, NoSquareException, SingularException{
+		return (inverse().equals(transpose()));
+	}
+	
+	/**
+	 * x<sub>i</sub>&#x22C5;x<sub>j</sub> =δ<sub>ij</sub>
+	 * @return true or false
+	 * @throws IllegalDimensionException
+	 */
+	public boolean isOrthonormalized() throws IllegalDimensionException{
+		int i, j;
+		double x;
+		for(i=0;i<getN();i++)
+			for(j=i;j<getN();j++){
+				x = getN(i).dot(getN(j));
+				if((i == j) && x != 1) return false;
+				else if((i != j) && x != 0) return false;
+			}
+		return true;
+	}
+	
+	/**
+	 * check if the matrix is diagonal
+	 * @return true or false
+	 */
+	public boolean isDiagonal(){
+		if(this instanceof Diagonal)
+			return true;
+		else{
+			int i, j;
+			for(i=0; i<m; i++)
+				for(j=0; j<n; j++)
+					if(i!=j && get(i,j)!=0)
+						return false;
+			return true;
+		}
+	}
+
+	/**
+	 * <b>A</b> is <b>B</b>?
+	 * @param B matrix
+	 * @return true or false
+	 * @throws IllegalDimensionException
+	 */
+    public boolean equals(Matrix B) throws IllegalDimensionException{
+    	int i,j;
+        if (B.m != m || B.n != n)
+        	throw new IllegalDimensionException("Illegal matrix dimensions.");
+        for(i=0; i<m; i++)
+            for(j=0; j<n; j++)
+                if(get(i,j) != B.get(i,j)) return false;
+        return true;
+    }
+    
+    /**
+     * column to vector
+     * @return vector
+     * @throws IllegalDimensionException
+     */
+    //TODO: add rows to vector
+    public Vector toVector() throws IllegalDimensionException{
+    	if(m!=1 && n!=1) throw new IllegalDimensionException();
+    	if(m==1) return getM(0);
+    	else{
+        	int i;
+        	double[] a = new double[m];
+        	for(i=0; i<m; i++)
+        		a[i] = get(i,0);
+        	return create(a);
+    	}
+    }
+    
+    public Diagonal toDiagonal(){
+    	if(!isDiagonal() || !isSquare()) throw new IllegalArgumentException();
+		Diagonal result = new Diagonal(m);
+		int i;
+        for(i=0; i<n; i++)
+        	result.set(i, get(i, i));
+		return result;
+    }
+	
+    /**
+     * print matrix to standard output
+     */
+    @Override
+	public String toString(){
+    	int i,j, log;
+    	StringBuffer outputMatrix = new StringBuffer();
+//    	double[][] mat = get();
+    	double[] mat = get();
+    	if(mat!=null){
+    		int pre = 3, post = 0, main = 0;
+    		double max = Maths.max(mat);
+    		log = Maths.maxOrderOfMagnitude(mat);
+    		if(max == max){
+    			// TODO: Spaltenweise pre ausrechenen
+        		pre = (log>=0?log:0)+1 + (Maths.min(mat)>=0.?0:1); // max/10+1 = Zehnerpotenz+1; 0/+1 = Vorzeichen
+        		String[] zahlaufteilung;
+                for(i=0; i<m; i++)
+                    for(j=0; j<n; j++)
+                    	if(main < 3){
+	    		    		zahlaufteilung = String.valueOf(get(i,j)).split("[.]"); // Am Komma trennen
+	    		    		if(zahlaufteilung.length>1 && !zahlaufteilung[1].equals("0") && Double.valueOf(zahlaufteilung[1]) > main)
+	    		        		main = zahlaufteilung[1].length(); // max. Nachstellen
+                    	}
+                if(main != 0) pre++;	  // false = +1 für Komma
+        		post = main>3 ? 3 : main; // 4 = max Nachstellen in der Ausgabe; 
+    		}
+//    		System.out.println(pre);
+    		// begin with new line, so the first line will not translate
+    		outputMatrix.append("\u250C"+WString.generateStringWithLength((pre+post+1)*n+1, ' ')+"\u2510\n");
+            for(i=0; i<m; i++) {
+            	outputMatrix.append("\u2502 ");
+                for(j=0; j<n; j++){
+                	set(i,j, (get(i,j) <= Maths.ε && get(i,j) >= -Maths.ε) ? 0.0 : get(i,j));
+                	outputMatrix.append(String.format("%"+(pre+post)+"."+post+"f", get(i,j)));
+                	outputMatrix.append(" ");
+                }
+                outputMatrix.append("\u2502\n");
+            }
+            outputMatrix.append("\u2514"+WString.generateStringWithLength((pre+post+1)*n+1, ' ')+"\u2518");
+    	} else outputMatrix.append(mat);
+    	StringBuffer output = outputMatrix;
+		if(getName() != null) {
+			i = 0;
+			StringBuffer outputBuffer = new StringBuffer();
+			Scanner scanner = new Scanner(output.toString());
+			String line; // = scanner.nextLine(); // cut off first row
+			while(scanner.hasNextLine()) {
+				line = scanner.nextLine();
+				if(i==((int)(getM()/2)+1))
+					outputBuffer.append(getName() + " = ");
+				else outputBuffer.append(WString.generateStringWithLength(getName().length(), ' ') + "   ");
+				if(scanner.hasNextLine()) outputBuffer.append(line + "\n");
+				else  outputBuffer.append(line);
+				i++;
+			}
+			scanner.close();
+			output = outputBuffer;
+		}
+        return output.toString();
+    }
+    
+    public static String toString(int[][] array){
+    	return new Matrix(array).toString();
+    }
+    
+    public static String toString(double[][] array){
+    	return new Matrix(array).toString();
+    }
+    
+    public void println(){
+		System.out.println(toString());
+    }
+
+    /**
+     * Print the matrix to the "standard" output stream.
+     * Line the elements up in columns.
+     */
+    public void print(){
+    	int significantDigits = 4;
+    	double[] minmax = minmax();
+    	int maxlog = (int) (Math.log10(minmax[1]));
+    	int numerator = maxlog<0 ? maxlog - 1 : maxlog + 1;
+    	int denominator = Maths.decimalPlaces(this);
+    	int number = numerator+denominator;
+    	int d = denominator;
+    	if(denominator>0 && number>significantDigits) d = number - significantDigits;
+    	int w = numerator + (minmax[0]<0?1:0) + d; // (+1 sign), +1 space
+    	if(numerator<0) w = 2 + (minmax[0]<0?1:0) + d; // +2 : 0,
+//    	System.out.println(numerator+" "+denominator+" "+w+" "+d);
+    	print(new PrintWriter(System.out,true),w,d);
+    }
+
+    /**
+     * Print the matrix to the "standard" output stream.
+     * Line the elements up in columns.
+     * @param w column width
+     * @param d number of digits after the decimal
+     */
+    public void print(int w, int d){
+    	print(new PrintWriter(System.out,true),w,d);
+    }
+
+    /**
+     * Print the matrix to the chosen output stream.
+     * Line the elements up in columns.
+     * @param output stream
+     * @param w column width
+     * @param d number of digits after the decimal
+     */
+    public void print(PrintWriter output, int w, int d){
+    	DecimalFormat format = new DecimalFormat();
+    	format.setDecimalFormatSymbols(new DecimalFormatSymbols(Locale.getDefault()));
+    	format.setMinimumIntegerDigits(1);
+    	format.setMaximumFractionDigits(d);
+    	format.setMinimumFractionDigits(d);
+    	format.setGroupingUsed(false);
+    	print(output,format,w+2);
+    }
+
+    /**
+     * Print the matrix to "standard" output stream. Line the elements up in columns.
+     * Use the format object, and right justify within columns of width characters.
+     * @param format a formatting object for individual elements
+     * @param width field width for each column
+     * @see java.text.DecimalFormat#setDecimalFormatSymbols
+     */
+    public void print(NumberFormat format, int width){
+    	print(new PrintWriter(System.out,true),format,width);
+    }
+
+    /**
+     * Print the matrix to the output stream. Line the elements up in columns.
+     * Use the format object, and right justify within columns of width characters.
+     * @param output the output stream
+     * @param format a formatting object to format the matrix elements 
+     * @param width column width
+     * @see java.text.DecimalFormat#setDecimalFormatSymbols
+     */
+    public void print(PrintWriter output, NumberFormat format, int width){
+    	int i, j, k, padding;
+    	String s;
+    	output.println(); // start on new line.
+    	for(i=0; i<m; i++){
+    		for(j=0; j<n; j++){
+    			s = format.format(get(i,j)); // format the number
+    			padding = Math.max(1,width-s.length()); // At _least_ 1 space
+    			for(k=0; k<padding; k++)
+    				output.print(' ');
+    			output.print(s);
+    		}
+    		output.println();
+    	}
+    	output.println(); // end with blank line.
+    }
+    
+    /**
+     * @param A double array
+     */
+    public static void print(double[][] A){
+    	int i, j;
+        int M = A.length;
+        int N = A[0].length;
+        for(i=0; i<M; i++){
+            for(j=0; j<N; j++)
+                System.out.printf("%8.5f ", A[i][j]);
+            System.out.println();
+        }
+    }
+
+    /**
+     * Read a matrix from a stream. The format is the same as the print method,
+     * so printed matrices can be read back in. Elements are separated by
+     * whitespace, all the elements for each row appear on a single line,
+     * the last row is followed by a blank line.
+     * @param input the input stream
+     */
+    public static Matrix read(BufferedReader input) throws IOException {
+    	// Although StreamTokenizer will parse numbers, it doesn't recognize
+    	// scientific notation (E or D); however, Double.valueOf does.
+    	// The strategy here is to disable StreamTokenizer's number parsing.
+    	// We'll only get whitespace delimited words, EOL's and EOF's.
+    	// These words should all be numbers, for Double.valueOf to parse.
+    	StreamTokenizer tokenizer = new StreamTokenizer(input);
+    	tokenizer.resetSyntax();
+    	tokenizer.wordChars(0,255);
+    	tokenizer.whitespaceChars(0, ' ');
+    	tokenizer.eolIsSignificant(true);
+    	java.util.Vector<Double> vD = new java.util.Vector<Double>();
+
+    	// Ignore initial empty lines
+    	while(tokenizer.nextToken() == StreamTokenizer.TT_EOL);
+    	if(tokenizer.ttype == StreamTokenizer.TT_EOF)
+    		throw new IOException("Unexpected EOF on matrix read.");
+    	do vD.addElement(Double.valueOf(tokenizer.sval)); // Read & store 1st row.
+    	while(tokenizer.nextToken() == StreamTokenizer.TT_WORD);
+
+    	int n = vD.size();  // Now we've got the number of columns!
+    	double row[] = new double[n];
+    	int j;
+    	for(j=0; j<n; j++)  // extract the elements of the 1st row.
+    		row[j] = vD.elementAt(j).doubleValue();
+    	java.util.Vector<double[]> v = new java.util.Vector<double[]>();
+    	v.addElement(row);  // Start storing rows instead of columns.
+    	while (tokenizer.nextToken() == StreamTokenizer.TT_WORD){
+    		// While non-empty lines
+    		v.addElement(row = new double[n]);
+    		j = 0;
+    		do {
+    			if(j >= n) throw new IOException("Row " + v.size() + " is too long.");
+    			row[j++] = Double.valueOf(tokenizer.sval).doubleValue();
+    		} while(tokenizer.nextToken() == StreamTokenizer.TT_WORD);
+    		if(j < n) throw new IOException("Row " + v.size() + " is too short.");
+    	}
+    	int m = v.size();  // Now we've got the number of rows.
+    	double[][] A = new double[m][];
+    	v.copyInto(A);  // copy the rows out of the vector
+    	return new Matrix(A);
+    }
+    
+    /**
+     * Demo
+     * @param args
+     */
+    public static void main(String[] args){
+    	try {
+	    	Matrix.setSupressErrorMessage(true);
+	    	
+	    	Matrix A, B, C, D;
+	
+	        A = Matrix.random(5, 5);
+	        A.setName("A").println();
+	
+	        A.swap(1, 2);
+	        System.out.println("A.swap(1, 2)");
+	        A.println();
+	
+	        B = A.transpose();
+	        System.out.println("A' = \n" + B.setName("B"));
+	
+	        C = Matrix.identity(5);
+	        System.out.println("identity(5) = \n"+C.setName("C"));
+	
+	        System.out.println(A);
+	
+	        // shouldn't be equal since AB != BA in general    
+	        System.out.println("A*B =? B*A : \n"+A.times(B).equals(B.times(A)));
+	        System.out.println(A);
+	
+	        Matrix f = Matrix.random(5, 1);
+	        System.out.println("random(5, 1) = \n"+f.setName("f"));
+	
+	        Matrix x = A.solve(f);
+	        System.out.println("A^-1*f = \n"+x.setName("x"));
+	        
+	    	Matrix2D a = new Matrix2D(0,1,1,0);
+	    	System.out.println(a.setName("a"));
+	    	System.out.println("rot90(a) = \n"+a.rotate90());
+	    	Matrix b = Matrix.random(2, 1).times(5);
+	    	System.out.println("Matrix.random(2, 1).times(5) = \n"+b.setName("b"));
+	    	Matrix c = a.sumDirect(b);
+	    	System.out.println("c = a.sumDirect(b) = \n"+c);
+	        System.out.println("m="+c.m+", n="+c.n);
+	        Maths.print(c.get());
+	        Maths.print(c.getD());
+
+    		D = new Matrix(new double[][]{{1,2,3,4},{5,6,7,8}});
+	    	System.out.println("D = \n"+D);
+	    	System.out.println("D' = \n"+D.transpose()+" : "+(D.transpose().equals(new Matrix(new double[][]{{1,5},{2,6},{3,7},{4,8}}))));
+	    	System.out.println("D = \n"+D);
+	    	System.out.println("D'' = \n"+D.transpose().transpose()+" : "+(D.transpose().transpose().equals(D)));
+
+    		Matrix S = new Matrix(new double[][]{{3,1,2},{2,0,5},{1,2,3}});
+	    	System.out.println("S = \n"+S);
+			System.out.println("det(S) = " + S.det(1) + " Rank " +S.rank());
+
+	        Matrix T = new Matrix(new double[][]{{3,1,2},{2,0,5},{5,1,7}});
+	    	System.out.println("det(T) = " + T.det(1) + " Rank " +T.rank());
+	
+	        Matrix U = new Matrix(new double[][]{{3,1},{2,0}});
+	    	System.out.println("det(U) = " + U.det(1) + " Rank " +U.rank());
+	
+	        Matrix V = new Matrix(new double[][]{
+	        		{1,-2,2,3,3},
+	        		{4,3,1,-1,0},
+	        		{5,-1,1,2,7},
+	        		{-1,0,0,1,6},
+	        		{1,10,3,4,3}});
+	    	System.out.println("det(V) = " + V.det(1) + " Rank " +V.rank());
+	
+	        Matrix W = new Matrix(new double[][]{
+	        		{1,-2,2,3,3,11},
+	        		{4,3,1,-1,0,6},
+	        		{5,-1,1,2,7,1},
+	        		{-1,0,0,1,6,5},
+	        		{1,10,3,4,3,10},
+	        		{1,12,14,1,2,6}});
+	    	System.out.println("det(W) = " + W.det(1) + " Rank " +W.rank());
+	
+	        Matrix W2 = new Matrix(new double[][]{
+	        		{1,-2,2,3,3,11},
+	        		{4,3,1,-1,0,6},
+	        		{5,-1,1,2,7,1},
+	        		{-1,0,0,1,6,5},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10}});
+	    	System.out.println("det(W2) = " + W2.det(1) + " Rank " +W2.rank());
+	
+	        Matrix W3 = new Matrix(new double[][]{
+	        		{1,-2,2,3,3,11},
+	        		{4,3,1,-1,0,6},
+	        		{5,-1,1,2,7,1},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10}});
+	    	System.out.println("det(W3) = " + W3.det(1) + " Rank " +W3.rank());
+	
+	        Matrix W4 = new Matrix(new double[][]{
+	        		{1,-2,2,3,3,11},
+	        		{4,3,1,-1,0,6},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10}});
+	    	System.out.println("det(W4) = " + W4.det(1) + " Rank " +W4.rank());
+	
+	        Matrix W5 = new Matrix(new double[][]{
+	        		{1,-2,2,3,3,11},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10}});
+	    	System.out.println("det(W5) = " + W5.det(1) + " Rank " +W5.rank());
+	
+	        Matrix W6 = new Matrix(new double[][]{
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10},
+	        		{1,10,3,4,3,10}});
+	    	System.out.println("det(W6) = " + W6.det(1) + " Rank " +W6.rank());
+	    	try{
+	    		System.out.println("inv(W6) = \n"+W6.inverse());
+	    	} catch(Exception e){
+	    		System.err.println("singulär");
+	    	}
+	    	
+	        Matrix E = new Matrix(new double[][]{ { 5 } });
+	        System.out.println(E.setName("E"));
+	    	System.out.println("det(E) = " + E.det(1) + " Rank " +E.rank());
+	        System.out.println("inv(E) = \n"+E.inverse());
+	    	
+	        Matrix E2 = new Matrix(new double[][]{ { 1, 2 }, { 3, 4 } });
+	        System.out.println(E2.setName("E2"));
+	    	System.out.println("det(E2) = " + E2.det(1) + " Rank " +E2.rank());
+	        System.out.println("inv(E2) = \n"+E2.inverse());
+	    	
+	        Matrix E3 = new Matrix(new double[][]{ { 1, 2, 3 }, { 4, 5, 6 }, { 9, 1, 3} });
+	        System.out.println(E3.setName("E3"));
+	    	System.out.println("det(E3) = " + E3.det(1) + " Rank " +E3.rank());
+	        System.out.println("inv(E3) = \n"+E3.inverse());
+	        System.out.println("sum(E3) = \n"+E3.sum());
+	        System.out.println(E3);
+	        System.out.println("eigenvalue(E3) = \n"+E3.eigenvalues());
+	        
+	        Matrix[] F = new Matrix[3];
+	        F[0] = new Matrix(1);
+	        F[1] = new Matrix(2);
+	        F[2] = new Matrix(3);
+	        System.out.println(first(F));
+	        
+	    	T = new Matrix(new double[][]{{1,0,0},{0,1,0},{0,0,1}});
+	        System.out.println(T.setName("T"));
+	    	T = T.inverse();
+	        System.out.println(T.setName("T^-1"));
+	        
+	        B = new Matrix(new double[][]{ { 2, -1 }, { -1, 1 } });
+	        System.out.println(B.setName("B"));
+	        Vector lambda = B.eigenvaluesRe();
+	        System.out.println("Eigenvalues(B) = \n"+lambda);
+	        Matrix X = B.eigenvectors();
+	        System.out.println("Eigenvectors(B) = \n"+X);
+//	        Matrix Lambda = new Matrix(2,2).set(0, 0, lambda.get(0)).set(1, 1, lambda.get(0));
+//	        C = B.minus(Lambda);
+//	        System.out.println("B-l1E = "+C);
+//	        System.out.println("(B-l1E)^-1 = "+C.inverse());
+//	        Vector X1 = Vector.zeros(2).set(0,1);
+//	        System.out.println("X1 = "+X1);
+//	        System.out.println("Eigenvectors(B) = "+C.times(X1));
+//	        System.out.println("Eigenvectors(B) = "+C.inverse().times(X1));
+//	        System.out.println("Eigenvectors(B) = "+C.inverse().times(X1));
+
+	        Matrix2D P = new Matrix2D(new double[][]{ { 1, 1 }, { 0, -1 } });
+	        System.out.println(P.setName("P"));
+	        System.out.println("P' = \n"+P.transpose()+" : "+(P.transpose().equals(new Matrix(new double[][]{ { 1, 0 }, { 1, -1 } }))));
+	        System.out.println("P^-1 = \n" + P.inverse());
+	        System.out.println("P'P = \n"+P.transpose().times(P));
+	        System.out.println("rot(P'P,45°) = \n" + ((Matrix2D) P.transpose().times(P)).rotate45());
+
+	        P = new Matrix2D(new double[][]{ { 1 }, { 1 } });
+	        System.out.println(P.setName("P"));
+	        System.out.println("rot(P,45°) = \n" + P.rotate(45));
+
+	        Vector2D v = new Vector2D(1,1);
+	        Vector2D w = new Vector2D(0,1);
+	        P = (Matrix2D) Matrix2D.mirror(v);
+	        System.out.println("Pw = ");
+	        P.times(w).println();
+	        
+	        A = new Matrix2D(new double[][]{ { 2, 1 }, { 0, 1 } });
+	        System.out.println(A.setName("A"));
+	        B = A.transpose().times(A);
+	        System.out.println("A'A = \n" + B);
+	        X = B.eigenvectors();
+	        System.out.println("X = \n" + X);
+	        System.out.println("X'A'A = \n" + X.transpose().times(B));
+	        Matrix L2 = X.transpose().times(B).times(X);
+	        System.out.println("L^2 = X'A'AX = \n" + L2);
+	        // TODO: Wurzeln für Diagonalmatrix hinzufügen
+	        Matrix L = L2.times(new Matrix(new double[][]{ { 1/Math.sqrt(L2.get(0, 0)), 0 }, { 0, 1/Math.sqrt(L2.get(1, 1)) } }));
+	        System.out.println("L = \n" + L);
+	        S = X.times(L).times(X);
+	        System.out.println("S = XLX' = \n" + S);
+	        Matrix2D Q = (Matrix2D) A.times(S.inverse());
+	        System.out.println("Q = AS^-1 = \n" + Q);
+	        Q.plot();
+	    	
+	        A = new Matrix(new double[][]{ { 3, 2, 6 }, { 2, 2, 5 }, { -2, -1, -4} });
+	        System.out.println(A.setName("A"));
+	        System.out.println("eigenvalue(A) = \n"+A.eigenvaluesRe());
+	    	
+	        A = new Matrix(new double[][]{
+	        		{   38,  -43./7, -63./2, -1149./14 },
+	        		{   14,  -19./7,     -7,   -181./7 },
+	        		{ 8./7, 122./49, -24./7,  -177./49 },
+	        		{   16,  -26./7,    -13,   -244./7 } });
+	        System.out.println(A.setName("A"));
+	        System.out.println("eigenvalue(A) = \n"+A.eigenvalues());
+	    	
+	        A = new Matrix(new double[][]{
+	        		{   7, -3, -11,   2 },
+	        		{  -7, -5, -11,   5 },
+	        		{ -11,  0,  -2, -15 },
+	        		{ -14, 15,  -5,   7 } });
+	        System.out.println(A.setName("A"));
+	        System.out.println("eigenvalue(A) = \n"+A.eigenvaluesIm());
+	        
+	        System.out.println("power of:");
+	        A = new Matrix(new double[][]{ { 0, 0, 0 }, { 0, 1, 0 }, {0, 0, 3 } });
+	        System.out.println(A.setName("A"));
+	        System.out.println("A^2 = \n"+A.pow(2));
+//	        A = new Matrix(new double[][]{ { 2, 2, 0 }, { 0, 1, 0 }, {0, 0, 3 } });
+//	        System.out.println(A.setName("A"));
+//	        System.out.println("A^1 = "+A.pow(1));
+	    	
+	        System.out.println("\nKronecker product:");
+	        A = new Matrix(new double[][]{ { 1, 2 }, { 3, 4 } });
+	        System.out.println(A.setName("A"));
+	        B = new Matrix(new double[][]{ { 0, 5 }, { 6, 7 } });
+	        System.out.println(B.setName("B"));
+	        C = A.kroneckerProduct(B);
+	        System.out.println("A \u24E7 B = \n"+C);
+	        D = new Matrix(new double[][]{ { 0, 5, 0, 10 }, { 6, 7, 12, 14 }, { 0, 15, 0, 20 }, { 18, 21, 24, 28 } });
+	        if(C.equals(D)) System.out.println("test successful");
+
+	        A.times(0.00001).print();
+	        A.times(0.001).print();
+	        A.times(10).print();
+	        A.times(25).print();
+	        A.times(250000).print();
+	        
+	        System.out.println(A);
+	        Maths.print(new Matrix(A).get());
+	        System.out.println(new Matrix(A.get()));
+	        System.out.println(new Matrix(A.n, A.get()));
+	        System.out.println(new Matrix(A.get(), A.m));
+	        
+	        D = new Diagonal(2.,3);
+	        System.out.println(A);
+	        System.out.println(D);
+	        System.out.println(A.plus(D));
+	        
+	        Matrix I = Matrix.identity(2);
+	        System.out.println(I);
+	        System.out.println(I instanceof Diagonal);
+	        I = I.toDiagonal();
+	        System.out.println(I instanceof Diagonal);
+	        
+	        Diagonal diag = new Diagonal(2.,3);
+	        System.out.println(diag);
+	        System.out.println(A.plus(diag));
+
+	        boolean speedtest = false;
+	        if(speedtest){
+		        A = Matrix.random(1000,500);
+		        B = Matrix.random(500,1000);
+		        System.out.println("C = A B");
+		        long startTime, endTime, duration;
+		        // V1
+		        System.out.println("V1");
+		        {
+			        startTime = System.nanoTime();
+			        C = A.times(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+		        {
+			        startTime = System.nanoTime();
+			        C = A.times(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+		        {
+			        startTime = System.nanoTime();
+			        C = A.times(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+		        {
+			        startTime = System.nanoTime();
+			        C = A.times(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+		        {
+			        startTime = System.nanoTime();
+			        C = A.times(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+	
+		        // V2
+		        System.out.println("V2");
+		        {
+			        startTime = System.nanoTime();
+			        C = A.timesV2(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+		        {
+			        startTime = System.nanoTime();
+			        C = A.timesV2(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+		        {
+			        startTime = System.nanoTime();
+			        C = A.timesV2(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+		        {
+			        startTime = System.nanoTime();
+			        C = A.timesV2(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+		        {
+			        startTime = System.nanoTime();
+			        C = A.timesV2(B);
+			        endTime = System.nanoTime();
+			        duration = (endTime - startTime)/1000000;
+			        System.out.println("That took " + duration + " milliseconds");
+		        }
+	        }
+
+	        A = new Matrix(new double[][]{ { 4, 1,  1 },
+	                         { 1, 5,  3 },
+	                         { 1, 3, 15 } });
+	        L = A.cholesky();
+	        L.println();
+	        /*
+	         * 2.00000  0.00000  0.00000 
+	         * 0.50000  2.17945  0.00000 
+	         * 0.50000  1.26179  3.62738 
+	         */
+
+	        
+		} catch (Exception e) {
+			e.printStackTrace();
+		}
+
+    }
+}