Rank-1 operation on real symmetric matrix.
420224
def cublasDsyr(handle, uplo, n, alpha, x, incx, A, lda):
    

    status = _libcublas.cublasDsyr_v2(handle,
                                      _CUBLAS_FILL_MODE[uplo], n, 
                                      ctypes.byref(ctypes.c_double(alpha)), 
                                      int(x), incx, int(A), lda)
    cublasCheckStatus(status)
1195500
static void compute_B(DMatrix2x2 B , DMatrixRMaj R_v , double v1 , double v2 ) {
		B.a11 = R_v.data[0]+R_v.data[6]*-v1;
		B.a12 = R_v.data[1]+R_v.data[7]*-v1;
		B.a21 = R_v.data[3]+R_v.data[6]*-v2;
		B.a22 = R_v.data[4]+R_v.data[7]*-v2;
	}
1046616
private void hqr2(Matrix H, Matrix V, Vector d, Vector e) {

        // This is derived from the Algol procedure hqr2,
        // by Martin and Wilkinson, Handbook for Auto. Comp.,
        // Vol.ii-Linear Algebra, and the corresponding
        // Fortran subroutine in EISPACK.

        // Initialize

        int nn = e.length();
        int n = nn - 1;
        int low = 0;
        int high = nn - 1;
        double eps = Math.pow(2.0, -52.0);
        double exshift = 0.0;
        double p = 0;
        double q = 0;
        double r = 0;
        double s = 0;
        double z = 0;
        double t;
        double w;
        double x;
        double y;

        // Store roots isolated by balanc and compute matrix norm

        double norm = 0.0;
        for (int i = 0; i < nn; i++) {
            if (i < low | i > high) {
                d.set(i, H.get(i, i));
                e.set(i, 0.0);
            }
            for (int j = Math.max(i - 1, 0); j < nn; j++) {
                norm = norm + Math.abs(H.get(i, j));
            }
        }

        // Outer loop over eigenvalue index

        int iter = 0;
        while (n >= low) {

            // Look for single small sub-diagonal element

            int l = n;
            while (l > low) {
                s = Math.abs(H.get(l - 1, l - 1)) 
                    + Math.abs(H.get(l, l));

                if (s == 0.0) {
                    s = norm;
                }
                if (Math.abs(H.get(l, l - 1)) < eps * s) {
                    break;
                }
                l--;
            }

            // Check for convergence
            // One root found

            if (l == n) {
                H.updateAt(n, n, Matrices.asPlusFunction(exshift));
                d.set(n, H.get(n, n));
                e.set(n, 0.0);
                n--;
                iter = 0;

                // Two roots found

            } else if (l == n - 1) {
                w = H.get(n, n - 1) * H.get(n - 1, n);
                p = (H.get(n - 1, n - 1) - H.get(n, n)) / 2.0;
                q = p * p + w;
                z = Math.sqrt(Math.abs(q));
                H.updateAt(n, n, Matrices.asPlusFunction(exshift));
                H.updateAt(n - 1, n - 1, Matrices.asPlusFunction(exshift));
                x = H.get(n, n);

                // Real pair

                if (q >= 0) {
                    if (p >= 0) {
                        z = p + z;
                    } else {
                        z = p - z;
                    }
                    d.set(n - 1, x + z);
                    d.set(n, d.get(n - 1));
                    if (z != 0.0) {
                        d.set(n, x - w / z);
                    }
                    e.set(n - 1, 0.0);
                    e.set(n, 0.0);
                    x = H.get(n, n - 1);
                    s = Math.abs(x) + Math.abs(z);
                    p = x / s;
                    q = z / s;
                    r = Math.sqrt(p * p + q * q);
                    p = p / r;
                    q = q / r;

                    // Row modification

                    for (int j = n - 1; j < nn; j++) {
                        z = H.get(n - 1, j);
                        H.set(n - 1, j, q * z + p * H.get(n, j));
                        H.set(n, j, q * H.get(n, j) - p * z);
                    }

                    // Column modification

                    for (int i = 0; i <= n; i++) {
                        z = H.get(i, n - 1);
                        H.set(i, n - 1, q * z + p * H.get(i, n));
                        H.set(i, n, q * H.get(i, n) - p * z);
                    }

                    // Accumulate transformations

                    for (int i = low; i <= high; i++) {
                        z = V.get(i, n - 1);
                        V.set(i, n - 1, q * z + p * V.get(i, n));
                        V.set(i, n, q * V.get(i, n) - p * z);
                    }

                    // Complex pair

                } else {
                    d.set(n - 1, x + p);
                    d.set(n, x + p);
                    e.set(n - 1, z);
                    e.set(n, -z);
                }
                n = n - 2;
                iter = 0;

                // No convergence yet

            } else {

                // Form shift

                x = H.get(n, n);
                y = 0.0;
                w = 0.0;
                if (l < n) {
                    y = H.get(n - 1, n - 1);
                    w = H.get(n, n - 1) * H.get(n - 1, n);
                }

                // Wilkinson's original ad hoc shift

                if (iter == 10) {
                    exshift += x;
                    for (int i = low; i <= n; i++) {
                        H.updateAt(i, i, Matrices.asMinusFunction(x));
                    }
                    s = Math.abs(H.get(n, n - 1)) 
                        + Math.abs(H.get(n - 1, n - 2));
                    x = y = 0.75 * s; // WTF ???
                    w = -0.4375 * s * s; // Are you kidding me???
                }

                // MATLAB's new ad hoc shift

                if (iter == 30) {
                    s = (y - x) / 2.0;
                    s = s * s + w;
                    if (s > 0) {
                        s = Math.sqrt(s);
                        if (y < x) {
                            s = -s;
                        }
                        s = x - w / ((y - x) / 2.0 + s);
                        for (int i = low; i <= n; i++) {
                            H.updateAt(i, i, Matrices.asMinusFunction(s));
                        }
                        exshift += s;
                        x = y = w = 0.964;
                    }
                }

                iter = iter + 1; // (Could check iteration count here.)

                // Look for two consecutive small sub-diagonal elements

                int m = n - 2;
                while (m >= l) {
                    z = H.get(m, m);
                    r = x - z;
                    s = y - z;
                    p = (r * s - w) / H.get(m + 1, m) 
                                      + H.get(m, m + 1);
                    q = H.get(m + 1, m + 1) - z - r - s;
                    r = H.get(m + 2, m + 1);
                    s = Math.abs(p) + Math.abs(q) + Math.abs(r);
                    p = p / s;
                    q = q / s;
                    r = r / s;
                    if (m == l) {
                        break;
                    }
                    if (Math.abs(H.get(m, m - 1)) * (Math.abs(q) + Math.abs(r)) < eps
                            * (Math.abs(p) * (Math.abs(H.get(m - 1, m - 1))
                                    + Math.abs(z) + Math.abs(H.get(m + 1, m + 1))))) {
                        break;
                    }
                    m--;
                }

                for (int i = m + 2; i <= n; i++) {
                    H.set(i, i - 2, 0.0);
                    if (i > m + 2) {
                        H.set(i, i - 3, 0.0);
                    }
                }

                // Double QR step involving rows l:n and columns m:n

                for (int k = m; k <= n - 1; k++) {
                    boolean notlast = (k != n - 1);
                    if (k != m) {
                        p = H.get(k, k - 1);
                        q = H.get(k + 1, k - 1);
                        r = (notlast ? H.get(k + 2, k - 1) : 0.0);
                        x = Math.abs(p) + Math.abs(q) + Math.abs(r);
                        if (x == 0.0) {
                            continue;
                        }
                        p = p / x;
                        q = q / x;
                        r = r / x;
                    }

                    s = Math.sqrt(p * p + q * q + r * r);
                    if (p < 0) {
                        s = -s;
                    }
                    if (s != 0) {
                        if (k != m) {
                            H.set(k, k - 1, -s * x);
                        } else if (l != m) {
                            H.updateAt(k, k - 1, Matrices.INV_FUNCTION);
                        }
                        p = p + s;
                        x = p / s;
                        y = q / s;
                        z = r / s;
                        q = q / p;
                        r = r / p;

                        // Row modification

                        for (int j = k; j < nn; j++) {
                            p = H.get(k, j) + q * H.get(k + 1, j);
                            if (notlast) {
                                p = p + r * H.get(k + 2, j);
                                H.updateAt(k + 2, j,
                                    Matrices.asMinusFunction(p * z));
                            }
                            H.updateAt(k, j, Matrices.asMinusFunction(p * x));
                            H.updateAt(k + 1, j, Matrices.asMinusFunction(p * y));
                        }

                        // Column modification

                        for (int i = 0; i <= Math.min(n, k + 3); i++) {
                            p = x * H.get(i, k) + y 
                                * H.get(i, k + 1);
                            if (notlast) {
                                p = p + z * H.get(i, k + 2);
                                H.updateAt(i, k + 2,
                                    Matrices.asMinusFunction(p * r));
                            }
                            H.updateAt(i, k, Matrices.asMinusFunction(p));
                            H.updateAt(i, k + 1, Matrices.asMinusFunction(p * q));
                        }

                        // Accumulate transformations

                        for (int i = low; i <= high; i++) {
                            p = x * V.get(i, k) + y 
                                * V.get(i, k + 1);
                            if (notlast) {
                                p = p + z * V.get(i, k + 2);
                                V.updateAt(i, k + 2,
                                    Matrices.asMinusFunction(p * r));
                            }
                            V.updateAt(i, k, Matrices.asMinusFunction(p));
                            V.updateAt(i, k + 1, Matrices.asMinusFunction(p * q));
                        }
                    } // (s != 0)
                } // k loop
            } // check convergence
        } // while (n >= low)

        // Backsubstitute to find vectors of upper triangular form

        if (norm == 0.0) {
            return;
        }

        for (n = nn - 1; n >= 0; n--) {
            p = d.get(n);
            q = e.get(n);

            // Real vector

            if (q == 0) {
                int l = n;
                H.set(n, n, 1.0);
                for (int i = n - 1; i >= 0; i--) {
                    w = H.get(i, i) - p;
                    r = 0.0;
                    for (int j = l; j <= n; j++) {
                        r = r + H.get(i, j) * H.get(j, n);
                    }
                    if (e.get(i) < 0.0) {
                        z = w;
                        s = r;
                    } else {
                        l = i;
                        if (e.get(i) == 0.0) {
                            if (w != 0.0) {
                                H.set(i, n, -r / w);
                            } else {
                                H.set(i, n, -r / (eps * norm));
                            }

                            // Solve real equations

                        } else {
                            x = H.get(i, i + 1);
                            y = H.get(i + 1, i);
                            q = (d.get(i) - p) * (d.get(i) - p) 
                                + e.get(i) * e.get(i);
                            t = (x * s - z * r) / q;
                            H.set(i, n, t);
                            if (Math.abs(x) > Math.abs(z)) {
                                H.set(i + 1, n, (-r - w * t) / x);
                            } else {
                                H.set(i + 1, n, (-s - y * t) / z);
                            }
                        }

                        // Overflow control

                        t = Math.abs(H.get(i, n));
                        if ((eps * t) * t > 1) {
                            for (int j = i; j <= n; j++) {
                                H.updateAt(j, n, Matrices.asDivFunction(t));
                            }
                        }
                    }
                }

                // Complex vector

            } else if (q < 0) {
                int l = n - 1;

                // Last vector component imaginary so matrix is triangular

                if (Math.abs(H.get(n, n - 1)) 
                        > Math.abs(H.get(n - 1, n))) {
                    H.set(n - 1, n - 1, q / H.get(n, n - 1));
                    H.set(n - 1, n, -(H.get(n, n) - p) 
                                            / H.get(n, n - 1));
                } else {
                    double[] cdiv = cdiv(0.0, -H.get(n - 1, n),
                         H.get(n - 1, n - 1) - p, q);
                    H.set(n - 1, n - 1, cdiv[0]);
                    H.set(n - 1, n, cdiv[1]);
                }
                H.set(n, n - 1, 0.0);
                H.set(n, n, 1.0);
                for (int i = n - 2; i >= 0; i--) {
                    double ra;
                    double sa;
                    double vr;
                    double vi;
                    ra = 0.0;
                    sa = 0.0;
                    for (int j = l; j <= n; j++) {
                        ra = ra + H.get(i, j) * H.get(j, n - 1);
                        sa = sa + H.get(i, j) * H.get(j, n);
                    }
                    w = H.get(i, i) - p;

                    if (e.get(i) < 0.0) {
                        z = w;
                        r = ra;
                        s = sa;
                    } else {
                        l = i;
                        if (e.get(i) == 0) {
                            double[] cdiv = cdiv(-ra, -sa, w, q);
                            H.set(i, n - 1, cdiv[0]);
                            H.set(i, n, cdiv[1]);
                        } else {

                            // Solve complex equations

                            x = H.get(i, i + 1);
                            y = H.get(i + 1, i);
                            vr = (d.get(i) - p) * (d.get(i) - p) 
                                 + e.get(i) * e.get(i) - q * q;
                            vi = (d.get(i) - p) * 2.0 * q;
                            if (vr == 0.0 & vi == 0.0) {
                                vr = eps
                                        * norm
                                        * (Math.abs(w) + Math.abs(q)
                                                + Math.abs(x) + Math.abs(y) + Math
                                                    .abs(z));
                            }
                            double[] cdiv = cdiv(x * r - z * ra + q * sa,
                                                 x * s - z * sa - q * ra, vr, vi);
                            H.set(i, n - 1, cdiv[0]);
                            H.set(i, n, cdiv[1]);
                            if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
                                H.set(i + 1, n - 1, (-ra - w 
                                             * H.get(i, n - 1) + q
                                             * H.get(i, n)) / x);
                                H.set(i + 1, n, (-sa - w 
                                             * H.get(i, n) - q 
                                             * H.get(i, n - 1)) / x);
                            } else {
                                cdiv = cdiv(-r - y 
                                     * H.get(i, n - 1), -s - y 
                                     * H.get(i, n), z, q);
                                H.set(i + 1, n - 1, cdiv[0]);
                                H.set(i + 1, n, cdiv[1]);
                            }
                        }

                        // Overflow control

                        t = Math.max(Math.abs(H.get(i, n - 1)), 
                                     Math.abs(H.get(i, n)));
                        if ((eps * t) * t > 1) {
                            for (int j = i; j <= n; j++) {
                                H.updateAt(j, n - 1, Matrices.asDivFunction(t));
                                H.updateAt(j, n, Matrices.asDivFunction(t));
                            }
                        }
                    }
                }
            }
        }

        // Vectors of isolated roots

        for (int i = 0; i < nn; i++) {
            if (i < low | i > high) {
                for (int j = i; j < nn; j++) {
                    V.set(i, j, H.get(i, j));
                }
            }
        }

        // Back transformation to get eigenvectors of original matrix

        for (int j = nn - 1; j >= low; j--) {
            for (int i = low; i <= high; i++) {
                z = 0.0;
                for (int k = low; k <= Math.min(j, high); k++) {
                    z = z + V.get(i, k) * H.get(k, j);
                }
                V.set(i, j, z);
            }
        }
    }
