Sparse QR factorization of an m-by-n matrix A, A= Q*R

@param A
column-compressed matrix
@param S
symbolic QR analysis
@return numeric QR factorization, null on error
435362
public static DZcsn cs_qr(DZcs A, DZcss S)
	{
		DZcsa Rx = new DZcsa(), Vx = new DZcsa(), Ax = new DZcsa(), x ;
		double Beta[] ;
		int i, k, p, n, vnz, p1, top, m2, len, col, rnz, s[], leftmost[], Ap[], Ai[],
			parent[], Rp[], Ri[], Vp[], Vi[], w[], pinv[], q[] ;
		DZcs R, V ;
		DZcsn N ;
		if (!CS_CSC (A) || S == null) return (null) ;
		n = A.n ; Ap = A.p ; Ai = A.i ; Ax.x = A.x ;
		q = S.q ; parent = S.parent ; pinv = S.pinv ; m2 = S.m2 ;
		vnz = S.lnz ; rnz = S.unz ; leftmost = S.leftmost ;
		w = new int [m2] ; 			/* get int workspace */
		x = new DZcsa (m2) ;			/* get double workspace */
		N = new DZcsn () ;			/* allocate result */
		s = new int [n] ; 			/* get int workspace, s is size n */
		//for (k = 0 ; k < m2 ; k++) x.set(k, cs_czero()) ; 	/* clear workspace x */
		N.L = V = cs_spalloc(m2, n, vnz, true, false) ;  	/* allocate result V */
		N.U = R = cs_spalloc(m2, n, rnz, true, false) ;  	/* allocate result R */
		N.B = Beta = new double [n] ;				/* allocate result Beta */
		if (R == null || V == null || Beta == null) return (cs_ndone (N, null, w, x, false)) ;
		Rp = R.p ; Ri = R.i ; Rx.x = R.x ;
		Vp = V.p ; Vi = V.i ; Vx.x = V.x ;
		for (i = 0 ; i < m2 ; i++) w [i] = -1 ;	/* clear w, to mark nodes */
		rnz = 0 ; vnz = 0 ;
		for (k = 0 ; k < n ; k++)		/* compute V and R */
		{
			Rp [k] = rnz ;			/* R(:,k) starts here */
			Vp [k] = p1 = vnz ;		/* V(:,k) starts here */
			w [k] = k ;			/* add V(k,k) to pattern of V */
			Vi [vnz++] = k ;
			top = n ;
			col = q != null ? q [k] : k ;
			for (p = Ap [col] ; p < Ap [col+1] ; p++)  /* find R(:,k) pattern */
			{
				i = leftmost [Ai [p]] ;	/* i = min(find(A(i,q))) */
				for (len = 0 ; w [i] != k ; i = parent [i])  /* traverse up to k */
				{
					s [len++] = i ;
					w [i] = k ;
				}
				while (len > 0)
					s [--top] = s [--len] ;  /* push path on stack */
				i = pinv [Ai [p]] ;	/* i = permuted row of A(:,col) */
				x.set(i, Ax.real(p), Ax.imag(p)) ;	/* x (i) = A(:,col) */
				if (i > k && w [i] < k)	/* pattern of V(:,k) = x (k+1:m) */
				{
					Vi [vnz++] = i ;  /* add i to pattern of V(:,k) */
					w [i] = k ;
				}
			}
			for (p = top ; p < n ; p++)	/* for each i in pattern of R(:,k) */
			{
				i = s [p] ;	/* R(i,k) is nonzero */
				cs_happly (V, i, Beta [i], x) ;  /* apply (V(i),Beta(i)) to x */
				Ri [rnz] = i ;		/* R(i,k) = x(i) */
				Rx.set(rnz++, x.real(i), x.imag(i)) ;
				x.set(i, 0.0, 0.0) ;
				if (parent [i] == k)
					vnz = cs_scatter (V, i, cs_czero(), w, null, k, V, vnz) ;
			}
			for (p = p1 ; p < vnz ; p++)	/* gather V(:,k) = x */
			{
				Vx.set(p, x.real(Vi [p]), x.imag(Vi [p])) ;
				x.set(Vi [p], 0.0, 0.0) ;
			}
			Ri [rnz] = k ;			/* R(k,k) = norm (x) */
			double[] beta = new double[] {Beta [k]} ;
			Rx.set(rnz++, cs_house(Vx, p1, beta, vnz - p1)) ;  /* [v,beta]=house(x) */
			Beta [k] = beta [0] ;
		}
		Rp [n] = rnz ; 				/* finalize R */
		Vp [n] = vnz ;				/* finalize V */
		return (N) ;
	}
