Add a factor times one column to another column

@param matrix
the matrix to modify
@param diag
coordinate on the diagonal
@param fromRow
first row to process
@param col
column to process
@param factor
factor to multiply
760022
public static void addColTimes(DenseDoubleMatrix2D matrix, long diag, long fromRow, long col,
			double factor) {
		long rows = matrix.getRowCount();
		for (long row = fromRow; row < rows; row++) {
			matrix.setDouble(matrix.getDouble(row, col) - factor * matrix.getDouble(row, diag),
					row, col);
		}
	}
521265
public static int[] factor (int n, int fromfactors[]){
    int factors[] = new int[64]; // Cant be more than 64 factors.
    int nf = 0;
    int ntest = n;
    int factor;

    if (n <= 0)			// Error case
      throw new Error("Number ("+n+") must be positive integer");

    /* deal with the preferred factors first */
    for(int i = 0; i < fromfactors.length && ntest != 1; i++){
      factor = fromfactors[i];
      while ((ntest % factor) == 0) {
	ntest /= factor;
	factors[nf++] = factor; }}

    /* deal with any other even prime factors (there is only one) */
    factor = 2;
    while ((ntest % factor) == 0 && (ntest != 1)) {
      ntest /= factor;
      factors[nf++] = factor; }

    /* deal with any other odd prime factors */
    factor = 3;
    while (ntest != 1) {
      while ((ntest % factor) != 0) {
	factor += 2; }
      ntest /= factor;
      factors[nf++] = factor; }

    /* check that the factorization is correct */
    int product = 1;
    for (int i = 0; i < nf; i++) {
      product *= factors[i]; }
    if (product != n)
      throw new Error("factorization failed for "+n);

    /* Now, make an array of the right length containing the factors... */
    int f[] = new int[nf];
    System.arraycopy(factors,0,f,0,nf);    
    return f; }
741450
public Matrix4x3f set(Vector3fc col0, Vector3fc col1, Vector3fc col2, Vector3fc col3) {
        if (col0 instanceof Vector3f &&
            col1 instanceof Vector3f &&
            col2 instanceof Vector3f &&
            col3 instanceof Vector3f) {
            MemUtil.INSTANCE.set(this, (Vector3f) col0, (Vector3f) col1, (Vector3f) col2, (Vector3f) col3);
        } else {
            setVector3fc(col0, col1, col2, col3);
        }
        this.properties = 0;
        return this;
    }
543395
public void addMatrix(double[][] m) {
    assert (m.length == dim);
    assert (m[0].length == dim);

    // reset inverse transformation - needs recomputation.
    inv = null;

    // extend the matrix with an extra row and column
    double[][] ht = new double[dim + 1][dim + 1];
    for(int i = 0; i < dim; i++) {
      for(int j = 0; j < dim; j++) {
        ht[i][j] = m[i][j];
      }
    }
    // the other cells default to identity matrix
    ht[dim][dim] = 1.0;
    // Multiply from left.
    trans = times(ht, trans);
  }
939324
public Matrix multiply(Matrix B)
    {
        Matrix C = new DenseMatrix(this.rows(), B.cols());
        multiply(B, C);
        return C;
    }
74443
public static Matrix createMatrix(int rowCount, int colCount) {
		if (rowCount<=0 || colCount<=0) throw new IllegalArgumentException();
		if (colCount==1)
			return createVector(rowCount);
		else if (rowCount==2 && colCount==2)
			return createMatrix2();
		else if (rowCount==3 && colCount==3)
			return createMatrix3();
		else
			return new MatrixImpl(rowCount, colCount);
	}
93982
public static void diag( DMatrix2x2 input , DMatrix2 out ) {
        out.a1 = input.a11;
        out.a2 = input.a22;
    }
73307
public static void addCol(Matrix A, int j, double c)
    {
        addCol(A, j, 0, A.rows(), c);
    }
56350
void add(DoubleVector dest, Vector src, int factor) {
        if (src instanceof SparseVector) {
            int[] nonZeros = ((SparseVector) src).getNonZeroIndices();
            for (int i : nonZeros)
                dest.add(i, factor * src.getValue(i).doubleValue());
        } else {
            for (int i = 0; i < src.length(); ++i)
                dest.add(i, factor * src.getValue(i).doubleValue());
        }
    }
944140
public void addConstantFactor(String factorName, Factor factor) {
    if (variables.containsAll(factor.getVars())) {
      int factorNum = factors.size();
      factors.add(factor);
      factorNames.add(factorName);

      VariableNumMap factorVars = factor.getVars();
      for (Integer varNum : factorVars.getVariableNums()) {
        variableFactorMap.put(varNum, factorNum);
        factorVariableMap.put(factorNum, varNum);
      }
    } else {
      Preconditions.checkState(!isRoot);
      crossingFactors.add(factor);
    }
  }
