# Direct solvers

    #!python

    """
    PyViennaCL provides direct solvers for dense triangular linear systems.
    The API is documented in ``help(pyviennacl.linalg.solve)``. In particular,
    the solver to use is determined by the tag instance supplied to the ``solve``
    function.

    LU factorisation is planned for a later release.

    Here, we demonstrate the use of the direct triangular solver for an
    upper triangular system; other forms and solvers are supported by the choice
    of a different solver tag class.
    """

    import pyviennacl as p
    import numpy as np
    import random

    # We want a square N x N system.
    N = 5

    # Create a NumPy matrix with float32 precision to hold the data on the host.
    # Firstly, we create an empty matrix, then fill the upper triangle with values.
    A = np.zeros((N, N), dtype = np.float32)

    for i in range(N):
        for j in range(N):
            if j >= i:
                A[i, j] = np.float32(random.randint(0,1000) / 100.0)

    # Transfer the system matrix to the compute device
    A = p.Matrix(A)

    print("A is\n%s" % A)

    # Create a right-hand-side vector on the host with random elements
    # and transfer it to the compute device
    b = p.Vector(np.random.rand(N).astype(np.float32))

    print("b is %s" % b)

    # Solve the system; note the choice of tag to denote an upper triangular system
    x = p.solve(A, b, p.upper_tag())

    # Copy the solution from the device to host and display it
    print("Solution of Ax = b for x:\n%s" % x)