diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
index d53b30e713e..0d11f5cd2c0 100644
--- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
+++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
@@ -17,803 +17,1146 @@
  * under the License.
  */
 
-package org.apache.sysds.runtime.matrix.data;
-
-import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.fft;
-import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.fft_linearized;
-import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.ifft;
-import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.ifft_linearized;
-
-import org.apache.commons.logging.Log;
-import org.apache.commons.logging.LogFactory;
-import org.apache.commons.math3.exception.MaxCountExceededException;
-import org.apache.commons.math3.linear.Array2DRowRealMatrix;
-import org.apache.commons.math3.linear.BlockRealMatrix;
-import org.apache.commons.math3.linear.CholeskyDecomposition;
-import org.apache.commons.math3.linear.DecompositionSolver;
-import org.apache.commons.math3.linear.EigenDecomposition;
-import org.apache.commons.math3.linear.LUDecomposition;
-import org.apache.commons.math3.linear.QRDecomposition;
-import org.apache.commons.math3.linear.RealMatrix;
-import org.apache.commons.math3.linear.SingularValueDecomposition;
-import org.apache.sysds.runtime.DMLRuntimeException;
-import org.apache.sysds.runtime.data.DenseBlock;
-import org.apache.sysds.runtime.functionobjects.Builtin;
-import org.apache.sysds.runtime.functionobjects.Divide;
-import org.apache.sysds.runtime.functionobjects.MinusMultiply;
-import org.apache.sysds.runtime.functionobjects.Multiply;
-import org.apache.sysds.runtime.functionobjects.SwapIndex;
-import org.apache.sysds.runtime.instructions.InstructionUtils;
-import org.apache.sysds.runtime.matrix.operators.AggregateBinaryOperator;
-import org.apache.sysds.runtime.matrix.operators.BinaryOperator;
-import org.apache.sysds.runtime.matrix.operators.LeftScalarOperator;
-import org.apache.sysds.runtime.matrix.operators.ReorgOperator;
-import org.apache.sysds.runtime.matrix.operators.RightScalarOperator;
-import org.apache.sysds.runtime.matrix.operators.ScalarOperator;
-import org.apache.sysds.runtime.matrix.operators.TernaryOperator;
-import org.apache.sysds.runtime.matrix.operators.UnaryOperator;
-import org.apache.sysds.runtime.util.DataConverter;
-
-/**
- * Library for matrix operations that need invocation of 
- * Apache Commons Math library. 
- * 
- * This library currently supports following operations:
- * matrix inverse, matrix decompositions (QR, LU, Eigen), solve 
- */
-public class LibCommonsMath 
-{
-	private static final Log LOG = LogFactory.getLog(LibCommonsMath.class.getName());
-	private static final double RELATIVE_SYMMETRY_THRESHOLD = 1e-6;
-	private static final double EIGEN_LAMBDA = 1e-8;
-
-	private LibCommonsMath() {
-		//prevent instantiation via private constructor
-	}
-	
-	public static boolean isSupportedUnaryOperation( String opcode ) {
-		return ( opcode.equals("inverse") || opcode.equals("cholesky") );
-	}
-	
-	public static boolean isSupportedMultiReturnOperation( String opcode ) {
-
-		switch (opcode) {
-			case "qr":
-			case "lu":
-			case "eigen":
-			case "fft":
-			case "ifft":
-			case "fft_linearized":
-			case "ifft_linearized":
-			case "stft":
-			case "svd": return true;
-			default: return false;
-		}
-
-	}
-	
-	public static boolean isSupportedMatrixMatrixOperation( String opcode ) {
-		return ( opcode.equals("solve") );
-	}
-		
-	public static MatrixBlock unaryOperations(MatrixBlock inj, String opcode) {
-		Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(inj);
-		if(opcode.equals("inverse"))
-			return computeMatrixInverse(matrixInput);
-		else if (opcode.equals("cholesky"))
-			return computeCholesky(matrixInput);
-		return null;
-	}
-
-	// public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode) {
-	// 	return multiReturnOperations(in, opcode, 1, 1);
-	// }
-
-	public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode, int threads) {
-		return multiReturnOperations(in, opcode, threads, 1);
-	}
-
-	public static MatrixBlock[] multiReturnOperations(MatrixBlock in1, MatrixBlock in2, String opcode) {
-		return multiReturnOperations(in1, in2, opcode, 1, 1);
-	}
-
-	public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode, int threads, int num_iterations, double tol) {
-		if(opcode.equals("eigen_qr"))
-			return computeEigenQR(in, num_iterations, tol, threads);
-		else
-			return multiReturnOperations(in, opcode, threads, 1);
-	}
-
-	public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode, int threads, long seed) {
-		if(opcode.equals("qr"))
-			return computeQR(in);
-		else if (opcode.equals("qr2"))
-			return computeQR2(in, threads);
-		else if (opcode.equals("lu"))
-			return computeLU(in);
-		else if (opcode.equals("eigen"))
-			return computeEigen(in);
-		else if (opcode.equals("eigen_lanczos"))
-			return computeEigenLanczos(in, threads, seed);
-		else if (opcode.equals("eigen_qr"))
-			return computeEigenQR(in, threads);
-		else if (opcode.equals("svd"))
-			return computeSvd(in);
-		else if (opcode.equals("fft"))
-			return computeFFT(in, threads);
-		else if (opcode.equals("ifft"))
-			return computeIFFT(in, threads);
-		else if (opcode.equals("fft_linearized"))
-			return computeFFT_LINEARIZED(in, threads);
-		else if (opcode.equals("ifft_linearized"))
-			return computeIFFT_LINEARIZED(in, threads);
-		return null;
-	}
-
-	public static MatrixBlock[] multiReturnOperations(MatrixBlock in1, MatrixBlock in2, String opcode, int threads,
-			long seed) {
-
-		switch (opcode) {
-			case "ifft":
-				return computeIFFT(in1, in2, threads);
-			case "ifft_linearized":
-				return computeIFFT_LINEARIZED(in1, in2, threads);
-			default:
-				return null;
-		}
-
-	}
-
-	public static MatrixBlock matrixMatrixOperations(MatrixBlock in1, MatrixBlock in2, String opcode) {
-		if(opcode.equals("solve")) {
-			if (in1.getNumRows() != in1.getNumColumns())
-				throw new DMLRuntimeException("The A matrix, in solve(A,b) should have squared dimensions.");
-			return computeSolve(in1, in2);
-		}
-		return null;
-	}
-	
-	/**
-	 * Function to solve a given system of equations.
-	 * 
-	 * @param in1 matrix object 1
-	 * @param in2 matrix object 2
-	 * @return matrix block
-	 */
-	private static MatrixBlock computeSolve(MatrixBlock in1, MatrixBlock in2) {
-		//convert to commons math BlockRealMatrix instead of Array2DRowRealMatrix
-		//to avoid unnecessary conversion as QR internally creates a BlockRealMatrix
-		BlockRealMatrix matrixInput = DataConverter.convertToBlockRealMatrix(in1);
-		BlockRealMatrix vectorInput = DataConverter.convertToBlockRealMatrix(in2);
-		
-		/*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput);
-		DecompositionSolver lusolver = ludecompose.getSolver();
-		RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/
-		
-		// Setup a solver based on QR Decomposition
-		QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
-		DecompositionSolver solver = qrdecompose.getSolver();
-		// Invoke solve
-		RealMatrix solutionMatrix = solver.solve(vectorInput);
-		
-		return DataConverter.convertToMatrixBlock(solutionMatrix);
-	}
-	
-	/**
-	 * Function to perform QR decomposition on a given matrix.
-	 * 
-	 * @param in matrix object
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeQR(MatrixBlock in) {
-		Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
-		
-		// Perform QR decomposition
-		QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
-		RealMatrix H = qrdecompose.getH();
-		RealMatrix R = qrdecompose.getR();
-		
-		// Read the results into native format
-		MatrixBlock mbH = DataConverter.convertToMatrixBlock(H.getData());
-		MatrixBlock mbR = DataConverter.convertToMatrixBlock(R.getData());
-
-		return new MatrixBlock[] { mbH, mbR };
-	}
-	
-	/**
-	 * Function to perform LU decomposition on a given matrix.
-	 * 
-	 * @param in matrix object
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeLU(MatrixBlock in) {
-		if(in.getNumRows() != in.getNumColumns()) {
-			throw new DMLRuntimeException(
-				"LU Decomposition can only be done on a square matrix. Input matrix is rectangular (rows="
-					+ in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
-		}
-		
-		Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
-		
-		// Perform LUP decomposition
-		LUDecomposition ludecompose = new LUDecomposition(matrixInput);
-		RealMatrix P = ludecompose.getP();
-		RealMatrix L = ludecompose.getL();
-		RealMatrix U = ludecompose.getU();
-		
-		// Read the results into native format
-		MatrixBlock mbP = DataConverter.convertToMatrixBlock(P.getData());
-		MatrixBlock mbL = DataConverter.convertToMatrixBlock(L.getData());
-		MatrixBlock mbU = DataConverter.convertToMatrixBlock(U.getData());
-
-		return new MatrixBlock[] { mbP, mbL, mbU };
-	}
-	
-	/**
-	 * Function to perform Eigen decomposition on a given matrix.
-	 * Input must be a symmetric matrix.
-	 * 
-	 * @param in matrix object
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeEigen(MatrixBlock in) {
-		if ( in.getNumRows() != in.getNumColumns() ) {
-			throw new DMLRuntimeException("Eigen Decomposition can only be done on a square matrix. "
-				+ "Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
-		}
-		
-		EigenDecomposition eigendecompose = null;
-		try {
-			Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
-			eigendecompose = new EigenDecomposition(matrixInput);
-		}
-		catch(MaxCountExceededException ex) {
-			LOG.warn("Eigen: "+ ex.getMessage()+". Falling back to regularized eigen factorization.");
-			eigendecompose = computeEigenRegularized(in);
-		}
-		
-		RealMatrix eVectorsMatrix = eigendecompose.getV();
-		double[][] eVectors = eVectorsMatrix.getData();
-		double[] eValues = eigendecompose.getRealEigenvalues();
-
-		return sortEVs(eValues, eVectors);
-	}
-
-	private static EigenDecomposition computeEigenRegularized(MatrixBlock in) {
-		if( in == null || in.isEmptyBlock(false) )
-			throw new DMLRuntimeException("Invalid empty block");
-		
-		//slightly modify input for regularization (pos/neg)
-		MatrixBlock in2 = new MatrixBlock(in, false);
-		DenseBlock a = in2.getDenseBlock();
-		for( int i=0; i<in2.rlen; i++ ) {
-			double[] avals = a.values(i);
-			int apos = a.pos(i);
-			for( int j=0; j<in2.clen; j++ ) {
-				double v = avals[apos+j];
-				avals[apos+j] += Math.signum(v) * EIGEN_LAMBDA;
-			}
-		}
-		
-		//run eigen decomposition
-		return new EigenDecomposition(
-			DataConverter.convertToArray2DRowRealMatrix(in2));
-	}
-
-	/**
-	 * Function to perform FFT on a given matrix.
-	 *
-	 * @param re      matrix object
-	 * @param threads number of threads
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeFFT(MatrixBlock re, int threads) {
-		if(re == null)
-			throw new DMLRuntimeException("Invalid empty block");
-		if(re.isEmptyBlock(false)) {
-			// Return the original matrix as the result
-			// Assuming you need to return two matrices: the real part and an imaginary part initialized to 0.
-			return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)};
-		}
-		// run fft
-		re.sparseToDense();
-		return fft(re, threads);
-	}
-
-	/**
-	 * Function to perform IFFT on a given matrix.
-	 *
-	 * @param re      matrix object
-	 * @param im      matrix object
-	 * @param threads number of threads
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeIFFT(MatrixBlock re, MatrixBlock im, int threads) {
-		if(re == null)
-			throw new DMLRuntimeException("Invalid empty block");
-
-		// run ifft
-		if(im != null && !im.isEmptyBlock(false)) {
-			re.sparseToDense();
-			im.sparseToDense();
-			return ifft(re, im, threads);
-		}
-		else {
-			if(re.isEmptyBlock(false)) {
-				// Return the original matrix as the result
-				// Assuming you need to return two matrices: the real part and an imaginary part initialized to 0.
-				return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)};
-			}
-			re.sparseToDense();
-			return ifft(re, threads);
-		}
-	}
-
-	/**
-	 * Function to perform IFFT on a given matrix.
-	 *
-	 * @param re      matrix object
-	 * @param threads number of threads
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeIFFT(MatrixBlock re, int threads) {
-		return computeIFFT(re, null, threads);
-	}
-
-	/**
-	 * Function to perform FFT_LINEARIZED on a given matrix.
-	 *
-	 * @param re      matrix object
-	 * @param threads number of threads
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeFFT_LINEARIZED(MatrixBlock re, int threads) {
-		if(re == null)
-			throw new DMLRuntimeException("Invalid empty block");
-		if(re.isEmptyBlock(false)) {
-			// Return the original matrix as the result
-			// Assuming you need to return two matrices: the real part and an imaginary part initialized to 0.
-			return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)};
-		}
-		// run fft
-		re.sparseToDense();
-		return fft_linearized(re, threads);
-	}
-
-	/**
-	 * Function to perform IFFT_LINEARIZED on a given matrix.
-	 *
-	 * @param re      matrix object
-	 * @param im      matrix object
-	 * @param threads number of threads
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeIFFT_LINEARIZED(MatrixBlock re, MatrixBlock im, int threads) {
-		if(re == null)
-			throw new DMLRuntimeException("Invalid empty block");
-
-		// run ifft
-		if(im != null && !im.isEmptyBlock(false)) {
-			re.sparseToDense();
-			im.sparseToDense();
-			return ifft_linearized(re, im, threads);
-		}
-		else {
-			if(re.isEmptyBlock(false)) {
-				// Return the original matrix as the result
-				// Assuming you need to return two matrices: the real part and an imaginary part initialized to 0.
-				return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)};
-			}
-			re.sparseToDense();
-			return ifft_linearized(re, threads);
-		}
-	}
-
-	/**
-	 * Function to perform IFFT_LINEARIZED on a given matrix
-	 * 
-	 * @param re      matrix object
-	 * @param threads number of threads
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeIFFT_LINEARIZED(MatrixBlock re, int threads) {
-		return computeIFFT_LINEARIZED(re, null, threads);
-	}
-
-
-	// /**
-	//  * Function to perform STFT on a given matrix.
-	//  *
-	//  * @param re matrix object
-	//  * @param im matrix object
-	//  * @param windowSize of stft
-	//  * @param overlap of stft
-	//  * @return array of matrix blocks
-	//  */
-	// private static MatrixBlock[] computeSTFT(MatrixBlock re, MatrixBlock im, int windowSize, int overlap, int threads) {
-	// 	if (re == null) {
-	// 		throw new DMLRuntimeException("Invalid empty block");
-	// 	} else if (im != null && !im.isEmptyBlock(false)) {
-	// 		re.sparseToDense();
-	// 		im.sparseToDense();
-	// 		return stft(re, im, windowSize, overlap, threads);
-	// 	} else {
-	// 		if (re.isEmptyBlock(false)) {
-	// 			// Return the original matrix as the result
-	// 			int rows = re.getNumRows();
-	// 			int cols = re.getNumColumns();
-
-	// 			int stepSize = windowSize - overlap;
-	// 			if (stepSize == 0) {
-	// 				throw new IllegalArgumentException("windowSize - overlap is zero");
-	// 			}
-
-	// 			int numberOfFramesPerRow = (cols - overlap + stepSize - 1) / stepSize;
-	// 			int rowLength= numberOfFramesPerRow * windowSize;
-	// 			int out_len = rowLength * rows;
-
-	// 			double[] out_zero = new double[out_len];
-
-	// 			return new MatrixBlock[]{new MatrixBlock(rows, rowLength, out_zero), new MatrixBlock(rows, rowLength, out_zero)};
-	// 			}
-	// 		re.sparseToDense();
-	// 		return stft(re, windowSize, overlap, threads);
-	// 	}
-	// }
-
-	// /**
-	//  * Function to perform STFT on a given matrix.
-	//  *
-	//  * @param re matrix object
-	//  * @return array of matrix blocks
-	//  */
-	// private static MatrixBlock[] computeSTFT(MatrixBlock re, int windowSize, int overlap, int threads) {
-	// 	return computeSTFT(re, null, windowSize, overlap, threads);
-	// }
-
-	/**
-	 * Performs Singular Value Decomposition. Calls Apache Commons Math SVD.
-	 * X = U * Sigma * Vt, where X is the input matrix,
-	 * U is the left singular matrix, Sigma is the singular values matrix returned as a
-	 * column matrix and Vt is the transpose of the right singular matrix V.
-	 * However, the returned array has  { U, Sigma, V}
-	 * 
-	 * @param in Input matrix
-	 * @return An array containing U, Sigma & V
-	 */
-	private static MatrixBlock[] computeSvd(MatrixBlock in) {
-		Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
-
-		SingularValueDecomposition svd = new SingularValueDecomposition(matrixInput);
-		double[] sigma = svd.getSingularValues();
-		RealMatrix u = svd.getU();
-		RealMatrix v = svd.getV();
-		MatrixBlock U = DataConverter.convertToMatrixBlock(u.getData());
-		MatrixBlock Sigma = DataConverter.convertToMatrixBlock(sigma, true);
-		Sigma = LibMatrixReorg.diag(Sigma, new MatrixBlock(Sigma.rlen, Sigma.rlen, true));
-		MatrixBlock V = DataConverter.convertToMatrixBlock(v.getData());
-
-		return new MatrixBlock[] { U, Sigma, V };
-	}
-	
-	/**
-	 * Function to compute matrix inverse via matrix decomposition.
-	 * 
-	 * @param in commons-math3 Array2DRowRealMatrix
-	 * @return matrix block
-	 */
-	private static MatrixBlock computeMatrixInverse(Array2DRowRealMatrix in) {
-		if(!in.isSquare())
-			throw new DMLRuntimeException("Input to inv() must be square matrix -- given: a " + in.getRowDimension()
-				+ "x" + in.getColumnDimension() + " matrix.");
-
-		QRDecomposition qrdecompose = new QRDecomposition(in);
-		DecompositionSolver solver = qrdecompose.getSolver();
-		RealMatrix inverseMatrix = solver.getInverse();
-
-		return DataConverter.convertToMatrixBlock(inverseMatrix.getData());
-	}
-
-	/**
-	 * Function to compute Cholesky decomposition of the given input matrix. 
-	 * The input must be a real symmetric positive-definite matrix.
-	 * 
-	 * @param in commons-math3 Array2DRowRealMatrix
-	 * @return matrix block
-	 */
-	private static MatrixBlock computeCholesky(Array2DRowRealMatrix in) {
-		if(!in.isSquare())
-			throw new DMLRuntimeException("Input to cholesky() must be square matrix -- given: a "
-				+ in.getRowDimension() + "x" + in.getColumnDimension() + " matrix.");
-		CholeskyDecomposition cholesky = new CholeskyDecomposition(in, RELATIVE_SYMMETRY_THRESHOLD,
-			CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
-		RealMatrix rmL = cholesky.getL();
-		return DataConverter.convertToMatrixBlock(rmL.getData());
-	}
-
-	/**
-	 * Creates a random normalized vector with dim elements.
-	 *
-	 * @param dim dimension of the created vector
-	 * @param threads number of threads
-	 * @param seed seed for the random MatrixBlock generation
-	 * @return random normalized vector
-	 */
-	private static MatrixBlock randNormalizedVect(int dim, int threads, long seed) {
-		MatrixBlock v1 = MatrixBlock.randOperations(dim, 1, 1.0, 0.0, 1.0, "UNIFORM", seed);
-
-		double v1_sum = v1.sum();
-		ScalarOperator op_div_scalar = new RightScalarOperator(Divide.getDivideFnObject(), v1_sum, threads);
-		v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
-		UnaryOperator op_sqrt = new UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT), threads, true);
-		v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
-
-		if(Math.abs(v1.sumSq() - 1.0) >= 1e-7)
-			throw new DMLRuntimeException("v1 not correctly normalized (maybe try changing the seed)");
-
-		return v1;
-	}
-
-	/**
-	 * Function to perform the Lanczos algorithm and then computes the Eigendecomposition.
-	 * Caution: Lanczos is not numerically stable (see https://en.wikipedia.org/wiki/Lanczos_algorithm)
-	 * Input must be a symmetric (and square) matrix.
-	 *
-	 * @param in matrix object
-	 * @param threads number of threads
-	 * @param seed seed for the random MatrixBlock generation
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeEigenLanczos(MatrixBlock in, int threads, long seed) {
-		if(in.getNumRows() != in.getNumColumns()) {
-			throw new DMLRuntimeException(
-				"Lanczos algorithm and Eigen Decomposition can only be done on a square matrix. "
-					+ "Input matrix is rectangular (rows=" + in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
-		}
-
-		int m = in.getNumRows();
-		MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
-		MatrixBlock v1 = randNormalizedVect(m, threads, seed);
-
-		MatrixBlock T = new MatrixBlock(m, m, 0.0);
-		MatrixBlock TV = new MatrixBlock(m, m, 0.0);
-		MatrixBlock w1;
-
-		ReorgOperator op_t = new ReorgOperator(SwapIndex.getSwapIndexFnObject(), threads);
-		TernaryOperator op_minus_mul = new TernaryOperator(MinusMultiply.getFnObject(), threads);
-		AggregateBinaryOperator op_mul_agg = InstructionUtils.getMatMultOperator(threads);
-		ScalarOperator op_div_scalar = new RightScalarOperator(Divide.getDivideFnObject(), 1, threads);
-
-		MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
-		for(int i = 0; i < m; i++) {
-			v1.putInto(TV, 0, i, false);
-			w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
-			MatrixBlock alpha = w1.aggregateBinaryOperations(v1.reorgOperations(op_t, new MatrixBlock(), 0, 0, m), w1, op_mul_agg);
-			if(i < m - 1) {
-				w1 = w1.ternaryOperations(op_minus_mul, v1, alpha, new MatrixBlock());
-				w1 = w1.ternaryOperations(op_minus_mul, v0, beta, new MatrixBlock());
-				beta.set(0, 0, Math.sqrt(w1.sumSq()));
-				v0.copy(v1);
-				op_div_scalar = op_div_scalar.setConstant(beta.getDouble(0, 0));
-				w1.scalarOperations(op_div_scalar, v1);
-
-				T.set(i + 1, i, beta.get(0, 0));
-				T.set(i, i + 1, beta.get(0, 0));
-			}
-			T.set(i, i, alpha.get(0, 0));
-		}
-
-		MatrixBlock[] e = computeEigen(T);
-		TV.setNonZeros((long) m*m);
-		e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
-		return e;
-	}
-
-	/**
-	 * Function to perform the QR decomposition.
-	 * Input must be a square matrix.
-	 * TODO: use Householder transformation and implicit shifts to further speed up QR decompositions
-	 *
-	 * @param in matrix object
-	 * @param threads number of threads
-	 * @return array of matrix blocks [Q, R]
-	 */
-	private static MatrixBlock[] computeQR2(MatrixBlock in, int threads) {
-		if(in.getNumRows() != in.getNumColumns()) {
-			throw new DMLRuntimeException("QR2 Decomposition can only be done on a square matrix. "
-				+ "Input matrix is rectangular (rows=" + in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
-		}
-
-		int m = in.rlen;
-
-		MatrixBlock A_n = new MatrixBlock();
-		A_n.copy(in);
-
-		MatrixBlock Q_n = new MatrixBlock(m, m, true);
-		for(int i = 0; i < m; i++) {
-			Q_n.set(i, i, 1.0);
-		}
-
-		ReorgOperator op_t = new ReorgOperator(SwapIndex.getSwapIndexFnObject(), threads);
-		AggregateBinaryOperator op_mul_agg = InstructionUtils.getMatMultOperator(threads);
-		BinaryOperator op_sub = InstructionUtils.parseExtendedBinaryOperator("-");
-		ScalarOperator op_div_scalar = new RightScalarOperator(Divide.getDivideFnObject(), 1, threads);
-		ScalarOperator op_mult_2 = new LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, threads);
-
-		for(int k = 0; k < m; k++) {
-			MatrixBlock z = A_n.slice(k, m - 1, k, k);
-			MatrixBlock uk = new MatrixBlock(m - k, 1, 0.0);
-			uk.copy(z);
-			uk.set(0, 0, uk.get(0, 0) + Math.signum(z.get(0, 0)) * Math.sqrt(z.sumSq()));
-			op_div_scalar = op_div_scalar.setConstant(Math.sqrt(uk.sumSq()));
-			uk = uk.scalarOperations(op_div_scalar, new MatrixBlock());
-
-			MatrixBlock vk = new MatrixBlock(m, 1, 0.0);
-			vk.copy(k, m - 1, 0, 0, uk, true);
-
-			MatrixBlock vkt = vk.reorgOperations(op_t, new MatrixBlock(), 0, 0, m);
-			MatrixBlock vkt2 = vkt.scalarOperations(op_mult_2, new MatrixBlock());
-			MatrixBlock vkvkt2 = vk.aggregateBinaryOperations(vk, vkt2, op_mul_agg);
-
-			A_n = A_n.binaryOperations(op_sub, A_n.aggregateBinaryOperations(vkvkt2, A_n, op_mul_agg));
-			Q_n = Q_n.binaryOperations(op_sub, Q_n.aggregateBinaryOperations(Q_n, vkvkt2, op_mul_agg));
-		}
-		// QR decomp: Q: Q_n; R: A_n
-		return new MatrixBlock[] {Q_n, A_n};
-	}
-
-	/**
-	 * Function that computes the Eigen Decomposition using the QR algorithm.
-	 * Caution: check if the QR algorithm is converged, if not increase iterations
-	 * Caution: if the input matrix has complex eigenvalues results will be incorrect
-	 *
-	 * @param in Input matrix
-	 * @param threads number of threads
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] computeEigenQR(MatrixBlock in, int threads) {
-		return computeEigenQR(in, 100, 1e-10, threads);
-	}
-
-	private static MatrixBlock[] computeEigenQR(MatrixBlock in, int num_iterations, double tol, int threads) {
-		if(in.getNumRows() != in.getNumColumns()) {
-			throw new DMLRuntimeException("Eigen Decomposition (QR) can only be done on a square matrix. "
-				+ "Input matrix is rectangular (rows=" + in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
-		}
-
-		int m = in.rlen;
-		AggregateBinaryOperator op_mul_agg = InstructionUtils.getMatMultOperator(threads);
-
-		MatrixBlock Q_prod = new MatrixBlock(m, m, 0.0);
-		for(int i = 0; i < m; i++) {
-			Q_prod.set(i, i, 1.0);
-		}
-
-		for(int i = 0; i < num_iterations; i++) {
-			MatrixBlock[] QR = computeQR2(in, threads);
-			Q_prod = Q_prod.aggregateBinaryOperations(Q_prod, QR[0], op_mul_agg);
-			in = QR[1].aggregateBinaryOperations(QR[1], QR[0], op_mul_agg);
-		}
-
-		// Is converged if all values are below tol and the there only is values on the diagonal.
-
-		double[] check = in.getDenseBlockValues();
-		double[] eval = new double[m];
-		for(int i = 0; i < m; i++)
-			eval[i] = check[i*m+i];
-		
-		double[] evec = Q_prod.getDenseBlockValues();
-		return sortEVs(eval, evec);
-	}
-
-	/**
-	 * Function to compute the Householder transformation of a Matrix.
-	 *
-	 * @param in Input Matrix
-	 * @param threads number of threads
-	 * @return transformed matrix
-	 */
-	@SuppressWarnings("unused")
-	private static MatrixBlock computeHouseholder(MatrixBlock in, int threads) {
-		int m = in.rlen;
-
-		MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
-		A_n.copy(in);
-
-		for(int k = 0; k < m - 2; k++) {
-			MatrixBlock ajk = A_n.slice(0, m - 1, k, k);
-			for(int i = 0; i <= k; i++) {
-				ajk.set(i, 0, 0.0);
-			}
-			double alpha = Math.sqrt(ajk.sumSq());
-			double ak1k = A_n.getDouble(k + 1, k);
-			if(ak1k > 0.0)
-				alpha *= -1;
-			double r = Math.sqrt(0.5 * (alpha * alpha - ak1k * alpha));
-			MatrixBlock v = new MatrixBlock(m, 1, 0.0);
-			v.copy(ajk);
-			v.set(k + 1, 0, ak1k - alpha);
-			ScalarOperator op_div_scalar = new RightScalarOperator(Divide.getDivideFnObject(), 2 * r, threads);
-			v = v.scalarOperations(op_div_scalar, new MatrixBlock());
-
-			MatrixBlock P = new MatrixBlock(m, m, 0.0);
-			for(int i = 0; i < m; i++) {
-				P.set(i, i, 1.0);
-			}
-
-			ReorgOperator op_t = new ReorgOperator(SwapIndex.getSwapIndexFnObject(), threads);
-			AggregateBinaryOperator op_mul_agg = InstructionUtils.getMatMultOperator(threads);
-			BinaryOperator op_add = InstructionUtils.parseExtendedBinaryOperator("+");
-			BinaryOperator op_sub = InstructionUtils.parseExtendedBinaryOperator("-");
-
-			MatrixBlock v_t = v.reorgOperations(op_t, new MatrixBlock(), 0, 0, m);
-			v_t = v_t.binaryOperations(op_add, v_t);
-			MatrixBlock v_v_t_2 = A_n.aggregateBinaryOperations(v, v_t, op_mul_agg);
-			P = P.binaryOperations(op_sub, v_v_t_2);
-			A_n = A_n.aggregateBinaryOperations(P, A_n.aggregateBinaryOperations(A_n, P, op_mul_agg), op_mul_agg);
-		}
-		return A_n;
-	}
-
-	/**
-	 * Sort the eigen values (and vectors) in increasing order (to be compatible w/ LAPACK.DSYEVR())
-	 *
-	 * @param eValues  Eigenvalues
-	 * @param eVectors Eigenvectors
-	 * @return array of matrix blocks
-	 */
-	private static MatrixBlock[] sortEVs(double[] eValues, double[][] eVectors) {
-		int n = eValues.length;
-		for(int i = 0; i < n; i++) {
-			int k = i;
-			double p = eValues[i];
-			for(int j = i + 1; j < n; j++) {
-				if(eValues[j] < p) {
-					k = j;
-					p = eValues[j];
-				}
-			}
-			if(k != i) {
-				eValues[k] = eValues[i];
-				eValues[i] = p;
-				for(int j = 0; j < n; j++) {
-					p = eVectors[j][i];
-					eVectors[j][i] = eVectors[j][k];
-					eVectors[j][k] = p;
-				}
-			}
-		}
-
-		MatrixBlock eval = DataConverter.convertToMatrixBlock(eValues, true);
-		MatrixBlock evec = DataConverter.convertToMatrixBlock(eVectors);
-		return new MatrixBlock[] {eval, evec};
-	}
-
-	private static MatrixBlock[] sortEVs(double[] eValues, double[] eVectors) {
-		int n = eValues.length;
-		for(int i = 0; i < n; i++) {
-			int k = i;
-			double p = eValues[i];
-			for(int j = i + 1; j < n; j++) {
-				if(eValues[j] < p) {
-					k = j;
-					p = eValues[j];
-				}
-			}
-			if(k != i) {
-				eValues[k] = eValues[i];
-				eValues[i] = p;
-				for(int j = 0; j < n; j++) {
-					p = eVectors[j*n+i];
-					eVectors[j*n+i] = eVectors[j*n+k];
-					eVectors[j*n+k] = p;
-				}
-			}
-		}
-
-		MatrixBlock eval = DataConverter.convertToMatrixBlock(eValues, true);
-		MatrixBlock evec = new MatrixBlock(n, n, false);
-		evec.init(eVectors, n, n);
-		return new MatrixBlock[] {eval, evec};
-	}
-}
+ package org.apache.sysds.runtime.matrix.data;
+
+ import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.fft;
+ import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.fft_linearized;
+ import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.ifft;
+ import static org.apache.sysds.runtime.matrix.data.LibMatrixFourier.ifft_linearized;
+ 
+ import java.util.Arrays;
+ 
+ import org.apache.commons.logging.Log;
+ import org.apache.commons.logging.LogFactory;
+ import org.apache.commons.math3.exception.MaxCountExceededException;
+ import org.apache.commons.math3.exception.util.LocalizedFormats;
+ import org.apache.commons.math3.linear.Array2DRowRealMatrix;
+ import org.apache.commons.math3.linear.BlockRealMatrix;
+ import org.apache.commons.math3.linear.CholeskyDecomposition;
+ import org.apache.commons.math3.linear.DecompositionSolver;
+ import org.apache.commons.math3.linear.EigenDecomposition;
+ import org.apache.commons.math3.linear.LUDecomposition;
+ import org.apache.commons.math3.linear.QRDecomposition;
+ import org.apache.commons.math3.linear.RealMatrix;
+ import org.apache.commons.math3.linear.SingularValueDecomposition;
+ import org.apache.commons.math3.util.FastMath;
+ import org.apache.sysds.runtime.DMLRuntimeException;
+ import org.apache.sysds.runtime.data.DenseBlock;
+ import org.apache.sysds.runtime.data.DenseBlockFactory;
+ import org.apache.sysds.runtime.functionobjects.Builtin;
+ import org.apache.sysds.runtime.functionobjects.Divide;
+ import org.apache.sysds.runtime.functionobjects.MinusMultiply;
+ import org.apache.sysds.runtime.functionobjects.Multiply;
+ import org.apache.sysds.runtime.functionobjects.SwapIndex;
+ import org.apache.sysds.runtime.instructions.InstructionUtils;
+ import org.apache.sysds.runtime.matrix.operators.AggregateBinaryOperator;
+ import org.apache.sysds.runtime.matrix.operators.BinaryOperator;
+ import org.apache.sysds.runtime.matrix.operators.LeftScalarOperator;
+ import org.apache.sysds.runtime.matrix.operators.ReorgOperator;
+ import org.apache.sysds.runtime.matrix.operators.RightScalarOperator;
+ import org.apache.sysds.runtime.matrix.operators.ScalarOperator;
+ import org.apache.sysds.runtime.matrix.operators.TernaryOperator;
+ import org.apache.sysds.runtime.matrix.operators.UnaryOperator;
+ import org.apache.sysds.runtime.util.DataConverter;
+ import org.apache.commons.math3.util.Precision;
+ 
+ /**
+  * Library for matrix operations that need invocation of 
+  * Apache Commons Math library. 
+  * 
+  * This library currently supports following operations:
+  * matrix inverse, matrix decompositions (QR, LU, Eigen), solve 
+  */
+ public class LibCommonsMath 
+ {
+	 private static final Log LOG = LogFactory.getLog(LibCommonsMath.class.getName());
+	 private static final double RELATIVE_SYMMETRY_THRESHOLD = 1e-6;
+	 private static final double EIGEN_LAMBDA = 1e-8;
+	 private static final double EIGEN_EPS = Precision.EPSILON;
+	 private static final int EIGEN_MAX_ITER = 30;
+ 
+	 private LibCommonsMath() {
+		 //prevent instantiation via private constructor
+	 }
+	 
+	 public static boolean isSupportedUnaryOperation( String opcode ) {
+		 return ( opcode.equals("inverse") || opcode.equals("cholesky") );
+	 }
+	 
+	 public static boolean isSupportedMultiReturnOperation( String opcode ) {
+ 
+		 switch (opcode) {
+			 case "qr":
+			 case "lu":
+			 case "eigen":
+			 case "fft":
+			 case "ifft":
+			 case "fft_linearized":
+			 case "ifft_linearized":
+			 case "stft":
+			 case "svd": return true;
+			 default: return false;
+		 }
+ 
+	 }
+	 
+	 public static boolean isSupportedMatrixMatrixOperation( String opcode ) {
+		 return ( opcode.equals("solve") );
+	 }
+		 
+	 public static MatrixBlock unaryOperations(MatrixBlock inj, String opcode) {
+		 Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(inj);
+		 if(opcode.equals("inverse"))
+			 return computeMatrixInverse(matrixInput);
+		 else if (opcode.equals("cholesky"))
+			 return computeCholesky(matrixInput);
+		 return null;
+	 }
+ 
+	 // public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode) {
+	 // 	return multiReturnOperations(in, opcode, 1, 1);
+	 // }
+ 
+	 public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode, int threads) {
+		 return multiReturnOperations(in, opcode, threads, 1);
+	 }
+ 
+	 public static MatrixBlock[] multiReturnOperations(MatrixBlock in1, MatrixBlock in2, String opcode) {
+		 return multiReturnOperations(in1, in2, opcode, 1, 1);
+	 }
+ 
+	 public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode, int threads, int num_iterations, double tol) {
+		 if(opcode.equals("eigen_qr"))
+			 return computeEigenQR(in, num_iterations, tol, threads);
+		 else
+			 return multiReturnOperations(in, opcode, threads, 1);
+	 }
+ 
+	 public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode, int threads, long seed) {
+		 if(opcode.equals("qr"))
+			 return computeQR(in);
+		 else if (opcode.equals("qr2"))
+			 return computeQR2(in, threads);
+		 else if (opcode.equals("lu"))
+			 return computeLU(in);
+		 else if (opcode.equals("eigen"))
+			 return computeEigen(in);
+		 else if (opcode.equals("eigen_lanczos"))
+			 return computeEigenLanczos(in, threads, seed);
+		 else if (opcode.equals("eigen_qr"))
+			 return computeEigenQR(in, threads);
+		 else if (opcode.equals("eigen_symm"))
+			 return computeEigenDecompositionSymm(in, threads);
+		 else if (opcode.equals("svd"))
+			 return computeSvd(in);
+		 else if (opcode.equals("fft"))
+			 return computeFFT(in, threads);
+		 else if (opcode.equals("ifft"))
+			 return computeIFFT(in, threads);
+		 else if (opcode.equals("fft_linearized"))
+			 return computeFFT_LINEARIZED(in, threads);
+		 else if (opcode.equals("ifft_linearized"))
+			 return computeIFFT_LINEARIZED(in, threads);
+		 return null;
+	 }
+ 
+	 public static MatrixBlock[] multiReturnOperations(MatrixBlock in1, MatrixBlock in2, String opcode, int threads,
+			 long seed) {
+ 
+		 switch (opcode) {
+			 case "ifft":
+				 return computeIFFT(in1, in2, threads);
+			 case "ifft_linearized":
+				 return computeIFFT_LINEARIZED(in1, in2, threads);
+			 default:
+				 return null;
+		 }
+ 
+	 }
+ 
+	 public static MatrixBlock matrixMatrixOperations(MatrixBlock in1, MatrixBlock in2, String opcode) {
+		 if(opcode.equals("solve")) {
+			 if (in1.getNumRows() != in1.getNumColumns())
+				 throw new DMLRuntimeException("The A matrix, in solve(A,b) should have squared dimensions.");
+			 return computeSolve(in1, in2);
+		 }
+		 return null;
+	 }
+	 
+	 /**
+	  * Function to solve a given system of equations.
+	  * 
+	  * @param in1 matrix object 1
+	  * @param in2 matrix object 2
+	  * @return matrix block
+	  */
+	 private static MatrixBlock computeSolve(MatrixBlock in1, MatrixBlock in2) {
+		 //convert to commons math BlockRealMatrix instead of Array2DRowRealMatrix
+		 //to avoid unnecessary conversion as QR internally creates a BlockRealMatrix
+		 BlockRealMatrix matrixInput = DataConverter.convertToBlockRealMatrix(in1);
+		 BlockRealMatrix vectorInput = DataConverter.convertToBlockRealMatrix(in2);
+		 
+		 /*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput);
+		 DecompositionSolver lusolver = ludecompose.getSolver();
+		 RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/
+		 
+		 // Setup a solver based on QR Decomposition
+		 QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
+		 DecompositionSolver solver = qrdecompose.getSolver();
+		 // Invoke solve
+		 RealMatrix solutionMatrix = solver.solve(vectorInput);
+		 
+		 return DataConverter.convertToMatrixBlock(solutionMatrix);
+	 }
+ 
+	 /**
+	  * Computes the eigen decomposition of a symmetric matrix.
+	  *
+	  * @param in The input matrix to compute the eigen decomposition on.
+	  * @param threads The number of threads to use for computation.
+	  * @return An array of MatrixBlock objects containing the real eigen values and eigen vectors.
+	  */
+	 public static MatrixBlock[] computeEigenDecompositionSymm(MatrixBlock in, int threads) {
+ 
+		 // TODO: Verify matrix is symmetric
+		 final double[] mainDiag = new double[in.rlen];
+		 final double[] secDiag = new double[in.rlen - 1];
+ 
+		 final MatrixBlock householderVectors = transformToTridiagonal(in, mainDiag, secDiag, threads);
+ 
+		 // TODO: Consider using sparse blocks
+		 final double[] hv = householderVectors.getDenseBlockValues();
+		 MatrixBlock houseHolderMatrix = getQ(hv, mainDiag, secDiag);
+ 
+		 MatrixBlock[] evResult = findEigenVectors(mainDiag, secDiag, houseHolderMatrix, EIGEN_MAX_ITER, threads);
+ 
+		 MatrixBlock realEigenValues = evResult[0];
+		 MatrixBlock eigenVectors = evResult[1];
+ 
+		 realEigenValues.setNonZeros(realEigenValues.rlen * realEigenValues.rlen);
+		 eigenVectors.setNonZeros(eigenVectors.rlen * eigenVectors.clen);
+ 
+		 return new MatrixBlock[] { realEigenValues, eigenVectors };
+	 }
+ 
+	 private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] main, double[] secondary, int threads) {
+		 final int m = matrix.rlen;
+ 
+		 // MatrixBlock householderVectors = ;
+		 MatrixBlock householderVectors = matrix.extractTriangular(new MatrixBlock(m, m, false), false, true, true);
+		 if (householderVectors.isInSparseFormat()) {
+			 householderVectors.sparseToDense(threads);      
+		 }
+ 
+		 final double[] z = new double[m];
+ 
+		 // TODO: Consider sparse block case
+		 final double[] hv = householderVectors.getDenseBlockValues();
+ 
+		 for (int k = 0; k < m - 1; k++) {
+			 final int k_kp1 = k * m + k + 1;
+			 final int km = k * m;
+ 
+			 // zero-out a row and a column simultaneously
+			 main[k] = hv[km + k];
+ 
+			 // TODO: may or may not improve, seems it does not
+			 // double xNormSqr = householderVectors.slice(k, k, k + 1, m - 1).sumSq();
+			 // double xNormSqr = LibMatrixMult.dotProduct(hv, hv, k_kp1, k_kp1, m - (k + 1));
+			 double xNormSqr = 0;
+			 for (int j = k + 1; j < m; ++j) {
+			 final double c = hv[k * m + j];
+			 xNormSqr += c * c;
+			 }
+ 
+			 final double a = (hv[k_kp1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
+ 
+			 secondary[k] = a;
+ 
+			 if (a != 0.0) {
+			 // apply Householder transform from left and right simultaneously
+			 hv[k_kp1] -= a;
+ 
+			 final double beta = -1 / (a * hv[k_kp1]);
+			 double gamma = 0;
+ 
+			 // compute a = beta A v, where v is the Householder vector
+			 // this loop is written in such a way
+			 // 1) only the upper triangular part of the matrix is accessed
+			 // 2) access is cache-friendly for a matrix stored in rows
+			 Arrays.fill(z, k + 1, m, 0);
+			 for (int i = k + 1; i < m; ++i) {
+				 final double hKI = hv[km + i];
+				 double zI = hv[i * m + i] * hKI;
+				 for (int j = i + 1; j < m; ++j) {
+				 final double hIJ = hv[i * m + j];
+				 zI += hIJ * hv[k * m + j];
+				 z[j] += hIJ * hKI;
+				 }
+				 z[i] = beta * (z[i] + zI);
+ 
+				 gamma += z[i] * hv[km + i];
+			 }
+ 
+			 gamma *= beta / 2;
+ 
+			 // compute z = z - gamma v
+			 // LibMatrixMult.vectMultiplyAdd(-gamma, hv, z, k_kp1, k + 1, m - (k + 1));
+			 for (int i = k + 1; i < m; ++i) {
+				 z[i] -= gamma * hv[km + i];
+			 }
+ 
+			 // update matrix: A = A - v zT - z vT
+			 // only the upper triangular part of the matrix is updated
+			 for (int i = k + 1; i < m; ++i) {
+				 final double hki = hv[km + i];
+				 for (int j = i; j < m; ++j) {
+				 final double hkj = hv[km + j];
+ 
+				 hv[i * m + j] -= (hki * z[j] + z[i] * hkj);
+				 }
+			 }
+			 }
+		 }
+ 
+		 main[m - 1] = hv[(m - 1) * m + m - 1];
+ 
+		 return householderVectors;
+	 }
+ 
+ 
+	 /**
+	  * Computes the orthogonal matrix Q using Householder transforms.
+	  * The matrix Q is built by applying Householder transforms to the input vectors.
+	  *
+	  * @param hv        The input vector containing the Householder vectors.
+	  * @param main      The main diagonal of the matrix.
+	  * @param secondary The secondary diagonal of the matrix.
+	  * @return The orthogonal matrix Q.
+	  */
+	 public static MatrixBlock getQ(final double[] hv, double[] main, double[] secondary) {
+		 final int m = main.length;
+ 
+		 // orthogonal matrix, most operations are column major (TODO)
+		 DenseBlock qaB = DenseBlockFactory.createDenseBlock(m, m);
+		 double[] qaV = qaB.valuesAt(0);
+ 
+		 // build up first part of the matrix by applying Householder transforms
+		 for (int k = m - 1; k >= 1; --k) {
+			 final int km = k * m;
+			 final int km1m = (k - 1) * m;
+ 
+			 qaV[km + k] = 1.0;
+			 if (hv[km1m + k] != 0.0) {
+			 final double inv = 1.0 / (secondary[k - 1] * hv[km1m + k]);
+ 
+			 double beta = 1.0 / secondary[k - 1];
+ 
+			 qaV[km + k] = 1 + beta * hv[km1m + k];
+ 
+			 // TODO: may speedup vector operations
+			 for (int i = k + 1; i < m; ++i) {
+				 qaV[i * m + k] = beta * hv[km1m + i];
+			 }
+ 
+			 for (int j = k + 1; j < m; ++j) {
+				 beta = 0;
+				 for (int i = k + 1; i < m; ++i) {
+				 beta += qaV[m * i + j] * hv[km1m + i];
+				 }
+				 beta *= inv;
+				 qaV[m * k + j] = hv[km1m + k] * beta;
+ 
+				 for (int i = k + 1; i < m; ++i) {
+				 qaV[m * i + j] += beta * hv[km1m + i];
+				 }
+			 }
+			 }
+		 }
+ 
+		 qaV[0] = 1.0;
+		 MatrixBlock res = new MatrixBlock(m, m, qaB);
+ 
+		 return res;
+	 }
+ 
+ 
+	 /**
+	  * Finds the eigen vectors corresponding to the given eigen values using the Householder transformation. 
+	  * (Dubrulle et al., 1971).
+	  *
+	  * @param main      The main diagonal of the tridiagonal matrix.
+	  * @param secondary The secondary diagonal of the tridiagonal matrix.
+	  * @param hhMatrix  The Householder matrix.
+	  * @param maxIter   The maximum number of iterations for convergence.
+	  * @param threads   The number of threads to use for computation.
+	  * @return  An array of two MatrixBlock objects: eigen values and eigen vectors.
+	  * @throws MaxCountExceededException If the maximum number of iterations is exceeded and convergence fails.
+	  */
+	 private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, final MatrixBlock hhMatrix,
+			 final int maxIter, int threads) {
+ 
+		 DenseBlock hhDense = hhMatrix.getDenseBlock();
+ 
+		 final int n = hhMatrix.rlen;
+		 MatrixBlock eigenValues = new MatrixBlock(n, 1, main);
+ 
+		 double[] ev = eigenValues.denseBlock.valuesAt(0);
+		 final double[] e = new double[n];
+ 
+		 System.arraycopy(secondary, 0, e, 0, n - 1);
+		 e[n - 1] = 0;
+ 
+		 // Determine the largest main and secondary value in absolute term.
+		 double maxAbsoluteValue = 0;
+		 for (int i = 0; i < n; i++) {
+			 maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(ev[i]));
+			 maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(e[i]));
+		 }
+		 // Make null any main and secondary value too small to be significant
+		 if (maxAbsoluteValue != 0) {
+			 for (int i = 0; i < n; i++) {
+			 if (FastMath.abs(ev[i]) <= EIGEN_EPS * maxAbsoluteValue && ev[i] != 0.0) {
+				 ev[i] = 0;
+			 }
+			 if (FastMath.abs(e[i]) <= EIGEN_EPS * maxAbsoluteValue && e[i] != 0.0) {
+				 e[i] = 0;
+			 }
+			 }
+		 }
+ 
+		 for (int j = 0; j < n; j++) {
+			 int its = 0;
+			 int m;
+			 do {
+			 for (m = j; m < n - 1; m++) {
+				 final double delta = FastMath.abs(ev[m]) +
+					 FastMath.abs(ev[m + 1]);
+				 if (FastMath.abs(e[m]) + delta == delta) {
+				 break;
+				 }
+			 }
+			 if (m != j) {
+				 if (its == maxIter) {
+				 throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, maxIter);
+				 }
+ 
+				 its++;
+				 double q = (ev[j + 1] - ev[j]) / (2 * e[j]);
+				 double t = FastMath.sqrt(1 + q * q);
+				 if (q < 0.0) {
+				 q = ev[m] - ev[j] + e[j] / (q - t);
+				 } else {
+				 q = ev[m] - ev[j] + e[j] / (q + t);
+				 }
+				 double u = 0.0;
+				 double s = 1.0;
+				 double c = 1.0;
+				 int i;
+				 for (i = m - 1; i >= j; i--) {
+				 double p = s * e[i];
+				 double h = c * e[i];
+				 if (FastMath.abs(p) >= FastMath.abs(q)) {
+					 c = q / p;
+					 t = FastMath.sqrt(c * c + 1.0);
+					 e[i + 1] = p * t;
+					 s = 1.0 / t;
+					 c *= s;
+				 } else {
+					 s = p / q;
+					 t = FastMath.sqrt(s * s + 1.0);
+					 e[i + 1] = q * t;
+					 c = 1.0 / t;
+					 s *= c;
+				 }
+				 if (e[i + 1] == 0.0) {
+					 ev[i + 1] -= u;
+					 e[m] = 0.0;
+					 break;
+				 }
+				 q = ev[i + 1] - u;
+				 t = (ev[i] - q) * s + 2.0 * c * h;
+				 u = s * t;
+				 ev[i + 1] = q + u;
+				 q = c * t - h;
+ 
+				 for (int ia = 0; ia < n; ++ia) {
+					 p = hhDense.get(ia, i + 1);
+					 hhDense.set(ia, i + 1, s * hhDense.get(ia, i) + c * p);
+					 hhDense.set(ia, i, c * hhDense.get(ia, i) - s * p);
+				 }
+				 }
+ 
+				 if (t == 0.0 && i >= j) {
+				 continue;
+				 }
+				 ev[j] -= u;
+				 e[j] = q;
+				 e[m] = 0.0;
+ 
+			 }
+			 } while (m != j);
+		 }
+ 
+		 // Sort the eigen values (and vectors) in increase order
+		 for (int i = 0; i < n; i++) {
+			 int k = i;
+			 double p = ev[i];
+			 for (int j = i + 1; j < n; j++) {
+			 // reversed order from original implementation
+			 if (ev[j] < p) {
+				 k = j;
+				 p = ev[j];
+			 }
+			 }
+			 if (k != i) {
+			 ev[k] = ev[i];
+			 ev[i] = p;
+			 for (int j = 0; j < n; j++) {
+				 // TODO: an operation like SwapIndex be faster?
+				 p = hhDense.get(j, i);
+				 hhDense.set(j, i, hhDense.get(j, k));
+				 hhDense.set(j, k, p);
+ 
+			 }
+			 }
+		 }
+ 
+		 // Determine the largest eigen value in absolute term.
+		 maxAbsoluteValue = 0;
+		 for (int i = 0; i < n; i++) {
+			 maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(ev[i]));
+		 }
+		 // Make null any eigen value too small to be significant
+		 if (maxAbsoluteValue != 0.0) {
+			 for (int i = 0; i < n; i++) {
+			 if (FastMath.abs(ev[i]) < EIGEN_EPS * maxAbsoluteValue) {
+				 ev[i] = 0;
+			 }
+			 }
+		 }
+ 
+		 // MatrixBlock realEigenValues = new MatrixBlock(z.rlen, 1, ev);
+ 
+		 return new MatrixBlock[] { eigenValues, hhMatrix };
+	 }
+ 
+	 /**
+	  * Function to perform QR decomposition on a given matrix.
+	  * 
+	  * @param in matrix object
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeQR(MatrixBlock in) {
+		 Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
+		 
+		 // Perform QR decomposition
+		 QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
+		 RealMatrix H = qrdecompose.getH();
+		 RealMatrix R = qrdecompose.getR();
+		 
+		 // Read the results into native format
+		 MatrixBlock mbH = DataConverter.convertToMatrixBlock(H.getData());
+		 MatrixBlock mbR = DataConverter.convertToMatrixBlock(R.getData());
+ 
+		 return new MatrixBlock[] { mbH, mbR };
+	 }
+	 
+	 /**
+	  * Function to perform LU decomposition on a given matrix.
+	  * 
+	  * @param in matrix object
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeLU(MatrixBlock in) {
+		 if(in.getNumRows() != in.getNumColumns()) {
+			 throw new DMLRuntimeException(
+				 "LU Decomposition can only be done on a square matrix. Input matrix is rectangular (rows="
+					 + in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
+		 }
+		 
+		 Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
+		 
+		 // Perform LUP decomposition
+		 LUDecomposition ludecompose = new LUDecomposition(matrixInput);
+		 RealMatrix P = ludecompose.getP();
+		 RealMatrix L = ludecompose.getL();
+		 RealMatrix U = ludecompose.getU();
+		 
+		 // Read the results into native format
+		 MatrixBlock mbP = DataConverter.convertToMatrixBlock(P.getData());
+		 MatrixBlock mbL = DataConverter.convertToMatrixBlock(L.getData());
+		 MatrixBlock mbU = DataConverter.convertToMatrixBlock(U.getData());
+ 
+		 return new MatrixBlock[] { mbP, mbL, mbU };
+	 }
+	 
+	 /**
+	  * Function to perform Eigen decomposition on a given matrix.
+	  * Input must be a symmetric matrix.
+	  * 
+	  * @param in matrix object
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeEigen(MatrixBlock in) {
+		 if ( in.getNumRows() != in.getNumColumns() ) {
+			 throw new DMLRuntimeException("Eigen Decomposition can only be done on a square matrix. "
+				 + "Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
+		 }
+		 
+		 EigenDecomposition eigendecompose = null;
+		 try {
+			 Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
+			 eigendecompose = new EigenDecomposition(matrixInput);
+		 }
+		 catch(MaxCountExceededException ex) {
+			 LOG.warn("Eigen: "+ ex.getMessage()+". Falling back to regularized eigen factorization.");
+			 eigendecompose = computeEigenRegularized(in);
+		 }
+		 
+		 RealMatrix eVectorsMatrix = eigendecompose.getV();
+		 double[][] eVectors = eVectorsMatrix.getData();
+		 double[] eValues = eigendecompose.getRealEigenvalues();
+ 
+		 return sortEVs(eValues, eVectors);
+	 }
+ 
+	 private static EigenDecomposition computeEigenRegularized(MatrixBlock in) {
+		 if( in == null || in.isEmptyBlock(false) )
+			 throw new DMLRuntimeException("Invalid empty block");
+		 
+		 //slightly modify input for regularization (pos/neg)
+		 MatrixBlock in2 = new MatrixBlock(in, false);
+		 DenseBlock a = in2.getDenseBlock();
+		 for( int i=0; i<in2.rlen; i++ ) {
+			 double[] avals = a.values(i);
+			 int apos = a.pos(i);
+			 for( int j=0; j<in2.clen; j++ ) {
+				 double v = avals[apos+j];
+				 avals[apos+j] += Math.signum(v) * EIGEN_LAMBDA;
+			 }
+		 }
+		 
+		 //run eigen decomposition
+		 return new EigenDecomposition(
+			 DataConverter.convertToArray2DRowRealMatrix(in2));
+	 }
+ 
+	 /**
+	  * Function to perform FFT on a given matrix.
+	  *
+	  * @param re      matrix object
+	  * @param threads number of threads
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeFFT(MatrixBlock re, int threads) {
+		 if(re == null)
+			 throw new DMLRuntimeException("Invalid empty block");
+		 if(re.isEmptyBlock(false)) {
+			 // Return the original matrix as the result
+			 // Assuming you need to return two matrices: the real part and an imaginary part initialized to 0.
+			 return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)};
+		 }
+		 // run fft
+		 re.sparseToDense();
+		 return fft(re, threads);
+	 }
+ 
+	 /**
+	  * Function to perform IFFT on a given matrix.
+	  *
+	  * @param re      matrix object
+	  * @param im      matrix object
+	  * @param threads number of threads
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeIFFT(MatrixBlock re, MatrixBlock im, int threads) {
+		 if(re == null)
+			 throw new DMLRuntimeException("Invalid empty block");
+ 
+		 // run ifft
+		 if(im != null && !im.isEmptyBlock(false)) {
+			 re.sparseToDense();
+			 im.sparseToDense();
+			 return ifft(re, im, threads);
+		 }
+		 else {
+			 if(re.isEmptyBlock(false)) {
+				 // Return the original matrix as the result
+				 // Assuming you need to return two matrices: the real part and an imaginary part initialized to 0.
+				 return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)};
+			 }
+			 re.sparseToDense();
+			 return ifft(re, threads);
+		 }
+	 }
+ 
+	 /**
+	  * Function to perform IFFT on a given matrix.
+	  *
+	  * @param re      matrix object
+	  * @param threads number of threads
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeIFFT(MatrixBlock re, int threads) {
+		 return computeIFFT(re, null, threads);
+	 }
+ 
+	 /**
+	  * Function to perform FFT_LINEARIZED on a given matrix.
+	  *
+	  * @param re      matrix object
+	  * @param threads number of threads
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeFFT_LINEARIZED(MatrixBlock re, int threads) {
+		 if(re == null)
+			 throw new DMLRuntimeException("Invalid empty block");
+		 if(re.isEmptyBlock(false)) {
+			 // Return the original matrix as the result
+			 // Assuming you need to return two matrices: the real part and an imaginary part initialized to 0.
+			 return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)};
+		 }
+		 // run fft
+		 re.sparseToDense();
+		 return fft_linearized(re, threads);
+	 }
+ 
+	 /**
+	  * Function to perform IFFT_LINEARIZED on a given matrix.
+	  *
+	  * @param re      matrix object
+	  * @param im      matrix object
+	  * @param threads number of threads
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeIFFT_LINEARIZED(MatrixBlock re, MatrixBlock im, int threads) {
+		 if(re == null)
+			 throw new DMLRuntimeException("Invalid empty block");
+ 
+		 // run ifft
+		 if(im != null && !im.isEmptyBlock(false)) {
+			 re.sparseToDense();
+			 im.sparseToDense();
+			 return ifft_linearized(re, im, threads);
+		 }
+		 else {
+			 if(re.isEmptyBlock(false)) {
+				 // Return the original matrix as the result
+				 // Assuming you need to return two matrices: the real part and an imaginary part initialized to 0.
+				 return new MatrixBlock[] {re, new MatrixBlock(re.getNumRows(), re.getNumColumns(), true)};
+			 }
+			 re.sparseToDense();
+			 return ifft_linearized(re, threads);
+		 }
+	 }
+ 
+	 /**
+	  * Function to perform IFFT_LINEARIZED on a given matrix
+	  * 
+	  * @param re      matrix object
+	  * @param threads number of threads
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeIFFT_LINEARIZED(MatrixBlock re, int threads) {
+		 return computeIFFT_LINEARIZED(re, null, threads);
+	 }
+ 
+ 
+	 // /**
+	 //  * Function to perform STFT on a given matrix.
+	 //  *
+	 //  * @param re matrix object
+	 //  * @param im matrix object
+	 //  * @param windowSize of stft
+	 //  * @param overlap of stft
+	 //  * @return array of matrix blocks
+	 //  */
+	 // private static MatrixBlock[] computeSTFT(MatrixBlock re, MatrixBlock im, int windowSize, int overlap, int threads) {
+	 // 	if (re == null) {
+	 // 		throw new DMLRuntimeException("Invalid empty block");
+	 // 	} else if (im != null && !im.isEmptyBlock(false)) {
+	 // 		re.sparseToDense();
+	 // 		im.sparseToDense();
+	 // 		return stft(re, im, windowSize, overlap, threads);
+	 // 	} else {
+	 // 		if (re.isEmptyBlock(false)) {
+	 // 			// Return the original matrix as the result
+	 // 			int rows = re.getNumRows();
+	 // 			int cols = re.getNumColumns();
+ 
+	 // 			int stepSize = windowSize - overlap;
+	 // 			if (stepSize == 0) {
+	 // 				throw new IllegalArgumentException("windowSize - overlap is zero");
+	 // 			}
+ 
+	 // 			int numberOfFramesPerRow = (cols - overlap + stepSize - 1) / stepSize;
+	 // 			int rowLength= numberOfFramesPerRow * windowSize;
+	 // 			int out_len = rowLength * rows;
+ 
+	 // 			double[] out_zero = new double[out_len];
+ 
+	 // 			return new MatrixBlock[]{new MatrixBlock(rows, rowLength, out_zero), new MatrixBlock(rows, rowLength, out_zero)};
+	 // 			}
+	 // 		re.sparseToDense();
+	 // 		return stft(re, windowSize, overlap, threads);
+	 // 	}
+	 // }
+ 
+	 // /**
+	 //  * Function to perform STFT on a given matrix.
+	 //  *
+	 //  * @param re matrix object
+	 //  * @return array of matrix blocks
+	 //  */
+	 // private static MatrixBlock[] computeSTFT(MatrixBlock re, int windowSize, int overlap, int threads) {
+	 // 	return computeSTFT(re, null, windowSize, overlap, threads);
+	 // }
+ 
+	 /**
+	  * Performs Singular Value Decomposition. Calls Apache Commons Math SVD.
+	  * X = U * Sigma * Vt, where X is the input matrix,
+	  * U is the left singular matrix, Sigma is the singular values matrix returned as a
+	  * column matrix and Vt is the transpose of the right singular matrix V.
+	  * However, the returned array has  { U, Sigma, V}
+	  * 
+	  * @param in Input matrix
+	  * @return An array containing U, Sigma & V
+	  */
+	 private static MatrixBlock[] computeSvd(MatrixBlock in) {
+		 Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
+ 
+		 SingularValueDecomposition svd = new SingularValueDecomposition(matrixInput);
+		 double[] sigma = svd.getSingularValues();
+		 RealMatrix u = svd.getU();
+		 RealMatrix v = svd.getV();
+		 MatrixBlock U = DataConverter.convertToMatrixBlock(u.getData());
+		 MatrixBlock Sigma = DataConverter.convertToMatrixBlock(sigma, true);
+		 Sigma = LibMatrixReorg.diag(Sigma, new MatrixBlock(Sigma.rlen, Sigma.rlen, true));
+		 MatrixBlock V = DataConverter.convertToMatrixBlock(v.getData());
+ 
+		 return new MatrixBlock[] { U, Sigma, V };
+	 }
+	 
+	 /**
+	  * Function to compute matrix inverse via matrix decomposition.
+	  * 
+	  * @param in commons-math3 Array2DRowRealMatrix
+	  * @return matrix block
+	  */
+	 private static MatrixBlock computeMatrixInverse(Array2DRowRealMatrix in) {
+		 if(!in.isSquare())
+			 throw new DMLRuntimeException("Input to inv() must be square matrix -- given: a " + in.getRowDimension()
+				 + "x" + in.getColumnDimension() + " matrix.");
+ 
+		 QRDecomposition qrdecompose = new QRDecomposition(in);
+		 DecompositionSolver solver = qrdecompose.getSolver();
+		 RealMatrix inverseMatrix = solver.getInverse();
+ 
+		 return DataConverter.convertToMatrixBlock(inverseMatrix.getData());
+	 }
+ 
+	 /**
+	  * Function to compute Cholesky decomposition of the given input matrix. 
+	  * The input must be a real symmetric positive-definite matrix.
+	  * 
+	  * @param in commons-math3 Array2DRowRealMatrix
+	  * @return matrix block
+	  */
+	 private static MatrixBlock computeCholesky(Array2DRowRealMatrix in) {
+		 if(!in.isSquare())
+			 throw new DMLRuntimeException("Input to cholesky() must be square matrix -- given: a "
+				 + in.getRowDimension() + "x" + in.getColumnDimension() + " matrix.");
+		 CholeskyDecomposition cholesky = new CholeskyDecomposition(in, RELATIVE_SYMMETRY_THRESHOLD,
+			 CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
+		 RealMatrix rmL = cholesky.getL();
+		 return DataConverter.convertToMatrixBlock(rmL.getData());
+	 }
+ 
+	 /**
+	  * Creates a random normalized vector with dim elements.
+	  *
+	  * @param dim dimension of the created vector
+	  * @param threads number of threads
+	  * @param seed seed for the random MatrixBlock generation
+	  * @return random normalized vector
+	  */
+	 private static MatrixBlock randNormalizedVect(int dim, int threads, long seed) {
+		 MatrixBlock v1 = MatrixBlock.randOperations(dim, 1, 1.0, 0.0, 1.0, "UNIFORM", seed);
+ 
+		 double v1_sum = v1.sum();
+		 ScalarOperator op_div_scalar = new RightScalarOperator(Divide.getDivideFnObject(), v1_sum, threads);
+		 v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+		 UnaryOperator op_sqrt = new UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT), threads, true);
+		 v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ 
+		 if(Math.abs(v1.sumSq() - 1.0) >= 1e-7)
+			 throw new DMLRuntimeException("v1 not correctly normalized (maybe try changing the seed)");
+ 
+		 return v1;
+	 }
+ 
+	 /**
+	  * Function to perform the Lanczos algorithm and then computes the Eigendecomposition.
+	  * Caution: Lanczos is not numerically stable (see https://en.wikipedia.org/wiki/Lanczos_algorithm)
+	  * Input must be a symmetric (and square) matrix.
+	  *
+	  * @param in matrix object
+	  * @param threads number of threads
+	  * @param seed seed for the random MatrixBlock generation
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeEigenLanczos(MatrixBlock in, int threads, long seed) {
+		 if(in.getNumRows() != in.getNumColumns()) {
+			 throw new DMLRuntimeException(
+				 "Lanczos algorithm and Eigen Decomposition can only be done on a square matrix. "
+					 + "Input matrix is rectangular (rows=" + in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
+		 }
+ 
+		 int m = in.getNumRows();
+		 MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+		 MatrixBlock v1 = randNormalizedVect(m, threads, seed);
+ 
+		 MatrixBlock T = new MatrixBlock(m, m, 0.0);
+		 MatrixBlock TV = new MatrixBlock(m, m, 0.0);
+		 MatrixBlock w1;
+ 
+		 ReorgOperator op_t = new ReorgOperator(SwapIndex.getSwapIndexFnObject(), threads);
+		 TernaryOperator op_minus_mul = new TernaryOperator(MinusMultiply.getFnObject(), threads);
+		 AggregateBinaryOperator op_mul_agg = InstructionUtils.getMatMultOperator(threads);
+		 ScalarOperator op_div_scalar = new RightScalarOperator(Divide.getDivideFnObject(), 1, threads);
+ 
+		 MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+		 for(int i = 0; i < m; i++) {
+			 v1.putInto(TV, 0, i, false);
+			 w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+			 MatrixBlock alpha = w1.aggregateBinaryOperations(v1.reorgOperations(op_t, new MatrixBlock(), 0, 0, m), w1, op_mul_agg);
+			 if(i < m - 1) {
+				 w1 = w1.ternaryOperations(op_minus_mul, v1, alpha, new MatrixBlock());
+				 w1 = w1.ternaryOperations(op_minus_mul, v0, beta, new MatrixBlock());
+				 beta.set(0, 0, Math.sqrt(w1.sumSq()));
+				 v0.copy(v1);
+				 op_div_scalar = op_div_scalar.setConstant(beta.getDouble(0, 0));
+				 w1.scalarOperations(op_div_scalar, v1);
+ 
+				 T.set(i + 1, i, beta.get(0, 0));
+				 T.set(i, i + 1, beta.get(0, 0));
+			 }
+			 T.set(i, i, alpha.get(0, 0));
+		 }
+ 
+		 MatrixBlock[] e = computeEigen(T);
+		 TV.setNonZeros((long) m*m);
+		 e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+		 return e;
+	 }
+ 
+	 /**
+	  * Function to perform the QR decomposition.
+	  * Input must be a square matrix.
+	  * TODO: use Householder transformation and implicit shifts to further speed up QR decompositions
+	  *
+	  * @param in matrix object
+	  * @param threads number of threads
+	  * @return array of matrix blocks [Q, R]
+	  */
+	 private static MatrixBlock[] computeQR2(MatrixBlock in, int threads) {
+		 if(in.getNumRows() != in.getNumColumns()) {
+			 throw new DMLRuntimeException("QR2 Decomposition can only be done on a square matrix. "
+				 + "Input matrix is rectangular (rows=" + in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
+		 }
+ 
+		 int m = in.rlen;
+ 
+		 MatrixBlock A_n = new MatrixBlock();
+		 A_n.copy(in);
+ 
+		 MatrixBlock Q_n = new MatrixBlock(m, m, true);
+		 for(int i = 0; i < m; i++) {
+			 Q_n.set(i, i, 1.0);
+		 }
+ 
+		 ReorgOperator op_t = new ReorgOperator(SwapIndex.getSwapIndexFnObject(), threads);
+		 AggregateBinaryOperator op_mul_agg = InstructionUtils.getMatMultOperator(threads);
+		 BinaryOperator op_sub = InstructionUtils.parseExtendedBinaryOperator("-");
+		 ScalarOperator op_div_scalar = new RightScalarOperator(Divide.getDivideFnObject(), 1, threads);
+		 ScalarOperator op_mult_2 = new LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, threads);
+ 
+		 for(int k = 0; k < m; k++) {
+			 MatrixBlock z = A_n.slice(k, m - 1, k, k);
+			 MatrixBlock uk = new MatrixBlock(m - k, 1, 0.0);
+			 uk.copy(z);
+			 uk.set(0, 0, uk.get(0, 0) + Math.signum(z.get(0, 0)) * Math.sqrt(z.sumSq()));
+			 op_div_scalar = op_div_scalar.setConstant(Math.sqrt(uk.sumSq()));
+			 uk = uk.scalarOperations(op_div_scalar, new MatrixBlock());
+ 
+			 MatrixBlock vk = new MatrixBlock(m, 1, 0.0);
+			 vk.copy(k, m - 1, 0, 0, uk, true);
+ 
+			 MatrixBlock vkt = vk.reorgOperations(op_t, new MatrixBlock(), 0, 0, m);
+			 MatrixBlock vkt2 = vkt.scalarOperations(op_mult_2, new MatrixBlock());
+			 MatrixBlock vkvkt2 = vk.aggregateBinaryOperations(vk, vkt2, op_mul_agg);
+ 
+			 A_n = A_n.binaryOperations(op_sub, A_n.aggregateBinaryOperations(vkvkt2, A_n, op_mul_agg));
+			 Q_n = Q_n.binaryOperations(op_sub, Q_n.aggregateBinaryOperations(Q_n, vkvkt2, op_mul_agg));
+		 }
+		 // QR decomp: Q: Q_n; R: A_n
+		 return new MatrixBlock[] {Q_n, A_n};
+	 }
+ 
+	 /**
+	  * Function that computes the Eigen Decomposition using the QR algorithm.
+	  * Caution: check if the QR algorithm is converged, if not increase iterations
+	  * Caution: if the input matrix has complex eigenvalues results will be incorrect
+	  *
+	  * @param in Input matrix
+	  * @param threads number of threads
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] computeEigenQR(MatrixBlock in, int threads) {
+		 return computeEigenQR(in, 100, 1e-10, threads);
+	 }
+ 
+	 private static MatrixBlock[] computeEigenQR(MatrixBlock in, int num_iterations, double tol, int threads) {
+		 if(in.getNumRows() != in.getNumColumns()) {
+			 throw new DMLRuntimeException("Eigen Decomposition (QR) can only be done on a square matrix. "
+				 + "Input matrix is rectangular (rows=" + in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
+		 }
+ 
+		 int m = in.rlen;
+		 AggregateBinaryOperator op_mul_agg = InstructionUtils.getMatMultOperator(threads);
+ 
+		 MatrixBlock Q_prod = new MatrixBlock(m, m, 0.0);
+		 for(int i = 0; i < m; i++) {
+			 Q_prod.set(i, i, 1.0);
+		 }
+ 
+		 for(int i = 0; i < num_iterations; i++) {
+			 MatrixBlock[] QR = computeQR2(in, threads);
+			 Q_prod = Q_prod.aggregateBinaryOperations(Q_prod, QR[0], op_mul_agg);
+			 in = QR[1].aggregateBinaryOperations(QR[1], QR[0], op_mul_agg);
+		 }
+ 
+		 // Is converged if all values are below tol and the there only is values on the diagonal.
+ 
+		 double[] check = in.getDenseBlockValues();
+		 double[] eval = new double[m];
+		 for(int i = 0; i < m; i++)
+			 eval[i] = check[i*m+i];
+		 
+		 double[] evec = Q_prod.getDenseBlockValues();
+		 return sortEVs(eval, evec);
+	 }
+ 
+	 /**
+	  * Function to compute the Householder transformation of a Matrix.
+	  *
+	  * @param in Input Matrix
+	  * @param threads number of threads
+	  * @return transformed matrix
+	  */
+	 @SuppressWarnings("unused")
+	 private static MatrixBlock computeHouseholder(MatrixBlock in, int threads) {
+		 int m = in.rlen;
+ 
+		 MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+		 A_n.copy(in);
+ 
+		 for(int k = 0; k < m - 2; k++) {
+			 MatrixBlock ajk = A_n.slice(0, m - 1, k, k);
+			 for(int i = 0; i <= k; i++) {
+				 ajk.set(i, 0, 0.0);
+			 }
+			 double alpha = Math.sqrt(ajk.sumSq());
+			 double ak1k = A_n.getDouble(k + 1, k);
+			 if(ak1k > 0.0)
+				 alpha *= -1;
+			 double r = Math.sqrt(0.5 * (alpha * alpha - ak1k * alpha));
+			 MatrixBlock v = new MatrixBlock(m, 1, 0.0);
+			 v.copy(ajk);
+			 v.set(k + 1, 0, ak1k - alpha);
+			 ScalarOperator op_div_scalar = new RightScalarOperator(Divide.getDivideFnObject(), 2 * r, threads);
+			 v = v.scalarOperations(op_div_scalar, new MatrixBlock());
+ 
+			 MatrixBlock P = new MatrixBlock(m, m, 0.0);
+			 for(int i = 0; i < m; i++) {
+				 P.set(i, i, 1.0);
+			 }
+ 
+			 ReorgOperator op_t = new ReorgOperator(SwapIndex.getSwapIndexFnObject(), threads);
+			 AggregateBinaryOperator op_mul_agg = InstructionUtils.getMatMultOperator(threads);
+			 BinaryOperator op_add = InstructionUtils.parseExtendedBinaryOperator("+");
+			 BinaryOperator op_sub = InstructionUtils.parseExtendedBinaryOperator("-");
+ 
+			 MatrixBlock v_t = v.reorgOperations(op_t, new MatrixBlock(), 0, 0, m);
+			 v_t = v_t.binaryOperations(op_add, v_t);
+			 MatrixBlock v_v_t_2 = A_n.aggregateBinaryOperations(v, v_t, op_mul_agg);
+			 P = P.binaryOperations(op_sub, v_v_t_2);
+			 A_n = A_n.aggregateBinaryOperations(P, A_n.aggregateBinaryOperations(A_n, P, op_mul_agg), op_mul_agg);
+		 }
+		 return A_n;
+	 }
+ 
+	 /**
+	  * Sort the eigen values (and vectors) in increasing order (to be compatible w/ LAPACK.DSYEVR())
+	  *
+	  * @param eValues  Eigenvalues
+	  * @param eVectors Eigenvectors
+	  * @return array of matrix blocks
+	  */
+	 private static MatrixBlock[] sortEVs(double[] eValues, double[][] eVectors) {
+		 int n = eValues.length;
+		 for(int i = 0; i < n; i++) {
+			 int k = i;
+			 double p = eValues[i];
+			 for(int j = i + 1; j < n; j++) {
+				 if(eValues[j] < p) {
+					 k = j;
+					 p = eValues[j];
+				 }
+			 }
+			 if(k != i) {
+				 eValues[k] = eValues[i];
+				 eValues[i] = p;
+				 for(int j = 0; j < n; j++) {
+					 p = eVectors[j][i];
+					 eVectors[j][i] = eVectors[j][k];
+					 eVectors[j][k] = p;
+				 }
+			 }
+		 }
+ 
+		 MatrixBlock eval = DataConverter.convertToMatrixBlock(eValues, true);
+		 MatrixBlock evec = DataConverter.convertToMatrixBlock(eVectors);
+		 return new MatrixBlock[] {eval, evec};
+	 }
+ 
+	 private static MatrixBlock[] sortEVs(double[] eValues, double[] eVectors) {
+		 int n = eValues.length;
+		 for(int i = 0; i < n; i++) {
+			 int k = i;
+			 double p = eValues[i];
+			 for(int j = i + 1; j < n; j++) {
+				 if(eValues[j] < p) {
+					 k = j;
+					 p = eValues[j];
+				 }
+			 }
+			 if(k != i) {
+				 eValues[k] = eValues[i];
+				 eValues[i] = p;
+				 for(int j = 0; j < n; j++) {
+					 p = eVectors[j*n+i];
+					 eVectors[j*n+i] = eVectors[j*n+k];
+					 eVectors[j*n+k] = p;
+				 }
+			 }
+		 }
+ 
+		 MatrixBlock eval = DataConverter.convertToMatrixBlock(eValues, true);
+		 MatrixBlock evec = new MatrixBlock(n, n, false);
+		 evec.init(eVectors, n, n);
+		 return new MatrixBlock[] {eval, evec};
+	 }
+ }
+ 
\ No newline at end of file
diff --git a/src/test/java/org/apache/sysds/test/TestUtils.java b/src/test/java/org/apache/sysds/test/TestUtils.java
index 0c407cc71a4..c96138e50df 100644
--- a/src/test/java/org/apache/sysds/test/TestUtils.java
+++ b/src/test/java/org/apache/sysds/test/TestUtils.java
@@ -2084,10 +2084,10 @@ public static MatrixBlock generateTestMatrixBlock(int rows, int cols, double min
 	 * @param seed seed
 	 * @return random symmetric MatrixBlock
 	 */
-	public static MatrixBlock generateTestMatrixBlockSym(int rows, int cols, double min, double max, double sparsity, long seed){
-		MatrixBlock m = MatrixBlock.randOperations(rows, cols, sparsity, min, max, "Uniform", seed);
-		for(int i = 0; i < rows; i++) {
-			for(int j = i+1; j < cols; j++) {
+	public static MatrixBlock generateTestMatrixBlockSym(int size, double min, double max, double sparsity, long seed){
+		MatrixBlock m = MatrixBlock.randOperations(size, size, sparsity, min, max, "Uniform", seed);
+		for(int i = 0; i < size; i++) {
+			for(int j = i+1; j < size; j++) {
 				m.set(i,j, m.get(j,i));
 			}
 		}
diff --git a/src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java b/src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java
index 6292a141389..8da74a4477a 100644
--- a/src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java
+++ b/src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java
@@ -37,31 +37,34 @@ public class EigenDecompTest {
 	protected static final Log LOG = LogFactory.getLog(EigenDecompTest.class.getName());
 
 	private enum type {
-		COMMONS, LANCZOS, QR,
+		COMMONS, LANCZOS, QR, COMMONS_NATIVE
 	}
 
 	@Test
+	@Ignore
 	public void testLanczosSimple() {
-		MatrixBlock in = new MatrixBlock(4, 4, new double[] {4, 1, -2, 2, 1, 2, 0, 1, -2, 0, 3, -2, 2, 1, -2, -1});
+		MatrixBlock in = new MatrixBlock(4, 4, new double[] { 4, 1, -2, 2, 1, 2, 0, 1, -2, 0, 3, -2, 2, 1, -2, -1 });
 		testEigen(in, 1e-4, 1, type.LANCZOS);
 	}
 
 	@Test
+	@Ignore
 	public void testLanczosSmall() {
-		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10, 0.0, 1.0, 1.0, 1);
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 0.0, 1.0, 1.0, 1);
 		testEigen(in, 1e-4, 1, type.LANCZOS);
 	}
 
 	@Test
+	@Ignore
 	public void testLanczosMedium() {
-		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(12, 12, 0.0, 1.0, 1.0, 1);
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(12, 0.0, 1.0, 1.0, 1);
 		testEigen(in, 1e-4, 1, type.LANCZOS);
 	}
 
 	@Test
 	@Ignore
 	public void testLanczosLarge() {
-		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100, 100, 0.0, 1.0, 1.0, 1);
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100, 0.0, 1.0, 1.0, 1);
 		testEigen(in, 1e-4, 1, type.LANCZOS);
 	}
 
@@ -69,34 +72,75 @@ public void testLanczosLarge() {
 	@Ignore
 	public void testLanczosLargeMT() {
 		int threads = 10;
-		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100, 100, 0.0, 1.0, 1.0, 1);
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100, 0.0, 1.0, 1.0, 1);
 		testEigen(in, 1e-4, threads, type.LANCZOS);
 	}
 
 	@Test
+	@Ignore
 	public void testQREigenSimple() {
 		MatrixBlock in = new MatrixBlock(4, 4,
-			new double[] {52, 30, 49, 28, 30, 50, 8, 44, 49, 8, 46, 16, 28, 44, 16, 22});
+				new double[] { 52, 30, 49, 28, 30, 50, 8, 44, 49, 8, 46, 16, 28, 44, 16, 22 });
 		testEigen(in, 1e-4, 1, type.QR);
 	}
 
 	@Test
