Create a permutation matrix from a LAPACK-style 'ipiv' vector.

@param ipiv row i was interchanged with row ipiv[i]
353571
public static DoubleMatrix permutationDoubleMatrixFromPivotIndices(int size, int[] ipiv) {
        int n = ipiv.length;
        //System.out.printf("size = %d n = %d\n", size, n);
        int indices[] = new int[size];
        for (int i = 0; i < size; i++)
            indices[i] = i;

        //for (int i = 0; i < n; i++)
        //    System.out.printf("ipiv[%d] = %d\n", i, ipiv[i]);

        for (int i = 0; i < n; i++) {
            int j = ipiv[i] - 1;
            int t = indices[i];
            indices[i] = indices[j];
            indices[j] = t;
        }
        DoubleMatrix result = new DoubleMatrix(size, size);
        for (int i = 0; i < size; i++)
            result.put(indices[i], i, 1.0);
        return result;
    }
311851
def get_matrix(self, indices):
        """"""
        new = numpy.empty(self.samples1.shape)
        for idx in range(len(indices)):
            if indices[idx]:
                new[idx] = self.samples1[idx]
            else:
                new[idx] = self.samples2[idx]

        if self.poly:
            new = self.poly(*new)
        return new
293772
public static void convert(DMatrixD1 input , ZMatrixD1 output ) {
        if( input.numCols != output.numCols || input.numRows != output.numRows ) {
            throw new IllegalArgumentException("The matrices are not all the same dimension.");
        }

        Arrays.fill(output.data, 0, output.getDataLength(), 0);

        final int length = output.getDataLength();

        for( int i = 0; i < length; i += 2 ) {
            output.data[i] = input.data[i/2];
        }
    }
1110825
private ColumnVector forwardSubstitution(ColumnVector b)
        throws MatrixException
    {
        ColumnVector y = new ColumnVector(nRows);

        // Do forward substitution.
        for (int r = 0; r < nRows; ++r) {
            int   pr  = permutation[r];     // permuted row index
            double dot = 0;
            for (int c = 0; c < r; ++c) {
                dot += LU.at(pr, c)*y.at(c);
            }
            y.set(r, b.at(pr) - dot);
        }

        return y;
    }
1005258
public static DMatrixRMaj convert(DMatrix2 input , DMatrixRMaj output ) {
        if( output == null)
            output = new DMatrixRMaj(2,1);

        if( output.getNumRows() != 1 && output.getNumCols() != 1 )
            throw new IllegalArgumentException("One row or column must have a length of 1 for it to be a vector");
        int length = Math.max(output.getNumRows(),output.getNumCols());
        if( length != 2 )
            throw new IllegalArgumentException("Length of input vector is not 2.  It is "+length);

        output.data[0] = input.a1;
        output.data[1] = input.a2;

        return output;
    }
384351
protected void decompose() throws MatrixException
    {
        // Return if the decomposition is valid.
        if (LU != null) return;

        // Create a new LU matrix and permutation vector.
        // LU is initially just a copy of the values of this system.
        LU = new SquareMatrix(this.copyValues2D());
        permutation = new int[nRows];

        double scales[] = new double[nRows];

        // Loop to initialize the permutation vector and scales.
        for (int r = 0; r < nRows; ++r) {
            permutation[r] = r;     // initially no row exchanges

            // Find the largest row element.
            double largestRowElmt = 0;
            for (int c = 0; c < nRows; ++c) {
                double elmt = Math.abs(LU.at(r, c));
                if (largestRowElmt < elmt) largestRowElmt = elmt;
            }

            // Set the scaling factor for row equilibration.
            if (largestRowElmt != 0) {
                scales[r] = 1/largestRowElmt;
            }
            else {
                throw new MatrixException(MatrixException.ZERO_ROW);
            }
        }

        // Do forward elimination with scaled partial row pivoting.
        forwardElimination(scales);

        // Check bottom right element of the permuted matrix.
        if (LU.at(permutation[nRows - 1], nRows - 1) == 0) {
            throw new MatrixException(MatrixException.SINGULAR);
        }
    }
