Convenience method for computing the full eigenvalue decomposition of the
given matrix

@param A
Matrix to factorize. Upper part extracted, and the matrix is
not modified
@return Newly allocated decomposition
@throws NotConvergedException
483340
public static SymmPackEVD factorize(Matrix A) throws NotConvergedException {
        return new SymmPackEVD(A.numRows(), true)
                .factor(new UpperSymmPackMatrix(A));
    }
27978
def factorize(self, A):
        
        
        A = csc_matrix(A)
        if self.prop == self.SYMMETRIC:
            A = (A + A.T) - triu(A)

        self.lu = self.umfpack.splu(A)
1320460
public static double elementDiagMaxMagnitude2(ZMatrixRMaj a) {
        final int size = Math.min(a.numRows,a.numCols);

        int rowStride = a.getRowStride();
        double max = 0;
        for( int i = 0; i < size; i++ ) {
            int index = i*rowStride + i*2;

            double real = a.data[index];
            double imaginary = a.data[index+1];

            double m = real*real + imaginary*imaginary;

            if( m > max ) {
                max = m;
            }
        }

        return max;
    }
1287543
public static FloatMatrix[] symmetricEigenvectors(FloatMatrix A) {
        A.assertSquare();
        FloatMatrix eigenvalues = new FloatMatrix(A.rows);
        FloatMatrix eigenvectors = A.dup();
        int isuppz[] = new int[2 * A.rows];
        SimpleBlas.syevr('V', 'A', 'U', A.dup(), 0, 0, 0, 0, 0, eigenvalues, eigenvectors, isuppz);
        return new FloatMatrix[]{eigenvectors, FloatMatrix.diag(eigenvalues)};
    }
876029
@Override
	public void computeHessian(DMatrixSparseCSC jacLeft , DMatrixSparseCSC jacRight) {
		A.reshape(jacLeft.numCols,jacLeft.numCols,1);
		B.reshape(jacLeft.numCols,jacRight.numCols,1);
		D.reshape(jacRight.numCols,jacRight.numCols,1);

		// take advantage of the inner product's symmetry when possible to reduce
		// the number of calculations
		CommonOps_DSCC.innerProductLower(jacLeft,tmp0,gw,gx);
		CommonOps_DSCC.symmLowerToFull(tmp0,A,gw);
		CommonOps_DSCC.multTransA(jacLeft,jacRight,B,gw,gx);
		CommonOps_DSCC.innerProductLower(jacRight,tmp0,gw,gx);
		CommonOps_DSCC.symmLowerToFull(tmp0,D,gw);
	}
258993
@Override
    public boolean setA(DMatrixRMaj A) {
        if( A.numRows > maxRows || A.numCols > maxCols ) {
            setMaxSize(A.numRows,A.numCols);
        }

        _setA(A);
        if( !decomposer.decompose(A) )
            return false;

        Q.reshape(numRows,numRows, false);
        R.reshape(numRows,numCols, false);
        decomposer.getQ(Q,false);
        decomposer.getR(R,false);

        return true;
    }
1302773
@Override
    public boolean decompose( ZMatrixRMaj A )
    {
        if( A.numRows != A.numCols )
            throw new IllegalArgumentException("A must be square.");
        if( A.numRows <= 0 )
            return false;

        QH = A;

        N = A.numCols;

        if( b.length < N*2 ) {
            b = new double[ N*2 ];
            gammas = new double[ N ];
            u = new double[ N*2 ];
        }
        return _decompose();
    }
1000736
public static void divide( DMatrix4 a , double alpha ) {
        a.a1 /= alpha;
        a.a2 /= alpha;
        a.a3 /= alpha;
        a.a4 /= alpha;
    }
1285072
@Override
	public DoubleMatrix1D[] getMatrixScalingFactors(DoubleMatrix2D A){
		DoubleFactory1D F1 = DoubleFactory1D.dense;
		Algebra ALG = Algebra.DEFAULT;
		int r = A.rows();
		int c = A.columns();
		DoubleMatrix1D D1 = F1.make(r, 1);
		DoubleMatrix1D D2 = F1.make(c, 1);
		DoubleMatrix2D AK = A.copy();
		DoubleMatrix1D DR = F1.make(r, 1);
		DoubleMatrix1D DC = F1.make(c, 1);
		DoubleMatrix1D DRInv = F1.make(r);
		DoubleMatrix1D DCInv = F1.make(c);
		log.debug("eps  : " + eps);
		int maxIteration = 50;
		for(int k=0; k<=maxIteration; k++){
			double normR = -Double.MAX_VALUE;
			double normC = -Double.MAX_VALUE;
			for(int i=0; i<r; i++){
				double dri = ALG.normInfinity(AK.viewRow(i));
				DR.setQuick(i, Math.sqrt(dri));
				DRInv.setQuick(i, 1./Math.sqrt(dri));
				normR = Math.max(normR, Math.abs(1-dri));
			}
			for(int j=0; j<c; j++){
				double dci = ALG.normInfinity(AK.viewColumn(j));
				DC.setQuick(j, Math.sqrt(dci));
				DCInv.setQuick(j, 1./Math.sqrt(dci));
				normC = Math.max(normC, Math.abs(1-dci));
			}
			
			log.debug("normR: " + normR);
			log.debug("normC: " + normC);
			if(normR < eps && normC < eps){
				break;
			}
			
			//D1 = ALG.mult(D1, DRInv);
			for(int i=0; i<r; i++){
				double prevD1I = D1.getQuick(i);
				double newD1I = prevD1I * DRInv.getQuick(i);
				D1.setQuick(i, newD1I);
			}
			//D2 = ALG.mult(D2, DCInv);
			for(int j=0; j<c; j++){
				double prevD2J = D2.getQuick(j);
				double newD2J = prevD2J * DCInv.getQuick(j);
				D2.setQuick(j, newD2J);
			}
			//log.debug("D1: " + ArrayUtils.toString(D1.toArray()));
			//log.debug("D2: " + ArrayUtils.toString(D2.toArray()));
			
			if(k==maxIteration){
				log.warn("max iteration reached");
			}
			
			//AK = ALG.mult(DRInv, ALG.mult(AK, DCInv));
			AK = ColtUtils.diagonalMatrixMult(DRInv, AK, DCInv);
		}
		
		return new DoubleMatrix1D[]{D1, D2};
	}
