SLH algorithm for increasing orthogonality of a matrix.
433559
def _proximity_to_association(proximity):
    """"""
    # pylint:disable=invalid-name
    # I'm afraid that the short names here are just a function of the
    # mathematical nature of the code.

    # Special case: zero-size matrix
    if proximity.shape[0] == 0 or proximity.shape[1] == 0:
        return proximity

    # Take SVD of proximity matrix
    U, sing_vals, Vh = np.linalg.svd(proximity)

    # Count the number of non-zero singular values
    nsv = np.count_nonzero(sing_vals)

    # Trim U and Vh to match
    U, Vh = U[:, :nsv], Vh[:nsv, :]

    # Compute new matrix as if singular values were unity
    return U.dot(Vh)
1072756
private function tql2(): void
    {
        for ($i = 1; $i < $this->n; ++$i) {
            $this->e[$i - 1] = $this->e[$i];
        }

        $this->e[$this->n - 1] = 0.0;
        $f = 0.0;
        $tst1 = 0.0;
        $eps = 2.0 ** -52.0;

        for ($l = 0; $l < $this->n; ++$l) {
            // Find small subdiagonal element
            $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
            $m = $l;
            while ($m < $this->n) {
                if (abs($this->e[$m]) <= $eps * $tst1) {
                    break;
                }

                ++$m;
            }

            // If m == l, $this->d[l] is an eigenvalue,
            // otherwise, iterate.
            if ($m > $l) {
                do {
                    // Compute implicit shift
                    $g = $this->d[$l];
                    $p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]);
                    $r = hypot($p, 1.0);
                    if ($p < 0) {
                        $r *= -1;
                    }

                    $this->d[$l] = $this->e[$l] / ($p + $r);
                    $this->d[$l + 1] = $this->e[$l] * ($p + $r);
                    $dl1 = $this->d[$l + 1];
                    $h = $g - $this->d[$l];
                    for ($i = $l + 2; $i < $this->n; ++$i) {
                        $this->d[$i] -= $h;
                    }

                    $f += $h;
                    // Implicit QL transformation.
                    $p = $this->d[$m];
                    $c = 1.0;
                    $c2 = $c3 = $c;
                    $el1 = $this->e[$l + 1];
                    $s = $s2 = 0.0;
                    for ($i = $m - 1; $i >= $l; --$i) {
                        $c3 = $c2;
                        $c2 = $c;
                        $s2 = $s;
                        $g = $c * $this->e[$i];
                        $h = $c * $p;
                        $r = hypot($p, $this->e[$i]);
                        $this->e[$i + 1] = $s * $r;
                        $s = $this->e[$i] / $r;
                        $c = $p / $r;
                        $p = $c * $this->d[$i] - $s * $g;
                        $this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]);
                        // Accumulate transformation.
                        for ($k = 0; $k < $this->n; ++$k) {
                            $h = $this->V[$k][$i + 1];
                            $this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h;
                            $this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
                        }
                    }

                    $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
                    $this->e[$l] = $s * $p;
                    $this->d[$l] = $c * $p;
                    // Check for convergence.
                } while (abs($this->e[$l]) > $eps * $tst1);
            }

            $this->d[$l] += $f;
            $this->e[$l] = 0.0;
        }

        // Sort eigenvalues and corresponding vectors.
        for ($i = 0; $i < $this->n - 1; ++$i) {
            $k = $i;
            $p = $this->d[$i];
            for ($j = $i + 1; $j < $this->n; ++$j) {
                if ($this->d[$j] < $p) {
                    $k = $j;
                    $p = $this->d[$j];
                }
            }

            if ($k != $i) {
                $this->d[$k] = $this->d[$i];
                $this->d[$i] = $p;
                for ($j = 0; $j < $this->n; ++$j) {
                    $p = $this->V[$j][$i];
                    $this->V[$j][$i] = $this->V[$j][$k];
                    $this->V[$j][$k] = $p;
                }
            }
        }
    }
