Given a symm matrix S that stores just its subdiagonal elements,
reconstructs the full symmetric matrix.
@FIXME: evitare il doppio setQuick
11940
public static final DoubleMatrix2D fillSubdiagonalSymmetricMatrix(DoubleMatrix2D S){

		if(S.rows() != S.columns()){
			throw new IllegalArgumentException("Not square matrix");
		}

		boolean isSparse = S instanceof SparseDoubleMatrix2D;
		DoubleFactory2D F2D = (isSparse)? DoubleFactory2D.sparse: DoubleFactory2D.dense;
		final DoubleMatrix2D SFull = F2D.make(S.rows(), S.rows());

		if (isSparse) {
			S.forEachNonZero(new IntIntDoubleFunction() {
				@Override
				public double apply(int i, int j, double hij) {
					SFull.setQuick(i, j, hij);
					SFull.setQuick(j, i, hij);
					return hij;
				}
			});
		} else {
			for (int i = 0; i < S.rows(); i++) {
				for (int j = 0; j < i + 1; j++) {
					double sij = S.getQuick(i, j);
					SFull.setQuick(i, j, sij);
					SFull.setQuick(j, i, sij);
				}
			}
		}

		return SFull;
	}
679722
public Matrix solve(final Matrix rhs)
    {
        if ((m_rows != m_columns) || (rhs.m_rows != m_columns) || (rhs.m_columns != 1))
        {
            throw new GeometryException("Illegal matrix dimensions");
        }
        // create copies of the data

        final Matrix A = new Matrix(this);

        final Matrix b = new Matrix(rhs);

        // Gaussian elimination with partial pivoting

        for (int i = 0; i < m_columns; i++)
        {
            // find pivot row and swap

            int max = i;

            for (int j = i + 1; j < m_columns; j++)
            {
                if (Math.abs(A.m_data[j][i]) > Math.abs(A.m_data[max][i]))
                {
                    max = j;
                }
            }
            A.swap(i, max);

            b.swap(i, max);

            // singular

            if (A.m_data[i][i] == 0.0)
            {
                throw new RuntimeException("Matrix is singular.");
            }
            // pivot within b

            for (int j = i + 1; j < m_columns; j++)
            {
                b.m_data[j][0] -= (b.m_data[i][0] * A.m_data[j][i]) / A.m_data[i][i];
            }
            // pivot within A

            for (int j = i + 1; j < m_columns; j++)
            {
                final double m = A.m_data[j][i] / A.m_data[i][i];

                for (int k = i + 1; k < m_columns; k++)
                {
                    A.m_data[j][k] -= A.m_data[i][k] * m;
                }
                A.m_data[j][i] = 0.0;
            }
        }
        // back substitution

        final Matrix x = new Matrix(m_columns, 1);

        for (int j = m_columns - 1; j >= 0; j--)
        {
            double t = 0.0;

            for (int k = j + 1; k < m_columns; k++)
            {
                t += A.m_data[j][k] * x.m_data[k][0];
            }
            x.m_data[j][0] = (b.m_data[j][0] - t) / A.m_data[j][j];
        }
        return x;
    }
19444
function(sbmA, sbmB) {
    // if not the same sizes
    if (sbmA.x != sbmB.x || sbmA.y != sbmB.y) {
      return
    }
    // create the matrix
    var mat = {}
    mat.x = sbmA.x
    mat.y = sbmA.y
    mat.data = []
    // populate
    for (var i = 0; i < mat.x; i++) {
      mat.data[i] = _.union(sbmA.data[i], sbmB.data[i]).sort()
    }
    // return matrix
    return mat
  }