609005
private ClusterResult fullCluster(Matrix matrix,
                                      final int depth) {
        verbose("Clustering at depth " + depth);

        // If the matrix has only one element or the depth is equal to the
        // maximum number of desired clusters then all items are in a single
        // cluster.
        if (matrix.rows() <= 1) 
            return new ClusterResult(new int[matrix.rows()], 1);

        // Get a fresh new eigen cutter and compute the specral cut of the
        // matrix.
        EigenCut eigenCutter = cutterGenerator.generate();
        eigenCutter.computeCut(matrix);

        final Matrix leftMatrix = eigenCutter.getLeftCut();
        final Matrix rightMatrix = eigenCutter.getRightCut();

        verbose(String.format("Splitting into two matricies %d-%d",
                              leftMatrix.rows(), rightMatrix.rows()));

        // If the compute decided that the matrix should not be split, short
        // circuit any attempts to further cut the matrix.
        if (leftMatrix.rows() == matrix.rows() ||
            rightMatrix.rows() == matrix.rows())
            return new ClusterResult(new int[matrix.rows()], 1);


        // Do clustering on the left and right branches.  We can do this in
        // parallel because all of the data members in this class are thread
        // safe, since each call to fullCluster uses a new instance of a
        // EigenCutter which has all of the state for a particular partition.
        final ClusterResult[] results = new ClusterResult[2];
        results[0] = fullCluster(leftMatrix, depth+1);
        results[1] = fullCluster(rightMatrix, depth+1);

        ClusterResult leftResult = results[0];
        ClusterResult rightResult = results[1];

        verbose("Merging at depth " + depth);

        // Compute the relaxed correlation objective function over the split
        // partitions found so far.
        double splitObjective = eigenCutter.getSplitObjective(
                alpha, beta,
                leftResult.numClusters, leftResult.assignments,
                rightResult.numClusters, rightResult.assignments);

        // Compute the objective when we merge the two branches together.
        double mergedObjective = eigenCutter.getMergedObjective(alpha, beta);

        // If the merged objective value is less than the split version, combine
        // all clusters into one.
        int[] assignments = new int[matrix.rows()];
        int numClusters = 1;
        if (mergedObjective < splitObjective) {
            verbose("Selecting to combine sub trees at depth " + depth);
            Arrays.fill(assignments, 0);
        } else  {
            verbose("Selecting to maintain sub trees at depth " + depth);

            // Copy over the left assignments and the right assignments, where
            // the cluster id's of the right assignments are incremented to
            // avoid duplicate cluster ids.
            numClusters = leftResult.numClusters + rightResult.numClusters;

            int[] leftReordering = eigenCutter.getLeftReordering();
            int[] rightReordering = eigenCutter.getRightReordering();

            for (int index = 0; index < leftReordering.length; ++index)
                assignments[leftReordering[index]] = 
                    leftResult.assignments[index];
            for (int index = 0; index < rightReordering.length; ++index)
                assignments[rightReordering[index]] =
                    rightResult.assignments[index] + leftResult.numClusters;
        }
        return new ClusterResult(assignments, numClusters);
    }
519789
public SymmBandEVD factor(UpperSymmBandMatrix A)
            throws NotConvergedException {
        if (uplo != UpLo.Upper)
            throw new IllegalArgumentException(
                    "Eigenvalue computer configured for upper-symmetrical matrices");

        return factor(A, A.getData(), A.numSuperDiagonals());
    }
443386
public static SymmDenseEVD factorize(Matrix A) throws NotConvergedException {
        return new SymmDenseEVD(A.numRows(), true)
                .factor(new UpperSymmDenseMatrix(A));
    }
987540
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.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");
				}
			}
		}

		switch (this.mode) {

		case DIAGONAL_PIVOLTING:
			pldltpt();
			break;

		case DIAGONAL_PIVOLTING_WITH_PARTIAL_PIVOTING:
			pldltptBK();
			break;

		default:
			throw new IllegalArgumentException("unknown mode " + mode);
		}
		
	}
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");
            }
        }
    }
957920
public static PackCholesky factorize(Matrix A) {
        return new PackCholesky(A.numRows(), true)
                .factor(new UpperSPDPackMatrix(A));
    }
1144392
@Override
    public DMatrixRMaj getA() {
        if( A.data.length < numRows*numCols ) {
            A = new DMatrixRMaj(numRows,numCols);
        }
        A.reshape(numRows,numCols, false);
        CommonOps_DDRM.mult(Q,R,A);

        return A;
    }
340054
public boolean process( DMatrixSparseCSC A ) {
        init(A);

        TriangularSolver_DSCC.eliminationTree(A,true,parent,gwork);

        countNonZeroInR(parent);
        countNonZeroInV(parent);

        // if more columns than rows it's possible that Q*R != A. That's because a householder
        // would need to be created that's outside the  m by m Q matrix. In reality it has
        // a partial solution. Column pivot are needed.
        if( m < n ) {
            for (int row = 0; row <m; row++) {
                if( gwork.data[head+row] < 0 ) {
                    return false;
                }
            }
        }
        return true;
    }