138715
public static Matrix create(int rows, int cols, Type matrixType) {
        switch (matrixType) {
            case SPARSE_IN_MEMORY:
                return new YaleSparseMatrix(rows, cols);
            case DENSE_IN_MEMORY:
                return new ArrayMatrix(rows, cols);
            case DIAGONAL_IN_MEMORY:
                return new DiagonalMatrix(rows);
            case SPARSE_ON_DISK:
                //return new SparseOnDiskMatrix(rows, cols);
                // REMDINER: implement me
                return new OnDiskMatrix(rows, cols);
            case DENSE_ON_DISK:
                return new OnDiskMatrix(rows, cols);
        }
        throw new IllegalArgumentException(
                "Unknown matrix type: " + matrixType);
    }
463781
public Assignments cluster(Matrix m,
                               int maxClusters,
                               Properties props) {
        int startSize = Integer.parseInt(props.getProperty(
                NUM_CLUSTERS_START, DEFAULT_NUM_CLUSTERS_START));
        int numGaps = Integer.parseInt(props.getProperty(
                NUM_REFERENCE_DATA_SETS, DEFAULT_NUM_REFERENCE_DATA_SETS));
        int numIterations = maxClusters - startSize;
        String criterion = props.getProperty(METHOD_PROPERTY, DEFAULT_METHOD);

        verbose("Transforming the original data set");
        Transform tfidf = new TfIdfDocStripedTransform();
        Transform rowMag = new RowMagnitudeTransform();
        m = rowMag.transform(tfidf.transform(m));

        verbose("Generating the reference data set");
        // Generate the reference data sets.
        ReferenceDataGenerator generator = new ReferenceDataGenerator(m);
        Matrix[] gapMatrices = new Matrix[numGaps];
        for (int i = 0; i < numGaps; ++i)
            gapMatrices[i] = rowMag.transform(tfidf.transform(
                    generator.generateTestData()));

        double[] gapResults = new double[numIterations];
        double[] gapStds = new double[numIterations];
        Assignments[] gapAssignments = new Assignments[numIterations];

        Assignments bestAssignments = null;
        double bestGap = Double.NEGATIVE_INFINITY;
        int bestK = 0;
        // Compute the gap statistic for each iteration.
        for (int i = 0; i < numIterations; ++i) {
            clusterIteration(i, startSize, criterion, m, gapMatrices,
                             gapResults, gapStds, gapAssignments);
            if (bestGap >= (gapResults[i] - gapStds[i])) {
                break;
            }

            // Otherwise, continue clustering with higher values of k.
            bestGap = gapResults[i];
            bestAssignments = gapAssignments[i];
            bestK = i + startSize;
        }

        return bestAssignments;
    }
