Compute A <- alpha * x * y^H + A (general rank-1 update)
509502
public static ComplexFloatMatrix gerc(ComplexFloat alpha, ComplexFloatMatrix x,
			ComplexFloatMatrix y, ComplexFloatMatrix a) {
		NativeBlas.cgerc(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
				0, a.rows);
		return a;
	}
35195
def cublasZhpr(handle, uplo, n, alpha, x, incx, AP):
    
    
    status = _libcublas.cublasZhpr_v2(handle,
                                      _CUBLAS_FILL_MODE[uplo], 
                                      n, ctypes.byref(ctypes.c_double(alpha)),
                                      int(x), incx, int(AP))
    cublasCheckStatus(status)
607444
def cublasDspr(handle, uplo, n, alpha, x, incx, AP):
    

    status = _libcublas.cublasDspr_v2(handle, 
                                      _CUBLAS_FILL_MODE[uplo], n,                                       
                                      ctypes.byref(ctypes.c_double(alpha)), 
                                      int(x), incx, int(AP))                                           
    cublasCheckStatus(status)
910674
def cublasZgerc(handle, m, n, alpha, x, incx, y, incy, A, lda):
    

    status = _libcublas.cublasZgerc_v2(handle,
                                       m, n, ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                               alpha.imag)),
                                       int(x), incx, int(y), incy, int(A), lda)
    cublasCheckStatus(status)
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;
	}
171276
def _h(self, X, tau):

        

        def cos_func(y):
            if isinstance(y, float) or isinstance(y, int):
                if y > 1:
                    return np.arccosh(y)
                else:
                    return np.arccos(y)
            else:
                values = np.ones_like(y)
                inds1 = np.where(y < 1)
                inds2 = np.where(y > 1)
                values[inds1] = np.arccos(y[inds1])
                # values[inds2] = np.arccosh(y[inds2])
                values[inds2] = np.arccos((1 - y[inds2]) ** .5)
                return values

        def nfw_func(x):
            if isinstance(x,float) or isinstance(x,int):
                if x<1:
                    return np.log(x / 2.) ** 2 - np.arccosh(1. / x) ** 2
                else:
                    return np.log(x / 2.)**2 - np.arccos(1. /x)**2
            else:
                inds1 = np.where(x < 1)
                inds2 = np.where(x >= 1)
                y = np.zeros_like(x)
                y[inds1] = np.log(x[inds1] / 2.) ** 2 - np.arccosh(1. / x[inds1]) ** 2
                y[inds2] = np.log(x[inds2] / 2.) ** 2 - np.arccosh(1. / x[inds2]) ** 2

                return y

        def tnfw_func(x,tau):

            u = x**2
            t2 = tau**2
            Lx = self.L(u**.5, tau)

            return (t2 + 1) ** -2 * (
                2 * t2 * np.pi * (tau - (t2 + u) ** .5 + tau * np.log(tau + (t2 + u) ** .5))
                +
                2 * (t2 - 1) * tau * (t2 + u) ** .5 * Lx
                +
                t2 * (t2 - 1) * Lx ** 2
                +
                4 * t2 * (u - 1) * self.F(u**.5)
                +
                t2 * (t2 - 1) * (cos_func(u ** -0.5)) ** 2
                +
                t2 * ((t2 - 1) * np.log(tau) - t2 - 1) * np.log(u)
                -
                t2 * (
                (t2 - 1) * np.log(tau) * np.log(4 * tau) + 2 * np.log(0.5 * tau) - 2 * tau * (tau - np.pi) * np.log(
                    tau * 2)))

        c = 1e-9
        rescale = 1

        if np.any(X-c<1):

            warnings.warn('Truncated NFW potential not yet implemented for x<1. Using the expression for the NFW '
                          'potential in this regime isntead.')
            rescale = tnfw_func(1.00001,tau)*nfw_func(1)**-1

        if isinstance(X,float) or isinstance(X,int):
            if X<1:
                X = min(X,0.001)
                return nfw_func(X)*rescale
            else:
                return tnfw_func(X,tau)
        else:

            X[np.where(X<0.001)] = 0.001

            values = np.zeros_like(X)
            inds1,inds2 = np.where(X<1),np.where(X>=1)
            values[inds1] = nfw_func(X[inds1])*rescale
            values[inds2] = tnfw_func(X[inds2],tau)

            return values
