Description: This routine finds the total number of non-zero elements and
the number of non-zero elements per row or column in the dense matrix A.
491598
public static int cusparseSnnz(
        cusparseHandle handle, 
        int dirA, 
        int m, 
        int n, 
        cusparseMatDescr descrA, 
        Pointer A, 
        int lda, 
        Pointer nnzPerRowCol, 
        Pointer nnzTotalDevHostPtr)
    {
        return checkResult(cusparseSnnzNative(handle, dirA, m, n, descrA, A, lda, nnzPerRowCol, nnzTotalDevHostPtr));
    }
370157
private static void hqr(DenseMatrix A, double[] d, double[] e) {
        int n = A.nrows();
        int nn, m, l, k, j, its, i, mmin;
        double z, y, x, w, v, u, t, s, r = 0.0, q = 0.0, p = 0.0, anorm = 0.0;

        for (i = 0; i < n; i++) {
            for (j = Math.max(i - 1, 0); j < n; j++) {
                anorm += Math.abs(A.get(i, j));
            }
        }
        nn = n - 1;
        t = 0.0;
        while (nn >= 0) {
            its = 0;
            do {
                for (l = nn; l > 0; l--) {
                    s = Math.abs(A.get(l - 1, l - 1) + Math.abs(A.get(l, l)));
                    if (s == 0.0) {
                        s = anorm;
                    }
                    if (Math.abs(A.get(l, l - 1)) <= Math.EPSILON * s) {
                        A.set(l, l - 1, 0.0);
                        break;
                    }
                }
                x = A.get(nn, nn);
                if (l == nn) {
                    d[nn--] = x + t;
                } else {
                    y = A.get(nn - 1, nn - 1);
                    w = A.get(nn, nn - 1) * A.get(nn - 1, nn);
                    if (l == nn - 1) {
                        p = 0.5 * (y - x);
                        q = p * p + w;
                        z = Math.sqrt(Math.abs(q));
                        x += t;
                        if (q >= 0.0) {
                            z = p + Math.copySign(z, p);
                            d[nn - 1] = d[nn] = x + z;
                            if (z != 0.0) {
                                d[nn] = x - w / z;
                            }
                        } else {
                            d[nn] = x + p;
                            e[nn] = -z;
                            d[nn - 1] = d[nn];
                            e[nn - 1] = -e[nn];
                        }
                        nn -= 2;
                    } else {
                        if (its == 30) {
                            throw new IllegalStateException("Too many iterations in hqr");
                        }
                        if (its == 10 || its == 20) {
                            t += x;
                            for (i = 0; i < nn + 1; i++) {
                                A.sub(i, i, x);
                            }
                            s = Math.abs(A.get(nn, nn - 1)) + Math.abs(A.get(nn - 1, nn - 2));
                            y = x = 0.75 * s;
                            w = -0.4375 * s * s;
                        }
                        ++its;
                        for (m = nn - 2; m >= l; m--) {
                            z = A.get(m, m);
                            r = x - z;
                            s = y - z;
                            p = (r * s - w) / A.get(m + 1, m) + A.get(m, m + 1);
                            q = A.get(m + 1, m + 1) - z - r - s;
                            r = A.get(m + 2, m + 1);
                            s = Math.abs(p) + Math.abs(q) + Math.abs(r);
                            p /= s;
                            q /= s;
                            r /= s;
                            if (m == l) {
                                break;
                            }
                            u = Math.abs(A.get(m, m - 1)) * (Math.abs(q) + Math.abs(r));
                            v = Math.abs(p) * (Math.abs(A.get(m - 1, m - 1)) + Math.abs(z) + Math.abs(A.get(m + 1, m + 1)));
                            if (u <= Math.EPSILON * v) {
                                break;
                            }
                        }
                        for (i = m; i < nn - 1; i++) {
                            A.set(i + 2, i, 0.0);
                            if (i != m) {
                                A.set(i + 2, i - 1, 0.0);
                            }
                        }
                        for (k = m; k < nn; k++) {
                            if (k != m) {
                                p = A.get(k, k - 1);
                                q = A.get(k + 1, k - 1);
                                r = 0.0;
                                if (k + 1 != nn) {
                                    r = A.get(k + 2, k - 1);
                                }
                                if ((x = Math.abs(p) + Math.abs(q) + Math.abs(r)) != 0.0) {
                                    p /= x;
                                    q /= x;
                                    r /= x;
                                }
                            }
                            if ((s = Math.copySign(Math.sqrt(p * p + q * q + r * r), p)) != 0.0) {
                                if (k == m) {
                                    if (l != m) {
                                        A.set(k, k - 1, -A.get(k, k - 1));
                                    }
                                } else {
                                    A.set(k, k - 1, -s * x);
                                }
                                p += s;
                                x = p / s;
                                y = q / s;
                                z = r / s;
                                q /= p;
                                r /= p;
                                for (j = k; j < nn + 1; j++) {
                                    p = A.get(k, j) + q * A.get(k + 1, j);
                                    if (k + 1 != nn) {
                                        p += r * A.get(k + 2, j);
                                        A.sub(k + 2, j, p * z);
                                    }
                                    A.sub(k + 1, j, p * y);
                                    A.sub(k, j, p * x);
                                }
                                mmin = nn < k + 3 ? nn : k + 3;
                                for (i = l; i < mmin + 1; i++) {
                                    p = x * A.get(i, k) + y * A.get(i, k + 1);
                                    if (k + 1 != nn) {
                                        p += z * A.get(i, k + 2);
                                        A.sub(i, k + 2, p * r);
                                    }
                                    A.sub(i, k + 1, p * q);
                                    A.sub(i, k, p);
                                }
                            }
                        }
                    }
                }
            } while (l + 1 < nn);
        }
    }