1018890
private  void pldltpt() throws Exception{
		//matrix S will be changed by the factorization, ad we do not want to change the matrix passed in by the client
		DoubleMatrix2D S = (rescaler==null)? this.Q.copy() : this.Q;
		int n = S.rows();
		DoubleMatrix2D A = S.copy();
		this.P = DoubleFactory2D.sparse.identity(n);
		this.D = DoubleFactory2D.sparse.make(n, n);
		this.L = DoubleFactory2D.sparse.make(n, n);
		DoubleMatrix2D LT = ALG.transpose(L);//remove and work only with L
		int s = 0;
		for(int j=0; j<n; j++){
			//log.debug("j: " + j);
			DoubleMatrix2D LPart = null;
			DoubleMatrix2D LTPart = null;
			DoubleMatrix2D APart = null;
			DoubleMatrix2D DPart = null;
			double ajj = A.getQuick(j, j);
			if(Math.abs(ajj) > 1.e-16){
				//1 x pivot
				D.setQuick(j, j, ajj);
				s = 1;
			}else{
				//ajj = 0, so the 2x2 matrix with ajj in its upper left position
				//is non singular if a(j+1, j)=a(j, j+1) !=0 
				int k=-1;
				for (int r = j + 1; r < n; r++) {
					if (Math.abs(A.getQuick(r, j)) > 1.e-16) {
						k = r;
						break;
					}
				}
				if(k<0){
					throw new Exception("singular matrix");
				}
				//Symmetrically permute row/column k to position j + 1.
				A = ColtUtils.symmPermutation(A, k, j+1);
				P.setQuick(k, k, 0);
				P.setQuick(j + 1, j + 1, 0);
				P.setQuick(k, j + 1, 1);
				P.setQuick(j + 1, k, 1);
				//Choose a 2 � 2 pivot,
				//D(j:j+1)(j:j+1) =  A(j:j+1)(j:j+1)
				D.setQuick(j, j, A.getQuick(j, j));
				D.setQuick(j, j+1, A.getQuick(j, j+1));
				D.setQuick(j+1, j, A.getQuick(j+1, j));
				D.setQuick(j+1, j+1, A.getQuick(j+1, j+1));
				s=2;
			}
			//log.debug("s: " + s);
			
			// L(j:n)(j:j+s-1) = A(j:n)(j:j+s-1).DInv(j:j+s-1)(j:j+s-1)
			APart = A.viewPart(j, j, n - j, 1 + s - 1);
			DPart = ALG.inverse(D.viewPart(j, j, 1 + s - 1, 1 + s - 1));
			LPart = L.viewPart(j, j, n - j, 1 + s - 1);
			DoubleMatrix2D AD = ALG.mult(APart, DPart);
			for (int r = 0; r < LPart.rows(); r++) {
				for (int c = 0; c < LPart.columns(); c++) {
					LPart.setQuick(r, c, AD.getQuick(r, c));
				}
			}
			
			// A(j+s-1:n)(j+s-1:n) = A(j+s-1:n)(j+s-1:n) - L(j+s-1:n)(j:j+s-1).D(j:j+s-1)(j:j+s-1).LT(j:j+s-1)(j+s-1:n)
			LPart = L.viewPart(j + s - 1, j, n - (j + s - 1), 1 + s - 1);
			DPart = D.viewPart(j, j, 1 + s - 1, 1 + s - 1);
			LTPart = LT.viewPart(j, j + s - 1, s, n - (j + s - 1));
			APart = A.viewPart(j + s - 1, j + s - 1, n - (j + s - 1), n - (j + s - 1));
			DoubleMatrix2D LDLT = ALG.mult(LPart, ALG.mult(DPart, LTPart));
			for (int r = 0; r < APart.rows(); r++) {
				for (int c = 0; c < APart.columns(); c++) {
					APart.setQuick(r, c, APart.getQuick(r, c) - LDLT.getQuick(r, c));
				}
			}
			//logger.debug("L: " + ArrayUtils.toString(L.toArray()));
			//logger.debug("A: " + ArrayUtils.toString(A.toArray()));
			
			j = j + s - 1;
		}
	}
336737
public void permute(final int[] arr) {
    checkArgument(arr.length == sources.length);
    final int[] tmp = new int[arr.length];
    for (int i = 0; i < tmp.length; ++i) {
      tmp[i] = arr[sources[i]];
    }
    System.arraycopy(tmp, 0, arr, 0, arr.length);
  }