96453
def _eig_complex_symmetric(M: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    
    if not np.allclose(M, M.transpose()):
        raise np.linalg.LinAlgError('Not a symmetric matrix')

    # The matrix of eigenvectors should be orthogonal.
    # But the standard 'eig' method will fail to return an orthogonal
    # eigenvector matrix when the eigenvalues are degenerate. However,
    # both the real and
    # imaginary part of M must be symmetric with the same orthogonal
    # matrix of eigenvectors. But either the real or imaginary part could
    # vanish. So we use a randomized algorithm where we diagonalize a
    # random linear combination of real and imaginary parts to find the
    # eigenvectors, taking advantage of the 'eigh' subroutine for
    # diagonalizing symmetric matrices.
    # This can fail if we're very unlucky with our random coefficient, so we
    # give the algorithm a few chances to succeed.

    # Empirically, never seems to fail on randomly sampled complex
    # symmetric 4x4 matrices.
    # If failure rate is less than 1 in a million, then 16 rounds
    # will have overall failure rate less than 1 in a googol.
    # However, cannot (yet) guarantee that there aren't special cases
    # which have much higher failure rates.

    # GEC 2018

    max_attempts = 16
    for _ in range(max_attempts):
        c = np.random.uniform(0, 1)
        matrix = c * M.real + (1-c) * M.imag
        _, eigvecs = np.linalg.eigh(matrix)
        eigvecs = np.array(eigvecs, dtype=complex)
        eigvals = np.diag(eigvecs.transpose() @ M @ eigvecs)

        # Finish if we got a correct answer.
        reconstructed = eigvecs @ np.diag(eigvals) @ eigvecs.transpose()
        if np.allclose(M, reconstructed):
            return eigvals, eigvecs

    # Should never happen. Hopefully.
    raise np.linalg.LinAlgError(
        'Cannot diagonalize complex symmetric matrix.')
1278575
def cublasZtrmm(handle, side, uplo, trans, diag, m, n, alpha, A, lda, B, ldb, C, ldc):
    
    
    status = _libcublas.cublasZtrmm_v2(handle, 
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       m, n, ctypes.byref(cuda.cuDoubleComplex(alpha.real,     
                                                                               alpha.imag)),
                                       int(A), lda, int(B), ldb, int(C), ldc)
    cublasCheckStatus(status)
1201810
def neighbors( open_list=[], closed_list=[] )
      i,j = free_cell      
      up = moved_matrix_with_parent(i, j, i-1, j, open_list, closed_list)
      down = moved_matrix_with_parent(i, j, i+1, j, open_list, closed_list)
      left = moved_matrix_with_parent(i, j, i, j-1, open_list, closed_list)
      right = moved_matrix_with_parent(i, j, i, j+1, open_list, closed_list)

      moved = []
      moved << up if !up.nil? 
      moved << down if !down.nil? 
      moved << left if !left.nil? 
      moved << right if !right.nil? 

      return moved
    end
721379
def cublasSsyr2k(handle, uplo, trans, n, k, alpha, A, lda, B, ldb, beta, C, ldc):
    
    
    status = _libcublas.cublasSsyr2k_v2(handle,
                                        _CUBLAS_FILL_MODE[uplo], 
                                        _CUBLAS_OP[trans], 
                                        n, k, ctypes.byref(ctypes.c_float(alpha)),
                                        int(A), lda, int(B), ldb, 
                                        ctypes.byref(ctypes.c_float(beta)), 
                                        int(C), ldc)
    cublasCheckStatus(status)
341621
def cublasSdgmm(handle, mode, m, n, A, lda, x, incx, C, ldc):
    

    status = _libcublas.cublasSdgmm(handle,
                                    _CUBLAS_SIDE[mode],
                                    m, n, 
                                    int(A), lda, 
                                    int(x), incx,
                                    int(C), ldc)
    cublasCheckStatus(status)
676362
private void tred2() {
    // This is derived from the Algol procedures tred2 by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.
    System.arraycopy(V[n - 1], 0, d, 0, n);

    // Householder reduction to tridiagonal form.
    for(int i = n - 1; i > 0; i--) {
      double[] V_i = V[i], V_im1 = V[i - 1];
      // Scale to avoid under/overflow.
      double scale = 0.;
      for(int k = 0; k < i; k++) {
        scale += abs(d[k]);
      }
      if(scale < Double.MIN_NORMAL) {
        e[i] = d[i - 1];
        for(int j = 0; j < i; j++) {
          d[j] = V_im1[j];
          V_i[j] = V[j][i] = 0.;
        }
        d[i] = 0;
        continue;
      }
      // Generate Householder vector.
      double h = 0.;
      for(int k = 0; k < i; k++) {
        double dk = d[k] /= scale;
        h += dk * dk;
      }
      {
        double f = d[i - 1], g = FastMath.sqrt(h);
        g = (f > 0) ? -g : g;
        e[i] = scale * g;
        h -= f * g;
        d[i - 1] = f - g;
        Arrays.fill(e, 0, i, 0.);
      }

      // Apply similarity transformation to remaining columns.
      for(int j = 0; j < i; j++) {
        final double[] Vj = V[j];
        double dj = Vj[i] = d[j], ej = e[j] + Vj[j] * dj;
        for(int k = j + 1; k <= i - 1; k++) {
          final double Vkj = V[k][j];
          ej += Vkj * d[k];
          e[k] += Vkj * dj;
        }
        e[j] = ej;
      }
      double sum = 0.;
      for(int j = 0; j < i; j++) {
        sum += (e[j] /= h) * d[j];
      }
      double hh = sum / (h + h);
      for(int j = 0; j < i; j++) {
        e[j] -= hh * d[j];
      }
      for(int j = 0; j < i; j++) {
        double dj = d[j], ej = e[j];
        for(int k = j; k <= i - 1; k++) {
          V[k][j] -= (dj * e[k] + ej * d[k]);
        }
        d[j] = V_im1[j];
        V_i[j] = 0.;
      }
      d[i] = h;
    }

    // Accumulate transformations.
    tred2AccumulateTransformations();
    System.arraycopy(V[n - 1], 0, d, 0, n);
    Arrays.fill(V[n - 1], 0.);
    V[n - 1][n - 1] = 1.;
    e[0] = 0.;
  }
929959
def cublasZsyrk(handle, uplo, trans, n, k, alpha, A, lda, beta, C, ldc):
    
    
    status = _libcublas.cublasZsyrk_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       n, k, ctypes.byref(cuda.cuDoubleComplex(alpha.real,    
                                                                               alpha.imag)),
                                       int(A), lda,
                                       ctypes.byref(cuda.cuDoubleComplex(beta.real,
                                                                         beta.imag)),
                                       int(C), ldc)
    cublasCheckStatus(status)