1327131
private void _symmetricTred2 ()
  {
    // 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.
    for (int j = 0; j < m_nDim; j++)
    {
      m_aEVd[j] = m_aEigenVector[m_nDim - 1][j];
    }

    // Householder reduction to tridiagonal form.

    for (int i = m_nDim - 1; i > 0; i--)
    {
      // Scale to avoid under/overflow.
      double scale = 0.0;
      double h = 0.0;
      for (int k = 0; k < i; k++)
        scale += MathHelper.abs (m_aEVd[k]);
      if (scale == 0.0)
      {
        m_aEVe[i] = m_aEVd[i - 1];
        for (int j = 0; j < i; j++)
        {
          m_aEVd[j] = m_aEigenVector[i - 1][j];
          m_aEigenVector[i][j] = 0.0;
          m_aEigenVector[j][i] = 0.0;
        }
      }
      else
      {
        // Generate Householder vector.

        for (int k = 0; k < i; k++)
        {
          m_aEVd[k] /= scale;
          h += m_aEVd[k] * m_aEVd[k];
        }
        double f = m_aEVd[i - 1];
        double g = Math.sqrt (h);
        if (f > 0)
        {
          g = -g;
        }
        m_aEVe[i] = scale * g;
        h -= f * g;
        m_aEVd[i - 1] = f - g;
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] = 0.0;
        }

        // Apply similarity transformation to remaining columns.
        for (int j = 0; j < i; j++)
        {
          f = m_aEVd[j];
          m_aEigenVector[j][i] = f;
          g = m_aEVe[j] + m_aEigenVector[j][j] * f;
          for (int k = j + 1; k <= i - 1; k++)
          {
            g += m_aEigenVector[k][j] * m_aEVd[k];
            m_aEVe[k] += m_aEigenVector[k][j] * f;
          }
          m_aEVe[j] = g;
        }
        f = 0.0;
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] /= h;
          f += m_aEVe[j] * m_aEVd[j];
        }
        final double hh = f / (h + h);
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] -= hh * m_aEVd[j];
        }
        for (int j = 0; j < i; j++)
        {
          f = m_aEVd[j];
          g = m_aEVe[j];
          for (int k = j; k <= i - 1; k++)
          {
            m_aEigenVector[k][j] -= (f * m_aEVe[k] + g * m_aEVd[k]);
          }
          m_aEVd[j] = m_aEigenVector[i - 1][j];
          m_aEigenVector[i][j] = 0.0;
        }
      }
      m_aEVd[i] = h;
    }

    // Accumulate transformations.
    for (int i = 0; i < m_nDim - 1; i++)
    {
      m_aEigenVector[m_nDim - 1][i] = m_aEigenVector[i][i];
      m_aEigenVector[i][i] = 1.0;
      final double h = m_aEVd[i + 1];
      if (h != 0.0)
      {
        for (int k = 0; k <= i; k++)
        {
          m_aEVd[k] = m_aEigenVector[k][i + 1] / h;
        }
        for (int j = 0; j <= i; j++)
        {
          double g = 0.0;
          for (int k = 0; k <= i; k++)
          {
            g += m_aEigenVector[k][i + 1] * m_aEigenVector[k][j];
          }
          for (int k = 0; k <= i; k++)
          {
            m_aEigenVector[k][j] -= g * m_aEVd[k];
          }
        }
      }
      for (int k = 0; k <= i; k++)
      {
        m_aEigenVector[k][i + 1] = 0.0;
      }
    }
    for (int j = 0; j < m_nDim; j++)
    {
      m_aEVd[j] = m_aEigenVector[m_nDim - 1][j];
      m_aEigenVector[m_nDim - 1][j] = 0.0;
    }
    m_aEigenVector[m_nDim - 1][m_nDim - 1] = 1.0;
    m_aEVe[0] = 0.0;
  }
766162
private void tql2() {
    // 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.

    System.arraycopy(e, 1, e, 0, n - 1);
    e[n - 1] = 0.;

    double f = 0., tst1 = 0.;
    for(int l = 0; l < n; l++) {
      // Find small subdiagonal element
      tst1 = max(tst1, abs(d[l]) + abs(e[l]));
      int m = l;
      for(; m < n; ++m) {
        if(abs(e[m]) <= 0x1P-52 * tst1) {
          break;
        }
      }

      // If m == l, d[l] is an eigenvalue,
      // otherwise, iterate.
      if(m > l) {
        do {
          // Compute implicit shift
          f += tql2ComputeImplicitShift(l);
          // Implicit QL transformation.
          tql2ImplicitQL(l, m, d[l + 1]);
          // Check for convergence.
          // TODO: Iteration limit?
        }
        while(abs(e[l]) > 0x1P-52 * tst1);
      }
      d[l] += f;
      e[l] = 0.;
    }
  }
320729
function svdPack() {
                var result = lapack$1.sgesvd('A', 'A', this.elements);

                return {
                    U: $M(result.U),
                    S: $M(result.S).column(1).toDiagonalMatrix(),
            	V: $M(result.VT).transpose()
                };
            }
977204
public static void inplaceAdjustCameraMatrix( double sx , double sy , double tx , double ty , DMatrixRMaj P ) {
		// multiply each column one at a time. Because of the zeros everything is decoupled
		for (int col = 0; col < 4; col++) {
			int row0 = col;
			int row1 = 4 + row0;
			int row2 = 4 + row1;
			P.data[row0] = P.data[row0]*sx + tx*P.data[row2];
			P.data[row1] = P.data[row1]*sy + ty*P.data[row2];
		}
	}