1340920
private static Matrix permuteMatrix(Matrix mat, int cpRows, int cpCols) {
		//fix up cp to be positive
		if(cpRows == 0 && cpCols == 0) {
			return mat.copy();
		}
		if(cpRows < 0) {
			cpRows = mat.getRowDimension()+cpRows;
		}
		if(cpRows < 0 || cpRows >= mat.getRowDimension()) {
			throw new ArrayIndexOutOfBoundsException( String.format(
					"Can't permute rows by %d: only %d rows.",
					cpRows, mat.getRowDimension() )
			);
		}

		if(cpCols < 0) {
			cpCols = mat.getColumnDimension()+cpCols;
		}
		if(cpCols < 0 || cpCols >= mat.getColumnDimension()) {
			throw new ArrayIndexOutOfBoundsException( String.format(
					"Can't permute cols by %d: only %d rows.",
					cpCols, mat.getColumnDimension() )
			);
		}

		int[] rows = new int[mat.getRowDimension()];
		for(int i=0;i<rows.length;i++) {
			rows[i] = (i+cpRows)%rows.length;
		}
		int[] cols = new int[mat.getColumnDimension()];
		for(int i=0;i<cols.length;i++) {
			cols[i] = (i+cpCols)%cols.length;
		}

		Matrix newMat = mat.getMatrix(rows, cols);
		assert(newMat.getRowDimension() == mat.getRowDimension());
		assert(newMat.getColumnDimension() == mat.getColumnDimension());
		assert(newMat.get(0, 0) ==
			mat.get(cpRows%mat.getRowDimension(), cpCols%mat.getColumnDimension()));


		return newMat;
	}
275311
private void pldltptBK() throws Exception {
		//matrix S will be changed by the factorization, ad we do not want to change the matrix passed in by the client
		DoubleMatrix2D S = (rescaler==null)? this.Q.copy() : this.Q;
		int n = S.rows();
		DoubleMatrix2D A = S.copy();
		this.P = DoubleFactory2D.sparse.identity(n);
		this.D = DoubleFactory2D.sparse.make(n, n);
		this.L = DoubleFactory2D.sparse.make(n, n);
		DoubleMatrix2D LT = ALG.transpose(L);//@TODO: remove and work only with L
		double mu = (1 + Math.sqrt(17)) / 8;
		int s = 0;
		for (int j = 0; j < n; j++) {
			//log.debug("j: " + j);
			DoubleMatrix2D LPart = null;
			DoubleMatrix2D LTPart = null;
			DoubleMatrix2D APart = null;
			DoubleMatrix2D DPart = null;
			
			//Determine largest off-diagonal element in column j
			double lambda = -1;
			int r = -1;
			for(int rr = j+1; rr<n; rr++){
				double d = Math.abs(A.getQuick(rr, j)); 
				if(d > lambda){
					lambda = d;
					r = rr;
				}
			}
			double ajj = A.getQuick(j, j);
			if (Math.abs(ajj) > mu*lambda) {
				// Choose a 1 x 1 pivot,
				D.setQuick(j, j, ajj);
				s = 1;
			} else {
				//Determine largest off-diagonal element in row/column r
				double sigma = -1;
				for(int rr = 0; rr<n; rr++){
					if(rr==r){
						continue;
					}
					double d = Math.abs(A.getQuick(rr, r)); 
					if(d > sigma){
						sigma = d;
					}
				}
				if(sigma * Math.abs(ajj) > mu*lambda*lambda){
					//Choose a 1 x 1 pivot,
					D.setQuick(j, j, A.getQuick(j, j));
					s = 1;
				}else if(Math.abs(A.getQuick(r, r)) > mu*sigma){
					//Symmetrically permute row/column r to position j
				  A = ColtUtils.symmPermutation(A, r, j);
					P.setQuick(r, r, 0);
					P.setQuick(j, j, 0);
					P.setQuick(r, j, 1);
					P.setQuick(j, r, 1);
					//Choose a 1 x 1 pivot,
					D.setQuick(j, j, A.getQuick(j, j));
					s = 1;
				}else{
					//Symmetrically permute rows/columns r and p to positions j and j + 1 respectively
					int t = j+1;
					//int v = j+1;
					A = ColtUtils.symmPermutation(A, r, t);
					P.setQuick(r, r, 0);
					P.setQuick(t, t, 0);
					P.setQuick(r, t, 1);
					P.setQuick(t, r, 1);
//					A = Utils.symmPermutation(A, p, v);
//					P.setQuick(p, p, 0);
//					P.setQuick(v, v, 0);
//					P.setQuick(p, v, 1);
//					P.setQuick(v, p, 1);
					//Choose a 2 x pivot,
					//D(j:j+1)(j:j+1) =  A(j:j+1)(j:j+1)
					D.setQuick(j, j, A.getQuick(j, j));
					D.setQuick(j, j+1, A.getQuick(j, j+1));
					D.setQuick(j+1, j, A.getQuick(j+1, j));
					D.setQuick(j+1, j+1, A.getQuick(j+1, j+1));
					s = 2;
				}
			}
			//log.debug("s: " + s);

			// L(j:n)(j:j+s-1) = A(j:n)(j:j+s-1).DInv(j:j+s-1)(j:j+s-1)
			APart = A.viewPart(j, j, n - j, 1 + s - 1);
//			if(s==2){
//				try{
//					ColtUtils.invert2x2Matrix(D.viewPart(j, j, 1 + s - 1, 1 + s - 1));
//				}catch(Exception e){
//					throw e;
//				}
//			}
			DPart = ALG.inverse(D.viewPart(j, j, 1 + s - 1, 1 + s - 1));
			LPart = L.viewPart(j, j, n - j, 1 + s - 1);
			DoubleMatrix2D AD = ALG.mult(APart, DPart);
			for (int rr = 0; rr < LPart.rows(); rr++) {
				for (int cc = 0; cc < LPart.columns(); cc++) {
					LPart.setQuick(rr, cc, AD.getQuick(rr, cc));
				}
			}

			// A(j+s-1:n)(j+s-1:n) = A(j+s-1:n)(j+s-1:n) -
			// L(j+s-1:n)(j:j+s-1).D(j:j+s-1)(j:j+s-1).LT(j:j+s-1)(j+s-1:n)
			LPart = L.viewPart(j + s - 1, j, n - (j + s - 1), 1 + s - 1);
			DPart = D.viewPart(j, j, 1 + s - 1, 1 + s - 1);
			LTPart = LT.viewPart(j, j + s - 1, s, n - (j + s - 1));
			APart = A.viewPart(j + s - 1, j + s - 1, n - (j + s - 1), n - (j + s - 1));
			DoubleMatrix2D LDLT = ALG.mult(LPart, ALG.mult(DPart, LTPart));
			for (int rr = 0; rr < APart.rows(); rr++) {
				for (int cc = 0; cc < APart.columns(); cc++) {
					APart.setQuick(rr, cc, APart.getQuick(rr, cc) - LDLT.getQuick(rr, cc));
				}
			}
			//logger.debug("L: " + ArrayUtils.toString(L.toArray()));
			//logger.debug("A: " + ArrayUtils.toString(A.toArray()));

			j = j + s - 1;
		}
	}
