Calculate the lower triangular matrix of the Cholesky decomposition of
    a symmetric, positive-definite matrix.
550875
def chol(A):
    
    A = np.array(A)
    assert A.shape[0] == A.shape[1], "Input matrix must be square"

    L = [[0.0] * len(A) for _ in range(len(A))]
    for i in range(len(A)):
        for j in range(i + 1):
            s = sum(L[i][k] * L[j][k] for k in range(j))
            L[i][j] = (
                (A[i][i] - s) ** 0.5 if (i == j) else (1.0 / L[j][j] * (A[i][j] - s))
            )

    return np.array(L)
371693
private void _symmetricTql2 ()
  {
    // 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 < m_nDim; i++)
      m_aEVe[i - 1] = m_aEVe[i];
    m_aEVe[m_nDim - 1] = 0.0;

    double f = 0.0;
    double tst1 = 0.0;
    for (int l = 0; l < m_nDim; l++)
    {
      // Find small subdiagonal element
      tst1 = Math.max (tst1, MathHelper.abs (m_aEVd[l]) + MathHelper.abs (m_aEVe[l]));
      int m = l;
      while (m < m_nDim)
      {
        if (MathHelper.abs (m_aEVe[m]) <= EPSILON * tst1)
          break;
        m++;
      }

      // If m == l, d[l] is an eigenvalue,
      // otherwise, iterate.

      if (m > l)
      {
        do
        {
          // Compute implicit shift
          double g = m_aEVd[l];
          double p = (m_aEVd[l + 1] - g) / (2.0 * m_aEVe[l]);
          double r = MathHelper.hypot (p, 1.0);
          if (p < 0)
          {
            r = -r;
          }
          m_aEVd[l] = m_aEVe[l] / (p + r);
          m_aEVd[l + 1] = m_aEVe[l] * (p + r);
          final double dl1 = m_aEVd[l + 1];
          double h = g - m_aEVd[l];
          for (int i = l + 2; i < m_nDim; i++)
          {
            m_aEVd[i] -= h;
          }
          f += h;

          // Implicit QL transformation.
          p = m_aEVd[m];
          double c = 1.0;
          double c2 = c;
          double c3 = c;
          final double el1 = m_aEVe[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 * m_aEVe[i];
            h = c * p;
            r = MathHelper.hypot (p, m_aEVe[i]);
            m_aEVe[i + 1] = s * r;
            s = m_aEVe[i] / r;
            c = p / r;
            p = c * m_aEVd[i] - s * g;
            m_aEVd[i + 1] = h + s * (c * g + s * m_aEVd[i]);

            // Accumulate transformation.

            for (int k = 0; k < m_nDim; k++)
            {
              h = m_aEigenVector[k][i + 1];
              m_aEigenVector[k][i + 1] = s * m_aEigenVector[k][i] + c * h;
              m_aEigenVector[k][i] = c * m_aEigenVector[k][i] - s * h;
            }
          }
          p = -s * s2 * c3 * el1 * m_aEVe[l] / dl1;
          m_aEVe[l] = s * p;
          m_aEVd[l] = c * p;

          // Check for convergence.

        } while (MathHelper.abs (m_aEVe[l]) > EPSILON * tst1);
      }
      m_aEVd[l] = m_aEVd[l] + f;
      m_aEVe[l] = 0.0;
    }

    // Sort eigenvalues and corresponding vectors.
    for (int i = 0; i < m_nDim - 1; i++)
    {
      int k = i;
      double p = m_aEVd[i];
      for (int j = i + 1; j < m_nDim; j++)
      {
        if (m_aEVd[j] < p)
        {
          k = j;
          p = m_aEVd[j];
        }
      }
      if (k != i)
      {
        m_aEVd[k] = m_aEVd[i];
        m_aEVd[i] = p;
        for (int j = 0; j < m_nDim; j++)
        {
          p = m_aEigenVector[j][i];
          m_aEigenVector[j][i] = m_aEigenVector[j][k];
          m_aEigenVector[j][k] = p;
        }
      }
    }
  }
241475
def d2Ibr_dV2(Ybr, V, lam):
    
    nb = len(V)
    diaginvVm = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))

    Haa = spdiag(mul(-(Ybr.T * lam), V))
    Hva = -1j * Haa * diaginvVm
    Hav = Hva
    Hvv = spmatrix([], [], [], (nb, nb))

    return Haa, Hav, Hva, Hvv
244980
def DSYR_blas(A, x, alpha=1.):
    
    At = blas.dsyr(lower=0, x=x, a=A, alpha=alpha, overwrite_a=False) #See https://github.com/scipy/scipy/issues/8155
    A[:] = At
    symmetrify(A, upper=True)
28120
def trisolve(U, b, uplo='U', transa='N', alpha=1., inplace=False):
    
    if inplace:
        x = b
    else:
        x = b.copy('F')
    if U.shape[0] == 0:
        raise ValueError('Attempted to solve zero-rank triangular system')
    dtrsm_wrap(a=U, b=x, side='L', uplo=uplo, transa=transa, alpha=alpha)
    return x
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;
            }
        }
    }
