python-sparse 1.1.1-1 / usr / lib / python2.7 / dist-packages / pysparse / pysparseSuperLU.py

"""
A framework for solving sparse linear systems of equations using an LU
factorization, by means of the supernodal sparse LU factorization package
SuperLU ([DEGLL99]_, [DGL99]_, [LD03]_).

This package is appropriate for factorizing sparse square unsymmetric or
rectangular matrices.

See [SLU]_ for more information.

**References:**

.. [DEGLL99] J. W. Demmel, S. C. Eisenstat, J. R. Gilbert, X. S. Li and
             J. W. H. Liu, *A supernodal approach to sparse partial pivoting*,
             SIAM Journal on Matrix Analysis and Applications **20**\ (3),
             pp. 720-755, 1999.
.. [DGL99] J. W. Demmel, J. R. Gilbert and X. S. Li,
           *An Asynchronous Parallel Supernodal Algorithm for Sparse Gaussian
           Elimination*, SIAM Journal on Matrix Analysis and Applications
           **20**\ (4), pp. 915-952, 1999.
.. [LD03] X. S. Li and J. W. Demmel, *SuperLU_DIST: A Scalable
          Distributed-Memory Sparse Direct Solver for Unsymmetric Linear
          Systems*, ACM Transactions on Mathematical Software **29**\ (2),
          pp. 110-140, 2003.
.. [SLU] http://crd.lbl.gov/~xiaoye/SuperLU

"""

# To look into:
#  - allow other data types

__docformat__ = 'restructuredtext'

import pysparseMatrix as psm
import numpy
import resource

from directSolver import PysparseDirectSolver
from pysparse  import superlu

def cputime():
    return resource.getrusage(resource.RUSAGE_SELF)[0]

class PysparseSuperLUSolver( PysparseDirectSolver ):
    """
    `PysparseSuperLUSolver` is a wrapper class around the SuperLu library for
    the factorization of full-rank n-by-m matrices. Only matrices with real
    coefficients are currently supported.

    :parameters:

       :A: The matrix to be factorized, supplied as a PysparseMatrix instance.

    :keywords:

       :symmetric: a boolean indicating that the user wishes to use symmetric
                   mode. In symmetric mode, ``permc_spec=2`` must be chosen and
                   ``diag_pivot_thresh`` must be small, e.g., 0.0 or 0.1. Since
                   the value of ``diag_pivot_thresh`` is up to the user, setting
                   ``symmetric`` to ``True`` does *not* automatically set
                   ``permc_spec`` and ``diag_pivot_thresh`` to appropriate
                   values.

       :diag_pivot_thresh: a float value between 0 and 1 representing the
                           threshold for partial pivoting (0 = no pivoting,
                           1 = always perform partial pivoting). Default: 1.0.

       :drop_tol: the value of a drop tolerance, between 0 and 1, if an
                  incomplete factorization is desired (0 = exact factorization).
                  This keyword does not exist if using SuperLU version 2.0 and
                  below. In more recent version of SuperLU, the keyword is
                  accepted but has no effect. Default: 0.0

       :relax: an integer controling the degree of relaxing supernodes.
               Default: 1.

       :panel_size: an integer specifying the maximum number of columns to form
                    a panel. Default: 10.

       :permc_spec: an integer specifying the ordering strategy used during the
                    factorization.

                    0. natural ordering,
                    1. MMD applied to the structure of
                       :math:`\mathbf{A}^T \mathbf{A}`
                    2. MMD applied to the structure of
                       :math:`\mathbf{A}^T + \mathbf{A}`
                    3. COLAMD.

                    Default: 2.

    .. attribute:: LU

       A :class:`superlu_context` object encapsulating the factorization.

    .. attribute:: sol

       The solution of the linear system after a call to :meth:`solve`.

    .. attribute:: factorizationTime

       The CPU time to perform the factorization.

    .. attribute:: solutionTime

       The CPU time to perform the forward and backward sweeps.

    .. attribute:: lunz

       The number of nonzero elements in the factors L and U together after a
       call to :meth:`fetch_lunz`.
    """
    def __init__(self, A, **kwargs):
        PysparseDirectSolver.__init__(self, A, **kwargs)

        self.type = numpy.float
        self.nrow, self.ncol = A.getShape()
        t = cputime()
        self.LU = superlu.factorize(A.matrix.to_csr(), **kwargs)
        self.factorizationTime = cputime() - t
        self.solutionTime = 0.0
        self.sol = None
        self.L = self.U = None
        return

    def solve(self, rhs, transpose = False):
        """
        Solve the linear system  ``A x = rhs``, where ``A`` is the input matrix
        and ``rhs`` is a Numpy vector of appropriate dimension. The result is
        placed in the :attr:`sol` member of the class instance.

        If the optional argument ``transpose`` is ``True``, the transpose system
        ``A^T x = rhs`` is solved.
        """
        if self.sol is None: self.sol = numpy.empty(self.ncol, self.type)
        transp = 'N'
        if transpose: transp = 'T'
        t = cputime()
        self.LU.solve(rhs, self.sol, transp)
        self.solutionTime = cputime() - t
        return

    def fetch_lunz(self):
        """
        Retrieve the number of nonzeros in the factors L and U together. The
        result is stored in the member :attr:`lunz` of the class instance.
        """
        self.lunz = self.LU.nnz

    def fetch_factors(self):
        """
        Not yet available.
        """
        raise NotImplementedError