494122
public static PermutationMatrix getInstance(int n)
    {
        PermutationMatrix pm = map.get(n);
        if (pm == null)
        {
            pm = new PermutationMatrix(n);
            map.put(n, pm);
        }
        return pm;
    }
176267
private ColumnVector backSubstitution(ColumnVector y)
        throws MatrixException
    {
        ColumnVector x = new ColumnVector(nRows);

        // Do back substitution.
        for (int r = nRows - 1; r >= 0; --r) {
            int   pr  = permutation[r];     // permuted row index
            double dot = 0;
            for (int c = r+1; c < nRows; ++c) {
                dot += LU.at(pr, c)*x.at(c);
            }
            x.set(r, (y.at(r) - dot)/LU.at(pr, r));
        }

        return x;
    }
867218
def permutation_t(vector)
      ByteVector.new(
        vector.each.with_index.inject([]) do |b_arr, (byte, index)|
          b_arr[T[index]] = byte
          b_arr
        end.map(&:to_i)
      )
    end
658174
def _generate_permutation(self, npoints):
        
        i = np.arange(0, npoints)
        # permutation
        p = np.random.permutation(npoints)
        ip = np.empty_like(p)
        # inverse permutation
        ip[p[i]] = i
        return p, ip
434175
private IntegerVector createSemanticVector() {
        IntegerVector v = (useSparseSemantics) 
            ? new CompactSparseIntegerVector(vectorLength)
            : new DenseIntVector(vectorLength);
        return v;
    }