56217
def rhalf_cauchy(alpha, beta, size=None):
    

    return abs(alpha + beta * np.tan(pi * random_number(size) - pi / 2.0))
115515
def low_rank_approximation(A, r):
    
    try:
        u, s, v = np.linalg.svd(A, full_matrices=False)
    except np.linalg.LinAlgError as e:
        print('Matrix:', A)
        print('Matrix rank:', np.linalg.matrix_rank(A))
        raise

    Ar = np.zeros((len(u), len(v)), dtype=u.dtype)
    buf_ = np.empty_like(Ar)
    sc_vec_ = np.empty((v.shape[1],), dtype=v.dtype)
    for i in range(r):
        np.multiply(v[i], s[i], out=sc_vec_)
        np.outer(u[:, i], sc_vec_, out=buf_)
        Ar += buf_
    return Ar
1044815
def cublasZsyr(handle, uplo, n, alpha, x, incx, A, lda):
    

    status = _libcublas.cublasZsyr_v2(handle,
                                        _CUBLAS_FILL_MODE[uplo], n, 
                                        ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                        alpha.imag)),
                                        int(x), incx, int(A), lda)
    cublasCheckStatus(status)
1333315
def relax_AX(self):
        """"""

        self.AXnr = self.X
        if self.rlx == 1.0:
            self.AX = self.X
        else:
            alpha = self.rlx
            self.AX = alpha*self.X + (1 - alpha)*self.Y[..., np.newaxis]
652284
public static FloatMatrix ger(float alpha, FloatMatrix x,
			FloatMatrix y, FloatMatrix a) {
		NativeBlas.sger(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
				0, a.rows);
		return a;
	}
675501
private void computeGradientAndHessian(DMatrixRMaj param  )
    {
        // residuals = f(x) - y
        function.compute(param, residuals);

        computeNumericalJacobian(param,jacobian);

        CommonOps_DDRM.multTransA(jacobian, residuals, g);
        CommonOps_DDRM.multTransA(jacobian, jacobian,  H);

        CommonOps_DDRM.extractDiag(H,Hdiag);
    }
1292473
def gradient_descent(alpha=0.01, num_iters=1500)
      m = @xs.length
      alpha_m = (alpha / m)

      (num_iters).times do

        @thetas = 
          @xs.map do |x_row|
            x_row.zip(@thetas).map do |x, theta|
             x * theta
            end.inject(:+)
          end.zip(@ys).map do |(predict_val, y)|  
            predict_val - y 
          end.zip(@xs).map do |(y, x_row)|
            x_row.map { |x| x * y }
          end.reduce([0.0] * 2) do |sum, val|
            sum.zip(val).map { |(a, b)| a + b}
          end.map do |value|
            alpha_m * value
          end.zip( @thetas ).map do |(a, b)|
            b - a
          end
      end

      @thetas
    end
888741
def newton_update(x, a, c, step=1.0):
    
    return x - step*np.linalg.inv(hessian(x, a)).dot(gradient(x, a, c))
174251
def par_u1step(i):
    
    """

    global mp_U1
    grpind = slice(mp_grp[i], mp_grp[i+1])
    mp_U1[grpind] += mp_alpha*(mp_X[grpind] - mp_Y1[grpind])
431795
def get_residual_norms(H, self_adjoint=False):
    ''''''
    H = H.copy()
    n_, n = H.shape
    y = numpy.eye(n_, 1, dtype=H.dtype)
    resnorms = [1.]
    for i in range(n_-1):
        G = Givens(H[i:i+2, [i]])
        if self_adjoint:
            H[i:i+2, i:i+3] = G.apply(H[i:i+2, i:i+3])
        else:
            H[i:i+2, i:] = G.apply(H[i:i+2, i:])
        y[i:i+2] = G.apply(y[i:i+2])
        resnorms.append(numpy.abs(y[i+1, 0]))
    if n_ == n:
        resnorms.append(0.)
    return numpy.array(resnorms)
362156
protected void upgradeSolution( DMatrixRMaj X ) {
        DMatrixRMaj z = Y; // recycle Y

        // compute the z which will minimize the 2-norm of X
        // because of the identity matrix tacked onto the end 'A' should never be singular
        if( !internalSolver.setA(W) )
            throw new RuntimeException("This should never happen.  Is input NaN?");
        z.reshape(numCols-rank,1);
        internalSolver.solve(X, z);

        // compute X by tweaking the original
        CommonOps_DDRM.multAdd(-1, W, z, X);
    }