494226
private static int scatter(SparseMatrix A, int j, double beta, int[] w, double[] x, int mark, SparseMatrix C, int nz) {
        int[] Ap = A.colIndex;
        int[] Ai = A.rowIndex;
        double[] Ax = A.x;

        int[] Ci = C.rowIndex;
        for (int p = Ap[j]; p < Ap[j + 1]; p++) {
            int i = Ai[p];                // A(i,j) is nonzero
            if (w[i] < mark) {
                w[i] = mark;              // i is new entry in column j
                Ci[nz++] = i;             // add i to pattern of C(:,j)
                x[i] = beta * Ax[p];      // x(i) = beta*A(i,j)
            } else {
                x[i] += beta * Ax[p];     // i exists in C(:,j) already
            }
        }

        return nz;
    }
714096
private void computePrincipleVectors(int numPrincipleVectors) {

        // Group all of the sspace's vectors into a Matrix so that we can
        // cluster them
        final int numTerms = sspace.getWords().size();
        final List<DoubleVector> termVectors = 
            new ArrayList<DoubleVector>(numTerms);
        String[] termAssignments = new String[numTerms];
        int k = 0;
        for (String term : sspace.getWords()) {
            termVectors.add(Vectors.asDouble(sspace.getVector(term)));
            termAssignments[k++] = term;
        }

        Random random = new Random();

        final DoubleVector[] principles = new DoubleVector[numPrincipleVectors];
        for (int i = 0; i < principles.length; ++i)
            principles[i] = new DenseVector(sspace.getVectorLength());

        // Randomly assign the data set to different intitial clusters 
        final MultiMap<Integer,Integer> clusterAssignment = 
            new HashMultiMap<Integer,Integer>();
        for (int i = 0; i < numTerms; ++i) 
            clusterAssignment.put(random.nextInt(numPrincipleVectors), i);
        
        int numIters = 1;
        for (int iter = 0; iter < numIters; ++iter) {
            verbose(LOGGER, "Computing principle vectors (round %d/%d)", iter+1, numIters);
            // Compute the centroids
            for (Map.Entry<Integer,Set<Integer>> e
                     : clusterAssignment.asMap().entrySet()) {
                int cluster = e.getKey();
                DoubleVector principle = new DenseVector(sspace.getVectorLength());
                principles[cluster] = principle;
                for (Integer row : e.getValue())
                    VectorMath.add(principle, termVectors.get(row));
            }

            // Reassign each element to the centroid to which it is closest
            clusterAssignment.clear();
            final int numThreads = workQueue.availableThreads();
            Object key = workQueue.registerTaskGroup(numThreads);

            for (int threadId_ = 0; threadId_ < numThreads; ++threadId_) {
                final int threadId = threadId_;
                workQueue.add(key, new Runnable() {
                        public void run() {
                            // Thread local cache of all the cluster assignments
                            MultiMap<Integer,Integer> clusterAssignment_ = 
                                new HashMultiMap<Integer,Integer>();
                            // For each of the vectors that this thread is
                            // responsible for, find the principle vector to
                            // which it is closest
                            for (int i = threadId; i < numTerms; i += numThreads) {
                                DoubleVector v = termVectors.get(i);
                                double highestSim = -Double.MAX_VALUE;
                                int pVec = -1;
                                for (int j = 0; j < principles.length; ++j) {
                                    DoubleVector principle = principles[j];
                                    double sim = Similarity.cosineSimilarity(
                                        v, principle);
                                    assert sim >= -1 && sim <= 1 : "similarity "
                                        + " to principle vector " + j + " is "
                                        + "outside the expected range: " + sim;
                                    if (sim > highestSim) {
                                        highestSim = sim;
                                        pVec = j;
                                    }
                                }
                                assert pVec != -1 : "Could not find match to "
                                    + "any of the " + principles.length + 
                                    " principle vectors";
                                clusterAssignment_.put(pVec, i);
                            }
                            // Once all the vectors for this thread have been
                            // mapped, update the thread-shared mapping.  Do
                            // this here, rather than for each comparison to
                            // minimize the locking
                            synchronized (clusterAssignment) {
                                clusterAssignment.putAll(clusterAssignment_);
                            }
                        }
                    });
            }
            workQueue.await(key);
        }

        double mean = numTerms / (double)numPrincipleVectors;
        double variance = 0d;
        for (Map.Entry<Integer,Set<Integer>> e
                 : clusterAssignment.asMap().entrySet()) {
            Set<Integer> rows = e.getValue();
            Set<String> terms = new HashSet<String>();
            for (Integer i : rows)
                terms.add(termAssignments[i]);
             verbose(LOGGER, "Principle vectod %d is closest to %d terms",
                     e.getKey(), terms.size());
            double diff = mean - terms.size();
            variance += diff * diff;
            principleVectorToNearestTerms.putMany(
                principles[e.getKey()], terms);
        }

        verbose(LOGGER, "Average number terms per principle vector: %f, "+
                "(%f stddev)", mean, Math.sqrt(variance / numPrincipleVectors));
    }
