def european_call(r, sigma, T, Smax, m, n, Smin=0.0, barrier=None):
    
    X = linspace(0.0, Smax, n+2)
    X = X[1:-1]
    
    Fp = clip(X-K, 0.0, 1e600)
    
    if barrier is None:
        Fu = Smax - K*exp(-r * linspace(0.0, T, m+1))
        Fl = zeros((m+1, ))
    elif barrier == 'up-and-out':
        Fu = Fl = zeros((m+1,))
    
    bss = BS_FDM_cn(r, sigma, T, Smin, Smax, Fl, Fu, Fp, m, n)
    return X, bss.solve()

 def test_vecself(self):
     """Ticket 844."""
     # Inner product of a vector with itself segfaults or give meaningless
     # result
     a = zeros(shape = (1, 80), dtype = float64)
     p = inner_(a, a)
     assert_almost_equal(p, 0, decimal = DECPREC)

def solve(a, b):
    """Return the solution of a*x = b
    """
    one_eq = len(b.shape) == 1
    if one_eq:
        b = b[:, newaxis]
    _assertRank2(a, b)
    _assertSquareness(a)
    n_eq = a.shape[0]
    n_rhs = b.shape[1]
    if n_eq != b.shape[0]:
        raise LinAlgError, 'Incompatible dimensions'
    t, result_t = _commonType(a, b)
#    lapack_routine = _findLapackRoutine('gesv', t)
    if isComplexType(t):
        lapack_routine = lapack_lite.zgesv
    else:
        lapack_routine = lapack_lite.dgesv
    a, b = _fastCopyAndTranspose(t, a, b)
    pivots = zeros(n_eq, fortran_int)
    results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
    if results['info'] > 0:
        raise LinAlgError, 'Singular matrix'
    if one_eq:
        return b.ravel().astype(result_t)
    else:
        return b.transpose().astype(result_t)

def european_put(r, sigma, T, Bu, m, n, Bl=0.0, barrier=None, method=None):
  """Compute prices for a European-style put option."""

  X = linspace(0.0, B, n+2)
  X = X[1:-1]

  Fp = clip(K-X, 0.0, 1e600)
  
  if barrier is None:
    Fu = zeros((m+1,))
    Fl = K*exp(-r * linspace(0.0, T, m+1))
  elif barrier == 'up-and-out':
    Fu = Fl = zeros((m+1,))

  bss = BlackScholesSolver(r, sigma, T, Bl, Bu, Fl, Fu, Fp, m, n)
  return X, bss.solve(method)

def det(a):
    """Compute the determinant of a matrix

    Parameters
    ----------
    a : array-like, shape (M, M)

    Returns
    -------
    det : float or complex
        Determinant of a

    Notes
    -----
    The determinant is computed via LU factorization, LAPACK routine z/dgetrf.
    """
    a = asarray(a)
    _assertRank2(a)
    _assertSquareness(a)
    t, result_t = _commonType(a)
    a = _fastCopyAndTranspose(t, a)
    n = a.shape[0]
    if isComplexType(t):
        lapack_routine = lapack_lite.zgetrf
    else:
        lapack_routine = lapack_lite.dgetrf
    pivots = zeros((n,), fortran_int)
    results = lapack_routine(n, n, a, n, pivots, 0)
    info = results['info']
    if (info < 0):
        raise TypeError, "Illegal input to Fortran routine"
    elif (info > 0):
        return 0.0
    sign = add.reduce(pivots != arange(1, n+1)) % 2
    return (1.-2.*sign)*multiply.reduce(diagonal(a), axis=-1)