1304808
private SimpleMatrix computeAffineH( List<AssociatedPair> observations ,
										 DMatrixRMaj H , DMatrixRMaj Hzero ) {

		SimpleMatrix A = new SimpleMatrix(observations.size(),3);
		SimpleMatrix b = new SimpleMatrix(A.numRows(),1);

		Point2D_F64 c = new Point2D_F64();
		Point2D_F64 k = new Point2D_F64();

		for( int i = 0; i < observations.size(); i++ ) {
			AssociatedPair a = observations.get(i);

			GeometryMath_F64.mult(Hzero, a.p1, k);
			GeometryMath_F64.mult(H,a.p2,c);

			A.setRow(i,0,k.x,k.y,1);
			b.set(i,0,c.x);
		}

		SimpleMatrix x = A.solve(b);

		SimpleMatrix Ha = SimpleMatrix.identity(3);
		Ha.setRow(0,0,x.getDDRM().data);

		return Ha;
	}
827941
public static LinearSolverDense<DMatrixRMaj> general(int numRows , int numCols ) {
        if( numRows == numCols )
            return linear(numRows);
        else
            return leastSquares(numRows,numCols);
    }
898374
def _coupling_matrix(self, lmax, nwin=None, weights=None):
        """"""
        if nwin is None:
            nwin = self.nwin

        if weights is None:
            weights = self.weights

        if weights is None:
            return _shtools.SHMTCouplingMatrix(lmax, self.tapers**2, k=nwin)
        else:
            return _shtools.SHMTCouplingMatrix(lmax, self.tapers**2, k=nwin,
                                               taper_wt=self.weights)