785859
private void tql2 (double[][] V) {

   //  This is derived from the Algol procedures tql2, by
   //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
   //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
   //  Fortran subroutine in EISPACK.
   
      for (int i = 1; i < n; i++) {
         e[i-1] = e[i];
      }
      e[n-1] = 0.0;
   
      double f = 0.0;
      double tst1 = 0.0;
      double eps = Math.pow(2.0,-52.0);
      for (int l = 0; l < n; l++) {

         // Find small subdiagonal element
   
         tst1 = Math.max(tst1,Math.abs(d[l]) + Math.abs(e[l]));
         int m = l;
         while (m < n) {
            if (Math.abs(e[m]) <= eps*tst1) {
               break;
            }
            m++;
         }
   
         // If m == l, d[l] is an eigenvalue,
         // otherwise, iterate.
   
         if (m > l) {
            int iter = 0;
            do {
               iter = iter + 1;  // (Could check iteration count here.)
   
               // Compute implicit shift
   
               double g = d[l];
               double p = (d[l+1] - g) / (2.0 * e[l]);
               double r = Math.hypot(p,1.0);
               if (p < 0) {
                  r = -r;
               }
               d[l] = e[l] / (p + r);
               d[l+1] = e[l] * (p + r);
               double dl1 = d[l+1];
               double h = g - d[l];
               for (int i = l+2; i < n; i++) {
                  d[i] -= h;
               }
               f = f + h;
   
               // Implicit QL transformation.
   
               p = d[m];
               double c = 1.0;
               double c2 = c;
               double c3 = c;
               double el1 = e[l+1];
               double s = 0.0;
               double s2 = 0.0;
               for (int i = m-1; i >= l; i--) {
                  c3 = c2;
                  c2 = c;
                  s2 = s;
                  g = c * e[i];
                  h = c * p;
                  r = Math.hypot(p,e[i]);
                  e[i+1] = s * r;
                  s = e[i] / r;
                  c = p / r;
                  p = c * d[i] - s * g;
                  d[i+1] = h + s * (c * g + s * d[i]);
   
                  // Accumulate transformation.
   
                  for (int k = 0; k < n; k++) {
                     h = V[k][i+1];
                     V[k][i+1] = s * V[k][i] + c * h;
                     V[k][i] = c * V[k][i] - s * h;
                  }
               }
               p = -s * s2 * c3 * el1 * e[l] / dl1;
               e[l] = s * p;
               d[l] = c * p;
   
               // Check for convergence.
   
            } while (Math.abs(e[l]) > eps*tst1);
         }
         d[l] = d[l] + f;
         e[l] = 0.0;
      }
     
      // Sort eigenvalues and corresponding vectors.
   
      for (int i = 0; i < n-1; i++) {
         int k = i;
         double p = d[i];
         for (int j = i+1; j < n; j++) {
            if (d[j] < p) {
               k = j;
               p = d[j];
            }
         }
         if (k != i) {
            d[k] = d[i];
            d[i] = p;
            for (int j = 0; j < n; j++) {
               p = V[j][i];
               V[j][i] = V[j][k];
               V[j][k] = p;
            }
         }
      }
   }
