Returns the result of QR decomposition of given matrix
<p>
See <a href="http://mathworld.wolfram.com/QRDecomposition.html">
http://mathworld.wolfram.com/QRDecomposition.html</a> for more details.
</p>

@return { Q, R }
777179
@Override
    public Matrix[] decompose() {

        Matrix[] qrr = super.decompose();
        Matrix qr = qrr[0];
        Matrix r = qrr[1];

        Matrix q = qr.blank();

        for (int k = q.columns() - 1; k >= 0; k--) {

            q.set(k, k, 1.0);

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

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

                    double acc = 0.0;

                    for (int i = k; i < q.rows(); i++) {
                        acc += qr.get(i, k) * q.get(i, j);
                    }

                    acc = -acc / qr.get(k, k);

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

        for (int i = 0; i < r.rows(); i++) {
            for (int j = i + 1; j < r.columns(); j++) {
                r.set(i, j, qr.get(i, j));
            }
        }

        return new Matrix[] { q, r };
    }
173164
function decomposeMatrix(matrix) {
    var m3d = [
      matrix.slice(0, 4),
      matrix.slice(4, 8),
      matrix.slice(8, 12),
      matrix.slice(12, 16)
    ];

    // skip normalization step as m3d[3][3] should always be 1
    if (m3d[3][3] !== 1) {
      throw 'attempt to decompose non-normalized matrix';
    }

    var perspectiveMatrix = m3d.concat(); // copy m3d
    for (var i = 0; i < 3; i++) {
      perspectiveMatrix[i][3] = 0;
    }

    if (determinant(perspectiveMatrix) === 0) {
      return false;
    }

    var rhs = [];

    var perspective;
    if (m3d[0][3] !== 0 || m3d[1][3] !== 0 || m3d[2][3] !== 0) {
      rhs.push(m3d[0][3]);
      rhs.push(m3d[1][3]);
      rhs.push(m3d[2][3]);
      rhs.push(m3d[3][3]);

      var inversePerspectiveMatrix = inverse(perspectiveMatrix);
      var transposedInversePerspectiveMatrix =
          transposeMatrix4(inversePerspectiveMatrix);
      perspective = multVecMatrix(rhs, transposedInversePerspectiveMatrix);
    } else {
      perspective = [0, 0, 0, 1];
    }

    var translate = m3d[3].slice(0, 3);

    var row = [];
    row.push(m3d[0].slice(0, 3));
    var scale = [];
    scale.push(length(row[0]));
    row[0] = normalize(row[0]);

    var skew = [];
    row.push(m3d[1].slice(0, 3));
    skew.push(dot(row[0], row[1]));
    row[1] = combine(row[1], row[0], 1.0, -skew[0]);

    scale.push(length(row[1]));
    row[1] = normalize(row[1]);
    skew[0] /= scale[1];

    row.push(m3d[2].slice(0, 3));
    skew.push(dot(row[0], row[2]));
    row[2] = combine(row[2], row[0], 1.0, -skew[1]);
    skew.push(dot(row[1], row[2]));
    row[2] = combine(row[2], row[1], 1.0, -skew[2]);

    scale.push(length(row[2]));
    row[2] = normalize(row[2]);
    skew[1] /= scale[2];
    skew[2] /= scale[2];

    var pdum3 = cross(row[1], row[2]);
    if (dot(row[0], pdum3) < 0) {
      for (var i = 0; i < 3; i++) {
        scale[i] *= -1;
        row[i][0] *= -1;
        row[i][1] *= -1;
        row[i][2] *= -1;
      }
    }

    var t = row[0][0] + row[1][1] + row[2][2] + 1;
    var s;
    var quaternion;

    if (t > 1e-4) {
      s = 0.5 / Math.sqrt(t);
      quaternion = [
        (row[2][1] - row[1][2]) * s,
        (row[0][2] - row[2][0]) * s,
        (row[1][0] - row[0][1]) * s,
        0.25 / s
      ];
    } else if (row[0][0] > row[1][1] && row[0][0] > row[2][2]) {
      s = Math.sqrt(1 + row[0][0] - row[1][1] - row[2][2]) * 2.0;
      quaternion = [
        0.25 * s,
        (row[0][1] + row[1][0]) / s,
        (row[0][2] + row[2][0]) / s,
        (row[2][1] - row[1][2]) / s
      ];
    } else if (row[1][1] > row[2][2]) {
      s = Math.sqrt(1.0 + row[1][1] - row[0][0] - row[2][2]) * 2.0;
      quaternion = [
        (row[0][1] + row[1][0]) / s,
        0.25 * s,
        (row[1][2] + row[2][1]) / s,
        (row[0][2] - row[2][0]) / s
      ];
    } else {
      s = Math.sqrt(1.0 + row[2][2] - row[0][0] - row[1][1]) * 2.0;
      quaternion = [
        (row[0][2] + row[2][0]) / s,
        (row[1][2] + row[2][1]) / s,
        0.25 * s,
        (row[1][0] - row[0][1]) / s
      ];
    }

    return {
      translate: translate, scale: scale, skew: skew,
      quaternion: quaternion, perspective: perspective
    };
  }
311724
function(matrix) {
            	var M = matrix.elements || matrix;
            	if (typeof(M[0][0]) == 'undefined') { M = Matrix$b.create(M).elements; }
            	var T = this.dup(), cols = T.elements[0].length;
            	var i = T.elements.length, nj = M[0].length, j;
            	if (i != M.length) { return null; }
            	while (i--) {
            	    j = nj;
            	    while (j--) {
            		T.elements[i][cols + j] = M[i][j];
            	    }
            	}
            	return T;
                }
358801
private void extractSolutionForQ( DMatrix4x4 Q ) {
		DMatrixRMaj nv = new DMatrixRMaj(10,1);
		SingularOps_DDRM.nullVector(svd,true,nv);

		// Convert the solution into a fixed sized matrix because it's easier to read
		encodeQ(Q,nv.data);

		// diagonal elements must be positive because Q = [K*K' .. ; ... ]
		// If they are a mix of positive and negative that's bad and can't be fixed
		// All 3 are checked because technically they could be zero. Not a physically useful solution but...
		if( Q.a11 < 0 || Q.a22 < 0 || Q.a33 < 0 ){
			CommonOps_DDF4.scale(-1,Q);
		}
	}
517864
public Matrix getD () {
		Matrix X = new Matrix(n,n);
		double[][] D = X.getArray();
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < n; j++) {
				D[i][j] = 0.0;
			}
			D[i][i] = d[i];
			if (e[i] > 0) {
				D[i][i+1] = e[i];
			} else if (e[i] < 0) {
				D[i][i-1] = e[i];
			}
		}
		return X;
	}