324694
public Matrix solve(Matrix rhs) {
        if (M != N || rhs.M != N || rhs.N != 1) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }

        // create copies of the data
        Matrix A = new Matrix(this);
        Matrix b = new Matrix(rhs);

        // 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.data[j][i]) > Math.abs(A.data[max][i])) {
                    max = j;
                }
            }
            A.swap(i, max);
            b.swap(i, max);

            // singular
            if (A.data[i][i] == 0.0) {
                throw new RuntimeException("Matrix is singular.");
            }

            // pivot within b
            for (int j = i + 1; j < N; j++) {
                b.data[j][0] -= b.data[i][0] * A.data[j][i] / A.data[i][i];
            }

            // pivot within A
            for (int j = i + 1; j < N; j++) {
                double m = A.data[j][i] / A.data[i][i];
                for (int k = i + 1; k < N; k++) {
                    A.data[j][k] -= A.data[i][k] * m;
                }
                A.data[j][i] = 0.0;
            }
        }

        // back substitution
        Matrix x = new Matrix(N, 1);
        for (int j = N - 1; j >= 0; j--) {
            double t = 0.0;
            for (int k = j + 1; k < N; k++) {
                t += A.data[j][k] * x.data[k][0];
            }
            x.data[j][0] = (b.data[j][0] - t) / A.data[j][j];
        }
        return x;

    }
343787
def inverse(self):
        
        invr = np.linalg.inv(self.affine_matrix)
        return SymmOp(invr)
860019
public void decompose( DMatrixRMaj A ) {

        this.QR = A.copy();

        int N = Math.min(A.numCols,A.numRows);

        gammas = new double[ A.numCols ];

        DMatrixRMaj A_small = new DMatrixRMaj(A.numRows,A.numCols);
        DMatrixRMaj A_mod = new DMatrixRMaj(A.numRows,A.numCols);
        DMatrixRMaj v = new DMatrixRMaj(A.numRows,1);
        DMatrixRMaj Q_k = new DMatrixRMaj(A.numRows,A.numRows);

        for( int i = 0; i < N; i++ ) {
            // reshape temporary variables
            A_small.reshape(QR.numRows-i,QR.numCols-i);
            v.reshape(A_small.numRows,1);
            Q_k.reshape(v.getNumElements(),v.getNumElements());

            // use extract matrix to get the column that is to be zeroed
            CommonOps_DDRM.extract(QR,i,QR.numRows,i,i+1,v,0,0);

            double max = CommonOps_DDRM.elementMaxAbs(v);

            if( max > 0 && v.getNumElements() > 1 ) {
                // normalize to reduce overflow issues
                CommonOps_DDRM.divide(v,max);

                // compute the magnitude of the vector
                double tau = NormOps_DDRM.normF(v);

                if( v.get(0) < 0 )
                    tau *= -1.0;

                double u_0 = v.get(0) + tau;
                double gamma = u_0/tau;

                CommonOps_DDRM.divide(v,u_0);
                v.set(0,1.0);

                // extract the submatrix of A which is being operated on
                CommonOps_DDRM.extract(QR,i,QR.numRows,i,QR.numCols,A_small,0,0);

                // A = (I - &gamma;*u*u<sup>T</sup>)A
                CommonOps_DDRM.setIdentity(Q_k);
                CommonOps_DDRM.multAddTransB(-gamma,v,v,Q_k);
                CommonOps_DDRM.mult(Q_k,A_small,A_mod);

                // save the results
                CommonOps_DDRM.insert(A_mod, QR, i,i);
                CommonOps_DDRM.insert(v, QR, i,i);
                QR.unsafe_set(i,i,-tau*max);

                // save gamma for recomputing Q later on
                gammas[i] = gamma;
            }
        }
    }
571815
def reconstruct_emds(edm, Om, all_points, method=None, **kwargs):
    
    from .point_set import dm_from_edm
    N = all_points.shape[0]
    d = all_points.shape[1]
    dm = dm_from_edm(edm)
    if method is None:
        from .mds import superMDS
        Xhat, __ = superMDS(all_points[0, :], N, d, Om=Om, dm=dm)
    else:
        C = kwargs.get('C', None)
        b = kwargs.get('b', None)
        if C is None or b is None:
            raise NameError(
                'Need constraints C and b for reconstruct_emds in iterative mode.')
        KE_noisy = np.multiply(np.outer(dm, dm), Om)
        if method == 'iterative':
            from .mds import iterativeEMDS
            Xhat, __ = iterativeEMDS(
                all_points[0, :], N, d, KE=KE_noisy, C=C, b=b)
        elif method == 'relaxed':
            from .mds import relaxedEMDS
            Xhat, __ = relaxedEMDS(
                all_points[0, :], N, d, KE=KE_noisy, C=C, b=b)
        else:
            raise NameError('Undefined method', method)
    Y, R, t, c = procrustes(all_points, Xhat, scale=False)
    return Y