545209
private static void eltran(DenseMatrix A, DenseMatrix V, int[] perm) {
        int n = A.nrows();
        for (int mp = n - 2; mp > 0; mp--) {
            for (int k = mp + 1; k < n; k++) {
                V.set(k, mp, A.get(k, mp - 1));
            }
            int i = perm[mp];
            if (i != mp) {
                for (int j = mp; j < n; j++) {
                    V.set(mp, j, V.get(i, j));
                    V.set(i, j, 0.0);
                }
                V.set(i, mp, 1.0);
            }
        }
    }
1132385
public static void symmRankNMinus_L(int blockLength ,
                                        DSubmatrixD1 A , DSubmatrixD1 B )
    {
        int widthB = B.col1-B.col0;
        if( widthB > blockLength )
            throw new IllegalArgumentException("Width of B cannot be greater than the block length");

        int N = B.row1-B.row0;

        if( A.col1-A.col0 != N )
            throw new IllegalArgumentException("A does not have the expected number of columns based on B's height");
        if( A.row1-A.row0 != N )
            throw new IllegalArgumentException("A does not have the expected number of rows based on B's height");

        for( int i = B.row0; i < B.row1; i += blockLength ) {


            int heightB_i = Math.min(blockLength,B.row1-i);
            int indexB_i = i*B.original.numCols + heightB_i*B.col0;

            int rowA = i-B.row0+A.row0;
            int heightA = Math.min( blockLength , A.row1 - rowA);

            for( int j = B.row0; j <= i; j += blockLength ) {
                
                int widthB_j = Math.min(blockLength,B.row1-j);

                int indexA = rowA * A.original.numCols + (j-B.row0+A.col0)*heightA;
                int indexB_j = j*B.original.numCols + widthB_j*B.col0;

                if( i == j ) {
                    multTransBBlockMinus_L( B.original.data,A.original.data,
                            indexB_i,indexB_j,indexA,widthB,heightB_i,widthB_j);
                } else {
                    multTransBBlockMinus( B.original.data,A.original.data,
                            indexB_i,indexB_j,indexA,widthB,heightB_i,widthB_j);
                }
            }
        }
    }