1361787
@Nonnull
  @ReturnsMutableCopy
  public Matrix getQ ()
  {
    final Matrix aNewMatrix = new Matrix (m_nRows, m_nCols);
    final double [] [] aNewArray = aNewMatrix.internalGetArray ();
    for (int k = m_nCols - 1; k >= 0; k--)
    {
      final double [] aQRk = m_aQR[k];
      for (int nRow = 0; nRow < m_nRows; nRow++)
        aNewArray[nRow][k] = 0.0;
      aNewArray[k][k] = 1.0;
      for (int j = k; j < m_nCols; j++)
      {
        if (aQRk[k] != 0)
        {
          double s = 0.0;
          for (int i = k; i < m_nRows; i++)
            s += m_aQR[i][k] * aNewArray[i][j];
          s = -s / aQRk[k];
          for (int i = k; i < m_nRows; i++)
          {
            aNewArray[i][j] += s * m_aQR[i][k];
          }
        }
      }
    }
    return aNewMatrix;
  }
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;
            }
        }
    }
1302699
public Matrix getMatrix(int[] r, int j0, int j1)
    {
        Matrix X = new Matrix(r.length, j1 - j0 + 1);
        double[][] B = X.getArray();
        try
        {
            for (int i = 0; i < r.length; i++)
            {
                for (int j = j0; j <= j1; j++)
                {
                    B[i][j - j0] = A[r[i]][j];
                }
            }
        }
        catch (ArrayIndexOutOfBoundsException e)
        {
            throw new ArrayIndexOutOfBoundsException("Submatrix indices");
        }
        return X;
    }