1029768
private void _nonsymetricHqr2 ()
  {
    // 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

    final int nn = m_nDim;
    int n = nn - 1;
    final int low = 0;
    final int high = nn - 1;
    double exshift = 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)
      {
        m_aEVd[i] = m_aHessenBerg[i][i];
        m_aEVe[i] = 0.0;
      }
      for (int j = Math.max (i - 1, 0); j < nn; j++)
      {
        norm += MathHelper.abs (m_aHessenBerg[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 = MathHelper.abs (m_aHessenBerg[l - 1][l - 1]) + MathHelper.abs (m_aHessenBerg[l][l]);
        if (s == 0.0)
          s = norm;
        if (MathHelper.abs (m_aHessenBerg[l][l - 1]) < EPSILON * s)
          break;
        l--;
      }

      // Check for convergence
      // One root found
      if (l == n)
      {
        m_aHessenBerg[n][n] = m_aHessenBerg[n][n] + exshift;
        m_aEVd[n] = m_aHessenBerg[n][n];
        m_aEVe[n] = 0.0;
        n--;
        iter = 0;
        // Two roots found
      }
      else
        if (l == n - 1)
        {
          w = m_aHessenBerg[n][n - 1] * m_aHessenBerg[n - 1][n];
          p = (m_aHessenBerg[n - 1][n - 1] - m_aHessenBerg[n][n]) / 2.0;
          q = p * p + w;
          z = Math.sqrt (MathHelper.abs (q));
          m_aHessenBerg[n][n] = m_aHessenBerg[n][n] + exshift;
          m_aHessenBerg[n - 1][n - 1] = m_aHessenBerg[n - 1][n - 1] + exshift;
          x = m_aHessenBerg[n][n];

          // Real pair

          if (q >= 0)
          {
            if (p >= 0)
              z = p + z;
            else
              z = p - z;
            m_aEVd[n - 1] = x + z;
            m_aEVd[n] = m_aEVd[n - 1];
            if (z != 0.0)
              m_aEVd[n] = x - w / z;
            m_aEVe[n - 1] = 0.0;
            m_aEVe[n] = 0.0;
            x = m_aHessenBerg[n][n - 1];
            s = MathHelper.abs (x) + MathHelper.abs (z);
            p = x / s;
            q = z / s;
            r = Math.sqrt (p * p + q * q);
            p /= r;
            q /= r;

            // Row modification
            for (int j = n - 1; j < nn; j++)
            {
              z = m_aHessenBerg[n - 1][j];
              m_aHessenBerg[n - 1][j] = q * z + p * m_aHessenBerg[n][j];
              m_aHessenBerg[n][j] = q * m_aHessenBerg[n][j] - p * z;
            }

            // Column modification
            for (int i = 0; i <= n; i++)
            {
              z = m_aHessenBerg[i][n - 1];
              m_aHessenBerg[i][n - 1] = q * z + p * m_aHessenBerg[i][n];
              m_aHessenBerg[i][n] = q * m_aHessenBerg[i][n] - p * z;
            }

            // Accumulate transformations
            for (int i = low; i <= high; i++)
            {
              z = m_aEigenVector[i][n - 1];
              m_aEigenVector[i][n - 1] = q * z + p * m_aEigenVector[i][n];
              m_aEigenVector[i][n] = q * m_aEigenVector[i][n] - p * z;
            }

            // Complex pair
          }
          else
          {
            m_aEVd[n - 1] = x + p;
            m_aEVd[n] = x + p;
            m_aEVe[n - 1] = z;
            m_aEVe[n] = -z;
          }
          n -= 2;
          iter = 0;

          // No convergence yet
        }
        else
        {
          // Form shift
          x = m_aHessenBerg[n][n];
          y = 0.0;
          w = 0.0;
          if (l < n)
          {
            y = m_aHessenBerg[n - 1][n - 1];
            w = m_aHessenBerg[n][n - 1] * m_aHessenBerg[n - 1][n];
          }

          // Wilkinson's original ad hoc shift
          if (iter == 10)
          {
            exshift += x;
            for (int i = low; i <= n; i++)
            {
              m_aHessenBerg[i][i] -= x;
            }
            s = MathHelper.abs (m_aHessenBerg[n][n - 1]) + MathHelper.abs (m_aHessenBerg[n - 1][n - 2]);
            x = 0.75 * s;
            y = x;
            w = -0.4375 * s * s;
          }

          // 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++)
              {
                m_aHessenBerg[i][i] -= s;
              }
              exshift += s;
              x = 0.964;
              y = x;
              w = x;
            }
          }

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

          // Look for two consecutive small sub-diagonal elements
          int m = n - 2;
          while (m >= l)
          {
            z = m_aHessenBerg[m][m];
            r = x - z;
            s = y - z;
            p = (r * s - w) / m_aHessenBerg[m + 1][m] + m_aHessenBerg[m][m + 1];
            q = m_aHessenBerg[m + 1][m + 1] - z - r - s;
            r = m_aHessenBerg[m + 2][m + 1];
            s = MathHelper.abs (p) + MathHelper.abs (q) + MathHelper.abs (r);
            p /= s;
            q /= s;
            r /= s;
            if (m == l)
              break;
            if (MathHelper.abs (m_aHessenBerg[m][m - 1]) *
                (MathHelper.abs (q) + MathHelper.abs (r)) < EPSILON *
                                                            (MathHelper.abs (p) *
                                                             (MathHelper.abs (m_aHessenBerg[m - 1][m - 1]) +
                                                              MathHelper.abs (z) +
                                                              MathHelper.abs (m_aHessenBerg[m + 1][m + 1]))))
            {
              break;
            }
            m--;
          }

          for (int i = m + 2; i <= n; i++)
          {
            m_aHessenBerg[i][i - 2] = 0.0;
            if (i > m + 2)
            {
              m_aHessenBerg[i][i - 3] = 0.0;
            }
          }

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

          for (int k = m; k <= n - 1; k++)
          {
            final boolean notlast = (k != n - 1);
            if (k != m)
            {
              p = m_aHessenBerg[k][k - 1];
              q = m_aHessenBerg[k + 1][k - 1];
              r = (notlast ? m_aHessenBerg[k + 2][k - 1] : 0.0);
              x = MathHelper.abs (p) + MathHelper.abs (q) + MathHelper.abs (r);
              if (x == 0.0)
              {
                continue;
              }
              p /= x;
              q /= x;
              r /= x;
            }

            s = Math.sqrt (p * p + q * q + r * r);
            if (p < 0)
            {
              s = -s;
            }
            if (s != 0)
            {
              if (k != m)
              {
                m_aHessenBerg[k][k - 1] = -s * x;
              }
              else
                if (l != m)
                {
                  m_aHessenBerg[k][k - 1] = -m_aHessenBerg[k][k - 1];
                }
              p += s;
              x = p / s;
              y = q / s;
              z = r / s;
              q /= p;
              r /= p;

              // Row modification

              for (int j = k; j < nn; j++)
              {
                p = m_aHessenBerg[k][j] + q * m_aHessenBerg[k + 1][j];
                if (notlast)
                {
                  p += r * m_aHessenBerg[k + 2][j];
                  m_aHessenBerg[k + 2][j] = m_aHessenBerg[k + 2][j] - p * z;
                }
                m_aHessenBerg[k][j] = m_aHessenBerg[k][j] - p * x;
                m_aHessenBerg[k + 1][j] = m_aHessenBerg[k + 1][j] - p * y;
              }

              // Column modification

              for (int i = 0; i <= Math.min (n, k + 3); i++)
              {
                p = x * m_aHessenBerg[i][k] + y * m_aHessenBerg[i][k + 1];
                if (notlast)
                {
                  p += z * m_aHessenBerg[i][k + 2];
                  m_aHessenBerg[i][k + 2] = m_aHessenBerg[i][k + 2] - p * r;
                }
                m_aHessenBerg[i][k] = m_aHessenBerg[i][k] - p;
                m_aHessenBerg[i][k + 1] = m_aHessenBerg[i][k + 1] - p * q;
              }

              // Accumulate transformations

              for (int i = low; i <= high; i++)
              {
                p = x * m_aEigenVector[i][k] + y * m_aEigenVector[i][k + 1];
                if (notlast)
                {
                  p += z * m_aEigenVector[i][k + 2];
                  m_aEigenVector[i][k + 2] = m_aEigenVector[i][k + 2] - p * r;
                }
                m_aEigenVector[i][k] = m_aEigenVector[i][k] - p;
                m_aEigenVector[i][k + 1] = m_aEigenVector[i][k + 1] - 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 = m_aEVd[n];
      q = m_aEVe[n];

      // Real vector

      if (q == 0)
      {
        int l = n;
        m_aHessenBerg[n][n] = 1.0;
        for (int i = n - 1; i >= 0; i--)
        {
          w = m_aHessenBerg[i][i] - p;
          r = 0.0;
          for (int j = l; j <= n; j++)
          {
            r += m_aHessenBerg[i][j] * m_aHessenBerg[j][n];
          }
          if (m_aEVe[i] < 0.0)
          {
            z = w;
            s = r;
          }
          else
          {
            l = i;
            if (m_aEVe[i] == 0.0)
            {
              if (w != 0.0)
              {
                m_aHessenBerg[i][n] = -r / w;
              }
              else
              {
                m_aHessenBerg[i][n] = -r / (EPSILON * norm);
              }

              // Solve real equations

            }
            else
            {
              x = m_aHessenBerg[i][i + 1];
              y = m_aHessenBerg[i + 1][i];
              q = (m_aEVd[i] - p) * (m_aEVd[i] - p) + m_aEVe[i] * m_aEVe[i];
              t = (x * s - z * r) / q;
              m_aHessenBerg[i][n] = t;
              if (MathHelper.abs (x) > MathHelper.abs (z))
              {
                m_aHessenBerg[i + 1][n] = (-r - w * t) / x;
              }
              else
              {
                m_aHessenBerg[i + 1][n] = (-s - y * t) / z;
              }
            }

            // Overflow control

            t = MathHelper.abs (m_aHessenBerg[i][n]);
            if ((EPSILON * t) * t > 1)
            {
              for (int j = i; j <= n; j++)
              {
                m_aHessenBerg[j][n] = m_aHessenBerg[j][n] / t;
              }
            }
          }
        }

        // Complex vector

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

          // Last vector component imaginary so matrix is triangular

          if (MathHelper.abs (m_aHessenBerg[n][n - 1]) > MathHelper.abs (m_aHessenBerg[n - 1][n]))
          {
            m_aHessenBerg[n - 1][n - 1] = q / m_aHessenBerg[n][n - 1];
            m_aHessenBerg[n - 1][n] = -(m_aHessenBerg[n][n] - p) / m_aHessenBerg[n][n - 1];
          }
          else
          {
            _cdiv (0.0, -m_aHessenBerg[n - 1][n], m_aHessenBerg[n - 1][n - 1] - p, q);
            m_aHessenBerg[n - 1][n - 1] = m_dCdivr;
            m_aHessenBerg[n - 1][n] = m_dCdivi;
          }
          m_aHessenBerg[n][n - 1] = 0.0;
          m_aHessenBerg[n][n] = 1.0;
          for (int i = n - 2; i >= 0; i--)
          {
            double ra = 0.0;
            double sa = 0.0;
            double vr;
            double vi;
            for (int j = l; j <= n; j++)
            {
              ra += m_aHessenBerg[i][j] * m_aHessenBerg[j][n - 1];
              sa += m_aHessenBerg[i][j] * m_aHessenBerg[j][n];
            }
            w = m_aHessenBerg[i][i] - p;

            if (m_aEVe[i] < 0.0)
            {
              z = w;
              r = ra;
              s = sa;
            }
            else
            {
              l = i;
              if (m_aEVe[i] == 0)
              {
                _cdiv (-ra, -sa, w, q);
                m_aHessenBerg[i][n - 1] = m_dCdivr;
                m_aHessenBerg[i][n] = m_dCdivi;
              }
              else
              {

                // Solve complex equations

                x = m_aHessenBerg[i][i + 1];
                y = m_aHessenBerg[i + 1][i];
                vr = (m_aEVd[i] - p) * (m_aEVd[i] - p) + m_aEVe[i] * m_aEVe[i] - q * q;
                vi = (m_aEVd[i] - p) * 2.0 * q;
                if (vr == 0.0 && vi == 0.0)
                {
                  vr = EPSILON *
                       norm *
                       (MathHelper.abs (w) +
                        MathHelper.abs (q) +
                        MathHelper.abs (x) +
                        MathHelper.abs (y) +
                        MathHelper.abs (z));
                }
                _cdiv (x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
                m_aHessenBerg[i][n - 1] = m_dCdivr;
                m_aHessenBerg[i][n] = m_dCdivi;
                if (MathHelper.abs (x) > (MathHelper.abs (z) + MathHelper.abs (q)))
                {
                  m_aHessenBerg[i + 1][n - 1] = (-ra - w * m_aHessenBerg[i][n - 1] + q * m_aHessenBerg[i][n]) / x;
                  m_aHessenBerg[i + 1][n] = (-sa - w * m_aHessenBerg[i][n] - q * m_aHessenBerg[i][n - 1]) / x;
                }
                else
                {
                  _cdiv (-r - y * m_aHessenBerg[i][n - 1], -s - y * m_aHessenBerg[i][n], z, q);
                  m_aHessenBerg[i + 1][n - 1] = m_dCdivr;
                  m_aHessenBerg[i + 1][n] = m_dCdivi;
                }
              }

              // Overflow control

              t = Math.max (MathHelper.abs (m_aHessenBerg[i][n - 1]), MathHelper.abs (m_aHessenBerg[i][n]));
              if ((EPSILON * t) * t > 1)
              {
                for (int j = i; j <= n; j++)
                {
                  m_aHessenBerg[j][n - 1] = m_aHessenBerg[j][n - 1] / t;
                  m_aHessenBerg[j][n] = m_aHessenBerg[j][n] / t;
                }
              }
            }
          }
        }
    }

    // Vectors of isolated roots

    for (int i = 0; i < nn; i++)
    {
      if (i < low || i > high)
      {
        for (int j = i; j < nn; j++)
        {
          m_aEigenVector[i][j] = m_aHessenBerg[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 += m_aEigenVector[i][k] * m_aHessenBerg[k][j];
        }
        m_aEigenVector[i][j] = z;
      }
    }
  }
124553
def cublasCsyr(handle, uplo, n, alpha, x, incx, A, lda):
    

    status = _libcublas.cublasCsyr_v2(handle,
                                        _CUBLAS_FILL_MODE[uplo], n, 
                                        ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                        alpha.imag)),
                                        int(x), incx, int(A), lda)
    cublasCheckStatus(status)
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)
1092191
private function hqr2()
    {
        //  Initialize
        $nn = $this->n;
        $n  = $nn - 1;
        $low = 0;
        $high = $nn - 1;
        $eps = pow(2.0, -52.0);
        $exshift = 0.0;
        $p = $q = $r = $s = $z = 0;
        // Store roots isolated by balanc and compute matrix norm
        $norm = 0.0;

        for ($i = 0; $i < $nn; ++$i) {
            if (($i < $low) or ($i > $high)) {
                $this->d[$i] = $this->H[$i][$i];
                $this->e[$i] = 0.0;
            }
            for ($j = max($i-1, 0); $j < $nn; ++$j) {
                $norm = $norm + abs($this->H[$i][$j]);
            }
        }

        // Outer loop over eigenvalue index
        $iter = 0;
        while ($n >= $low) {
            // Look for single small sub-diagonal element
            $l = $n;
            while ($l > $low) {
                $s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]);
                if ($s == 0.0) {
                    $s = $norm;
                }
                if (abs($this->H[$l][$l-1]) < $eps * $s) {
                    break;
                }
                --$l;
            }
            // Check for convergence
            // One root found
            if ($l == $n) {
                $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
                $this->d[$n] = $this->H[$n][$n];
                $this->e[$n] = 0.0;
                --$n;
                $iter = 0;
            // Two roots found
            } elseif ($l == $n-1) {
                $w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
                $p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0;
                $q = $p * $p + $w;
                $z = sqrt(abs($q));
                $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
                $this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift;
                $x = $this->H[$n][$n];
                // Real pair
                if ($q >= 0) {
                    if ($p >= 0) {
                        $z = $p + $z;
                    } else {
                        $z = $p - $z;
                    }
                    $this->d[$n-1] = $x + $z;
                    $this->d[$n] = $this->d[$n-1];
                    if ($z != 0.0) {
                        $this->d[$n] = $x - $w / $z;
                    }
                    $this->e[$n-1] = 0.0;
                    $this->e[$n] = 0.0;
                    $x = $this->H[$n][$n-1];
                    $s = abs($x) + abs($z);
                    $p = $x / $s;
                    $q = $z / $s;
                    $r = sqrt($p * $p + $q * $q);
                    $p = $p / $r;
                    $q = $q / $r;
                    // Row modification
                    for ($j = $n-1; $j < $nn; ++$j) {
                        $z = $this->H[$n-1][$j];
                        $this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j];
                        $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
                    }
                    // Column modification
                    for ($i = 0; $i <= $n; ++$i) {
                        $z = $this->H[$i][$n-1];
                        $this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n];
                        $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
                    }
                    // Accumulate transformations
                    for ($i = $low; $i <= $high; ++$i) {
                        $z = $this->V[$i][$n-1];
                        $this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n];
                        $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
                    }
                // Complex pair
                } else {
                    $this->d[$n-1] = $x + $p;
                    $this->d[$n] = $x + $p;
                    $this->e[$n-1] = $z;
                    $this->e[$n] = -$z;
                }
                $n = $n - 2;
                $iter = 0;
            // No convergence yet
            } else {
                // Form shift
                $x = $this->H[$n][$n];
                $y = 0.0;
                $w = 0.0;
                if ($l < $n) {
                    $y = $this->H[$n-1][$n-1];
                    $w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
                }
                // Wilkinson's original ad hoc shift
                if ($iter == 10) {
                    $exshift += $x;
                    for ($i = $low; $i <= $n; ++$i) {
                        $this->H[$i][$i] -= $x;
                    }
                    $s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]);
                    $x = $y = 0.75 * $s;
                    $w = -0.4375 * $s * $s;
                }
                // MATLAB's new ad hoc shift
                if ($iter == 30) {
                    $s = ($y - $x) / 2.0;
                    $s = $s * $s + $w;
                    if ($s > 0) {
                        $s = sqrt($s);
                        if ($y < $x) {
                            $s = -$s;
                        }
                        $s = $x - $w / (($y - $x) / 2.0 + $s);
                        for ($i = $low; $i <= $n; ++$i) {
                            $this->H[$i][$i] -= $s;
                        }
                        $exshift += $s;
                        $x = $y = $w = 0.964;
                    }
                }
                // Could check iteration count here.
                $iter = $iter + 1;
                // Look for two consecutive small sub-diagonal elements
                $m = $n - 2;
                while ($m >= $l) {
                    $z = $this->H[$m][$m];
                    $r = $x - $z;
                    $s = $y - $z;
                    $p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
                    $q = $this->H[$m+1][$m+1] - $z - $r - $s;
                    $r = $this->H[$m+2][$m+1];
                    $s = abs($p) + abs($q) + abs($r);
                    $p = $p / $s;
                    $q = $q / $s;
                    $r = $r / $s;
                    if ($m == $l) {
                        break;
                    }
                    if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) <
                        $eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
                        break;
                    }
                    --$m;
                }
                for ($i = $m + 2; $i <= $n; ++$i) {
                    $this->H[$i][$i-2] = 0.0;
                    if ($i > $m+2) {
                        $this->H[$i][$i-3] = 0.0;
                    }
                }
                // Double QR step involving rows l:n and columns m:n
                for ($k = $m; $k <= $n-1; ++$k) {
                    $notlast = ($k != $n-1);
                    if ($k != $m) {
                        $p = $this->H[$k][$k-1];
                        $q = $this->H[$k+1][$k-1];
                        $r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
                        $x = abs($p) + abs($q) + abs($r);
                        if ($x != 0.0) {
                            $p = $p / $x;
                            $q = $q / $x;
                            $r = $r / $x;
                        }
                    }
                    if ($x == 0.0) {
                        break;
                    }
                    $s = sqrt($p * $p + $q * $q + $r * $r);
                    if ($p < 0) {
                        $s = -$s;
                    }
                    if ($s != 0) {
                        if ($k != $m) {
                            $this->H[$k][$k-1] = -$s * $x;
                        } elseif ($l != $m) {
                            $this->H[$k][$k-1] = -$this->H[$k][$k-1];
                        }
                        $p = $p + $s;
                        $x = $p / $s;
                        $y = $q / $s;
                        $z = $r / $s;
                        $q = $q / $p;
                        $r = $r / $p;
                        // Row modification
                        for ($j = $k; $j < $nn; ++$j) {
                            $p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
                            if ($notlast) {
                                $p = $p + $r * $this->H[$k+2][$j];
                                $this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
                            }
                            $this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
                            $this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y;
                        }
                        // Column modification
                        for ($i = 0; $i <= min($n, $k+3); ++$i) {
                            $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1];
                            if ($notlast) {
                                $p = $p + $z * $this->H[$i][$k+2];
                                $this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r;
                            }
                            $this->H[$i][$k] = $this->H[$i][$k] - $p;
                            $this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
                        }
                        // Accumulate transformations
                        for ($i = $low; $i <= $high; ++$i) {
                            $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1];
                            if ($notlast) {
                                $p = $p + $z * $this->V[$i][$k+2];
                                $this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r;
                            }
                            $this->V[$i][$k] = $this->V[$i][$k] - $p;
                            $this->V[$i][$k+1] = $this->V[$i][$k+1] - $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 = $this->d[$n];
            $q = $this->e[$n];
            // Real vector
            if ($q == 0) {
                $l = $n;
                $this->H[$n][$n] = 1.0;
                for ($i = $n-1; $i >= 0; --$i) {
                    $w = $this->H[$i][$i] - $p;
                    $r = 0.0;
                    for ($j = $l; $j <= $n; ++$j) {
                        $r = $r + $this->H[$i][$j] * $this->H[$j][$n];
                    }
                    if ($this->e[$i] < 0.0) {
                        $z = $w;
                        $s = $r;
                    } else {
                        $l = $i;
                        if ($this->e[$i] == 0.0) {
                            if ($w != 0.0) {
                                $this->H[$i][$n] = -$r / $w;
                            } else {
                                $this->H[$i][$n] = -$r / ($eps * $norm);
                            }
                        // Solve real equations
                        } else {
                            $x = $this->H[$i][$i+1];
                            $y = $this->H[$i+1][$i];
                            $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
                            $t = ($x * $s - $z * $r) / $q;
                            $this->H[$i][$n] = $t;
                            if (abs($x) > abs($z)) {
                                $this->H[$i+1][$n] = (-$r - $w * $t) / $x;
                            } else {
                                $this->H[$i+1][$n] = (-$s - $y * $t) / $z;
                            }
                        }
                        // Overflow control
                        $t = abs($this->H[$i][$n]);
                        if (($eps * $t) * $t > 1) {
                            for ($j = $i; $j <= $n; ++$j) {
                                $this->H[$j][$n] = $this->H[$j][$n] / $t;
                            }
                        }
                    }
                }
            // Complex vector
            } elseif ($q < 0) {
                $l = $n-1;
                // Last vector component imaginary so matrix is triangular
                if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) {
                    $this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1];
                    $this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1];
                } else {
                    $this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q);
                    $this->H[$n-1][$n-1] = $this->cdivr;
                    $this->H[$n-1][$n]   = $this->cdivi;
                }
                $this->H[$n][$n-1] = 0.0;
                $this->H[$n][$n] = 1.0;
                for ($i = $n-2; $i >= 0; --$i) {
                    // double ra,sa,vr,vi;
                    $ra = 0.0;
                    $sa = 0.0;
                    for ($j = $l; $j <= $n; ++$j) {
                        $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1];
                        $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
                    }
                    $w = $this->H[$i][$i] - $p;
                    if ($this->e[$i] < 0.0) {
                        $z = $w;
                        $r = $ra;
                        $s = $sa;
                    } else {
                        $l = $i;
                        if ($this->e[$i] == 0) {
                            $this->cdiv(-$ra, -$sa, $w, $q);
                            $this->H[$i][$n-1] = $this->cdivr;
                            $this->H[$i][$n]   = $this->cdivi;
                        } else {
                            // Solve complex equations
                            $x = $this->H[$i][$i+1];
                            $y = $this->H[$i+1][$i];
                            $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
                            $vi = ($this->d[$i] - $p) * 2.0 * $q;
                            if ($vr == 0.0 & $vi == 0.0) {
                                $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
                            }
                            $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
                            $this->H[$i][$n-1] = $this->cdivr;
                            $this->H[$i][$n]   = $this->cdivi;
                            if (abs($x) > (abs($z) + abs($q))) {
                                $this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x;
                                $this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x;
                            } else {
                                $this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q);
                                $this->H[$i+1][$n-1] = $this->cdivr;
                                $this->H[$i+1][$n]   = $this->cdivi;
                            }
                        }
                        // Overflow control
                        $t = max(abs($this->H[$i][$n-1]), abs($this->H[$i][$n]));
                        if (($eps * $t) * $t > 1) {
                            for ($j = $i; $j <= $n; ++$j) {
                                $this->H[$j][$n-1] = $this->H[$j][$n-1] / $t;
                                $this->H[$j][$n]   = $this->H[$j][$n] / $t;
                            }
                        }
                    } // end else
                } // end for
            } // end else for complex case
        } // end for

        // Vectors of isolated roots
        for ($i = 0; $i < $nn; ++$i) {
            if ($i < $low | $i > $high) {
                for ($j = $i; $j < $nn; ++$j) {
                    $this->V[$i][$j] = $this->H[$i][$j];
                }
            }
        }

        // Back transformation to get eigenvectors of original matrix
        for ($j = $nn-1; $j >= $low; --$j) {
            for ($i = $low; $i <= $high; ++$i) {
                $z = 0.0;
                for ($k = $low; $k <= min($j, $high); ++$k) {
                    $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
                }
                $this->V[$i][$j] = $z;
            }
        }
    }