314543
public static void update(DMatrixRMaj H ,
							  DMatrixRMaj s ,
							  DMatrixRMaj y ,
							  DMatrixRMaj tempV0, DMatrixRMaj tempV1) {
		double p = VectorVectorMult_DDRM.innerProd(y,s);
		if( p == 0 )
			return;

		p = 1.0/p;

		double sBs = VectorVectorMult_DDRM.innerProdA(s,H,s);
		if( sBs == 0 )
			return;

		CommonOps_DDRM.mult(H,s,tempV0);
		CommonOps_DDRM.multTransA(s,H,tempV1);

		VectorVectorMult_DDRM.rank1Update(-p, H , tempV0, y);
		VectorVectorMult_DDRM.rank1Update(-p, H , y, tempV1);
		VectorVectorMult_DDRM.rank1Update(p*(p*sBs+1), H , y, y);
	}
114136
private long[] refine(long[] invariants, boolean[] hydrogens) {

        int ord = g.length;

        InvariantRanker ranker = new InvariantRanker(ord);

        // current/next vertices, these only hold the vertices which are
        // equivalent
        int[] currVs = new int[ord];
        int[] nextVs = new int[ord];

        // fill with identity (also set number of non-unique)
        int nnu = ord;
        for (int i = 0; i < ord; i++)
            currVs[i] = i;

        long[] prev = invariants;
        long[] curr = Arrays.copyOf(invariants, ord);

        // initially all labels are 1, the input invariants are then used to
        // refine this coarse partition
        Arrays.fill(prev, 1L);

        // number of ranks
        int n = 0, m = 0;

        // storage of symmetry classes
        long[] symmetry = null;

        while (n < ord) {

            // refine the initial invariants using product of primes from
            // adjacent ranks
            while ((n = ranker.rank(currVs, nextVs, nnu, curr, prev)) > m && n < ord) {
                nnu = 0;
                for (int i = 0; i < ord && nextVs[i] >= 0; i++) {
                    int v = nextVs[i];
                    currVs[nnu++] = v;
                    curr[v] = hydrogens[v] ? prev[v] : primeProduct(g[v], prev, hydrogens);
                }
                m = n;
            }

            if (symmetry == null) {

                // After symmetry classes have been found without hydrogens we add
                // back in the hydrogens and assign ranks. We don't refine the
                // partition until the next time round the while loop to avoid
                // artificially splitting due to hydrogen representation, for example
                // the two hydrogens are equivalent in this SMILES for ethane '[H]CC'
                for (int i = 0; i < g.length; i++) {
                    if (hydrogens[i]) {
                        curr[i] = prev[g[i][0]];
                        hydrogens[i] = false;
                    }
                }
                n = ranker.rank(currVs, nextVs, nnu, curr, prev);
                symmetry = Arrays.copyOf(prev, ord);

                // Update the buffer of non-unique vertices as hydrogens next
                // to discrete heavy atoms are also discrete (and removed from
                // 'nextVs' during ranking.
                nnu = 0;
                for (int i = 0; i < ord && nextVs[i] >= 0; i++) {
                    currVs[nnu++] = nextVs[i];
                }
            }

            // partition is discrete or only symmetry classes are needed
            if (symOnly || n == ord) return symmetry;

            // artificially split the lowest cell, we perturb the value
            // of all vertices with equivalent rank to the lowest non-unique
            // vertex
            int lo = nextVs[0];
            for (int i = 1; i < ord && nextVs[i] >= 0 && prev[nextVs[i]] == prev[lo]; i++)
                prev[nextVs[i]]++;

            // could also swap but this is cleaner
            System.arraycopy(nextVs, 0, currVs, 0, nnu);
        }

        return symmetry;
    }