45219
public static double[][] getMatrix(final double[][] m1, final int r0, final int r1, final int[] c) {
    assert r0 <= r1 : ERR_INVALID_RANGE;
    assert r1 <= m1.length : ERR_MATRIX_DIMENSIONS;
    final int rowdim = r1 - r0, coldim = c.length;
    final double[][] X = new double[rowdim][coldim];
    for(int i = r0; i < r1; i++) {
      final double[] row = m1[i];
      final double[] Xi = X[i - r0];
      for(int j = 0; j < coldim; j++) {
        Xi[j] = row[c[j]];
      }
    }
    return X;
  }
1068159
@Nonnull
  @ReturnsMutableCopy
  public Matrix getR ()
  {
    final Matrix aNewMatrix = new Matrix (m_nCols, m_nCols);
    final double [] [] aNewArray = aNewMatrix.internalGetArray ();
    for (int nRow = 0; nRow < m_nCols; nRow++)
    {
      final double [] aDstRow = aNewArray[nRow];
      for (int j = 0; j < m_nCols; j++)
      {
        if (nRow < j)
          aDstRow[j] = m_aQR[nRow][j];
        else
          if (nRow == j)
            aDstRow[j] = m_aRdiag[nRow];
          else
            aDstRow[j] = 0.0;
      }
    }
    return aNewMatrix;
  }
779188
@Override
    public DMatrixRMaj getQ(DMatrixRMaj Q , boolean compact ) {
        if( compact ) {
            if( Q == null ) {
                Q = CommonOps_DDRM.identity(numRows,minLength);
            } else {
                if( Q.numRows != numRows || Q.numCols != minLength ) {
                    throw new IllegalArgumentException("Unexpected matrix dimension.");
                } else {
                    CommonOps_DDRM.setIdentity(Q);
                }
            }
        } else {
            if( Q == null ) {
                Q = CommonOps_DDRM.identity(numRows);
            } else {
                if( Q.numRows != numRows || Q.numCols != numRows ) {
                    throw new IllegalArgumentException("Unexpected matrix dimension.");
                } else {
                    CommonOps_DDRM.setIdentity(Q);
                }
            }
        }

        for( int j = minLength-1; j >= 0; j-- ) {
            u[j] = 1;
            for( int i = j+1; i < numRows; i++ ) {
                u[i] = QR.get(i,j);
            }
            QrHelperFunctions_DDRM.rank1UpdateMultR(Q, u, gammas[j], j, j, numRows, v);
        }

        return Q;
    }
43026
public double[][] getR() {
    double[][] R = new double[n][n];
    for(int i = 0; i < n; i++) {
      System.arraycopy(QR[i], i + 1, R[i], i + 1, n - i - 1);
      R[i][i] = Rdiag[i];
    }
    return R;
  }
377988
def _q_matrix(self):
        
        return np.array([
            [self.q[0], -self.q[1], -self.q[2], -self.q[3]],
            [self.q[1],  self.q[0], -self.q[3],  self.q[2]],
            [self.q[2],  self.q[3],  self.q[0], -self.q[1]],
            [self.q[3], -self.q[2],  self.q[1],  self.q[0]]])
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 };
    }
583711
public static QRP factorize(Matrix A) {
        return new QRP(A.numRows(), A.numColumns()).factor(A);
    }
654610
function qr(A) {
    if (A.shape.length !== 2) {
        throw new Error('Input is not matrix');
    }
    var _a = A.shape, m = _a[0], n = _a[1];
    if (m === 0 || n === 0) {
        throw new Error('Empty matrix');
    }
    var copyA = A.clone();
    var minmn = Math.min(m, n);
    copyA.swapOrder();
    var tau = new common_1.NDArray({ shape: [minmn] });
    lapack.geqrf(copyA.data, m, n, tau.data);
    var r = copyA.clone();
    r.swapOrder();
    r = triu(r.get(':', ':' + minmn));
    var q = copyA.get(':' + n);
    lapack.orgqr(q.data, m, n, minmn, tau.data);
    q.swapOrder();
    return [q, r];
}
419013
def _q_to_dcm(self, q):
        
        assert(len(q) == 4)
        arr = super(Quaternion, self)._q_to_dcm(q)
        return self._dcm_array_to_matrix3(arr)