def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti,
             work_function=fftpack.cfftf, fft_cache = _fft_cache ):
    a = asarray(a)

    if n == None: n = a.shape[axis]

    if n < 1: raise ValueError("Invalid number of FFT data points (%d) specified." % n)

    try:
        wsave = fft_cache[n]
    except(KeyError):
        wsave = init_function(n)
        fft_cache[n] = wsave

    if a.shape[axis] != n:
        s = list(a.shape)
        if s[axis] > n:
            index = [slice(None)]*len(s)
            index[axis] = slice(0,n)
            a = a[index]
        else:
            index = [slice(None)]*len(s)
            index[axis] = slice(0,s[axis])
            s[axis] = n
            z = zeros(s, a.dtype.char)
            z[index] = a
            a = z

    if axis != -1:
        a = swapaxes(a, axis, -1)
    r = work_function(a, wsave)
    if axis != -1:
        r = swapaxes(r, axis, -1)
    return r

 def test_log2(self):
     a = nx.array([4.5, 2.3, 6.5])
     out = nx.zeros(a.shape, float)
     tgt = nx.array([2.169925, 1.20163386, 2.70043972])
     res = ufl.log2(a)
     assert_almost_equal(res, tgt)
     res = ufl.log2(a, out)
     assert_almost_equal(res, tgt)
     assert_almost_equal(out, tgt)

    def test_isneginf(self):
        a = nx.array([nx.inf, -nx.inf, nx.nan, 0.0, 3.0, -3.0])
        out = nx.zeros(a.shape, bool)
        tgt = nx.array([False, True, False, False, False, False])

        res = ufl.isneginf(a)
        assert_equal(res, tgt)
        res = ufl.isneginf(a, out)
        assert_equal(res, tgt)
        assert_equal(out, tgt)

    def test_fix(self):
        a = nx.array([[1.0, 1.1, 1.5, 1.8], [-1.0, -1.1, -1.5, -1.8]])
        out = nx.zeros(a.shape, float)
        tgt = nx.array([[1., 1., 1., 1.], [-1., -1., -1., -1.]])

        res = ufl.fix(a)
        assert_equal(res, tgt)
        res = ufl.fix(a, out)
        assert_equal(res, tgt)
        assert_equal(out, tgt)
        assert_equal(ufl.fix(3.14), 3)

def rand(*args):
    """Returns an array of random numbers with the given shape.

    This only uses the standard library, so it is useful for testing purposes.
    """
    import random
    from numpy.core import zeros, float64
    results = zeros(args, float64)
    f = results.flat
    for i in range(len(f)):
        f[i] = random.random()
    return results

    def test_isneginf(self):
        a = nx.array([nx.inf, -nx.inf, nx.nan, 0.0, 3.0, -3.0])
        out = nx.zeros(a.shape, bool)
        tgt = nx.array([False, True, False, False, False, False])

        res = ufl.isneginf(a)
        assert_equal(res, tgt)
        res = ufl.isneginf(a, out)
        assert_equal(res, tgt)
        assert_equal(out, tgt)

        a = a.astype(np.complex)
        with assert_raises(TypeError):
            ufl.isneginf(a)

def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti,
             work_function=fftpack.cfftf, fft_cache=_fft_cache):
    a = asarray(a)
    axis = normalize_axis_index(axis, a.ndim)

    if n is None:
        n = a.shape[axis]

    if n < 1:
        raise ValueError("Invalid number of FFT data points (%d) specified."
                         % n)

    # We have to ensure that only a single thread can access a wsave array
    # at any given time. Thus we remove it from the cache and insert it
    # again after it has been used. Multiple threads might create multiple
    # copies of the wsave array. This is intentional and a limitation of
    # the current C code.
    wsave = fft_cache.pop_twiddle_factors(n)
    if wsave is None:
        wsave = init_function(n)

    if a.shape[axis] != n:
        s = list(a.shape)
        if s[axis] > n:
            index = [slice(None)]*len(s)
            index[axis] = slice(0, n)
            a = a[tuple(index)]
        else:
            index = [slice(None)]*len(s)
            index[axis] = slice(0, s[axis])
            s[axis] = n
            z = zeros(s, a.dtype.char)
            z[tuple(index)] = a
            a = z

    if axis != a.ndim - 1:
        a = swapaxes(a, axis, -1)
    r = work_function(a, wsave)
    if axis != a.ndim - 1:
        r = swapaxes(r, axis, -1)

    # As soon as we put wsave back into the cache, another thread could pick it
    # up and start using it, so we must not do this until after we're
    # completely done using it ourselves.
    fft_cache.put_twiddle_factors(n, wsave)

    return r