1035019
private static int[] elmhes(DenseMatrix A) {
        int n = A.nrows();
        int[] perm = new int[n];

        for (int m = 1; m < n - 1; m++) {
            double x = 0.0;
            int i = m;
            for (int j = m; j < n; j++) {
                if (Math.abs(A.get(j, m - 1)) > Math.abs(x)) {
                    x = A.get(j, m - 1);
                    i = j;
                }
            }
            perm[m] = i;
            if (i != m) {
                for (int j = m - 1; j < n; j++) {
                    double swap = A.get(i, j);
                    A.set(i, j, A.get(m, j));
                    A.set(m, j, swap);
                }
                for (int j = 0; j < n; j++) {
                    double swap = A.get(j, i);
                    A.set(j, i, A.get(j, m));
                    A.set(j, m, swap);
                }
            }
            if (x != 0.0) {
                for (i = m + 1; i < n; i++) {
                    double y = A.get(i, m - 1);
                    if (y != 0.0) {
                        y /= x;
                        A.set(i, m - 1, y);
                        for (int j = m; j < n; j++) {
                            A.sub(i, j, y * A.get(m, j));
                        }
                        for (int j = 0; j < n; j++) {
                            A.add(j, m, y * A.get(j, i));
                        }
                    }
                }
            }
        }

        return perm;
    }
1261167
def eigenvectors_from_displacements(disp,masses):
    
    nphonons,natoms,ndirections = disp.shape
    sqrt_masses = np.sqrt(masses)
    return np.einsum("nax,a->nax",disp,sqrt_masses)
322757
def full_s(self):
        
        x = np.zeros((self.shape),dtype=np.float32)

        x[:self.s.shape[0],:self.s.shape[0]] = self.s.as_2d
        s = Matrix(x=x, row_names=self.row_names,
                          col_names=self.col_names, isdiagonal=False,
                          autoalign=False)
        return s
111486
function (projectionMatrix) {
        Matrix4.copy(this.projectionMatrix, projectionMatrix);
        Matrix4.invert(this.invProjectionMatrix, projectionMatrix);
        this.decomposeProjectionMatrix();
    }
898706
function(sbm) {
    // if the matrix ain't square, symmetry is off
    if (sbm.x !== sbm.y) {
      return false
    }

    // symmetry status
    var symmetry = true
    // iterate
    for (var i = 0; i < sbm.x; i++) {
      for (var j = 0; j <= i; j++) {
        symmetry = symmetry && ((_.includes(sbm.data[i], j) ^ _.includes(sbm.data[j], i)) === 0)
      }
    }
    // return
    return symmetry
  }
732322
def inversion(origin=(0, 0, 0)):
        
        mat = -np.eye(4)
        mat[3, 3] = 1
        mat[0:3, 3] = 2 * np.array(origin)
        return SymmOp(mat)
320729
function svdPack() {
                var result = lapack$1.sgesvd('A', 'A', this.elements);

                return {
                    U: $M(result.U),
                    S: $M(result.S).column(1).toDiagonalMatrix(),
            	V: $M(result.VT).transpose()
                };
            }
473232
def contained_segments_matrix(segments):
    
    x1, y1 = segments[:, 0], segments[:, 1]
    x2, y2 = x1 + segments[:, 2], y1 + segments[:, 3]
    n = len(segments)

    x1so, x2so, y1so, y2so = list(map(numpy.argsort, (x1, x2, y1, y2)))
    x1soi, x2soi, y1soi, y2soi = list(map(numpy.argsort, (x1so, x2so, y1so, y2so)))  # inverse transformations
    # let rows be x1 and collumns be x2. this array represents where x1<x2
    o1 = numpy.triu(numpy.ones((n, n)), k=1).astype(bool)
    # let rows be x1 and collumns be x2. this array represents where x1>x2
    o2 = numpy.tril(numpy.ones((n, n)), k=0).astype(bool)
    a_inside_b_x = o2[x1soi][:, x1soi] * o1[x2soi][:, x2soi]  # (x1[a]>x1[b] and x2[a]<x2[b])
    a_inside_b_y = o2[y1soi][:, y1soi] * o1[y2soi][:, y2soi]  # (y1[a]>y1[b] and y2[a]<y2[b])
    a_inside_b = a_inside_b_x * a_inside_b_y
    return a_inside_b