1233496
def matrix_chain_mult(M):
    
    opt, arg = matrix_mult_opt_order(M)
    return _apply_order(M, arg, 0, len(M)-1)
221722
private void svdImpute(double[][] raw, double[][] data) {
        SVD svd = Matrix.newInstance(data).svd();

        int d = data[0].length;

        for (int i = 0; i < raw.length; i++) {
            int missing = 0;
            for (int j = 0; j < d; j++) {
                if (Double.isNaN(raw[i][j])) {
                    missing++;
                } else {
                    data[i][j] = raw[i][j];
                }
            }

            if (missing == 0) {
                continue;
            }

            DenseMatrix A = Matrix.zeros(d - missing, k);
            double[] b = new double[d - missing];

            for (int j = 0, m = 0; j < d; j++) {
                if (!Double.isNaN(raw[i][j])) {
                    for (int l = 0; l < k; l++) {
                        A.set(m, l, svd.getV().get(j, l));
                    }
                    b[m++] = raw[i][j];
                }
            }

            double[] s = new double[k];
            QR qr = A.qr();
            qr.solve(b, s);

            for (int j = 0; j < d; j++) {
                if (Double.isNaN(raw[i][j])) {
                    data[i][j] = 0;
                    for (int l = 0; l < k; l++) {
                        data[i][j] += s[l] * svd.getV().get(j, l);
                    }
                }
            }
        }
    }