(npts, spacing=1.0, beta=0.1)
return (1.0/spacing**2)*A
Construct a Gauge Laplacian from Quantum Chromodynamics for regular 2D grids.

    Note that this function is not written efficiently, but should be
    fine for N x N grids where N is in the low hundreds.

    Parameters
    ----------
    npts : int
        number of pts in x and y directions

    spacing : float
        grid spacing between points

    beta : float
        temperature
        Note that if beta=0, then we get the typical 5pt Laplacian stencil

    Returns
    -------
    A : csr matrix
        A is Hermitian positive definite for beta > 0.0
        A is Symmetric semi-definite for beta = 0.0

    Examples
    --------
    >>> from pyamg.gallery import gauge_laplacian
    >>> A = gauge_laplacian(10)

    References
    ----------
    .. [1] MacLachlan, S. and Oosterlee, C.,
       "Algebraic Multigrid Solvers for Complex-Valued Matrices",
       Vol. 30, SIAM J. Sci. Comp, 2008
Construct a Gauge Laplacian from Quantum Chromodynamics for regular 2D grids.
Construct a Gauge Laplacian from Quantum Chromodynamics for regular 2D grids .
def gauge_laplacian(npts, spacing=1.0, beta=0.1):
    """Construct a Gauge Laplacian from Quantum Chromodynamics for regular 2D grids.

    Note that this function is not written efficiently, but should be
    fine for N x N grids where N is in the low hundreds.

    Parameters
    ----------
    npts : int
        number of pts in x and y directions

    spacing : float
        grid spacing between points

    beta : float
        temperature
        Note that if beta=0, then we get the typical 5pt Laplacian stencil

    Returns
    -------
    A : csr matrix
        A is Hermitian positive definite for beta > 0.0
        A is Symmetric semi-definite for beta = 0.0

    Examples
    --------
    >>> from pyamg.gallery import gauge_laplacian
    >>> A = gauge_laplacian(10)

    References
    ----------
    .. [1] MacLachlan, S. and Oosterlee, C.,
       "Algebraic Multigrid Solvers for Complex-Valued Matrices",
       Vol. 30, SIAM J. Sci. Comp, 2008

    """
    # The gauge Laplacian has the same sparsity structure as a normal
    # Laplacian, so we start out with a Poisson Operator
    N = npts
    A = poisson((N, N), format='coo', dtype=complex)

    # alpha is a random function of a point's integer position
    # on a 1-D grid along the x or y direction.  e.g. the first
    # point at (0,0) would be evaluate at alpha_*[0], while the
    # last point at (N*spacing, N*spacing) would evaluate at alpha_*[-1]
    alpha_x = 1.0j * 2.0 * np.pi * beta * np.random.randn(N*N)
    alpha_y = 1.0j * 2.0 * np.pi * beta * np.random.randn(N*N)

    # Replace off diagonals of A
    for i in range(A.nnz):
        r = A.row[i]
        c = A.col[i]
        diff = np.abs(r - c)
        index = min(r, c)
        if r > c:
            s = -1.0
        else:
            s = 1.0
        if diff == 1:
            # differencing in the x-direction
            A.data[i] = -1.0 * np.exp(s * alpha_x[index])
        if diff == N:
            # differencing in the y-direction
            A.data[i] = -1.0 * np.exp(s * alpha_y[index])

    # Handle periodic BCs
    alpha_x = 1.0j * 2.0 * np.pi * beta * np.random.randn(N*N)
    alpha_y = 1.0j * 2.0 * np.pi * beta * np.random.randn(N*N)
    new_r = []
    new_c = []
    new_data = []
    new_diff = []
    for i in range(0, N):
        new_r.append(i)
        new_c.append(i + N*N - N)
        new_diff.append(N)

    for i in range(N*N - N, N*N):
        new_r.append(i)
        new_c.append(i - N*N + N)
        new_diff.append(N)

    for i in range(0, N*N-1, N):
        new_r.append(i)
        new_c.append(i + N - 1)
        new_diff.append(1)

    for i in range(N-1, N*N, N):
        new_r.append(i)
        new_c.append(i - N + 1)
        new_diff.append(1)

    for i, (r, c) in enumerate(zip(new_r, new_c)):
        diff = new_diff[i]
        index = min(r, c)
        if r > c:
            s = -1.0
        else:
            s = 1.0
        if diff == 1:
            # differencing in the x-direction
            new_data.append(-1.0 * np.exp(s * alpha_x[index]))
        if diff == N:
            # differencing in the y-direction
            new_data.append(-1.0 * np.exp(s * alpha_y[index]))

    # Construct Final Matrix
    data = np.hstack((A.data, np.array(new_data)))
    row = np.hstack((A.row, np.array(new_r)))
    col = np.hstack((A.col, np.array(new_c)))
    A = sp.sparse.coo_matrix((data, (row, col)), shape=(N*N, N*N)).tocsr()

    return (1.0/spacing**2)*A