1250225
public static ColumnMajorSparseMatrix zero(int rows, int columns, int capacity) {
        return CCSMatrix.zero(rows, columns, capacity);
    }
864754
public static void permutationVector( DMatrixSparseCSC P , int[] vector) {
        if( P.numCols != P.numRows ) {
            throw new MatrixDimensionException("Expected a square matrix");
        } else if( P.nz_length != P.numCols ) {
            throw new IllegalArgumentException("Expected N non-zero elements in permutation matrix");
        } else if( vector.length < P.numCols ) {
            throw new IllegalArgumentException("vector is too short");
        }

        int M = P.numCols;

        for (int i = 0; i < M; i++) {
            if( P.col_idx[i+1] != i+1 )
                throw new IllegalArgumentException("Unexpected number of elements in a column");

            vector[P.nz_rows[i]] = i;
        }
    }
883012
def p_matrix(p):
    
    if len(p) == 3:
        p[0] = node.matrix()
    else:
        p[0] = node.matrix(p[2])
582198
def +(matrix)
    result = @v
    if have_the_same_dimensions matrix
      @lines.times do |x|
        @columns.times do |y|
          result[x][y] += matrix.v[x][y]
        end
      end
    elsif matrix.is_a?Numeric
      @lines.times do |x|
        @columns.times do |y|
          result[x][y] += matrix
        end
      end
    else
      Error.call "Matrix::+ : Impossible de calculer cela", Error::ERR_HIGH
    end
    Matrix.new result
  end
230565
def fill
      @matrix[0][0] = 0
      # Set up the first column on each row.
      1.upto(@rows-1) {|i| @matrix[i][0] = @matrix[i-1][0] + @scoring.score_delete(@seq1[i])}
      # Set up the first row 
      1.upto(@cols-1) {|j| @matrix[0][j] = @matrix[0][j-1] + @scoring.score_insert(@seq2[j])}

      1.upto(@rows-1) do |i|
        prv_row = @matrix[i-1]
        cur_row = @matrix[i]

        1.upto(@cols-1) do |j|
          
          seq1_obj = @seq1[i-1]
          seq2_obj = @seq2[j-1]

          # Calculate the score.
          score_align = prv_row[j-1] + @scoring.score_align(seq1_obj, seq2_obj)
          score_delete = prv_row[j] + @scoring.score_delete(seq1_obj)
          score_insert = cur_row[j-1] + @scoring.score_insert(seq2_obj)
          max = max3(score_align, score_delete, score_insert)

          @matrix[i][j] = max
        end
      end
    end
