diff --git a/.gitignore b/.gitignore index c2e9891..4f76999 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,6 @@ *.pyc +.coverage +coverage.xml *.so build/ dist/ diff --git a/Makefile b/Makefile index 08a96fe..88a3dd5 100644 --- a/Makefile +++ b/Makefile @@ -4,8 +4,7 @@ build: python setup.py build_ext --inplace coverage: - nosetests --logging-level=INFO --with-coverage --cover-package=pymatsolver --cover-html - open cover/index.html + pytest --cov-config=.coveragerc --cov-report=xml --cov=pymatsolver -s -v lint: pylint --output-format=html pymatsolver > pylint.html @@ -14,7 +13,7 @@ graphs: pyreverse -my -A -o pdf -p pymatsolver pymatsolver/**.py pymatsolver/**/**.py tests: - nosetests --logging-level=INFO + pytest docs: cd docs;make html diff --git a/README.rst b/README.rst index b321ee4..74557b3 100644 --- a/README.rst +++ b/README.rst @@ -9,7 +9,7 @@ pymatsolver :target: https://github.com/simpeg/pymatsolver/blob/master/LICENSE :alt: MIT license. -.. image:: https://img.shields.io/travis/rowanc1/pymatsolver.svg +.. image:: https://img.shields.io/travis/simpeg/pymatsolver.svg :target: https://travis-ci.org/simpeg/pymatsolver :alt: Travis CI build status @@ -64,9 +64,9 @@ From a clean install on Ubuntu: ./miniconda.sh -b export PATH=/root/anaconda/bin:/root/miniconda/bin:$PATH conda update --yes conda - conda install --yes numpy scipy matplotlib cython ipython nose + conda install --yes numpy scipy matplotlib cython ipython pytest coverage - git clone https://github.com/rowanc1/pymatsolver.git + git clone https://github.com/simpeg/pymatsolver.git cd pymatsolver make mumps @@ -79,7 +79,7 @@ This assumes that you have Brew and some python installed (numpy, scipy): brew install mumps --with-scotch5 --without-mpi - git clone https://github.com/rowanc1/pymatsolver.git + git clone https://github.com/simpeg/pymatsolver.git cd pymatsolver make mumps_mac diff --git a/pymatsolver/__init__.py b/pymatsolver/__init__.py index ead34b0..a155254 100644 --- a/pymatsolver/__init__.py +++ b/pymatsolver/__init__.py @@ -71,31 +71,31 @@ AvailableSolvers['Pardiso'] = True PardisoSolver = Pardiso # backwards compatibility except ImportError: - SolverHelp['Pardiso'] = """Pardiso is not working + SolverHelp['Pardiso'] = """Pardiso is not working. Ensure that you have pydiso installed, which may also require Python to be installed through conda. """ -# try: -# from pymatsolver.mumps import Mumps -# AvailableSolvers['Mumps'] = True -# MumpsSolver = Mumps # backwards compatibility -# except Exception: -# SolverHelp['Mumps'] = """Mumps is not working. -# -# Ensure that you have Mumps installed, and know where the path to it is. -# Try something like: -# cd pymatsolver/mumps -# make -# -# When that doesn't work, you may need to edit the make file, to modify the -# path to the mumps installation, or any other compiler flags. -# If you find a good way of doing it, please share. -# -# brew install mumps --with-scotch5 --without-mpi -# mpicc --showme -# """ +try: + from pymatsolver.mumps import Mumps + AvailableSolvers['Mumps'] = True + MumpsSolver = Mumps # backwards compatibility +except ImportError: + SolverHelp['Mumps'] = """Mumps is not working. + +Ensure that you have Mumps installed, and know where the path to it is. +Try something like: + cd pymatsolver/mumps + make + +When that doesn't work, you may need to edit the make file, to modify the +path to the mumps installation, or any other compiler flags. +If you find a good way of doing it, please share. + +brew install mumps --with-scotch5 --without-mpi +mpicc --showme +""" __version__ = '0.2.0' __author__ = 'Rowan Cockett' diff --git a/pymatsolver/mumps/Makefile b/pymatsolver/mumps/Makefile index e6501e4..c29b41a 100644 --- a/pymatsolver/mumps/Makefile +++ b/pymatsolver/mumps/Makefile @@ -43,4 +43,4 @@ clean: rm -rf *.dSYM test: - cd ..;nosetests + cd ../..;pytest diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index 6826e25..524f888 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -4,162 +4,254 @@ from __future__ import absolute_import import gc +import warnings + +import numpy as np + from pymatsolver.solvers import Base from . import MumpsInterface as _MUMPSINT -_mumpsErrors = { - -1: "An error occurred on processor INFO(2).", - -2: "NZ is out of range. INFO(2)=NZ.", - -3: "MUMPS was called with an invalid value for JOB. This may happen for example if the analysis (JOB=1) was not performed before the factorization (JOB=2), or the factorization was not performed before the solve (JOB=3), or the initialization phase (JOB=-1) was performed a second time on an instance not freed (JOB=-2). See description of JOB in Section 3. This error also occurs if JOB does not contain the same value on all processes on entry to MUMPS.", - -4: "Error in user-provided permutation array PERM IN at position INFO(2). This error may only occuron the host.", - -5: "Problem of REAL or COMPLEX workspace allocation of size INFO(2) during analysis.", - -6: "Matrix is singular in structure. INFO(2) holds the structural rank.", - -7: "Problem of INTEGER workspace allocation of size INFO(2) during analysis.", - -8: "Main internal integer workarray IS too small for factorization. This may happen, for example, if numerical pivoting leads to significantly more fill-in than was predicted by the analysis. The user should increase the value of ICNTL(14) before calling the factorization again (JOB=2).", - -9: "Main internal real/complex workarray S too small. If INFO(2) is positive, then the number of entries that are missing in S at the moment when the error is raised is available in INFO(2). If INFO(2) is negative, then its absolute value should be multiplied by 1 million. If an error -9 occurs, the user should increase the value of ICNTL(14) before calling the factorization (JOB=2) again, except if ICNTL(23) is provided, in which case ICNTL(23) should be increased.", - -10: "Numerically singular matrix.", - -11: "Internal real/complex workarray S too small for solution. Please contact us. If INFO(2) is positive, then the number of entries that are missing in S at the moment when the error is raised is available in INFO(2).", - -12: "Internal real/complex workarray S too small for iterative refinement. Please contact us.", - -13: "An error occurred in a Fortran ALLOCATE statement. The size that the package requested is available in INFO(2). If INFO(2) is negative, then the size that the package requested is obtained by multiplying the absolute value of INFO(2) by 1 million.", - -14: "Internal integer workarray IS too small for solution. See error INFO(1) = -8.", - -15: "Integer workarray IS too small for iterative refinement and/or error analysis. See error INFO(1) = -8.", - -16: "N is out of range. INFO(2)=N.", - -17: "The internal send buffer that was allocated dynamically by MUMPS on the processor is too small. The user should increase the value of ICNTL(14) before calling MUMPS again.", - -20: "The internal reception buffer that was allocated dynamically by MUMPS is too small. INFO(2) holds the minimum size of the reception buffer required (in bytes). The user should increase the value of ICNTL(14) before calling MUMPS again.", - -21: "Value of PAR=0 is not allowed because only one processor is available; Running node mode (the host is not a slave processor itself) requires at least two processors. The user should either set PAR to 1 or increase the number of processors.", - -22: "A pointer array is provided by the user that is either not associated, or has insufficient size, or is associated and should not be associated (for example, RHS on non-host processors).", - -23: "MPI was not initialized by the user prior to a call to MUMPS with JOB = -1.", - -24: "NELT is out of range. INFO(2)=NELT.", - -25: "A problem has occurred in the initialization of the BLACS. This may be because you are using a vendor's BLACS. Try using a BLACS version from netlib instead.", - -26: "LRHS is out of range. INFO(2)=LRHS.", - -27: "NZ RHS and IRHS PTR(NRHS+1) do not match. INFO(2) = IRHS PTR(NRHS+1).", - -28: "IRHS PTR(1) is not equal to 1. INFO(2) = IRHS PTR(1).", - -29: "LSOL loc is smaller than INFO(23). INFO(2)=LSOL loc.", - -30: "SCHUR LLD is out of range. INFO(2) = SCHUR LLD.", - -31: "A 2D block cyclic symmetric (SYM=1 or 2) Schur complement is required with the option ICNTL(19)=3, but the user has provided a process grid that does not satisfy the constraint MBLOCK=NBLOCK. INFO(2)=MBLOCK-NBLOCK.", - -32: "Incompatible values of NRHS and ICNTL(25). Either ICNTL(25) was set to -1 and NRHS is different from INFOG(28); or ICNTL(25) was set to i, 1 <= i <= INFOG(28) and NRHS is different from 1. Value of NRHS is stored in INFO(2).", - -33: "ICNTL(26) was asked for during solve phase but the Schur complement was not asked for at the analysis phase (ICNTL(19)). INFO(2)=ICNTL(26).", - -34: "LREDRHS is out of range. INFO(2)=LREDRHS.", - -35: "This error is raised when the expansion phase is called (ICNTL(26) = 2) but reduction phase (ICNTL(26)=1) was not called before. INFO(2) contains the value of ICNTL(26).", - -36: "Incompatible values of ICNTL(25) and INFOG(28). The value of ICNTL(25) is stored in INFO(2).", - -37: "Value of ICNTL(25) incompatible with some other parameter. with SYM or ICNTL(xx). If INFO(2)=0 then ICNTL(25) is incompatible with SYM: in current version, the null space basis functionality is not available for unsymmetric matrices (SYM=0). Otherwise, ICNTL(25) is incompatible with ICNTL(xx), and the index xx is stored in INFO(2).", - -38: "Parallel analysis was set (i.e., ICNTL(28)=2) but PT-SCOTCH or ParMetis were not provided.", - -39: "Incompatible values for ICNTL(28) and ICNTL(5) and/or ICNTL(19) and/or ICNTL(6). Parallel analysis is not possible in the cases where the matrix is unassembled and/or a Schur complement is requested and/or a maximum transversal is requested on the matrix.", - -40: "The matrix was indicated to be positive definite (SYM=1) by the user but a negative or null pivot was encountered during the processing of the root by ScaLAPACK. SYM=2 should be used.", - -44: "The solve phase (JOB=3) cannot be performed because the factors or part of the factors are not available. INFO(2) contains the value of ICNTL(31).", - -45: "NRHS <= 0. INFO(2) contains the value of NRHS.", - -46: "NZ RHS <= 0. This is currently not allowed in case of reduced right-hand-side (ICNTL(26)=1) and in case entries of A-1 are requested (ICNTL(30)=1). INFO(2) contains the value of NZ RHS.", - -47: "Entries of A-1 were requested during the solve phase (JOB=3, ICNTL(30)=1) but the constraint NRHS=N is not respected. The value of NRHS is provided in INFO(2).", - -48: "A-1 Incompatible values of ICNTL(30) and ICNTL(xx). xx is stored in INFO(2).", - -49: "SIZE SCHUR has an incorrect value (SIZE SCHUR < 0 or SIZE SCHUR >=N, or SIZE SCHUR was modified on the host since the analysis phase. The value of SIZE SCHUR is provided in INFO(2).", - -90: "Error in out-of-core management. See the error message returned on output unit ICNTL(1) for more information.", - +1: "Index (in IRN or JCN) out of range. Action taken by subroutine is to ignore any such entries and continue. INFO(2) is set to the number of faulty entries. Details of the first ten are printed on unit ICNTL(2).", - +2: "During error analysis the max-norm of the computed solution was found to be zero.", - +4: "User data JCN has been modified (internally) by the solver.", - +8: "Warning return from the iterative refinement routine. More than ICNTL(10) iterations are required.", -} - - -class _Pointer(object): - """Gets an int and a destroy call that gets called on garbage collection. - - There can be multiple Solvers around the place that are pointing to the same factor in memory. - - This ensures that when all references are removed, there is automatic garbage collection. - - You can always clean the Solver, and it will explicitly call the destroy method. + +def _mumps_message_from_exit_code(ierr): + """Convert a MUMPS routine's exit code to a message. + + See section 8 of the MUMPS user guide for a full catalog of these messages. + """ + error_messages = { + -1: "An error occurred on processor INFO(2).", + -2: "NNZ/NZ, NNZ loc/NZ loc or P NNZ loc/P NZ loc are out of range. INFO(2)=NNZ/NZ, NNZ loc/NZ loc or P NNZ loc/P NZ loc.", + -3: "MUMPS was called with an invalid value for JOB. This may happen if the analysis (JOB=1) was not performed (or failed) before the factorization (JOB=2), or the factorization was not performed (or failed) before the solve (JOB=3), or the initialization phase (JOB=-1) was performed a second time on an instance not freed (JOB=-2). See description of JOB in Section 4. This error also occurs if JOB does not contain the same value on all processes on entry to MUMPS. INFO(2) is then set to the local value of JOB.", + -4: "Error in user-provided permutation array PERM IN at position INFO(2). This error may only occur on the host.", + -5: 'Problem of real workspace allocation of size INFO(2) during analysis. The unit for INFO(2) is the number of real values (single precision for SMUMPS/CMUMPS, double precision for DMUMPS/ZMUMPS), in the Fortran "ALLOCATE" statement that did not succeed. If INFO(2) is negative, then its absolute value should be multiplied by 1 million.', + -6: "Matrix is singular in structure. INFO(2) holds the structural rank.", + -7: "Problem of integer workspace allocation of size INFO(2) during analysis. The unit for INFO(2) is the number of integer values that MUMPS tried to allocate in the Fortran ALLOCATE statement that did not succeed. If INFO(2) is negative, then its absolute value should be multiplied by 1 million.", + -8: "Main internal integer workarray IS too small for factorization. This may happen, for example, if numerical pivoting leads to significantly more fill-in than was predicted by the analysis. The user should increase the value of ICNTL(14) before calling the factorization again (JOB=2).", + -9: "The main internal real/complex workarray S is too small. If INFO(2) is positive, then the number of entries that are missing in S at the moment when the error is raised is available in INFO(2). If INFO(2) is negative, then its absolute value should be multiplied by 1 million. If an error -9 occurs, the user should increase the value of ICNTL(14) before calling the factorization (JOB=2) again, except if LWK USER is provided LWK USER should be increased.", + -10: "Numerically singular matrix. INFO(2) holds the number of eliminated pivots.", + -11: "Internal real/complex workarray S or LWK USER too small for solution. If INFO(2) is positive, then the number of entries that are missing in S/LWK USER at the moment when the error is raised is available in INFO(2). If the numerical phases are out-of-core and LWK USER is provided for the solution phase and is smaller than the value provided for the factorization, it should be increased by at least INFO(2). In other cases, please contact us.", + -12: "Internal real/complex workarray S too small for iterative refinement. Please contact us.", + -13: "Problem of workspace allocation of size INFO(2) during the factorization or solve steps. The size that the package tried to allocate with a Fortran ALLOCATE statement is available in INFO(2). If INFO(2) is negative, then the size that the package requested is obtained by multiplying the absolute value of INFO(2) by 1 million. In general, the unit for INFO(2) is the number of scalar entries of the type of the input matrix (real, complex, single or double precision).", + -14: "Internal integer workarray IS too small for solution. See error INFO(1) = -8.", + -15: "Integer workarray IS too small for iterative refinement and/or error analysis. See error INFO(1) = -8.", + -16: "N is out of range. INFO(2)=N.", + -17: "The internal send buffer that was allocated dynamically by MUMPS on the processor is too small. The user should increase the value of ICNTL(14) before calling MUMPS again.", + -18: "The blocking size for multiple RHS (ICNTL(27)) is too large and may lead to an integer overflow. This error may only occurs for very large matrices with large values of ICNTL(27) (e.g., several thousands). INFO(2) provides an estimate of the maximum value of ICNTL(27) that should be used.", + -19: "The maximum allowed size of working memory ICNTL(23) is too small to run the factorization phase and should be increased. If INFO(2) is positive, then the number of entries that are missing at the moment when the error is raised is available in INFO(2). If INFO(2) is negative, then its absolute value should be multiplied by 1 million.", + -20: "The internal reception buffer that was allocated dynamically by MUMPS is too small. Normally, this error is raised on the sender side when detecting that the message to be sent is too large for the reception buffer on the receiver. INFO(2) holds the minimum size of the reception buffer required (in bytes). The user should increase the value of ICNTL(14) before calling MUMPS again.", + -21: "Value of PAR=0 is not allowed because only one processor is available; Running MUMPS in host-node mode (the host is not a slave processor itself) requires at least two processors. The user should either set PAR to 1 or increase the number of processors.", + -22: "A pointer array is provided by the user that is either not associated, or has insufficient size, or is associated and should not be associated (for example, RHS on non-host processors).", + -23: "MPI was not initialized by the user prior to a call to MUMPS with JOB = -1.", + -24: "NELT is out of range. INFO(2)=NELT.", + -25: "A problem has occurred in the initialization of the BLACS. This may be because you are using a vendor's BLACS. Try using a BLACS version from netlib instead.", + -26: "LRHS is out of range. INFO(2)=LRHS.", + -27: "NZ RHS and IRHS PTR(NRHS+1) do not match. INFO(2) = IRHS PTR(NRHS+1).", + -28: "IRHS PTR(1) is not equal to 1. INFO(2) = IRHS PTR(1).", + -29: "LSOL loc is smaller than INFO(23). INFO(2)=LSOL loc.", + -30: "SCHUR LLD is out of range. INFO(2) = SCHUR LLD.", + -31: "A 2D block cyclic symmetric (SYM=1 or 2) Schur complement is required with the option ICNTL(19)=3, but the user has provided a process grid that does not satisfy the constraint MBLOCK=NBLOCK. INFO(2)=MBLOCK-NBLOCK.", + -32: "Incompatible values of NRHS and ICNTL(25). Either ICNTL(25) was set to -1 and NRHS is different from INFOG(28); or ICNTL(25) was set to i, 1 <= i <= INFOG(28) and NRHS is different from 1. Value of NRHS is stored in INFO(2).", + -33: "ICNTL(26) was asked for during solve phase (or during the factorization - see ICNTL(32)) but the Schur complement was not asked for at the analysis phase (ICNTL(19)). INFO(2)=ICNTL(26).", + -34: "LREDRHS is out of range. INFO(2)=LREDRHS.", + -35: "This error is raised when the expansion phase is called (ICNTL(26) = 2) but reduction phase (ICNTL(26)=1) was not called before. This error also occurs in case the reduction phase (ICNTL(26) = 1) is asked for at the solution phase (JOB=3) but the forward elimination was already performed during the factorization phase (JOB=2 and ICNTL(32)=1). INFO(2) contains the value of ICNTL(26).", + -36: "Incompatible values of ICNTL(25) and INFOG(28). The value of ICNTL(25) is stored in INFO(2).", + -37: "Value of ICNTL(25) incompatible with some other parameter. If ICNTL(25) is incompatible with ICNTL(xx), the index xx is stored in INFO(2).", + -38: "Parallel analysis was set (i.e., ICNTL(28)=2) but PT-SCOTCH or ParMetis were not provided.", + -39: "Incompatible values for ICNTL(28) and ICNTL(5) and/or ICNTL(19) and/or ICNTL(6). Parallel analysis is not possible in the cases where the matrix is unassembled and/or a Schur complement is requested and/or a maximum transversal is requested on the matrix.", + -40: "The matrix was indicated to be positive definite (SYM=1) by the user but a negative or null pivot was encountered during the processing of the root by ScaLAPACK. SYM=2 should be used.", + -41: "Incompatible value of LWK USER between factorization and solution phases. This error may only occur when the factorization is in-core (ICNTL(22)=1), in which case both the contents of WK USER and LWK USER should be passed unchanged between the factorization (JOB=2) and solution (JOB=3) phases.", + -42: "ICNTL(32) was set to 1 (forward during factorization), but the value of NRHS on the host processor is incorrect: either the value of NRHS provided at analysis is negative or zero, or the value provided at factorization or solve is different from the value provided at analysis. INFO(2) holds the value of id%NRHS that was provided at analysis.", + -43: "Incompatible values of ICNTL(32) and ICNTL(xx). The index xx is stored in INFO(2).", + -44: "The solve phase (JOB=3) cannot be performed because the factors or part of the factors are not available. INFO(2) contains the value of ICNTL(31).", + -45: "NRHS <= 0. INFO(2) contains the value of NRHS.", + -46: "NZ RHS <= 0. This is currently not allowed in case of reduced right-hand-side (ICNTL(26)=1) and in case entries of A-1 are requested (ICNTL(30)=1). INFO(2) contains the value of NZ RHS.", + -47: "Entries of A-1 were requested during the solve phase (JOB=3, ICNTL(30)=1) but the constraint NRHS=N is not respected. The value of NRHS is provided in INFO(2).", + -48: "A-1 Incompatible values of ICNTL(30) and ICNTL(xx). xx is stored in INFO(2).", + -49: "SIZE SCHUR has an incorrect value: SIZE SCHUR < 0 or SIZE SCHUR >= N, or SIZE SCHUR was modified on the host since the analysis phase. The value of SIZE SCHUR is provided in INFO(2).", + -50: "An error occurred while computing the fill-reducing ordering during the analysis phase. This commonly happens when an (external) ordering tool returns an error code or a wrong result.", + -51: "An external ordering (Metis/ParMetis, SCOTCH/PT-SCOTCH, PORD), with 32-bit default integers, is invoked to processing a graph of size larger than 231 - 1. INFO(2) holds the size required to store the graph as a number of integer values; it is negative and its absolute value should be multiplied by 1 million.", + -52: "When default Fortran integers are 64 bit (e.g. Fortran compiler flag -i8 -fdefault-integer-8 or something equivalent depending on your compiler) then external ordering libraries (Metis/ParMetis, SCOTCH/PT-SCOTCH, PORD) should also have 64-bit default integers. INFO(2) = 1, 2, 3 means that respectively Metis/ParMetis, SCOTCH/PT-SCOTCH or PORD were invoked and were not generated with 64-bit default integers.", + -53: "Internal error that could be due to inconsistent input data between two consecutive calls.", + -54: "The analysis phase (JOB=1) was called with ICNTL(35)=0 but the factorization phase was called with ICNTL(35)=1, 2 or 3. In order to perform the factorization with BLR compression, please perform the analysis phase again using ICNTL(35)=1, 2 or 3 (see the documentation of ICNTL(35)).", + -55: "During a call to MUMPS including the solve phase with distributed right-hand side, LRHS loc was detected to be smaller than Nloc RHS. INFO(2)=LRHS loc.", + -56: "During a call to MUMPS including the solve phase with distributed right-hand side and distributed solution, RHS loc and SOL loc point to the same workarray but LRHS loc < LSOL loc. INFO(2)=LRHS loc.", + -57: "During a call to MUMPS analysis phase with a block format (ICNTL(15) != 0), an error in the interface provided by the user was detected. INFO(2) holds additional information about the issue.", + -70: "During a call to MUMPS with JOB=7, the file specified to save the current instance, as derived from SAVE DIR and/or SAVE PREFIX, already exists. Before saving an instance into this file, it should be first suppressed (see JOB=-3). Otherwise, a different file should be specified by changing the values of SAVE DIR and/or SAVE PREFIX.", + -71: "An error has occured during the creation of one of the files needed to save MUMPS data (JOB=7).", + -72: "Error while saving data (JOB=7); a write operation did not succeed (e.g., disk full, I/O error, . . . ). INFO(2) is the size that should have been written during that operation. If INFO(2) is negative, then its absolute value should be multiplied by 1 million.", + -73: "During a call to MUMPS with JOB=8, one parameter of the current instance is not compatible with the corresponding one in the saved instance.", + -74: "The file resulting from the setting of SAVE DIR and SAVE PREFIX could not be opened for restoring data (JOB=8). INFO(2) is the rank of the process (in the communicator COMM) on which the error was detected.", + -75: "Error while restoring data (JOB=8); a read operation did not succeed (e.g., end of file reached, I/O error, . . . ). INFO(2) is the size still to be read. If INFO(2) is negative, then the size that the package requested is obtained by multiplying the absolute value of INFO(2) by 1 million.", + -76: "Error while deleting the files (JOB=-3); some files to be erased were not found or could not be suppressed. INFO(2) is the rank of the process (in the communicator COMM) on which the error was detected.", + -77: "Neither SAVE DIR nor the environment variable MUMPS SAVE DIR are defined.", + -78: "Problem of workspace allocation during the restore step. The size still to be allocated is available in INFO(2). If INFO(2) is negative, then the size that the package requested is obtained by multiplying the absolute value of INFO(2) by 1 million.", + -79: "MUMPS could not find a Fortran file unit to perform I/O's.", + -89: "Internal error during SCOTCH kway-partitioning in SCOTCHFGRAPHPART. METIS package should be made available to MUMPS.", + -90: "Error in out-of-core management. See the error message returned on output unit ICNTL(1) for more information.", + -800: "Temporary error associated to the current MUMPS release, subject to change or disappearance in the future. If INFO(2)=5, then this error is due to the fact that the elemental matrix format (ICNTL(5)=1) is currently incompatible with a BLR factorization (ICNTL(35) != 0).", + } + + warning_messages = { + 1: "Index (in IRN or JCN) out of range. Action taken by subroutine is to ignore any such entries and continue. INFO(2) is set to the number of faulty entries. Details of the first ten are printed on unit ICNTL(2).", + 2: "During error analysis the max-norm of the computed solution was found to be zero.", + 4: "User data JCN has been modified (internally) by the solver.", + 8: "Warning return from the iterative refinement routine. More than ICNTL(10) iterations are required.", + } + + if ierr == 0: + return "No error occurred." + elif ierr < 0: + return error_messages.get(ierr, "An unknown error occurred.") + else: + # The warning values are actually bitwise flags. + messages = [] + for bit in warning_messages.keys(): + if ierr & bit != 0: + messages.append(warning_messages[bit]) + + if not messages: + messages.append("An unknown warning occurred.") + + return "\n".join(messages) + + +class _Pointer: + """Store a pointer with a destroy function that gets called on garbage collection. + + Multiple Mumps solver instances can point to the same factor in memory. + This class ensures that the factor is garbage collected only when all references are removed. + + You can always clean the Solver, and it will explicitly call the destroy method. """ - def __init__(self, pointerINT, destroyCall): + def __init__(self, pointerINT, dtype): self.INT = pointerINT - self.destroyCall = destroyCall + # The destroy function should only be called from this class upon deallocation. + if dtype == float: + self._destroy_func = _MUMPSINT.destroy_mumps + elif dtype == complex: + self._destroy_func = _MUMPSINT.destroy_mumps_cmplx + else: + raise ValueError(f"Attempted to use an invalid data type ({dtype})") def __del__(self): - self.destroyCall(self.INT) + self._destroy_func(self.INT) + class Mumps(Base): """ documentation:: - http://mumps.enseeiht.fr/doc/userguide_4.10.0.pdf + https://mumps-solver.org/doc/userguide_5.6.1.pdf """ - transpose = False - symmetric = False + _transpose = False + _conjugate = False @property def T(self): - newMS = self.__class__( + """Transpose operator. + + Allows solving A^T * x = b without needing to factor again. + """ + return self._make_copy(transpose=True) + + def conjugate(self): + """Conjugate operator. + + Allows solving \\bar(A) * x = b without needing to factor again. + """ + return self._make_copy(conjugate=True) + + def _make_copy(self, *, transpose=False, conjugate=False): + properties_with_setters = [] + for a in dir(self.__class__): + attr = getattr(self.__class__, a) + if hasattr(attr, "fset") and attr.fset is not None: + properties_with_setters.append(a) + kwargs = {attr: getattr(self, attr) for attr in properties_with_setters} + + copy = self.__class__( self.A, - symmetric=self.symmetric, - fromPointer=self.pointer + from_pointer=self.pointer, + **kwargs, ) - newMS.transpose = not self.transpose - return newMS + copy._transpose = (not self._transpose) if transpose else self._transpose + copy._conjugate = (not self._conjugate) if conjugate else self._conjugate + return copy - def __init__(self, A, symmetric=False, fromPointer=None): + def __init__(self, A, from_pointer=None, **kwargs): self.A = A.tocsc() - self.symmetric = symmetric + # As of 5.6.1, MUMPS has no optimizations for Hermitian matrices + self.set_kwargs(**kwargs, ignore="is_hermitian") - if fromPointer is None: + if from_pointer is None: self.factor() - elif isinstance(fromPointer, _Pointer): - self.pointer = fromPointer + elif isinstance(from_pointer, _Pointer): + self.pointer = from_pointer else: raise Exception('Unknown pointer for construction.') @property - def isfactored(self): + def _is_factored(self): return getattr(self, 'pointer', None) is not None - def _funhandle(self, ftype): - """ - switches the function handle between real and complex. - - ftype in ['F','S','D'] + @property + def _matrix_type(self): + if self.is_symmetric: + if self.is_positive_definite: + return 1 # symmetric, positive definite + return 2 # general symmetric + return 0 # unsymmetric - means factor, solve, destroy - """ + def _funhandle(self, function): + """Switch the function handle between real and complex.""" if self.A.dtype == float: - return {'F': _MUMPSINT.factor_mumps, - 'S': _MUMPSINT.solve_mumps, - 'D': _MUMPSINT.destroy_mumps}[ftype] + return { + 'factor': _MUMPSINT.factor_mumps, + 'solve': _MUMPSINT.solve_mumps + }[function] elif self.A.dtype == complex: - return {'F': _MUMPSINT.factor_mumps_cmplx, - 'S': _MUMPSINT.solve_mumps_cmplx, - 'D': _MUMPSINT.destroy_mumps_cmplx}[ftype] + return { + 'factor': _MUMPSINT.factor_mumps_cmplx, + 'solve': _MUMPSINT.solve_mumps_cmplx + }[function] + + raise ValueError(f"Attempted to use an invalid data type ({self.A.dtype})") def factor(self): - if self.isfactored: + if self._is_factored: return - sym = 1 if self.symmetric else 0 - ierr, p = self._funhandle('F')( - sym, + ierr, p = self._funhandle('factor')( + self._matrix_type, self.A.data, self.A.indices+1, self.A.indptr+1 ) if ierr < 0: - raise Exception("Mumps Exception [{}] - {}".format(ierr, _mumpsErrors[ierr])) + raise Exception(f"Mumps Exception [{ierr}] - {_mumps_message_from_exit_code(ierr)}") elif ierr > 0: - print("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) + warnings.warn(f"Mumps Warning [{ierr}] - {_mumps_message_from_exit_code(ierr)}") - self.pointer = _Pointer(p, self._funhandle('D')) + self.pointer = _Pointer(p, self.A.dtype) def _solveM(self, rhs): self.factor() rhs = rhs.flatten(order='F') n = self.A.shape[0] nrhs = rhs.size // n - T = 1 if self.transpose else 0 - sol = self._funhandle('S')(self.pointer.INT, nrhs, rhs, T) - return sol + T = 1 if self._transpose else 0 + sol = self._funhandle('solve')( + self.pointer.INT, + nrhs, + np.conjugate(rhs) if self._conjugate else rhs, + T + ) + return np.conjugate(sol) if self._conjugate else sol _solve1 = _solveM def clean(self): - self._funhandle('D')(self.pointer.INT) del self.pointer gc.collect() diff --git a/pymatsolver/mumps/mumps_cmplx_p.f90 b/pymatsolver/mumps/mumps_cmplx_p.f90 index 4fb056d..bdd0502 100644 --- a/pymatsolver/mumps/mumps_cmplx_p.f90 +++ b/pymatsolver/mumps/mumps_cmplx_p.f90 @@ -229,17 +229,11 @@ subroutine solve( mumps_par, nrhs, rhs, x, transpose ) mumps_par%icntl(9) = 1 ! for solving Ax = b end if - if (transpose) then - mumps_par%RHS = conjg(mumps_par%RHS) - end if mumps_par%JOB = 3 ! Solve the system. CALL ZMUMPS(mumps_par) ! At this point mumps_par%RHS (rhs) contains the solution. - if (transpose) then - mumps_par%RHS = conjg(mumps_par%RHS) - end if return end subroutine solve diff --git a/pymatsolver/mumps/mumps_p.f90 b/pymatsolver/mumps/mumps_p.f90 old mode 100755 new mode 100644 diff --git a/pymatsolver/solvers.py b/pymatsolver/solvers.py index 890e508..1111300 100644 --- a/pymatsolver/solvers.py +++ b/pymatsolver/solvers.py @@ -65,6 +65,18 @@ def T(self): newS = self._transposeClass(self.A.T) return newS + @property + def conjugate(self): + "The element-wise complex conjugate for this class" + if self._transposeClass is None: + raise Exception( + 'The complex conjugate for the {} class is not possible.'.format( + self.__name__ + ) + ) + newS = self._transposeClass(self.A.conjugate()) + return newS + def _compute_accuracy(self, rhs, x): nrm = np.linalg.norm(np.ravel(self.A*x - rhs), np.inf) nrm_rhs = np.linalg.norm(np.ravel(rhs), np.inf) @@ -87,13 +99,16 @@ def _solve(self, rhs): else: x = self._solveM(rhs) + if nrhs == 1: + x = x.flatten() + elif nrhs > 1: + x = x.reshape((n, nrhs), order='F') + if self.check_accuracy: self._compute_accuracy(rhs, x) - if nrhs == 1: - return x.flatten() - elif nrhs > 1: - return x.reshape((n, nrhs), order='F') + return x + def clean(self): pass diff --git a/setup.py b/setup.py index f7681ef..977e738 100644 --- a/setup.py +++ b/setup.py @@ -5,8 +5,7 @@ """ -from distutils.core import setup -from setuptools import find_packages +from setuptools import find_packages, setup CLASSIFIERS = [ 'Development Status :: 4 - Beta', @@ -42,7 +41,7 @@ license="MIT", keywords="matrix solver", url="http://simpeg.xyz/", - download_url="http://github.com/rowanc1/pymatsolver", + download_url="http://github.com/simpeg/pymatsolver", classifiers=CLASSIFIERS, platforms=["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"], use_2to3=False diff --git a/tests/test_Basic.py b/tests/test_Basic.py index adcd4c5..bd01a6d 100644 --- a/tests/test_Basic.py +++ b/tests/test_Basic.py @@ -1,4 +1,4 @@ -import unittest +import pytest import numpy as np import scipy.sparse as sp from pymatsolver import Diagonal @@ -6,10 +6,9 @@ TOL = 1e-12 -class TestBasic(unittest.TestCase): +class TestBasic: def test_DiagonalSolver(self): - A = sp.identity(5)*2.0 rhs = np.c_[np.arange(1, 6), np.arange(2, 11, 2)] X = Diagonal(A) * rhs @@ -17,14 +16,10 @@ def test_DiagonalSolver(self): sol = rhs/2.0 - with self.assertRaises(TypeError): + with pytest.raises(TypeError): Diagonal(A, check_accuracy=np.array([1, 2, 3])) - with self.assertRaises(TypeError): + with pytest.raises(TypeError): Diagonal(A, accuracy_tol=0) - self.assertLess(np.linalg.norm(sol-X, np.inf), TOL) - self.assertLess(np.linalg.norm(sol[:, 0]-x, np.inf), TOL) - - -if __name__ == '__main__': - unittest.main() + assert np.linalg.norm(sol-X, np.inf) < TOL + assert np.linalg.norm(sol[:, 0]-x, np.inf) < TOL diff --git a/tests/test_BicgJacobi.py b/tests/test_BicgJacobi.py index 24e41a5..cd01b3f 100644 --- a/tests/test_BicgJacobi.py +++ b/tests/test_BicgJacobi.py @@ -1,4 +1,3 @@ -import unittest from pymatsolver import BicgJacobi import numpy as np import scipy.sparse as sp @@ -6,9 +5,10 @@ TOL = 1e-6 -class TestBicgJacobi(unittest.TestCase): +class TestBicgJacobi: - def setUp(self): + @classmethod + def setup_class(cls): nSize = 100 A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) @@ -19,9 +19,9 @@ def setUp(self): sol = np.random.rand(nSize, 4) rhs = A.dot(sol) - self.A = A - self.rhs = rhs - self.sol = sol + cls.A = A + cls.rhs = rhs + cls.sol = sol def test(self): rhs = self.rhs @@ -30,7 +30,7 @@ def test(self): for i in range(3): err = np.linalg.norm( self.A*solb[:, i] - rhs[:, i]) / np.linalg.norm(rhs[:, i]) - self.assertLess(err, TOL) + assert err < TOL Ainv.clean() def test_T(self): @@ -42,13 +42,14 @@ def test_T(self): for i in range(3): err = np.linalg.norm( self.A.T*solb[:, i] - rhs[:, i]) / np.linalg.norm(rhs[:, i]) - self.assertLess(err, TOL) + assert err < TOL Ainv.clean() -class TestPardisoComplex(unittest.TestCase): +class TestBicgJacobiComplex: - def setUp(self): + @classmethod + def setup_class(cls): nSize = 100 A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) A.data = A.data + 1j*np.random.rand(A.nnz) @@ -59,9 +60,9 @@ def setUp(self): sol = np.random.rand(nSize, 5) + 1j*np.random.rand(nSize, 5) rhs = A.dot(sol) - self.A = A - self.rhs = rhs - self.sol = sol + cls.A = A + cls.rhs = rhs + cls.sol = sol def test(self): rhs = self.rhs @@ -70,7 +71,7 @@ def test(self): for i in range(3): err = np.linalg.norm( self.A*solb[:, i] - rhs[:, i]) / np.linalg.norm(rhs[:, i]) - self.assertLess(err, TOL) + assert err < TOL Ainv.clean() def test_T(self): @@ -82,9 +83,5 @@ def test_T(self): for i in range(3): err = np.linalg.norm( self.A.T*solb[:, i] - rhs[:, i]) / np.linalg.norm(rhs[:, i]) - self.assertLess(err, TOL) + assert err < TOL Ainv.clean() - - -if __name__ == '__main__': - unittest.main() diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index 1bf4afc..b3c505a 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -1,16 +1,192 @@ -import unittest +import os +import warnings + import numpy as np +import pytest import scipy.sparse as sp -import pymatsolver + +try: + from pymatsolver import Mumps + should_run = True +except ImportError: + should_run = False TOL = 1e-11 -""" -class TestMumps(unittest.TestCase): - if pymatsolver.AvailableSolvers['Mumps']: +if should_run: + + class TestMumps: + + @classmethod + def setup_class(cls): + + nSize = 100 + A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) + A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) + A = A.T*A + A = A.tocsr() + np.random.seed(1) + sol = np.random.rand(nSize, 5) + rhs = A.dot(sol) + + cls.A = A + cls.rhs = rhs + cls.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True) + for i in range(3): + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL + + def test_T(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True) + AinvT = Ainv.T + x = AinvT * rhs + + for i in range(3): + assert np.linalg.norm(x[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(x - sol, np.inf) < TOL + + class TestMumpsNotSymmetric: + + @classmethod + def setup_class(cls): + + nSize = 100 + A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) + A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) + A = A.tocsr() + np.random.seed(1) + sol = np.random.rand(nSize, 5) + rhs = A.dot(sol) + + cls.A = A + cls.rhs = rhs + cls.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True, check_accuracy=True) + with pytest.raises(Exception): + Ainv * rhs + Ainv.clean() + + Ainv = Mumps(self.A, check_accuracy=True) + for i in range(3): + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL + Ainv.clean() + + + class TestMumpsFDEM: + + @classmethod + def setup_class(cls): + + base_path = os.path.join(os.path.split(os.path.abspath(__file__))[0], 'fdem') + + data = np.load(os.path.join(base_path, 'A_data.npy')) + indices = np.load(os.path.join(base_path, 'A_indices.npy')) + indptr = np.load(os.path.join(base_path, 'A_indptr.npy')) + + cls.A = sp.csr_matrix((data, indices, indptr), shape=(13872, 13872)) + cls.rhs = np.load(os.path.join(base_path, 'RHS.npy')) + + def test(self): + rhs = self.rhs + Ainv = Mumps(self.A, check_accuracy=True) + sol = Ainv * rhs + sol = Ainv * rhs.real + + + class TestMumpsComplex: + + @classmethod + def setup_class(cls): + nSize = 100 + A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) + A.data = A.data + 1j*np.random.rand(A.nnz) + A = A.T.dot(A) + sp.spdiags(np.ones(nSize), 0, nSize, nSize) + A = A.tocsr() + + np.random.seed(1) + sol = np.random.rand(nSize, 5) + 1j*np.random.rand(nSize, 5) + rhs = A.dot(sol) + + cls.A = A + cls.rhs = rhs + cls.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True) + for i in range(3): + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL + Ainv.clean() + + def test_T(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True) + AinvT = Ainv.T + x = AinvT * rhs + for i in range(3): + assert np.linalg.norm(x[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(x - sol, np.inf) < TOL + + def test_T_congruent(self): + rhs = self.rhs + sol = self.sol + + AinvT1 = Mumps(self.A, is_symmetric=True).T + AinvT2 = Mumps(self.A.T, is_symmetric=True) + x1 = AinvT2 * rhs + x2 = AinvT1 * rhs + + for i in range(3): + assert np.linalg.norm(x1[:, i] - x2[:, i]) < TOL + assert np.linalg.norm(x1 - x2, np.inf) < TOL + + def test_conjugate_congruent(self): + rhs = self.rhs + sol = self.sol - def setUp(self): + Ainv_conj1 = Mumps(self.A).conjugate() + Ainv_conj2 = Mumps(self.A.conjugate()) + x1 = Ainv_conj1 * rhs + x2 = Ainv_conj2 * rhs + + for i in range(3): + assert np.linalg.norm(x1[:, i] - x2[:, i]) < TOL + assert np.linalg.norm(x1 - x2, np.inf) < TOL + + def test_complex_conjugate_congruent(self): + rhs = self.rhs + sol = self.sol + + Ainv_conj1 = Mumps(self.A).T.conjugate() + Ainv_conj2 = Mumps(self.A.T.conjugate()) + x1 = Ainv_conj1 * rhs + x2 = Ainv_conj2 * rhs + + for i in range(3): + assert np.linalg.norm(x1[:, i] - x2[:, i]) < TOL + assert np.linalg.norm(x1 - x2, np.inf) < TOL + + + class TestMumps1to5: + + @classmethod + def setup_class(cls): n = 5 irn = np.r_[0, 1, 3, 4, 1, 0, 4, 2, 1, 2, 0, 2] jcn = np.r_[1, 2, 2, 4, 0, 0, 1, 3, 4, 1, 2, 2] @@ -23,64 +199,58 @@ def setUp(self): sol = np.r_[1., 2., 3., 4., 5.] sol = np.c_[sol, 10*sol, 100*sol] A = sp.coo_matrix((a, (irn, jcn)), shape=(n, n)) - self.A = A - self.rhs = rhs - self.sol = sol + cls.A = A + cls.rhs = rhs + cls.sol = sol - def test_1to5(self): + def test(self): rhs = self.rhs sol = self.sol - Ainv = pymatsolver.Mumps(self.A) + Ainv = Mumps(self.A) for i in range(3): - self.assertLess( - np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL - ) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL - def test_1to5_cmplx(self): + def test_cmplx(self): rhs = self.rhs.astype(complex) sol = self.sol.astype(complex) self.A = self.A.astype(complex) - Ainv = pymatsolver.Mumps(self.A) + Ainv = Mumps(self.A) for i in range(3): - self.assertLess( - np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL - ) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL - def test_1to5_T(self): + def test_T(self): rhs = self.rhs sol = self.sol - Ainv = pymatsolver.Mumps(self.A) + Ainv = Mumps(self.A) AinvT = Ainv.T for i in range(3): - self.assertLess( - np.linalg.norm(AinvT.T * rhs[:, i] - sol[:, i]), TOL - ) - self.assertLess(np.linalg.norm(AinvT.T * rhs - sol, np.inf), TOL) + assert np.linalg.norm(AinvT.T * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(AinvT.T * rhs - sol, np.inf) < TOL - # def test_singular(self): - # A = sp.identity(5).tocsr() - # A[-1, -1] = 0 - # self.assertRaises(Exception, pymatsolver.Mumps, A) + def test_singular(self): + A = sp.identity(5).tocsr() + A[-1, -1] = 0 + with pytest.raises(Exception): + Mumps(A) - def test_multiFactorsInMem(self): + def test_multi_factors_in_mem(self): n = 100 A = sp.rand(n, n, 0.7)+sp.identity(n) x = np.ones((n, 10)) rhs = A * x - solvers = [pymatsolver.Mumps(A) for _ in range(20)] + solvers = [Mumps(A) for _ in range(20)] for Ainv in solvers: - self.assertLess( - np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs), TOL) + assert np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs) < TOL Ainv.clean() for Ainv in solvers: - self.assertLess( - np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs), TOL - ) + assert np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs) < TOL + +else: -if __name__ == '__main__': - unittest.main() -""" + def test_mumps_not_available(): + """Run only when Mumps is not available.""" + warnings.warn("NOTE: Mumps is not available so we're skipping its tests.") diff --git a/tests/test_Pardiso.py b/tests/test_Pardiso.py index da10975..53d5e3f 100644 --- a/tests/test_Pardiso.py +++ b/tests/test_Pardiso.py @@ -1,180 +1,195 @@ -import unittest -from pymatsolver import Pardiso -from pydiso.mkl_solver import ( - get_mkl_pardiso_max_threads, - PardisoTypeConversionWarning -) +import os +import warnings + import numpy as np +import pytest import scipy.sparse as sp -import os + +try: + from pymatsolver import Pardiso + from pydiso.mkl_solver import ( + get_mkl_pardiso_max_threads, + PardisoTypeConversionWarning + ) + should_run = True +except ImportError: + should_run = False + TOL = 1e-10 +if should_run: + + class TestPardiso: -class TestPardiso(unittest.TestCase): - - def setUp(self): - - nSize = 100 - A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) - A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) - A = A.T*A - A = A.tocsr() - np.random.seed(1) - sol = np.random.rand(nSize, 5) - rhs = A.dot(sol) - - self.A = A - self.rhs = rhs - self.sol = sol - - def test(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - - def test_refactor(self): - rhs = self.rhs - sol = self.sol - A = self.A - Ainv = Pardiso(A, is_symmetric=True) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - - # scale rows and collumns - D = sp.diags(np.random.rand(A.shape[0]) + 1.0) - A2 = D.T @ A @ D - - rhs2 = A2 @ sol - Ainv.factor(A2) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs2[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs2 - sol, np.inf), TOL) - - def test_T(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True) - - with self.assertWarns(PardisoTypeConversionWarning): - AinvT = Ainv.T - x = AinvT * rhs + @classmethod + def setup_class(cls): + nSize = 100 + A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) + A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) + A = A.T*A + A = A.tocsr() + np.random.seed(1) + sol = np.random.rand(nSize, 5) + rhs = A.dot(sol) + + cls.A = A + cls.rhs = rhs + cls.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True) + for i in range(3): + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL + + def test_refactor(self): + rhs = self.rhs + sol = self.sol + A = self.A + Ainv = Pardiso(A, is_symmetric=True) for i in range(3): - self.assertLess(np.linalg.norm(x[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(x - sol, np.inf), TOL) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL - def test_n_threads(self): - max_threads = get_mkl_pardiso_max_threads() - print(f'testing setting n_threads to 1 and {max_threads}') - Ainv = Pardiso(self.A, is_symmetric=True, n_threads=1) - self.assertEqual(Ainv.n_threads, 1) - - Ainv2 = Pardiso(self.A, is_symmetric=True, n_threads=max_threads) - self.assertEqual(Ainv2.n_threads, max_threads) - self.assertEqual(Ainv2.n_threads, Ainv.n_threads) - - Ainv.n_threads = 1 - self.assertEqual(Ainv.n_threads, 1) - self.assertEqual(Ainv2.n_threads, Ainv.n_threads) - - with self.assertRaises(TypeError): - Ainv.n_threads = "2" - - -class TestPardisoNotSymmetric(unittest.TestCase): - - def setUp(self): - - nSize = 100 - A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) - A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) - A = A.tocsr() - np.random.seed(1) - sol = np.random.rand(nSize, 5) - rhs = A.dot(sol) - - self.A = A - self.rhs = rhs - self.sol = sol - - def test(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True, check_accuracy=True) - self.assertRaises(Exception, lambda: Ainv * rhs) - Ainv.clean() - - Ainv = Pardiso(self.A) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - Ainv.clean() - - -class TestPardisoFDEM(unittest.TestCase): - - def setUp(self): - - base_path = os.path.join(os.path.split(os.path.abspath(__file__))[0], 'fdem') - - data = np.load(os.path.join(base_path, 'A_data.npy')) - indices = np.load(os.path.join(base_path, 'A_indices.npy')) - indptr = np.load(os.path.join(base_path, 'A_indptr.npy')) - - self.A = sp.csr_matrix((data, indices, indptr), shape=(13872, 13872)) - self.rhs = np.load(os.path.join(base_path, 'RHS.npy')) - - def test(self): - rhs = self.rhs - Ainv = Pardiso(self.A, check_accuracy=True) - sol = Ainv * rhs - with self.assertWarns(PardisoTypeConversionWarning): - sol = Ainv * rhs.real - - -class TestPardisoComplex(unittest.TestCase): - - def setUp(self): - nSize = 100 - A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) - A.data = A.data + 1j*np.random.rand(A.nnz) - A = A.T.dot(A) + sp.spdiags(np.ones(nSize), 0, nSize, nSize) - A = A.tocsr() - - np.random.seed(1) - sol = np.random.rand(nSize, 5) + 1j*np.random.rand(nSize, 5) - rhs = A.dot(sol) - - self.A = A - self.rhs = rhs - self.sol = sol - - def test(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - Ainv.clean() - - def test_T(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True) - with self.assertWarns(PardisoTypeConversionWarning): - AinvT = Ainv.T - x = AinvT * rhs + # scale rows and columns + D = sp.diags(np.random.rand(A.shape[0]) + 1.0) + A2 = D.T @ A @ D + + rhs2 = A2 @ sol + Ainv.factor(A2) + for i in range(3): + assert np.linalg.norm(Ainv * rhs2[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs2 - sol, np.inf) < TOL + + def test_T(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True) + + with pytest.warns(PardisoTypeConversionWarning): + AinvT = Ainv.T + x = AinvT * rhs + + for i in range(3): + assert np.linalg.norm(x[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(x - sol, np.inf) < TOL + + def test_n_threads(self): + max_threads = get_mkl_pardiso_max_threads() + print(f'testing setting n_threads to 1 and {max_threads}') + Ainv = Pardiso(self.A, is_symmetric=True, n_threads=1) + assert Ainv.n_threads == 1 + + Ainv2 = Pardiso(self.A, is_symmetric=True, n_threads=max_threads) + assert Ainv2.n_threads == max_threads + assert Ainv2.n_threads == Ainv.n_threads + + Ainv.n_threads = 1 + assert Ainv.n_threads == 1 + assert Ainv2.n_threads == Ainv.n_threads + + with pytest.raises(TypeError): + Ainv.n_threads = "2" + + + class TestPardisoNotSymmetric: + + @classmethod + def setup_class(cls): + + nSize = 100 + A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) + A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) + A = A.tocsr() + np.random.seed(1) + sol = np.random.rand(nSize, 5) + rhs = A.dot(sol) + + cls.A = A + cls.rhs = rhs + cls.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True, check_accuracy=True) + with pytest.raises(Exception): + Ainv * rhs + Ainv.clean() + + Ainv = Pardiso(self.A) for i in range(3): - self.assertLess( - np.linalg.norm(x[:, i] - sol[:, i]), TOL - ) - self.assertLess(np.linalg.norm(x - sol, np.inf), TOL) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL + Ainv.clean() + -if __name__ == '__main__': - unittest.main() + class TestPardisoFDEM: + + @classmethod + def setup_class(cls): + + base_path = os.path.join(os.path.split(os.path.abspath(__file__))[0], 'fdem') + + data = np.load(os.path.join(base_path, 'A_data.npy')) + indices = np.load(os.path.join(base_path, 'A_indices.npy')) + indptr = np.load(os.path.join(base_path, 'A_indptr.npy')) + + cls.A = sp.csr_matrix((data, indices, indptr), shape=(13872, 13872)) + cls.rhs = np.load(os.path.join(base_path, 'RHS.npy')) + + def test(self): + rhs = self.rhs + Ainv = Pardiso(self.A, check_accuracy=True) + sol = Ainv * rhs + with pytest.warns(PardisoTypeConversionWarning): + sol = Ainv * rhs.real + + + class TestPardisoComplex: + + @classmethod + def setup_class(cls): + nSize = 100 + A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) + A.data = A.data + 1j*np.random.rand(A.nnz) + A = A.T.dot(A) + sp.spdiags(np.ones(nSize), 0, nSize, nSize) + A = A.tocsr() + + np.random.seed(1) + sol = np.random.rand(nSize, 5) + 1j*np.random.rand(nSize, 5) + rhs = A.dot(sol) + + cls.A = A + cls.rhs = rhs + cls.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True) + for i in range(3): + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL + Ainv.clean() + + def test_T(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True) + with pytest.warns(PardisoTypeConversionWarning): + AinvT = Ainv.T + x = AinvT * rhs + for i in range(3): + assert np.linalg.norm(x[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(x - sol, np.inf) < TOL + +else: + + def test_pardiso_not_available(): + """Run only when Pardiso is not available.""" + warnings.warn("NOTE: Pardiso is not available so we're skipping its tests.") diff --git a/tests/test_Scipy.py b/tests/test_Scipy.py index eb77135..8969ba7 100644 --- a/tests/test_Scipy.py +++ b/tests/test_Scipy.py @@ -1,4 +1,3 @@ -import unittest from pymatsolver import Solver, Diagonal, SolverCG, SolverLU import scipy.sparse as sp import numpy as np @@ -50,38 +49,28 @@ def dotest(MYSOLVER, multi=False, A=None, **solverOpts): return np.linalg.norm(e-x, np.inf) -class TestSolver(unittest.TestCase): +class TestSolver: def test_direct_spsolve_1(self): - self.assertLess(dotest(Solver, False), TOLD) + assert dotest(Solver, False) < TOLD def test_direct_spsolve_M(self): - self.assertLess(dotest(Solver, True), TOLD) + assert dotest(Solver, True) < TOLD def test_direct_splu_1(self): - self.assertLess(dotest(SolverLU, False), TOLD) + assert dotest(SolverLU, False) < TOLD def test_direct_splu_M(self): - self.assertLess(dotest(SolverLU, True), TOLD) + assert dotest(SolverLU, True) < TOLD def test_iterative_diag_1(self): - self.assertLess(dotest( - Diagonal, False, - A=sp.diags(np.random.rand(10)+1.0) - ), TOLI) + assert dotest(Diagonal, False, A=sp.diags(np.random.rand(10)+1.0)) < TOLI def test_iterative_diag_M(self): - self.assertLess(dotest( - Diagonal, True, - A=sp.diags(np.random.rand(10)+1.0) - ), TOLI) + assert dotest(Diagonal, True, A=sp.diags(np.random.rand(10)+1.0)) < TOLI def test_iterative_cg_1(self): - self.assertLess(dotest(SolverCG, False), TOLI) + assert dotest(SolverCG, False) < TOLI def test_iterative_cg_M(self): - self.assertLess(dotest(SolverCG, True), TOLI) - - -if __name__ == '__main__': - unittest.main() + assert dotest(SolverCG, True) < TOLI diff --git a/tests/test_Triangle.py b/tests/test_Triangle.py index 241d1e0..98ee279 100644 --- a/tests/test_Triangle.py +++ b/tests/test_Triangle.py @@ -1,4 +1,3 @@ -import unittest import numpy as np import scipy.sparse as sp import pymatsolver @@ -6,29 +5,27 @@ TOL = 1e-12 -class TestTriangle(unittest.TestCase): +class TestTriangle: - def setUp(self): + @classmethod + def setup_class(cls): n = 50 nrhs = 20 - self.A = sp.rand(n, n, 0.4) + sp.identity(n) - self.sol = np.ones((n, nrhs)) - self.rhsU = sp.triu(self.A) * self.sol - self.rhsL = sp.tril(self.A) * self.sol + cls.A = sp.rand(n, n, 0.4) + sp.identity(n) + cls.sol = np.ones((n, nrhs)) + cls.rhsU = sp.triu(cls.A) * cls.sol + cls.rhsL = sp.tril(cls.A) * cls.sol def test_directLower(self): ALinv = pymatsolver.Forward(sp.tril(self.A)) X = ALinv * self.rhsL x = ALinv * self.rhsL[:, 0] - self.assertLess(np.linalg.norm(self.sol-X, np.inf), TOL) - self.assertLess(np.linalg.norm(self.sol[:, 0]-x, np.inf), TOL) + assert np.linalg.norm(self.sol-X, np.inf) < TOL + assert np.linalg.norm(self.sol[:, 0]-x, np.inf) < TOL def test_directLower_1(self): AUinv = pymatsolver.Backward(sp.triu(self.A)) X = AUinv * self.rhsU x = AUinv * self.rhsU[:, 0] - self.assertLess(np.linalg.norm(self.sol-X, np.inf), TOL) - self.assertLess(np.linalg.norm(self.sol[:, 0]-x, np.inf), TOL) - -if __name__ == '__main__': - unittest.main() + assert np.linalg.norm(self.sol-X, np.inf) < TOL + assert np.linalg.norm(self.sol[:, 0]-x, np.inf) < TOL