1247897
public static void symmLowerToFull( DMatrixRMaj A )
    {
        if( A.numRows != A.numCols )
            throw new MatrixDimensionException("Must be a square matrix");

        final int cols = A.numCols;

        for (int row = 0; row < A.numRows; row++) {
            for (int col = row+1; col < cols; col++) {
                A.data[row*cols+col] = A.data[col*cols+row];
            }
        }
    }
876006
public void decompose( DMatrixRMaj A ) {
        init(A);

        for( int k = 1; k < N; k++ ) {
            similarTransform(k);
//            System.out.println("k=="+k);
//            QT.print();
        }
    }
259553
protected void convertToColumnMajor(DMatrixRMaj A) {
        for( int x = 0; x < numCols; x++ ) {
            double colQ[] = dataQR[x];
            for( int y = 0; y < numRows; y++ ) {
                colQ[y] = A.data[y*A.numCols+x];
            }
        }
    }
311030
public static <T extends SparseDoubleVector> SparseMatrix asSparseMatrix(
            List<T> vectors, int columns) {
        return new SparseListMatrix<T>(vectors, columns);
    }
409352
public DMatrixRMaj getQ() {
        DMatrixRMaj Q = CommonOps_DDRM.identity(QR.numRows);
        DMatrixRMaj Q_k = new DMatrixRMaj(QR.numRows,QR.numRows);
        DMatrixRMaj u = new DMatrixRMaj(QR.numRows,1);

        DMatrixRMaj temp = new DMatrixRMaj(1,1);

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

        // compute Q by first extracting the householder vectors from the columns of QR and then applying it to Q
        for( int j = N-1; j>= 0; j-- ) {
            CommonOps_DDRM.extract(QR,j, QR.numRows,j,j+1,u,j,0);
            u.set(j,1.0);

            // A = (I - &gamma;*u*u<sup>T</sup>)*A<br>
            CommonOps_DDRM.setIdentity(Q_k);
            CommonOps_DDRM.multAddTransB(-gammas[j],u,u,Q_k);
            CommonOps_DDRM.mult(Q_k,Q,temp);
            Q.set(temp);
        }

        return Q;
    }
291077
public static EVD factorize(Matrix A) throws NotConvergedException {
        return new EVD(A.numRows()).factor(new DenseMatrix(A));
    }
583711
public static QRP factorize(Matrix A) {
        return new QRP(A.numRows(), A.numColumns()).factor(A);
    }
417852
public static DMatrixRMaj nullspaceQR( DMatrixRMaj A , int totalSingular ) {
        SolveNullSpaceQR_DDRM solver = new SolveNullSpaceQR_DDRM();

        DMatrixRMaj nullspace = new DMatrixRMaj(1,1);

        if( !solver.process(A,totalSingular,nullspace))
            throw new RuntimeException("Solver failed. try SVD based method instead?");

        return nullspace;
    }
625279
public Matrix timesEquals(double s)
    {
        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < n; j++)
            {
                A[i][j] = s * A[i][j];
            }
        }
        return this;
    }
1293364
public static int cusolverSpScsrlsqvqrHost(
        cusolverSpHandle handle, 
        int m, 
        int n, 
        int nnz, 
        cusparseMatDescr descrA, 
        Pointer csrValA, 
        Pointer csrRowPtrA, 
        Pointer csrColIndA, 
        Pointer b, 
        float tol, 
        Pointer rankA, 
        Pointer x, 
        Pointer p, 
        Pointer min_norm)
    {
        return checkResult(cusolverSpScsrlsqvqrHostNative(handle, m, n, nnz, descrA, csrValA, csrRowPtrA, csrColIndA, b, tol, rankA, x, p, min_norm));
    }
340962
protected void convertToColumnMajor(ZMatrixRMaj A) {
        for( int x = 0; x < numCols; x++ ) {
            double colQ[] = dataQR[x];
            int indexCol = 0;
            for( int y = 0; y < numRows; y++ ) {
                int index = (y*numCols+x)*2;
                colQ[indexCol++] = A.data[index];
                colQ[indexCol++] = A.data[index+1];
            }
        }
    }