105926
def coarsen_matrix(Z, xlevel=0, ylevel=0, method='average'):
    

    # coarsen x
    if not ylevel:
        Z_coarsened = Z
    else:
        temp = []
        for z in Z: temp.append(coarsen_array(z, ylevel, method))
        Z_coarsened = _n.array(temp)

    # coarsen y
    if xlevel:
        Z_coarsened = Z_coarsened.transpose()
        temp = []
        for z in Z_coarsened: temp.append(coarsen_array(z, xlevel, method))
        Z_coarsened = _n.array(temp).transpose()

    return Z_coarsened



    # first coarsen the columns (if necessary)
    if ylevel:
        Z_ycoarsened = []
        for c in Z: Z_ycoarsened.append(coarsen_array(c, ylevel, method))
        Z_ycoarsened = _n.array(Z_ycoarsened)

    # now coarsen the rows
    if xlevel: return coarsen_array(Z_ycoarsened, xlevel, method)
    else:      return _n.array(Z_ycoarsened)
958530
public static double innerProdA(DMatrixD1 x, DMatrixD1 A , DMatrixD1 y )
    {
        int n = A.numRows;
        int m = A.numCols;

        if( x.getNumElements() != n )
            throw new IllegalArgumentException("Unexpected number of elements in x");
        if( y.getNumElements() != m )
            throw new IllegalArgumentException("Unexpected number of elements in y");

        double result = 0;

        for( int i = 0; i < m; i++ ) {
            double total = 0;

            for( int j = 0; j < n; j++ ) {
                total += x.get(j)*A.unsafe_get(j,i);
            }

            result += total*y.get(i);
        }

        return result;
    }
411255
public Matrix getD()
    {
        Matrix X = new DenseMatrix(n, n);
        for (int i = 0; i < n; i++)
        {
            X.set(i, i, d[i]);
            
            if (e[i] > 0)
                X.set(i, i+1, e[i]);
            else if (e[i] < 0)
                X.set(i, i-1, e[i]);
        }
        return X;
    }
599517
public AbstractMatrix transpose()
    {
        int m = rows();
        int n = columns();
        AbstractMatrix tr = getInstance(n, m, cls);
        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < n; j++)
            {
                Array.set(tr.array, M.at(j, i), Array.get(array, M.at(i, j)));
            }
        }
        return tr;
    }
1151377
function() {
    var matrix_rows = [];
    var n = this.elements.length, k = n, i;
    do { i = k - n;
      matrix_rows.push(Vector.create(this.elements[i]).inspect());
    } while (--n);
    return matrix_rows.join('\n');
  }
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;
            }
        }
    }
54798
function() {
    if (!this.isSquare() || this.isSingular()) { return null; }
    var ni = this.elements.length, ki = ni, i, j;
    var M = this.augment(Matrix.I(ni)).toRightTriangular();
    var np, kp = M.elements[0].length, p, els, divisor;
    var inverse_elements = [], new_element;
    // Matrix is non-singular so there will be no zeros on the diagonal
    // Cycle through rows from last to first
    do { i = ni - 1;
      // First, normalise diagonal elements to 1
      els = []; np = kp;
      inverse_elements[i] = [];
      divisor = M.elements[i][i];
      do { p = kp - np;
        new_element = M.elements[i][p] / divisor;
        els.push(new_element);
        // Shuffle of the current row of the right hand side into the results
        // array as it will not be modified by later runs through this loop
        if (p >= ki) { inverse_elements[i].push(new_element); }
      } while (--np);
      M.elements[i] = els;
      // Then, subtract this row from those above it to
      // give the identity matrix on the left hand side
      for (j = 0; j < i; j++) {
        els = []; np = kp;
        do { p = kp - np;
          els.push(M.elements[j][p] - M.elements[i][p] * M.elements[j][i]);
        } while (--np);
        M.elements[j] = els;
      }
    } while (--ni);
    return Matrix.create(inverse_elements);
  }
67059
function dfFromMatrixHelper(matrix, startRow, colNames) {
  var nCol = colNames.size();
  var nRow = matrix.length - startRow;
  var columns = allocArray(nCol);
  var j;
  for (j = 0; j < nCol; j++) {
    columns[j] = allocArray(nRow);
  }
  for (var i = 0; i < nRow; i++) {
    var rowArray = matrix[i + startRow];
    if (rowArray.length !== nCol) {
      throw new Error('all row arrays must be of the same size');
    }
    for (j = 0; j < nCol; j++) {
      columns[j][i] = rowArray[j];
    }
  }
  return newDataFrame(columns, colNames);
}