721457
public void factorize(boolean checkSymmetry) throws Exception{
		if (checkSymmetry && !Property.TWELVE.isSymmetric(Q)) {
			throw new Exception("Matrix is not symmetric");
		}
		
		if(this.rescaler != null){
			double[] cn_00_original = null;
			double[] cn_2_original = null;
			double[] cn_00_scaled = null;
			double[] cn_2_scaled = null;
			if(log.isDebugEnabled()){
				cn_00_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
				log.debug("cn_00_original Q before scaling: " + ArrayUtils.toString(cn_00_original));
				cn_2_original = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
				log.debug("cn_2_original Q before scaling : " + ArrayUtils.toString(cn_2_original));
			}
			//scaling the Q matrix, we have:
			//Q1 = U.Q.U[T] = U.L.L[T].U[T] = (U.L).(U.L)[T] 
			//and because U is diagonal it preserves the triangular form of U.L, so
			//Q1 = U.Q.U[T] = L1.L1[T] is the new Cholesky decomposition  
			DoubleMatrix1D Uv = rescaler.getMatrixScalingFactorsSymm(Q);
			if(log.isDebugEnabled()){
				boolean checkOK = rescaler.checkScaling(ColtUtils.fillSubdiagonalSymmetricMatrix(Q), Uv, Uv);
				if(!checkOK){
					log.warn("Scaling failed (checkScaling = false)");
				}
			}
			this.U = Uv;
			this.Q = ColtUtils.diagonalMatrixMult(Uv, Q, Uv);
			if(log.isDebugEnabled()){
				cn_00_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), Integer.MAX_VALUE);
				log.debug("cn_00_scaled Q after scaling : " + ArrayUtils.toString(cn_00_scaled));
				cn_2_scaled = ColtUtils.getConditionNumberRange(new Array2DRowRealMatrix(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()), 2);
				log.debug("cn_2_scaled Q after scaling  : " + ArrayUtils.toString(cn_2_scaled));
				
				if(cn_00_original[0] < cn_00_scaled[0] || cn_2_original[0] < cn_2_scaled[0]){
					//log.info("Q: " + ArrayUtils.toString(ColtUtils.fillSubdiagonalSymmetricMatrix(Q).toArray()));
					log.warn("Problematic scaling");
					//throw new RuntimeException("Scaling failed");
				}
			}
		}
		
		//double threshold = Math.pow(Utils.getDoubleMachineEpsilon(), 2);
		double threshold = Utils.getDoubleMachineEpsilon();
		this.LData = new double[dim][];

		for (int i = 0; i < dim; i++) {
			LData[i] = new double[i+1];
			//DoubleMatrix1D QI = Q.viewRow(i);
			double[] LDataI = LData[i];
			//j < i
			for (int j = 0; j < i; j++) {
				double[] LDataJ = LData[j];
				double sum = 0.0;
				for (int k = 0; k < j; k++) {
					sum += LDataI[k] * LDataJ[k];
				}
				LDataI[j] = 1.0 / LDataJ[j] * (Q.getQuick(i, j) - sum);
			}
			//j==i
			double sum = 0.0;
			for (int k = 0; k < i; k++) {
				sum += Math.pow(LDataI[k], 2);
			}
			double d = Q.getQuick(i, i) - sum;
			if(!(d > threshold)){
				throw new Exception("not positive definite matrix");
			}
			LDataI[i] = Math.sqrt(d);
		}
	}