776553
def cublasDsyrk(handle, uplo, trans, n, k, alpha, A, lda, beta, C, ldc):
    
    
    status = _libcublas.cublasDsyrk_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       n, k, ctypes.byref(cuda.cuFloatComplex(alpha.real,     
                                                                              alpha.imag)),
                                       int(A), lda, 
                                       ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                        beta.imag)),
                                       int(C), ldc)
    cublasCheckStatus(status)
976081
def _normalize_eigensystem(u, lv, rv):
    
    # first normalize the stationary distribution separately
    lv[:, 0] = lv[:, 0] / np.sum(lv[:, 0])

    for i in range(1, lv.shape[1]):
        # the remaining left eigenvectors to satisfy
        # <\phi_i, \phi_i>_{\mu^{-1}} = 1
        lv[:, i] = lv[:, i] / np.sqrt(np.dot(lv[:, i], lv[:, i] / lv[:, 0]))

    for i in range(rv.shape[1]):
        # the right eigenvectors to satisfy <\phi_i, \psi_j> = \delta_{ij}
        rv[:, i] = rv[:, i] / np.dot(lv[:, i], rv[:, i])

    return u, lv, rv
1035851
private Matrix calculateV(IBasis basis) {
        int size = basis.getSize();
        Matrix V = new Matrix(size, size);
        int i, j;
        for (i = 0; i < size; i++)
            for (j = 0; j < size; j++)
                V.matrix[i][j] = basis.calcV(i, j);

        return V;
    }
268773
def _computePartialLikelihoods(self):
        """"""
        for n in range(self.ntips, self.nnodes):
            ni = n - self.ntips # internal node number
            nright = self.rdescend[ni]
            nleft = self.ldescend[ni]
            if nright < self.ntips:
                istipr = True
            else:
                istipr = False
            if nleft < self.ntips:
                istipl = True
            else:
                istipl = False
            tright = self.t[nright]
            tleft = self.t[nleft]
            self.L[n] = scipy.ndarray(self._Lshape, dtype='float')
            if self.dparamscurrent:
                for param in self._paramlist_PartialLikelihoods:
                    self.dL[param][n] = scipy.ndarray(self._dLshape[param],
                            dtype='float')
            if self.dtcurrent:
                for n2 in self.dL_dt.keys():
                    self.dL_dt[n2][n] = scipy.zeros(self._Lshape, dtype='float')
            for k in self._catindices:
                if istipr:
                    Mright = MLright = self._M(k, tright,
                            self.tips[nright], self.gaps[nright])
                else:
                    Mright = self._M(k, tright)
                    MLright = broadcastMatrixVectorMultiply(Mright,
                            self.L[nright][k])
                if istipl:
                    Mleft = MLleft = self._M(k, tleft,
                            self.tips[nleft], self.gaps[nleft])
                else:
                    Mleft = self._M(k, tleft)
                    MLleft = broadcastMatrixVectorMultiply(Mleft,
                            self.L[nleft][k])
                self.L[n][k] = MLright * MLleft

                if self.dtcurrent:
                    for (tx, Mx, nx, MLxother, istipx) in [
                            (tright, Mright, nright, MLleft, istipr),
                            (tleft, Mleft, nleft, MLright, istipl)]:
                        if istipx:
                            tipsx = self.tips[nx]
                            gapsx = self.gaps[nx]
                        else:
                            tipsx = gapsx = None
                        dM_dt = self._dM(k, tx, 't', Mx, tipsx, gapsx)
                        if istipx:
                            LdM_dt = dM_dt
                        else:
                            LdM_dt = broadcastMatrixVectorMultiply(
                                    dM_dt, self.L[nx][k])
                        self.dL_dt[nx][n][k] = LdM_dt * MLxother
                        for ndx in self.descendants[nx]:
                            self.dL_dt[ndx][n][k] = broadcastMatrixVectorMultiply(
                                    Mx, self.dL_dt[ndx][nx][k]) * MLxother

                if self.dparamscurrent:
                    for param in self._paramlist_PartialLikelihoods:
                        if istipr:
                            dMright = self._dM(k, tright, param, Mright,
                                    self.tips[nright], self.gaps[nright])
                        else:
                            dMright = self._dM(k, tright, param, Mright)
                        if istipl:
                            dMleft = self._dM(k, tleft, param, Mleft,
                                    self.tips[nleft], self.gaps[nleft])
                        else:
                            dMleft = self._dM(k, tleft, param, Mleft)
                        for j in self._sub_index_param(param):
                            if istipr:
                                dMLright = dMright[j]
                                MdLright = 0
                            else:
                                dMLright = broadcastMatrixVectorMultiply(
                                        dMright[j], self.L[nright][k])
                                MdLright = broadcastMatrixVectorMultiply(
                                        Mright, self.dL[param][nright][k][j])
                            if istipl:
                                dMLleft = dMleft[j]
                                MdLleft = 0
                            else:
                                dMLleft = broadcastMatrixVectorMultiply(
                                        dMleft[j], self.L[nleft][k])
                                MdLleft = broadcastMatrixVectorMultiply(
                                        Mleft, self.dL[param][nleft][k][j])
                            self.dL[param][n][k][j] = ((dMLright + MdLright)
                                    * MLleft + MLright * (dMLleft + MdLleft))

            if ni > 0 and ni % self.underflowfreq == 0:
                # rescale by same amount for each category k
                scale = scipy.amax(scipy.array([scipy.amax(self.L[n][k],
                        axis=1) for k in self._catindices]), axis=0)
                assert scale.shape == (self.nsites,)
                self.underflowlogscale += scipy.log(scale)
                for k in self._catindices:
                    self.L[n][k] /= scale[:, scipy.newaxis]
                    if self.dtcurrent:
                        for n2 in self.dL_dt.keys():
                            self.dL_dt[n2][n][k] /= scale[:, scipy.newaxis]
                    if self.dparamscurrent:
                        for param in self._paramlist_PartialLikelihoods:
                            for j in self._sub_index_param(param):
                                self.dL[param][n][k][j] /= scale[:, scipy.newaxis]

            # free unneeded memory by deleting already used values
            for ntodel in [nright, nleft]:
                if ntodel in self.L:
                    del self.L[ntodel]
                if self.dparamscurrent:
                    for param in self._paramlist_PartialLikelihoods:
                        if ntodel in self.dL[param]:
                            del self.dL[param][ntodel]
                if self.dtcurrent:
                    for n2 in self.dL_dt.keys():
                        if ntodel in self.dL_dt[n2]:
                            del self.dL_dt[n2][ntodel]