+	@Ignore
 	public void testQREigenSymSmall() {
-		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10, 0.0, 1.0, 1.0, 1);
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 0.0, 1.0, 1.0, 1);
 		testEigen(in, 1e-3, 1, type.QR);
 	}
 
 	@Test
+	@Ignore
 	public void testQREigenSymSmallMT() {
 		int threads = 10;
-		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10, 0.0, 1.0, 1.0, 1);
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 0.0, 1.0, 1.0, 1);
 		testEigen(in, 1e-3, threads, type.QR);
 	}
 
+	@Test
+	public void testEigenSymSmall() {
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 0.0, 1.0, 1.0, 1);
+		testEigen(in, 1e-4, 10, type.COMMONS);
+	}
+
+	@Test
+	public void testEigenSymMedium() {
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(200, 0.0, 1.0, 1.0, 1);
+		testEigen(in, 1e-4, 10, type.COMMONS);
+	}
+
+	@Test
+	public void testEigenSymLarge() {
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(1000, 0.0, 1.0, 1.0, 1);
+		testEigen(in, 1e-4, 10, type.COMMONS);
+	}
+
+	@Test
+	public void testNativeEigenSymSmall() {
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 0.0, 1.0, 1.0, 1);
+		testEigen(in, 1e-4, 10, type.COMMONS_NATIVE);
+	}
+
+	@Test
+	public void testNativeEigenSymMedium() {
+		int threads = 10;
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(200, 0.0, 1.0, 1.0, 1);
+		testEigen(in, 1e-4, threads, type.COMMONS_NATIVE);
+	}
+
+	@Test
+	public void testNativeEigenSymLarge() {
+		int threads = 10;
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(1000, 0.0, 1.0, 1.0, 1);
+		testEigen(in, 1e-4, threads, type.COMMONS_NATIVE);
+	}
+
 	@Test
 	@Ignore
 	public void testQREigenSymLarge() {
-		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 50, 0.0, 1.0, 1.0, 1);
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 0.0, 1.0, 1.0, 1);
 		testEigen(in, 1e-4, 1, type.QR);
 	}
 