108386
public static BandLU factorize(BandMatrix A) {
        return new BandLU(A.numRows(), A.kl, A.ku).factor(A, false);
    }
698269
public SVD factor(DenseMatrix A) throws NotConvergedException {
        if (A.numRows() != m)
            throw new IllegalArgumentException("A.numRows() != m");
        else if (A.numColumns() != n)
            throw new IllegalArgumentException("A.numColumns() != n");

        intW info = new intW(0);
        LAPACK.getInstance().dgesdd(job.netlib(), m, n, A.getData(),
                Matrices.ld(m), S, vectors ? U.getData() : new double[0],
                Matrices.ld(m), vectors ? Vt.getData() : new double[0],
                Matrices.ld(n), work, work.length, iwork, info);

        if (info.val > 0)
            throw new NotConvergedException(
                    NotConvergedException.Reason.Iterations);
        else if (info.val < 0)
            throw new IllegalArgumentException();

        return this;
    }
21674
public static Scsn cs_qr(Scs A, Scss S) {
        float Rx[], Vx[], Ax[], x[], Beta[];
        int i, k, p, n, vnz, p1, top, m2, len, col, rnz, s[], leftmost[], Ap[], Ai[], parent[], Rp[], Ri[], Vp[], Vi[], w[], pinv[], q[];
        Scs R, V;
        Scsn N;
        if (!Scs_util.CS_CSC(A) || S == null)
            return (null);
        n = A.n;
        Ap = A.p;
        Ai = A.i;
        Ax = A.x;
        q = S.q;
        parent = S.parent;
        pinv = S.pinv;
        m2 = S.m2;
        vnz = S.lnz;
        rnz = S.unz;
        leftmost = S.leftmost;
        w = new int[m2 + n]; /* get int workspace */
        x = new float[m2]; /* get float workspace */
        N = new Scsn(); /* allocate result */
        s = w;
        int s_offset = m2; /* s is size n */
        for (k = 0; k < m2; k++)
            x[k] = 0; /* clear workspace x */
        N.L = V = Scs_util.cs_spalloc(m2, n, vnz, true, false); /* allocate result V */
        N.U = R = Scs_util.cs_spalloc(m2, n, rnz, true, false); /* allocate result R */
        N.B = Beta = new float[n]; /* allocate result Beta */
        Rp = R.p;
        Ri = R.i;
        Rx = R.x;
        Vp = V.p;
        Vi = V.i;
        Vx = V.x;
        for (i = 0; i < m2; i++)
            w[i] = -1; /* clear w, to mark nodes */
        rnz = 0;
        vnz = 0;
        for (k = 0; k < n; k++) /* compute V and R */
        {
            Rp[k] = rnz; /* R(:,k) starts here */
            Vp[k] = p1 = vnz; /* V(:,k) starts here */
            w[k] = k; /* add V(k,k) to pattern of V */
            Vi[vnz++] = k;
            top = n;
            col = q != null ? q[k] : k;
            for (p = Ap[col]; p < Ap[col + 1]; p++) /* find R(:,k) pattern */
            {
                i = leftmost[Ai[p]]; /* i = min(find(A(i,q))) */
                for (len = 0; w[i] != k; i = parent[i]) /* traverse up to k */
                {
                    s[s_offset + (len++)] = i;
                    w[i] = k;
                }
                while (len > 0)
                    s[s_offset + (--top)] = s[s_offset + (--len)]; /* push path on stack */
                i = pinv[Ai[p]]; /* i = permuted row of A(:,col) */
                x[i] = Ax[p]; /* x (i) = A(:,col) */
                if (i > k && w[i] < k) /* pattern of V(:,k) = x (k+1:m) */
                {
                    Vi[vnz++] = i; /* add i to pattern of V(:,k) */
                    w[i] = k;
                }
            }
            for (p = top; p < n; p++) /* for each i in pattern of R(:,k) */
            {
                i = s[s_offset + p]; /* R(i,k) is nonzero */
                Scs_happly.cs_happly(V, i, Beta[i], x); /* apply (V(i),Beta(i)) to x */
                Ri[rnz] = i; /* R(i,k) = x(i) */
                Rx[rnz++] = x[i];
                x[i] = 0;
                if (parent[i] == k)
                    vnz = Scs_scatter.cs_scatter(V, i, 0, w, null, k, V, vnz);
            }
            for (p = p1; p < vnz; p++) /* gather V(:,k) = x */
            {
                Vx[p] = x[Vi[p]];
                x[Vi[p]] = 0;
            }
            Ri[rnz] = k; /* R(k,k) = norm (x) */
            float[] beta = new float[1];
            beta[0] = Beta[k];
            Rx[rnz++] = Scs_house.cs_house(Vx, p1, beta, vnz - p1); /* [v,beta]=house(x) */
            Beta[k] = beta[0];
        }
        Rp[n] = rnz; /* finalize R */
        Vp[n] = vnz; /* finalize V */
        return N;
    }