def det(a):
    """
    Compute the determinant of an array.

    Parameters
    ----------
    a : array_like, shape (M, M)
        Input array.

    Returns
    -------
    det : ndarray
        Determinant of `a`.

    Notes
    -----
    The determinant is computed via LU factorization using the LAPACK
    routine z/dgetrf.

    Examples
    --------
    The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:

    >>> a = np.array([[1, 2], [3, 4]])
    >>> np.linalg.det(a)
    -2.0

    """
    a = asarray(a)
    _assertRank2(a)
    _assertSquareness(a)
    t, result_t = _commonType(a)
    a = _fastCopyAndTranspose(t, a)
    n = a.shape[0]
    if isComplexType(t):
        lapack_routine = lapack_lite.zgetrf
    else:
        lapack_routine = lapack_lite.dgetrf
    pivots = zeros((n,), fortran_int)
    results = lapack_routine(n, n, a, n, pivots, 0)
    info = results['info']
    if (info < 0):
        raise TypeError, "Illegal input to Fortran routine"
    elif (info > 0):
        return 0.0
    sign = add.reduce(pivots != arange(1, n+1)) % 2
    return (1.-2.*sign)*multiply.reduce(diagonal(a), axis=-1)

def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti,
             work_function=fftpack.cfftf, fft_cache=_fft_cache):
    a = asarray(a)

    if n is None:
        n = a.shape[axis]

    if n < 1:
        raise ValueError("Invalid number of FFT data points (%d) specified."
                         % n)

    try:
        # Thread-safety note: We rely on list.pop() here to atomically
        # retrieve-and-remove a wsave from the cache.  This ensures that no
        # other thread can get the same wsave while we're using it.
        wsave = fft_cache.setdefault(n, []).pop()
    except (IndexError):
        wsave = init_function(n)

    if a.shape[axis] != n:
        s = list(a.shape)
        if s[axis] > n:
            index = [slice(None)]*len(s)
            index[axis] = slice(0, n)
            a = a[index]
        else:
            index = [slice(None)]*len(s)
            index[axis] = slice(0, s[axis])
            s[axis] = n
            z = zeros(s, a.dtype.char)
            z[index] = a
            a = z

    if axis != -1:
        a = swapaxes(a, axis, -1)
    r = work_function(a, wsave)
    if axis != -1:
        r = swapaxes(r, axis, -1)

    # As soon as we put wsave back into the cache, another thread could pick it
    # up and start using it, so we must not do this until after we're
    # completely done using it ourselves.
    fft_cache[n].append(wsave)

    return r

