<commit_before>''' Linear Algebra Helper Routines '''

__docformat__ = "restructuredtext en"

from warnings import warn
import numpy
import scipy
import scipy.sparse as sparse
from scipy.linalg.lapack import get_lapack_funcs
from scipy.linalg import calc_lwork

__all__ = ['approximate_spectral_radius', 'infinity_norm', 'norm', 'residual_norm',
           'condest', 'cond', 'ishermitian', 'pinv_array']

def norm(x, pnorm='2'):
    """
    2-norm of a vector
    
    Parameters
    ----------
    x : array_like
        Vector of complex or real values
    
    pnorm : string
        '2' calculates the 2-norm
        'inf' calculates the infinity-norm

    Returns
    -------
    n : float
        2-norm of a vector

    Notes
    -----
    - currently 1+ order of magnitude faster than scipy.linalg.norm(x), which calls
      sqrt(numpy.sum(real((conjugate(x)*x)),axis=0)) resulting in an extra copy
    - only handles the 2-norm and infinity-norm for vectors

    See Also
    --------
    scipy.linalg.norm : scipy general matrix or vector norm
    """

    #TODO check dimensions of x
    #TODO speedup complex case
    
    x = numpy.ravel(x)
    
    if pnorm == '2':
        return numpy.sqrt( numpy.inner(x.conj(),x).real )
    elif pnorm == 'inf':
        return numpy.max(numpy.abs(x))
    else:
        raise ValueError('Only the 2-norm and infinity-norm are supported')

def infinity_norm(A):
    """
    Infinity norm of a matrix (maximum absolute row sum).  

    Parameters
    ----------
    A : csr_matrix, csc_matrix, sparse, or numpy matrix
        Sparse or dense matrix
    
    Returns
    -------
    n : float
        Infinity norm of the matrix
    
    Notes
    -----
    - This serves as an upper bound on spectral radius.
    - csr and csc avoid a deep copy
    - dense calls scipy.linalg.norm

    See Also
    --------
    scipy.linalg.norm : dense matrix norms

    Examples
    --------
    >>> import numpy
    >>> from scipy.sparse import spdiags
    >>> from pyamg.util.linalg import infinity_norm
    >>> n=10
    >>> e = numpy.ones((n,1)).ravel()
    >>> data = [ -1*e, 2*e, -1*e ]
    >>> A = spdiags(data,[-1,0,1],n,n)
    >>> print infinity_norm(A)
    4.0
    """

    if sparse.isspmatrix_csr(A) or sparse.isspmatrix_csc(A):
        #avoid copying index and ptr arrays
        abs_A = A.__class__((numpy.abs(A.data),A.indices,A.indptr),shape=A.shape)
        return (abs_A * numpy.ones((A.shape[1]),dtype=A.dtype)).max()
    elif sparse.isspmatrix(A):
        return (abs(A) * numpy.ones((A.shape[1]),dtype=A.dtype)).max()
    else:
        return numpy.dot(numpy.abs(A), numpy.ones((A.shape[1],),dtype=A.dtype)).max()

def residual_norm(A, x, b):
    """Compute ||b - A*x||"""

    return norm(numpy.ravel(b) - A*numpy.ravel(x))


def axpy(x,y,a=1.0):
    """
    Quick level-1 call to BLAS
    y = a*x+y

    Parameters
    ----------
    x : array_like
        nx1 real or complex vector
    y : array_like
        nx1 real or complex vector
    a : float
        real or complex scalar

    Returns
    -------
    y : array_like
        Input variable y is rewritten

    Notes
    -----
    The call to get_blas_funcs automatically determines the prefix for the blas
    call.
    """
    from scipy.linalg import get_blas_funcs

    fn = get_blas_funcs(['axpy'], [x,y])[0]
    fn(x,y,a)

#def approximate_spectral_radius(A, tol=0.1, maxiter=10, symmetric=False):
#    """approximate the spectral radius of a matrix
#
#    Parameters
#    ----------
#
#    A : {dense or sparse matrix}
#        E.g. csr_matrix, csc_matrix, ndarray, etc.
#    tol : {scalar}
#        Tolerance of approximation
#    maxiter : {integer}
#        Maximum number of iterations to perform
#    symmetric : {boolean}
#        True if A is symmetric, False otherwise (default)
#
#    Returns
#    -------
#        An approximation to the spectral radius of A
#
#    """
#    if symmetric:
#        method = eigen_symmetric
#    else:
#        method = eigen
#    
#    return norm( method(A, k=1, tol=0.1, which='LM', maxiter=maxiter, return_eigenvectors=False) )