241819
public static void symmLowerToFull( DMatrixSparseCSC A , DMatrixSparseCSC B , @Nullable IGrowArray gw )
    {
        if( A.numCols != A.numRows )
            throw new IllegalArgumentException("Must be a lower triangular square matrix");

        int N = A.numCols;
        int w[] = adjust(gw,N,N);
        B.reshape(N,N,A.nz_length*2);
        B.indicesSorted = false;

        //=== determine the row counts of the full matrix
        for (int col = 0; col < N; col++) {
            int idx0 = A.col_idx[col];
            int idx1 = A.col_idx[col+1];

            // We know the length of the lower part of this column already
            w[col] += idx1 - idx0;

            // add elements to the top of the other columns along row with index 'col'
            for (int i = idx0; i < idx1; i++) {
                int row = A.nz_rows[i];
                if( row > col ) {
                    w[row]++;
                }
            }
        }

        // Update the structure of B
        B.histogramToStructure(w);

        // Zero W again. It's being used to keep track of how many elements have been added to a column already
        Arrays.fill(w,0,N,0);
        // Fill in matrix
        for (int col = 0; col < N; col++) {

            int idx0 = A.col_idx[col];
            int idx1 = A.col_idx[col+1];

            int lengthA = idx1 - idx0;
            int lengthB = B.col_idx[col+1] - B.col_idx[col];

            // Copy the non-zero values from A into B along the columns while taking in account the upper
            // elements already copied
            System.arraycopy(A.nz_values,idx0,B.nz_values,B.col_idx[col]+lengthB-lengthA,lengthA);
            System.arraycopy(A.nz_rows,idx0,B.nz_rows,B.col_idx[col]+lengthB-lengthA,lengthA);

            // Copy this column into the upper portion of B
            for (int i = idx0; i < idx1; i++) {
                int row = A.nz_rows[i];
                if( row > col ) {
                    int indexB = B.col_idx[row] + w[row]++;
                    B.nz_rows[indexB] = col;
                    B.nz_values[indexB] = A.nz_values[i];
                }
            }
        }

    }
896318
public void decompose( DMatrixRMaj A ) {

        Equation eq = new Equation();

        this.QR = A.copy();

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

        gammas = new double[ A.numCols ];

        for( int i = 0; i < N; i++ ) {
            // update temporary variables
            eq.alias(QR.numRows-i,"Ni",QR,"QR",i,"i");

            // Place the column that should be zeroed into v
            eq.process("v=QR(i:,i)");
            eq.process("maxV=max(abs(v))");

            // Note that v is lazily created above.  Need direct access to it, which is done below.
            DMatrixRMaj v = eq.lookupDDRM("v");

            double maxV = eq.lookupDouble("maxV");
            if( maxV > 0 && v.getNumElements() > 1 ) {
                // normalize to reduce overflow issues
                eq.process("v=v/maxV");

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

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

                eq.alias(tau,"tau");
                eq.process("u_0 = v(0,0)+tau");
                eq.process("gamma = u_0/tau");
                eq.process("v=v/u_0");
                eq.process("v(0,0)=1");
                eq.process("QR(i:,i:) = (eye(Ni) - gamma*v*v')*QR(i:,i:)");
                eq.process("QR(i:,i) = v");
                eq.process("QR(i,i) = -1*tau*maxV");

                // save gamma for recomputing Q later on
                gammas[i] = eq.lookupDouble("gamma");
            }
        }
    }