def ifftpad(a, n, scale=True):
    """
    Pad the spectrum at high frequencies.

    The padding done by the `ifft` function appends zeros to the end of the
    spectrum which shifts the frequencies and can make the resulting signal
    differ from the non-padded version. This function pads the spectrum
    by putting the zeros in the middle where the highest frequencies are.
    Taking the `ifft` of this padded version result in a signal that is
    an interpolated version of the unpadded signal, which is what is expected.

    Parameters
    ----------
    a : array_like
        Input array, can be complex.
    n : int
        Length of the padded spectrum.
        `n` should be larger than the length of the input.
    scale : bool, optional
        Whether to scale the spectrum or not. The `ifft` function
        divides by the input length which will be the incorrect length
        for a padded spectrum. Setting this parameter will pre-scale
        the spectrum so that dividing by the padded length will be correct.

    Returns
    -------
    out : ndarray
        The spectrum padded to length `n`. Possibly scaled as well.

    Examples
    --------
    >>> spectrum = np.array([0, 1, 2, -3, -2, -1])
    >>> np.fft.ifftpad(spectrum, 10, scale=False)
    array([ 0.,  1.,  2.,  0.,  0.,  0.,  0., -3., -2., -1.])
    """
    spectrum = concatenate((a[:len(a) // 2],
                            zeros(n - len(a)),
                            a[len(a) // 2:]))
    if scale:
        spectrum *= (n / len(a))
    return spectrum

def solve(a, b):
    """Solve the equation a x = b

    Parameters
    ----------
    a : array-like, shape (M, M)
    b : array-like, shape (M,)

    Returns
    -------
    x : array, shape (M,)

    Raises LinAlgError if a is singular or not square

    """
    a, _ = _makearray(a)
    b, wrap = _makearray(b)
    one_eq = len(b.shape) == 1
    if one_eq:
        b = b[:, newaxis]
    _assertRank2(a, b)
    _assertSquareness(a)
    n_eq = a.shape[0]
    n_rhs = b.shape[1]
    if n_eq != b.shape[0]:
        raise LinAlgError, 'Incompatible dimensions'
    t, result_t = _commonType(a, b)
#    lapack_routine = _findLapackRoutine('gesv', t)
    if isComplexType(t):
        lapack_routine = lapack_lite.zgesv
    else:
        lapack_routine = lapack_lite.dgesv
    a, b = _fastCopyAndTranspose(t, a, b)
    pivots = zeros(n_eq, fortran_int)
    results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
    if results['info'] > 0:
        raise LinAlgError, 'Singular matrix'
    if one_eq:
        return wrap(b.ravel().astype(result_t))
    else:
        return wrap(b.transpose().astype(result_t))

def det(a):
    "The determinant of the 2-d array a"
    a = asarray(a)
    _assertRank2(a)
    _assertSquareness(a)
    t, result_t = _commonType(a)
    a = _fastCopyAndTranspose(t, a)
    n = a.shape[0]
    if isComplexType(t):
        lapack_routine = lapack_lite.zgetrf
    else:
        lapack_routine = lapack_lite.dgetrf
    pivots = zeros((n,), fortran_int)
    results = lapack_routine(n, n, a, n, pivots, 0)
    info = results['info']
    if (info < 0):
        raise TypeError, "Illegal input to Fortran routine"
    elif (info > 0):
        return 0.0
    sign = add.reduce(pivots != arange(1, n+1)) % 2
    return (1.-2.*sign)*multiply.reduce(diagonal(a), axis=-1)

def eigh(a, UPLO='L'):
    """
    Eigenvalues and eigenvectors of a Hermitian or real symmetric matrix.

    Parameters
    ----------
    a : array_like, shape (M, M)
        A complex Hermitian or symmetric real matrix.
    UPLO : {'L', 'U'}, optional
        Specifies whether the calculation is done with data from the
        lower triangular part of `a` ('L', default) or the upper triangular
        part ('U').

    Returns
    -------
    w : ndarray, shape (M,)
        The eigenvalues. The eigenvalues are not necessarily ordered.
    v : ndarray, shape (M, M)
        The normalized eigenvector corresponding to the eigenvalue w[i] is
        the column v[:,i].

    Raises
    ------
    LinAlgError
        If the eigenvalue computation does not converge.

    See Also
    --------
    eigvalsh : eigenvalues of symmetric or Hemitiean arrays.
    eig : eigenvalues and right eigenvectors for non-symmetric arrays
    eigvals : eigenvalues of non-symmetric array.

    Notes
    -----
    A simple interface to the LAPACK routines dsyevd and zheevd that compute
    the eigenvalues and eigenvectors of real symmetric and complex Hermitian
    arrays respectively.

    The number w is an eigenvalue of a if there exists a vector v
    satisfying the equation dot(a,v) = w*v. Alternately, if w is a root of
    the characteristic equation det(a - w[i]*I) = 0, where det is the
    determinant and I is the identity matrix. The eigenvalues of real
    symmetric or complex Hermitean matrices are always real. The array v
    of eigenvectors is unitary and a, w, and v satisfy the equation
    dot(a,v[i]) = w[i]*v[:,i].

    """
    a, wrap = _makearray(a)
    _assertRank2(a)
    _assertSquareness(a)
    t, result_t = _commonType(a)
    real_t = _linalgRealType(t)
    a = _fastCopyAndTranspose(t, a)
    n = a.shape[0]
    liwork = 5*n+3
    iwork = zeros((liwork,), fortran_int)
    if isComplexType(t):
        lapack_routine = lapack_lite.zheevd
        w = zeros((n,), real_t)
        lwork = 1
        work = zeros((lwork,), t)
        lrwork = 1
        rwork = zeros((lrwork,), real_t)
        results = lapack_routine('V', UPLO, n, a, n, w, work, -1,
                                 rwork, -1, iwork, liwork,  0)
        lwork = int(abs(work[0]))
        work = zeros((lwork,), t)
        lrwork = int(rwork[0])
        rwork = zeros((lrwork,), real_t)
        results = lapack_routine('V', UPLO, n, a, n, w, work, lwork,
                                 rwork, lrwork, iwork, liwork,  0)
    else:
        lapack_routine = lapack_lite.dsyevd
        w = zeros((n,), t)
        lwork = 1
        work = zeros((lwork,), t)
        results = lapack_routine('V', UPLO, n, a, n, w, work, -1,
                iwork, liwork, 0)
        lwork = int(work[0])
        work = zeros((lwork,), t)
        results = lapack_routine('V', UPLO, n, a, n, w, work, lwork,
                iwork, liwork, 0)
    if results['info'] > 0:
        raise LinAlgError, 'Eigenvalues did not converge'
    at = a.transpose().astype(result_t)
    return w.astype(_realType(result_t)), wrap(at)

def eig(a):
    """
    Compute eigenvalues and right eigenvectors of an array.

    Parameters
    ----------
    a : array_like, shape (M, M)
        A complex or real 2-D array.

    Returns
    -------
    w : ndarray, shape (M,)
        The eigenvalues, each repeated according to its multiplicity.
        The eigenvalues are not necessarily ordered, nor are they
        necessarily real for real matrices.
    v : ndarray, shape (M, M)
        The normalized eigenvector corresponding to the eigenvalue ``w[i]`` is
        the column ``v[:,i]``.

    Raises
    ------
    LinAlgError
        If the eigenvalue computation does not converge.

    See Also
    --------
    eigvalsh : eigenvalues of symmetric or Hemitiean arrays.
    eig : eigenvalues and right eigenvectors for non-symmetric arrays
    eigvals : eigenvalues of non-symmetric array.

    Notes
    -----
    This is a simple interface to the LAPACK routines dgeev and zgeev
    that compute the eigenvalues and eigenvectors of general real and
    complex arrays respectively.

    The number `w` is an eigenvalue of a if there exists a vector `v`
    satisfying the equation ``dot(a,v) = w*v``. Alternately, if `w` is
    a root of the characteristic equation ``det(a - w[i]*I) = 0``, where
    `det` is the determinant and `I` is the identity matrix. The arrays
    `a`, `w`, and `v` satisfy the equation ``dot(a,v[i]) = w[i]*v[:,i]``.

    The array `v` of eigenvectors may not be of maximum rank, that is, some
    of the columns may be dependent, although roundoff error may
    obscure that fact. If the eigenvalues are all different, then theoretically
    the eigenvectors are independent. Likewise, the matrix of eigenvectors
    is unitary if the matrix `a` is normal, i.e., if
    ``dot(a, a.H) = dot(a.H, a)``.

    The left and right eigenvectors are not necessarily the (Hermitian)
    transposes of each other.

    """
    a, wrap = _makearray(a)
    _assertRank2(a)
    _assertSquareness(a)
    _assertFinite(a)
    a, t, result_t = _convertarray(a) # convert to double or cdouble type
    real_t = _linalgRealType(t)
    n = a.shape[0]
    dummy = zeros((1,), t)
    if isComplexType(t):
        # Complex routines take different arguments
        lapack_routine = lapack_lite.zgeev
        w = zeros((n,), t)
        v = zeros((n, n), t)
        lwork = 1
        work = zeros((lwork,), t)
        rwork = zeros((2*n,), real_t)
        results = lapack_routine('N', 'V', n, a, n, w,
                                 dummy, 1, v, n, work, -1, rwork, 0)
        lwork = int(abs(work[0]))
        work = zeros((lwork,), t)
        results = lapack_routine('N', 'V', n, a, n, w,
                                 dummy, 1, v, n, work, lwork, rwork, 0)
    else:
        lapack_routine = lapack_lite.dgeev
        wr = zeros((n,), t)
        wi = zeros((n,), t)
        vr = zeros((n, n), t)
        lwork = 1
        work = zeros((lwork,), t)
        results = lapack_routine('N', 'V', n, a, n, wr, wi,
                                  dummy, 1, vr, n, work, -1, 0)
        lwork = int(work[0])
        work = zeros((lwork,), t)
        results = lapack_routine('N', 'V', n, a, n, wr, wi,
                                  dummy, 1, vr, n, work, lwork, 0)
        if all(wi == 0.0):
            w = wr
            v = vr
            result_t = _realType(result_t)
        else:
            w = wr+1j*wi
            v = array(vr, w.dtype)
            ind = flatnonzero(wi != 0.0)      # indices of complex e-vals
            for i in range(len(ind)/2):
                v[ind[2*i]] = vr[ind[2*i]] + 1j*vr[ind[2*i+1]]
                v[ind[2*i+1]] = vr[ind[2*i]] - 1j*vr[ind[2*i+1]]
            result_t = _complexType(result_t)
#.........这里部分代码省略.........

def eigvalsh(a, UPLO='L'):
    """
    Compute the eigenvalues of a Hermitean or real symmetric matrix.

    Parameters
    ----------
    a : array_like, shape (M, M)
        A complex or real matrix whose eigenvalues and eigenvectors
        will be computed.
    UPLO : {'L', 'U'}, optional
        Specifies whether the calculation is done with data from the
        lower triangular part of `a` ('L', default) or the upper triangular
        part ('U').

    Returns
    -------
    w : ndarray, shape (M,)
        The eigenvalues, each repeated according to its multiplicity.
        They are not necessarily ordered.

    Raises
    ------
    LinAlgError
        If the eigenvalue computation does not converge.

    See Also
    --------
    eigh : eigenvalues and eigenvectors of symmetric/Hermitean arrays.
    eigvals : eigenvalues of general real or complex arrays.
    eig : eigenvalues and eigenvectors of general real or complex arrays.

    Notes
    -----
    This is a simple interface to the LAPACK routines dsyevd and
    zheevd that sets the flags to return only the eigenvalues of real
    symmetric and complex Hermetian arrays respectively.

    The number w is an eigenvalue of a if there exists a vector v
    satisfying the equation dot(a,v) = w*v. Alternately, if w is a root of
    the characteristic equation det(a - w[i]*I) = 0, where det is the
    determinant and I is the identity matrix.

    """
    a, wrap = _makearray(a)
    _assertRank2(a)
    _assertSquareness(a)
    t, result_t = _commonType(a)
    real_t = _linalgRealType(t)
    a = _fastCopyAndTranspose(t, a)
    n = a.shape[0]
    liwork = 5*n+3
    iwork = zeros((liwork,), fortran_int)
    if isComplexType(t):
        lapack_routine = lapack_lite.zheevd
        w = zeros((n,), real_t)
        lwork = 1
        work = zeros((lwork,), t)
        lrwork = 1
        rwork = zeros((lrwork,), real_t)
        results = lapack_routine('N', UPLO, n, a, n, w, work, -1,
                                 rwork, -1, iwork, liwork,  0)
        lwork = int(abs(work[0]))
        work = zeros((lwork,), t)
        lrwork = int(rwork[0])
        rwork = zeros((lrwork,), real_t)
        results = lapack_routine('N', UPLO, n, a, n, w, work, lwork,
                                rwork, lrwork, iwork, liwork,  0)
    else:
        lapack_routine = lapack_lite.dsyevd
        w = zeros((n,), t)
        lwork = 1
        work = zeros((lwork,), t)
        results = lapack_routine('N', UPLO, n, a, n, w, work, -1,
                                 iwork, liwork, 0)
        lwork = int(work[0])
        work = zeros((lwork,), t)
        results = lapack_routine('N', UPLO, n, a, n, w, work, lwork,
                                 iwork, liwork, 0)
    if results['info'] > 0:
        raise LinAlgError, 'Eigenvalues did not converge'
    return w.astype(result_t)