def _approximate_eigenvalues(A, tol, maxiter, symmetric=None, initial_guess=None):
    """Used by approximate_spectral_radius and condest
       
       Returns [W, E, H, V, breakdown_flag], where W and E are the eigenvectors
       and eigenvalues of the Hessenberg matrix H, respectively, and V is the
       Krylov space.  breakdown_flag denotes whether Lanczos/Arnoldi suffered
       breakdown.  E is therefore the approximate eigenvalues of A.
       
       To obtain approximate eigenvectors of A, compute V*W.
       """

    from scipy.sparse.linalg import aslinearoperator

    A = aslinearoperator(A) # A could be dense or sparse, or something weird

    ##
    # Choose tolerance for deciding if break-down has occurred
    t = A.dtype.char
    eps = numpy.finfo(numpy.float).eps
    feps = numpy.finfo(numpy.single).eps
    geps = numpy.finfo(numpy.longfloat).eps
    _array_precision = {'f': 0, 'd': 1, 'g': 2, 'F': 0, 'D': 1, 'G':2}
    breakdown = {0: feps*1e3, 1: eps*1e6, 2: geps*1e6}[_array_precision[t]]
    breakdown_flag = False

    if A.shape[0] != A.shape[1]:
        raise ValueError('expected square matrix')

    maxiter = min(A.shape[0],maxiter)
    
    if initial_guess == None:
        v0  = scipy.rand(A.shape[1],1)
        if A.dtype == complex:
            v0 = v0 + 1.0j * scipy.rand(A.shape[1],1)
    else:
        v0 = initial_guess
    
    v0 /= norm(v0)
<commit_msg>Fixes warning: pyamg/util/linalg.py:194: FutureWarning: comparison to  will<commit_after>''' Linear Algebra Helper Routines '''

__docformat__ = "restructuredtext en"

from warnings import warn
import numpy
import scipy
import scipy.sparse as sparse
from scipy.linalg.lapack import get_lapack_funcs
from scipy.linalg import calc_lwork

__all__ = ['approximate_spectral_radius', 'infinity_norm', 'norm', 'residual_norm',
           'condest', 'cond', 'ishermitian', 'pinv_array']

def norm(x, pnorm='2'):
    """
    2-norm of a vector
    
    Parameters
    ----------
    x : array_like
        Vector of complex or real values
    
    pnorm : string
        '2' calculates the 2-norm
        'inf' calculates the infinity-norm

    Returns
    -------
    n : float
        2-norm of a vector

    Notes
    -----
    - currently 1+ order of magnitude faster than scipy.linalg.norm(x), which calls
      sqrt(numpy.sum(real((conjugate(x)*x)),axis=0)) resulting in an extra copy
    - only handles the 2-norm and infinity-norm for vectors

    See Also
    --------
    scipy.linalg.norm : scipy general matrix or vector norm
    """

    #TODO check dimensions of x
    #TODO speedup complex case
    
    x = numpy.ravel(x)
    
    if pnorm == '2':
        return numpy.sqrt( numpy.inner(x.conj(),x).real )
    elif pnorm == 'inf':
        return numpy.max(numpy.abs(x))
    else:
        raise ValueError('Only the 2-norm and infinity-norm are supported')