def gauge_laplacian ( npts , spacing = 1.0 , beta = 0.1 ) : # The gauge Laplacian has the same sparsity structure as a normal # Laplacian, so we start out with a Poisson Operator N = npts A = poisson ( ( N , N ) , format = 'coo' , dtype = complex ) # alpha is a random function of a point's integer position # on a 1-D grid along the x or y direction.  e.g. the first # point at (0,0) would be evaluate at alpha_*[0], while the # last point at (N*spacing, N*spacing) would evaluate at alpha_*[-1] alpha_x = 1.0j * 2.0 * np . pi * beta * np . random . randn ( N * N ) alpha_y = 1.0j * 2.0 * np . pi * beta * np . random . randn ( N * N ) # Replace off diagonals of A for i in range ( A . nnz ) : r = A . row [ i ] c = A . col [ i ] diff = np . abs ( r - c ) index = min ( r , c ) if r > c : s = - 1.0 else : s = 1.0 if diff == 1 : # differencing in the x-direction A . data [ i ] = - 1.0 * np . exp ( s * alpha_x [ index ] ) if diff == N : # differencing in the y-direction A . data [ i ] = - 1.0 * np . exp ( s * alpha_y [ index ] ) # Handle periodic BCs alpha_x = 1.0j * 2.0 * np . pi * beta * np . random . randn ( N * N ) alpha_y = 1.0j * 2.0 * np . pi * beta * np . random . randn ( N * N ) new_r = [ ] new_c = [ ] new_data = [ ] new_diff = [ ] for i in range ( 0 , N ) : new_r . append ( i ) new_c . append ( i + N * N - N ) new_diff . append ( N ) for i in range ( N * N - N , N * N ) : new_r . append ( i ) new_c . append ( i - N * N + N ) new_diff . append ( N ) for i in range ( 0 , N * N - 1 , N ) : new_r . append ( i ) new_c . append ( i + N - 1 ) new_diff . append ( 1 ) for i in range ( N - 1 , N * N , N ) : new_r . append ( i ) new_c . append ( i - N + 1 ) new_diff . append ( 1 ) for i , ( r , c ) in enumerate ( zip ( new_r , new_c ) ) : diff = new_diff [ i ] index = min ( r , c ) if r > c : s = - 1.0 else : s = 1.0 if diff == 1 : # differencing in the x-direction new_data . append ( - 1.0 * np . exp ( s * alpha_x [ index ] ) ) if diff == N : # differencing in the y-direction new_data . append ( - 1.0 * np . exp ( s * alpha_y [ index ] ) ) # Construct Final Matrix data = np . hstack ( ( A . data , np . array ( new_data ) ) ) row = np . hstack ( ( A . row , np . array ( new_r ) ) ) col = np . hstack ( ( A . col , np . array ( new_c ) ) ) A = sp . sparse . coo_matrix ( ( data , ( row , col ) ) , shape = ( N * N , N * N ) ) . tocsr ( ) return ( 1.0 / spacing ** 2 ) * A