1269882
def _coupling_matrix(self, lmax, nwin=None, weights=None):
        """"""
        if nwin is None:
            nwin = self.nwin

        if weights is None:
            weights = self.weights

        tapers_power = _np.zeros((self.lwin+1, nwin))
        for i in range(nwin):
            tapers_power[:, i] = _spectrum(self.to_array(i),
                                           normalization='4pi',
                                           convention='power', unit='per_l')

        if weights is None:
            return _shtools.SHMTCouplingMatrix(lmax, tapers_power, k=nwin)
        else:
            return _shtools.SHMTCouplingMatrix(lmax, tapers_power, k=nwin,
                                               taper_wt=self.weights)
237334
public Matrix3x2f mulLocal(Matrix3x2fc left, Matrix3x2f dest) {
        float nm00 = left.m00() * m00 + left.m10() * m01;
        float nm01 = left.m01() * m00 + left.m11() * m01;
        float nm10 = left.m00() * m10 + left.m10() * m11;
        float nm11 = left.m01() * m10 + left.m11() * m11;
        float nm20 = left.m00() * m20 + left.m10() * m21 + left.m20();
        float nm21 = left.m01() * m20 + left.m11() * m21 + left.m21();
        dest.m00 = nm00;
        dest.m01 = nm01;
        dest.m10 = nm10;
        dest.m11 = nm11;
        dest.m20 = nm20;
        dest.m21 = nm21;
        return dest;
    }
27487
def solve_chol(A,B):
    
    # X = linalg.solve(A,linalg.solve(A.transpose(),B))
    # much faster version
    X = linalg.cho_solve((A, True), B)
    return X
979000
def MatrixTriangularSolve(a, rhs, lower, adj):
    
    trans = 0 if not adj else 2

    r = np.empty(rhs.shape).astype(a.dtype)
    for coord in np.ndindex(a.shape[:-2]):
        pos = coord + (Ellipsis,)
        r[pos] = sp.linalg.solve_triangular(a[pos] if not adj else np.conj(a[pos]), rhs[pos],
                                            trans=trans, lower=lower)

    return r,
776373
@Override
	public void computeHessian(DMatrixRMaj jacLeft , DMatrixRMaj jacRight) {
		A.reshape(jacLeft.numCols,jacLeft.numCols);
		B.reshape(jacLeft.numCols,jacRight.numCols);
		D.reshape(jacRight.numCols,jacRight.numCols);

		// take advantage of the inner product's symmetry when possible to reduce
		// the number of calculations
		MatrixMultProduct_DDRM.inner_reorder_lower(jacLeft,A);
		CommonOps_DDRM.symmLowerToFull(A);
		CommonOps_DDRM.multTransA(jacLeft,jacRight,B);
		MatrixMultProduct_DDRM.inner_reorder_lower(jacRight,D);
		CommonOps_DDRM.symmLowerToFull(D);
	}
205648
def MatrixSolve(a, rhs, adj):
    
    return np.linalg.solve(a if not adj else _adjoint(a), rhs),
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)