1088348
@Override
    public Matrix[] decompose() {

        // AHTIUNG: this code derived from Jama

        Matrix a = matrix.copy();

        int n = Math.min(a.rows(), a.columns());

        Matrix u = matrix.blankOfShape(a.rows(), n);
        Matrix s = matrix.blankOfShape(a.columns(), a.columns());
        Matrix v = matrix.blankOfShape(a.columns(), a.columns());

        Vector e = DenseVector.zero(a.columns());
        Vector work = DenseVector.zero(a.rows());

        int nct = Math.min(a.rows() - 1, a.columns());
        int nrt = Math.max(0, Math.min(a.columns() - 2, a.rows()));

        for (int k = 0; k < Math.max(nct, nrt); k++) {

            if (k < nct) {

                s.set(k, k, 0.0);

                for (int i = k; i < a.rows(); i++) {
                    s.set(k, k, hypot(s.get(k, k), a.get(i, k)));
                }

                if (Math.abs(s.get(k, k)) > Matrices.EPS) {

                    if (a.get(k, k) < 0.0) {
                        s.updateAt(k, k, Matrices.INV_FUNCTION);
                    }

                    double skk = s.get(k, k);
                    for (int i = k; i < a.rows(); i++) {
                        a.updateAt(i, k, Matrices.asDivFunction(skk));
                    }

                    a.updateAt(k, k, Matrices.INC_FUNCTION);
                }

                s.updateAt(k, k, Matrices.INV_FUNCTION);
            }

            for (int j = k + 1; j < a.columns(); j++) {

                if ((k < nct) & (Math.abs(s.get(k, k)) > Matrices.EPS)) {

                    double t = 0;

                    for (int i = k; i < a.rows(); i++) {
                        t += a.get(i, k) * a.get(i, j);
                    }

                    t = -t / a.get(k, k);

                    for (int i = k; i < a.rows(); i++) {
                        a.updateAt(i, j, Matrices.asPlusFunction(t * a.get(i, k)));
                    }
                }

                e.set(j, a.get(k, j));
            }

            if (k < nct) {

                for (int i = k; i < a.rows(); i++) {
                    u.set(i, k, a.get(i, k));
                }

            }

            if (k < nrt) {

                e.set(k, 0);

                for (int i = k + 1; i < a.columns(); i++) {
                    e.set(k, hypot(e.get(k), e.get(i)));
                }

                if (Math.abs(e.get(k)) > Matrices.EPS) {

                    if (e.get(k + 1) < 0.0) {

                        e.updateAt(k, Vectors.INV_FUNCTION);
                    }

                    double ek = e.get(k);
                    for (int i = k + 1; i < a.columns(); i++) {
                        e.updateAt(i, Vectors.asDivFunction(ek));
                    }

                    e.updateAt(k + 1, Vectors.INC_FUNCTION);
                }

                e.updateAt(k, Vectors.INV_FUNCTION);

                if ((k + 1 < a.rows()) && (Math.abs(e.get(k)) > Matrices.EPS)) {

                    for (int i = k + 1; i < a.rows(); i++) {
                        work.set(i, 0.0);
                    }

                    for (int j = k + 1; j < a.columns(); j++) {
                        for (int i = k + 1; i < a.rows(); i++) {
                            work.updateAt(i, Vectors.asPlusFunction(e.get(j) *
                                a.get(i, j)));
                        }
                    }

                    for (int j = k + 1; j < a.columns(); j++) {

                        double t = -e.get(j) / e.get(k + 1);

                        for (int i = k + 1; i < a.rows(); i++) {
                            a.updateAt(i, j, Matrices.asPlusFunction(t *
                                work.get(i)));
                        }
                    }
                }

                for (int i = k + 1; i < a.columns(); i++) {
                    v.set(i, k, e.get(i));
                }
            }
        }

        int p = Math.min(a.columns(), a.rows() + 1);

        if (nct < a.columns()) {
            s.set(nct, nct, a.get(nct, nct));
        }

        if (a.rows() < p) {
            s.set(p - 1, p - 1, 0.0);
        }

        if (nrt + 1 < p) {
            e.set(nrt, a.get(nrt, p - 1));
        }

        e.set(p - 1, 0.0);

        for (int j = nct; j < n; j++) {

            for (int i = 0; i < a.rows(); i++) {
                u.set(i, j, 0.0);
            }

            u.set(j, j, 1.0);
        }

        for (int k = nct - 1; k >= 0; k--) {

            if (Math.abs(s.get(k, k)) > Matrices.EPS) {

                for (int j = k + 1; j < n; j++) {

                    double t = 0;
                    for (int i = k; i < a.rows(); i++) {
                        t += u.get(i, k) * u.get(i, j);
                    }

                    t = -t / u.get(k, k);

                    for (int i = k; i < a.rows(); i++) {
                        u.updateAt(i, j, Matrices.asPlusFunction(t * u.get(i, k)));
                    }
                }

                for (int i = k; i < a.rows(); i++) {
                    u.updateAt(i, k, Matrices.INV_FUNCTION);
                }

                u.updateAt(k, k, Matrices.INC_FUNCTION);

                for (int i = 0; i < k - 1; i++) {
                    u.set(i, k, 0.0);
                }

            } else {

                for (int i = 0; i < a.rows(); i++) {
                    u.set(i, k, 0.0);
                }

                u.set(k, k, 1.0);
            }
        }

        for (int k = n - 1; k >= 0; k--) {

            if ((k < nrt) & (Math.abs(e.get(k)) > Matrices.EPS)) {

                for (int j = k + 1; j < n; j++) {

                    double t = 0;

                    for (int i = k + 1; i < a.columns(); i++) {
                        t += v.get(i, k) * v.get(i, j);
                    }

                    t = -t / v.get(k + 1, k);

                    for (int i = k + 1; i < a.columns(); i++) {
                        v.updateAt(i, j, Matrices.asPlusFunction(t * v.get(i, k)));
                    }
                }
            }

            for (int i = 0; i < a.columns(); i++) {
                v.set(i, k, 0.0);
            }

            v.set(k, k, 1.0);
        }

        int pp = p - 1;
        int iter = 0;
        double eps = Math.pow(2.0, -52.0);
        double tiny = Math.pow(2.0, -966.0);

        while (p > 0) {

            int k;
            int kase;

            for (k = p - 2; k >= -1; k--) {
                if (k == -1)
                    break;

                if (Math.abs(e.get(k)) <= tiny + eps * 
                        (Math.abs(s.get(k, k)) + Math.abs(s.get(k + 1, k + 1)))) {

                    e.set(k, 0.0);
                    break;
                }
            }

            if (k == p - 2) {

                kase = 4;

            } else {

                int ks;

                for (ks = p - 1; ks >= k; ks--) {

                    if (ks == k)
                        break;

                    double t = (ks != p ? Math.abs(e.get(ks)) : 0.)
                            + (ks != k + 1 ? Math.abs(e.get(ks - 1)) 
                            			   : 0.);

                    if (Math.abs(s.get(ks, ks)) <= tiny + eps * t) {
                        s.set(ks, ks, 0.0);
                        break;
                    }
                }

                if (ks == k) {
                    kase = 3;
                } else if (ks == p - 1) {
                    kase = 1;
                } else {
                    kase = 2;
                    k = ks;
                }
            }

            k++;

            switch (kase) {

            case 1: {
                double f = e.get(p - 2);
                e.set(p - 2, 0.0);

                for (int j = p - 2; j >= k; j--) {

                    double t = hypot(s.get(j, j), f);
                    double cs = s.get(j, j) / t;
                    double sn = f / t;

                    s.set(j, j, t);

                    if (j != k) {
                        f = -sn * e.get(j - 1);
                        e.set(j - 1, cs * e.get(j - 1));
                    }

                    for (int i = 0; i < a.columns(); i++) {
                        t = cs * v.get(i, j) + sn 
                        		* v.get(i, p - 1);
                        v.set(i, p - 1,
                                -sn * v.get(i, j) 
                                + cs * v.get(i, p - 1));
                        v.set(i, j, t);
                    }
                }
            }
                break;

            case 2: {
                double f = e.get(k - 1);
                e.set(k - 1, 0.0);

                for (int j = k; j < p; j++) {

                    double t = hypot(s.get(j, j), f);
                    double cs = s.get(j, j) / t;
                    double sn = f / t;

                    s.set(j, j, t);
                    f = -sn * e.get(j);
                    e.set(j, cs * e.get(j));

                    for (int i = 0; i < a.rows(); i++) {
                        t = cs * u.get(i, j) 
                        		+ sn * u.get(i, k - 1);
                        u.set(i, k - 1,
                                -sn * u.get(i, j) 
                                + cs * u.get(i, k - 1));
                        u.set(i, j, t);
                    }
                }
            }
                break;

            case 3: {

                double scale = Math.max(
                        Math.max(Math.max(
                                Math.max(Math.abs(s.get(p - 1, p - 1)),
                                        Math.abs(s.get(p - 2, p - 2))),
                                    Math.abs(e.get(p - 2))), 
                                Math.abs(s.get(k, k))),
                        Math.abs(e.get(k)));

                double sp = s.get(p - 1, p - 1) / scale;
                double spm1 = s.get(p - 2, p - 2) / scale;
                double epm1 = e.get(p - 2) / scale;
                double sk = s.get(k, k) / scale;
                double ek = e.get(k) / scale;
                double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
                double c = (sp * epm1) * (sp * epm1);
                double shift = 0.0;

                if ((b != 0.0) | (c != 0.0)) {
                    shift = Math.sqrt(b * b + c);
                    if (b < 0.0) {
                        shift = -shift;
                    }
                    shift = c / (b + shift);
                }

                double f = (sk + sp) * (sk - sp) + shift;
                double g = sk * ek;

                for (int j = k; j < p - 1; j++) {
                    double t = hypot(f, g);
                    double cs = f / t;
                    double sn = g / t;

                    if (j != k) {
                        e.set(j - 1, t);
                    }

                    f = cs * s.get(j, j) + sn * e.get(j);
                    e.set(j, 
                    		cs * e.get(j) - sn * s.get(j, j));
                    g = sn * s.get(j + 1, j + 1);
                    s.set(j + 1, j + 1, 
                    		cs * s.get(j + 1, j + 1));

                    for (int i = 0; i < a.columns(); i++) {
                        t = cs * v.get(i, j) 
                        		+ sn * v.get(i, j + 1);
                        v.set(i, j + 1,
                                -sn * v.get(i, j) 
                                + cs * v.get(i, j + 1));
                        v.set(i, j, t);
                    }

                    t = hypot(f, g);
                    cs = f / t;
                    sn = g / t;
                    s.set(j, j, t);
                    f = cs * e.get(j) + sn * s.get(j + 1, j + 1);
                    s.set(j + 1, j + 1, 
                            -sn * e.get(j) 
                            + cs * s.get(j + 1, j + 1));
                    g = sn * e.get(j + 1);
                    e.updateAt(j + 1, Vectors.asMulFunction(cs));

                    if (j < a.rows() - 1) {
                        for (int i = 0; i < a.rows(); i++) {
                            t = cs * u.get(i, j) 
                                    + sn * u.get(i, j + 1);
                            u.set(i, j + 1,
                                    -sn * u.get(i, j) 
                                    + cs * u.get(i, j + 1));
                            u.set(i, j, t);
                        }
                    }
                }

                e.set(p - 2, f);
                iter = iter + 1;
            }
                break;

            case 4: {

                if (s.get(k, k) <= 0.0) {
                    s.set(k, k, s.get(k, k) < 0.0 ? -s.get(k, k) : 0.0);
                    for (int i = 0; i <= pp; i++) {
                        v.updateAt(i, k, Matrices.INV_FUNCTION);
                    }
                }

                while (k < pp) {

                    if (s.get(k, k) >= s.get(k + 1, k + 1)) {
                        break;
                    }

                    double t = s.get(k, k);
                    s.set(k, k, s.get(k + 1, k + 1));
                    s.set(k + 1, k + 1, t);

                    if (k < a.columns() - 1) {
                        for (int i = 0; i < a.columns(); i++) {
                            t = v.get(i, k + 1);
                            v.set(i, k + 1, v.get(i, k));
                            v.set(i, k, t);
                        }
                    }

                    if (k < a.rows() - 1) {
                        for (int i = 0; i < a.rows(); i++) {
                            t = u.get(i, k + 1);
                            u.set(i, k + 1, u.get(i, k));
                            u.set(i, k, t);
                        }
                    }

                    k++;
                }

                iter = 0;
                p--;
            }
                break;
            }
        }

        // TODO:
        // We need to change the logic of this code in order do not use resize
        // It is also not a good idea to use SVN for [m < n] matrices
        // We need to change it in further releases
        return new Matrix[] { u, s.copyOfShape(n, a.columns()), v };
    }