123512
public static Matrix diag(Vec a)
    {
        DenseMatrix A = new DenseMatrix(a.length(), a.length());
        for(Iterator<IndexValue> iter = a.getNonZeroIterator(); iter.hasNext();)
        {
            IndexValue iv = iter.next();
            A.set(iv.getIndex(), iv.getIndex(), iv.getValue());
        }
            
        return A;
    }
1014128
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 = (SparseDoubleMatrix2D) 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 = Utils.getDoubleMachineEpsilon();
		//elements of L are stored in a single array for i<diagonalLength, 
		//and with arrays of dimension i+1 for i>=diagonalLength
		this.LData = new double[dim - diagonalLength + 1][];
		double[] LData0 = new double[diagonalLength];
		LData[0] = LData0;
		
		if(dim - diagonalLength == 1 && Double.compare(Q.getQuick(dim-1, dim-1), 1.)==0){
			//TODO: the second condition can always be true if the matrix is normalized on the last element
			//closed form decomposition
			//@see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 671"
			for (int i = 0; i < dim-1; i++) {
				double d = Q.getQuick(i, i);
				if (d < Utils.getDoubleMachineEpsilon()) {
					throw new Exception("not positive definite matrix");
				}
				LData0[i] = Math.sqrt(d);
			}
			//now the last row of L
			double[] LData1 = new double[dim];
			for(int j=0; j<dim - 1; j++){
				LData1[j] = Q.getQuick(dim-1, j) / LData0[j];
			}
			double d = 0;
			for(int k=0; k<diagonalLength; k++){
				d += Math.pow(Q.getQuick(dim-1, k), 2) / Q.getQuick(k, k);
			}
			LData1[dim-1] = Math.sqrt(1 - d);
			LData[1] = LData1;
		}else{
			for (int i = 0; i < dim; i++) {
			    if (i < diagonalLength) {
					double d = Q.getQuick(i, i);
					if(!(d > threshold)){
						throw new Exception("not positive definite matrix");
					}
					LData0[i] = Math.sqrt(d);
				} else {

					LData[i - diagonalLength + 1] = new double[i + 1];
					double[] LDataI = LData[i - diagonalLength + 1];

					// j < i
					for (int j = 0; j < diagonalLength; j++) {
						// here LData is not null only in its diagonal
						LDataI[j] = 1.0 / LData[0][j] * (Q.getQuick(i, j));
					}
					for (int j = diagonalLength; j < i; j++) {
						double[] LDataJ = LData[j - diagonalLength + 1];
						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 < Utils.getDoubleMachineEpsilon()) {
						throw new Exception("not positive definite matrix");
					}
					LDataI[i] = Math.sqrt(d);
				}
			}
		}
	}
1130849
public static double meanDoublematrixColumn( double[][] matrix, int column ) {

        double mean;

        mean = 0;

        int length = matrix.length;
        for( int i = 0; i < length; i++ ) {
            mean += matrix[i][column];
        }

        return mean / length;

    }
1181769
public static void swapCols(double[][] matrix, int col1, int col2) {
		double temp = 0;
		int rows = matrix.length;
		double[] r = null;
		for (int row = 0; row < rows; row++) {
			r = matrix[row];
			temp = r[col1];
			r[col1] = r[col2];
			r[col2] = temp;
		}
	}
401939
public static DMatrixRMaj concatRowsMulti(DMatrixRMaj ...m ) {
        int rows = 0;
        int cols = 0;

        for (int i = 0; i < m.length; i++) {
            rows += m[i].numRows;
            cols = Math.max(cols,m[i].numCols);
        }
        DMatrixRMaj R = new DMatrixRMaj(rows,cols);

        int row = 0;
        for (int i = 0; i < m.length; i++) {
            insert(m[i],R,row,0);
            row += m[i].numRows;
        }

        return R;
    }
227425
protected static void initializeDistanceMatrix(MatrixParadigm mat, DistanceQuery<?> dq, Linkage linkage) {
    final DBIDArrayIter ix = mat.ix, iy = mat.iy;
    final double[] matrix = mat.matrix;
    final boolean issquare = dq.getDistanceFunction().isSquared();
    FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Distance matrix computation", matrix.length, LOG) : null;
    int pos = 0;
    for(ix.seek(0); ix.valid(); ix.advance()) {
      final int x = ix.getOffset();
      assert (pos == MatrixParadigm.triangleSize(x));
      for(iy.seek(0); iy.getOffset() < x; iy.advance()) {
        matrix[pos++] = linkage.initial(dq.distance(ix, iy), issquare);
      }
      if(prog != null) {
        prog.setProcessed(pos, LOG);
      }
    }
    // Avoid logging errors in case scratch space was too large:
    if(prog != null) {
      prog.setProcessed(matrix.length, LOG);
    }
    LOG.ensureCompleted(prog);
  }
163446
public void apply(Matrix H, int column, int i1, int i2) {
        double temp = c * H.get(i1, column) + s * H.get(i2, column);
        H.set(i2, column, -s * H.get(i1, column) + c * H.get(i2, column));
        H.set(i1, column, temp);
    }