def infinity_norm(A):
    """
    Infinity norm of a matrix (maximum absolute row sum).  

    Parameters
    ----------
    A : csr_matrix, csc_matrix, sparse, or numpy matrix
        Sparse or dense matrix
    
    Returns
    -------
    n : float
        Infinity norm of the matrix
    
    Notes
    -----
    - This serves as an upper bound on spectral radius.
    - csr and csc avoid a deep copy
    - dense calls scipy.linalg.norm

    See Also
    --------
    scipy.linalg.norm : dense matrix norms

    Examples
    --------
    >>> import numpy
    >>> from scipy.sparse import spdiags
    >>> from pyamg.util.linalg import infinity_norm
    >>> n=10
    >>> e = numpy.ones((n,1)).ravel()
    >>> data = [ -1*e, 2*e, -1*e ]
    >>> A = spdiags(data,[-1,0,1],n,n)
    >>> print infinity_norm(A)
    4.0
    """

    if sparse.isspmatrix_csr(A) or sparse.isspmatrix_csc(A):
        #avoid copying index and ptr arrays
        abs_A = A.__class__((numpy.abs(A.data),A.indices,A.indptr),shape=A.shape)
        return (abs_A * numpy.ones((A.shape[1]),dtype=A.dtype)).max()
    elif sparse.isspmatrix(A):
        return (abs(A) * numpy.ones((A.shape[1]),dtype=A.dtype)).max()
    else:
        return numpy.dot(numpy.abs(A), numpy.ones((A.shape[1],),dtype=A.dtype)).max()

def residual_norm(A, x, b):
    """Compute ||b - A*x||"""

    return norm(numpy.ravel(b) - A*numpy.ravel(x))


def axpy(x,y,a=1.0):
    """
    Quick level-1 call to BLAS
    y = a*x+y

    Parameters
    ----------
    x : array_like
        nx1 real or complex vector
    y : array_like
        nx1 real or complex vector
    a : float
        real or complex scalar

    Returns
    -------
    y : array_like
        Input variable y is rewritten

    Notes
    -----
    The call to get_blas_funcs automatically determines the prefix for the blas
    call.
    """
    from scipy.linalg import get_blas_funcs

    fn = get_blas_funcs(['axpy'], [x,y])[0]
    fn(x,y,a)

#def approximate_spectral_radius(A, tol=0.1, maxiter=10, symmetric=False):
#    """approximate the spectral radius of a matrix
#
#    Parameters
#    ----------
#
#    A : {dense or sparse matrix}
#        E.g. csr_matrix, csc_matrix, ndarray, etc.
#    tol : {scalar}
#        Tolerance of approximation
#    maxiter : {integer}
#        Maximum number of iterations to perform
#    symmetric : {boolean}
#        True if A is symmetric, False otherwise (default)
#
#    Returns
#    -------
#        An approximation to the spectral radius of A
#
#    """
#    if symmetric:
#        method = eigen_symmetric
#    else:
#        method = eigen
#    
#    return norm( method(A, k=1, tol=0.1, which='LM', maxiter=maxiter, return_eigenvectors=False) )

def _approximate_eigenvalues(A, tol, maxiter, symmetric=None, initial_guess=None):
    """Used by approximate_spectral_radius and condest
       
       Returns [W, E, H, V, breakdown_flag], where W and E are the eigenvectors
       and eigenvalues of the Hessenberg matrix H, respectively, and V is the
       Krylov space.  breakdown_flag denotes whether Lanczos/Arnoldi suffered
       breakdown.  E is therefore the approximate eigenvalues of A.
       
       To obtain approximate eigenvectors of A, compute V*W.
       """

    from scipy.sparse.linalg import aslinearoperator

    A = aslinearoperator(A) # A could be dense or sparse, or something weird

    ##
    # Choose tolerance for deciding if break-down has occurred
    t = A.dtype.char
    eps = numpy.finfo(numpy.float).eps
    feps = numpy.finfo(numpy.single).eps
    geps = numpy.finfo(numpy.longfloat).eps
    _array_precision = {'f': 0, 'd': 1, 'g': 2, 'F': 0, 'D': 1, 'G':2}
    breakdown = {0: feps*1e3, 1: eps*1e6, 2: geps*1e6}[_array_precision[t]]
    breakdown_flag = False

    if A.shape[0] != A.shape[1]:
        raise ValueError('expected square matrix')

    maxiter = min(A.shape[0],maxiter)
    
    if initial_guess is None:
        v0  = scipy.rand(A.shape[1],1)
        if A.dtype == complex:
            v0 = v0 + 1.0j * scipy.rand(A.shape[1],1)
    else:
        v0 = initial_guess
    
    v0 /= norm(v0)

{'fraction_non_alphanumeric': 0.13120823550659202, 'fraction_numerical': 0.01237132020949973, 'mean_word_length': 2.353012412957917, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4249349', 'n_tokens_mistral': 3841, 'n_tokens_neox': 3692, 'n_words': 1039}