@@ -104,14 +148,14 @@ public void testQREigenSymLarge() {
 	@Ignore
 	public void testQREigenSymLargeMT() {
 		int threads = 10;
-		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 50, 0.0, 1.0, 1.0, 1);
+		MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 0.0, 1.0, 1.0, 1);
 		testEigen(in, 1e-4, threads, type.QR);
 	}
 
-	private void testEigen(MatrixBlock in, double tol, int threads, type t) {
+	private void testEigen(MatrixBlock in, double tol, int threads, type t, boolean validate) {
 		try {
 			MatrixBlock[] m = null;
-			switch(t) {
+			switch (t) {
 				case COMMONS:
 					m = LibCommonsMath.multiReturnOperations(in, "eigen", threads, 1);
 					break;
@@ -121,19 +165,26 @@ private void testEigen(MatrixBlock in, double tol, int threads, type t) {
 				case QR:
 					m = LibCommonsMath.multiReturnOperations(in, "eigen_qr", threads, 1);
 					break;
+				case COMMONS_NATIVE:
+					m = LibCommonsMath.multiReturnOperations(in, "eigen_symm", threads, 1);
+					break;
 				default:
 					throw new NotImplementedException();
 			}
 
-			isValidDecomposition(in, m[1], m[0], tol, t.toString());
+			if (validate)
+				isValidDecomposition(in, m[1], m[0], tol, t.toString());
 
-		}
-		catch(Exception e) {
+		} catch (Exception e) {
 			e.printStackTrace();
 			fail(e.getMessage());
 		}
 	}
 
+	private void testEigen(MatrixBlock in, double tol, int threads, type t) {
+		testEigen(in, tol, threads, t, true);
+	}
+
 	private void isValidDecomposition(MatrixBlock A, MatrixBlock V, MatrixBlock vD, double tol, String message) {
 		// Any eigen decomposition is valid if A = V %*% D %*% t(V)
 		// A is the input of the eigen decomposition