570764
protected static void adjustUpLeft(byte[][] matrix, byte[][] tempMatrix, int maxX, int maxY) {
		int currentMinX = Integer.MAX_VALUE, currentMinY = Integer.MAX_VALUE;
		for(int i = 0; i < maxX; i++) {
			for(int j = 0; j < maxY; j++) {
				if(matrix[i][j] == 1) {
					if(i < currentMinX) {
						currentMinX = i;
					}
					if(j < currentMinY) {
						currentMinY = j;
					}
				}
			}
		}
		for(int i = currentMinX; i < maxX; i++) {
			for(int j = currentMinY; j < maxY; j++) {
				tempMatrix[i - currentMinX][j - currentMinY] = matrix[i][j];
			}
		}
		for(int i = 0; i < maxX; i++) {
			for(int j = 0; j < maxY; j++) {
				matrix[i][j] = tempMatrix[i][j];
			}
		}
	}
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);
            }
        }
    }
446790
def concatenate_slh(self, other):
        """"""
        selfS = self.S
        otherS = other.S
        new_S = block_matrix(
                selfS, zerosm((selfS.shape[0], otherS.shape[1]), dtype=int),
                zerosm((otherS.shape[0], selfS.shape[1]), dtype=int), otherS)
        new_L = vstackm((self.L, other.L))
        new_H = self.H + other.H

        return SLH(new_S, new_L, new_H)
490009
def precompute_sharp_round(nxk, nyk, xc, yc):
    

    # Create arrays for the two- and four-fold symmetry computations:
    s4m = np.ones((nyk,nxk),dtype=np.int16)
    s4m[yc, xc] = 0

    s2m = np.ones((nyk,nxk),dtype=np.int16)
    s2m[yc, xc] = 0
    s2m[yc:nyk, 0:xc]     = -1;
    s2m[0:yc+1, xc+1:nxk] = -1;

    return s2m, s4m