From a6b0d64a99b62a461cc94218ececfe6b6f2f1205 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 7 Jan 2017 20:36:13 -0800 Subject: [PATCH 1/6] Bring back mumps --- Makefile | 29 +++ pymatsolver/__init__.py | 29 ++- pymatsolver/mumps/Makefile | 45 +++++ pymatsolver/mumps/__init__.py | 171 ++++++++++++++++ pymatsolver/mumps/mumps_cmplx_p.f90 | 272 ++++++++++++++++++++++++++ pymatsolver/mumps/mumps_interface.f90 | 192 ++++++++++++++++++ pymatsolver/mumps/mumps_p.f90 | 266 +++++++++++++++++++++++++ tests/test_Mumps.py | 85 ++++++++ 8 files changed, 1087 insertions(+), 2 deletions(-) create mode 100644 Makefile create mode 100644 pymatsolver/mumps/Makefile create mode 100644 pymatsolver/mumps/__init__.py create mode 100644 pymatsolver/mumps/mumps_cmplx_p.f90 create mode 100644 pymatsolver/mumps/mumps_interface.f90 create mode 100755 pymatsolver/mumps/mumps_p.f90 create mode 100644 tests/test_Mumps.py diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..08a96fe --- /dev/null +++ b/Makefile @@ -0,0 +1,29 @@ +.PHONY: build coverage lint graphs tests docs mumps mumps_mac mumps_install_mac + +build: + python setup.py build_ext --inplace + +coverage: + nosetests --logging-level=INFO --with-coverage --cover-package=pymatsolver --cover-html + open cover/index.html + +lint: + pylint --output-format=html pymatsolver > pylint.html + +graphs: + pyreverse -my -A -o pdf -p pymatsolver pymatsolver/**.py pymatsolver/**/**.py + +tests: + nosetests --logging-level=INFO + +docs: + cd docs;make html + +mumps: + cd pymatsolver/mumps;make build + +mumps_mac: + cd pymatsolver/mumps;make build_mac + +mumps_install_mac: + brew install mumps --with-scotch5 --without-mpi diff --git a/pymatsolver/__init__.py b/pymatsolver/__init__.py index eb9cdeb..db3bc41 100644 --- a/pymatsolver/__init__.py +++ b/pymatsolver/__init__.py @@ -13,20 +13,45 @@ from pymatsolver.iterative import BicgJacobi +SolverHelp = {} AvailableSolvers = { "Diagonal": True, "Solver": True, "SolverLU": True, "SolverCG": True, "Triangle": True, - "Pardiso": False + "Pardiso": False, + "Mumps": False } try: from pymatsolver.direct import Pardiso AvailableSolvers['Pardiso'] = True except ImportError: - pass + SolverHelp['Pardiso'] = """Pardiso is not working + +Ensure that you have pyMKL installed, which may also require Python +to be installed through conda. +""" + +try: + from pymatsolver.mumps import Mumps + AvailableSolvers['Mumps'] = True +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.0.3' __author__ = 'Rowan Cockett' diff --git a/pymatsolver/mumps/Makefile b/pymatsolver/mumps/Makefile new file mode 100644 index 0000000..f6a809b --- /dev/null +++ b/pymatsolver/mumps/Makefile @@ -0,0 +1,45 @@ +INCMUMPS = -I/usr/include -L/usr/lib \ + -L/usr/local/opt/libevent/lib +INCMAC = -L/usr/local/opt/scotch5/lib \ + -L/usr/local/Cellar/open-mpi/2.0.1/lib \ + -L/usr/local/opt/libevent/lib \ + -L/usr/local/opt/mumps/libexec/lib \ + -I/usr/local/opt/mumps/libexec/include \ + -I/usr/local/opt/scotch5/include \ + -I/usr/local/Cellar/open-mpi/2.0.1/include + +# brew install mumps --with-scotch5 --without-mpi +# Helpful: mpicc --showme + +LIB = -lmpi -lblas -lpord -lcmumps +LIBSEQ = -lmpiseq_seq -lzmumps_seq -ldmumps_seq -lsmumps_seq -lmumps_common_seq +LIBPAR = -lmpiseq -lzmumps -ldmumps -lsmumps -lmumps_common +LIBSCOTCH = -lesmumps -lptscotch -lptscotcherrexit -lscotch -lscotcherrexit -lptesmumps -lptscotcherr -lptscotchparmetis -lscotcherr -lscotchmetis + +all: clean build build_mac + +build_mac: + gfortran -c $(INCMAC) mumps_p.f90 -fPIC + gfortran -c $(INCMAC) mumps_cmplx_p.f90 -fPIC + f2py -c mumps_interface.f90 -m MumpsInterface \ + $(INCMAC) \ + --f90flags='-fcray-pointer' \ + $(LIB) $(LIBPAR) $(LIBSCOTCH)\ + mumps_p.o mumps_cmplx_p.o + rm -f *.o *.mod + +build: + gfortran -c $(INCMUMPS) mumps_p.f90 -fPIC + gfortran -c $(INCMUMPS) mumps_cmplx_p.f90 -fPIC + f2py -c mumps_interface.f90 -m MumpsInterface \ + $(INCMUMPS) \ + --f90flags='-fcray-pointer' \ + $(LIB) $(LIBSEQ)\ + mumps_p.o mumps_cmplx_p.o + rm -f *.o *.mod + +clean: + rm -f *.o *.mod *.so + +test: + cd ..;nosetests diff --git a/pymatsolver/mumps/__init__.py b/pymatsolver/mumps/__init__.py new file mode 100644 index 0000000..a554a0e --- /dev/null +++ b/pymatsolver/mumps/__init__.py @@ -0,0 +1,171 @@ +from __future__ import unicode_literals +from __future__ import print_function +from __future__ import division +# from __future__ import absolute_import + +import MumpsInterface +import gc +from pymatsolver.solvers import Base + +_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.", +} + + +# if _mumpsExists: + +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 __init__(self, pointerINT, destroyCall): + self.INT = pointerINT + self.destroyCall = destroyCall + + def __del__(self): + self.destroyCall(self.INT) + + +class Mumps(Base): + """ + + documentation:: + + http://mumps.enseeiht.fr/doc/userguide_4.10.0.pdf + + """ + + transpose = False + symmetric = False + + @property + def T(self): + newMS = self.__class__( + self.A, + symmetric=self.symmetric, + fromPointer=self.pointer + ) + newMS.transpose = not self.transpose + return newMS + + def __init__(self, A, symmetric=False, fromPointer=None): + self.A = A.tocsc() + self.symmetric = symmetric + + if fromPointer is None: + self.factor() + elif isinstance(fromPointer, _Pointer): + self.pointer = fromPointer + else: + raise Exception('Unknown pointer for construction.') + + @property + def isfactored(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'] + + means factor, solve, destroy + """ + if self.A.dtype == float: + return {'F': MumpsInterface.factor_mumps, + 'S': MumpsInterface.solve_mumps, + 'D': MumpsInterface.destroy_mumps}[ftype] + elif self.A.dtype == complex: + return {'F': MumpsInterface.factor_mumps_cmplx, + 'S': MumpsInterface.solve_mumps_cmplx, + 'D': MumpsInterface.destroy_mumps_cmplx}[ftype] + + def factor(self): + if self.isfactored: + return + + sym = 1 if self.symmetric else 0 + ierr, p = self._funhandle('F')( + sym, + self.A.data, + self.A.indices+1, + self.A.indptr+1 + ) + if ierr < 0: + raise Exception("Mumps Exception [{}] - {}".format(ierr, _mumpsErrors[ierr])) + elif ierr > 0: + print("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) + + self.pointer = _Pointer(p, self._funhandle('D')) + + 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 + + _solve1 = _solveM + + def clean(self): + self._funhandle('D')(self.pointer.INT) + del self.pointer + gc.collect() + + +# del _mumpsExists diff --git a/pymatsolver/mumps/mumps_cmplx_p.f90 b/pymatsolver/mumps/mumps_cmplx_p.f90 new file mode 100644 index 0000000..4fb056d --- /dev/null +++ b/pymatsolver/mumps/mumps_cmplx_p.f90 @@ -0,0 +1,272 @@ + + + +module mumps_cmplx_mod + + save + ! integer,parameter:: ui_out = 10 ! unit for mumps log file + !integer,parameter:: st_unit=11 ! unit for output statistics file + + + INCLUDE 'zmumps_struc.h' + + !logical,private:: mumps_initialized = .false. + + +contains + +!--------------------------------------------------------------- + +subroutine init( mumps_par, sym ) + +implicit none + +TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par +integer(kind=8),intent(in):: sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. + + !mumps_par%COMM = MPI_COMM_WORLD + mumps_par%SYM = sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. + mumps_par%PAR = 1 + + mumps_par%JOB = -1 ! initialize + CALL ZMUMPS(mumps_par) + + + !if (mumps_par%MYID == 0) then + ! open(unit=ui_out, file='mumps.log', action='write') + !end if + + !mumps_initialized = .true. + +return +end subroutine init + +!--------------------------------------------------------------- + +subroutine convert_to_mumps_format( mumps_par, n, A,jA,iA, ierr ) +! A is actually a transpose. +implicit none + +TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par +integer(kind=8),intent(in):: n ! # of rows in A +complex(kind=8),intent(in):: A(*) +integer(kind=8),intent(in):: jA(*), iA(n+1) +integer(kind=8),intent(out):: ierr + +integer nonz, i,j, j1,j2, ind, jcol, istat + + +nonz = iA(n+1) - 1 +!n = size(iA) - 1 + +mumps_par%N = n + +if (mumps_par%SYM == 0) then ! not symmetric matrix + mumps_par%NZ = nonz + allocate( mumps_par%IRN(nonz), mumps_par%JCN(nonz), mumps_par%A(nonz), stat=istat) + if ( istat /= 0 ) then + ierr = -13 ! allocation error + return + end if + + + do i = 1, n + mumps_par%JCN( iA(i) : iA(i+1)-1 ) = i + end do ! i + + mumps_par%IRN = jA(1:nonz) + mumps_par%A = A(1:nonz) + +else ! symmetric matrix + ! Keep only the lower half of the matrix (row >= column). + + nonz = nonz/2 + n ! should be +n/2, but I'm using n just in case. + allocate( mumps_par%IRN(nonz), mumps_par%JCN(nonz), mumps_par%A(nonz), stat=istat) + if ( istat /= 0 ) then + ierr = -13 ! allocation error + return + end if + + ind = 0 + do i = 1, n + + j1 = iA(i) + j2 = iA(i+1) - 1 + do j = j1, j2 + jcol = jA(j) + + if (i >= jcol) then + ind = ind + 1 + mumps_par%A(ind) = A(j) + mumps_par%JCN(ind) = jcol + mumps_par%IRN(ind) = i + end if + + end do ! j + end do ! i + + mumps_par%NZ = ind + !if (nonz < ind) call errorstop('nonz < ind') ! debug +end if + +ierr = 0 +return +end subroutine convert_to_mumps_format + +!--------------------------------------------------------------- + +subroutine factor_matrix( mumps_par, ierr ) + +implicit none + +TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par +integer(kind=8),intent(out):: ierr + +!mumps_par%icntl(2) = 6 ! output stream for diagnostics +mumps_par%icntl(4) = 0 ! 1 ! amount of output + + +mumps_par%icntl(2) = 0 ! ui_out ! output stream for diagnostic printing +mumps_par%icntl(3) = 0 ! ui_out ! output stream for global information + +!mumps_par%icntl(14) = 40 ! % increase in working space + +!mumps_par%icntl(11) = 0 ! return statistics +! mumps_par%icntl(7) = 5 ! ordering type +!mumps_par%cntl(1) = 0.d0 ! rel. threshold for numerical pivoting + + +mumps_factorization: do + + mumps_par%JOB = 1 ! analysis + CALL ZMUMPS(mumps_par) + + mumps_par%JOB = 2 ! factorization + CALL ZMUMPS(mumps_par) + + ierr = mumps_par%INFOG(1) + + + if ( mumps_par%INFOG(1) == -9 .or. mumps_par%INFOG(1) == -8 ) then + ! Main internal real/complex workarray S too small. + mumps_par%icntl(14) = mumps_par%icntl(14) + 10 + if ( mumps_par%MYID == 0 ) then + write(*,30) mumps_par%icntl(14) + 30 format(/'MUMPS percentage increase in the estimated working space', & + /'increased to',i4) + end if + + !else if (mumps_par%INFOG(1) == -13) then + ! if ( mumps_par%MYID == 0 ) then + ! write(*,40) + ! 40 format(/'MUMPS memory allocation error.'/) + ! end if + ! stop + + !else if (mumps_par%INFOG(1) == -10) then + ! if ( mumps_par%MYID == 0 ) then + ! write(*,45) + ! 45 format(/'MUMPS ERROR: Numerically singular matrix.'/) + ! end if + ! stop + + !else if (mumps_par%INFOG(1) == -40) then + ! if ( mumps_par%MYID == 0 ) then + ! write(*,46) + ! 46 format(/'MUMPS ERROR: matrix is not positive definite.'/) + ! end if + ! stop + + !else if ( mumps_par%INFOG(1) < 0 ) then + ! if ( mumps_par%MYID == 0 ) then + ! write(*,20) mumps_par%INFOG(1), mumps_par%INFOG(2) + ! 20 format(/'ERROR occured in MUMPS!',/,'INFOG(1), INFOG(2) ', 2i6,/) + ! end if + ! stop + + else + exit mumps_factorization ! factorization successful + end if + +end do mumps_factorization + + +if ( mumps_par%MYID == 0 ) then + !flush(ui_out) + ! Turn off mumps output. + mumps_par%icntl(2) = 0 ! output stream for diagnostic printing + mumps_par%icntl(3) = 0 ! output stream for global information +end if + + +return +end subroutine factor_matrix + +!-------------------------------------------------------------- + +subroutine solve( mumps_par, nrhs, rhs, x, transpose ) +! Solve A*x = rhs + +implicit none + +TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par +integer(kind=8),intent(in):: nrhs ! # of right-hand-sides +complex(kind=8),intent(in):: rhs(nrhs * mumps_par%N) +complex(kind=8),intent(out),target:: x(nrhs * mumps_par%N) ! solution +logical,intent(in):: transpose ! if .true. take the transpose + + x = rhs + + mumps_par%RHS => x + + ! The following is significant only on the host cpu. + mumps_par%NRHS = nrhs ! # of right-hand-sides + mumps_par%LRHS = mumps_par%N ! size of system + + if (transpose) then + mumps_par%icntl(9) = 0 ! for solving A'x = b + else + 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 + +!-------------------------------------------------------------- + +subroutine destroy(mumps_par) + +implicit none +TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par + +!if (.not. mumps_initialized) return + +! Destroy the instance (deallocate internal data structures) +mumps_par%JOB = -2 +CALL ZMUMPS(mumps_par) + +if (associated(mumps_par%A)) then + deallocate(mumps_par%IRN, mumps_par%JCN, mumps_par%A) +end if + +!if (mumps_par%MYID == 0) then +! close(ui_out) +!end if + +return +end subroutine destroy + + +end module mumps_cmplx_mod diff --git a/pymatsolver/mumps/mumps_interface.f90 b/pymatsolver/mumps/mumps_interface.f90 new file mode 100644 index 0000000..943124e --- /dev/null +++ b/pymatsolver/mumps/mumps_interface.f90 @@ -0,0 +1,192 @@ + +!------------------------------------------------- +!------------- REAL ------------ +!------------------------------------------------- + +subroutine factor_mumps( n, nnz, sym, A, jA, iA, ierr, pm_out ) + + +use mumps_mod, only: init, convert_to_mumps_format, factor_matrix, destroy +implicit none + +INCLUDE 'dmumps_struc.h' + +integer(kind=8),intent(out):: pm_out ! mumps pointer +integer(kind=8),intent(in):: n ! # of rows in A +integer(kind=8),intent(in):: nnz ! # of non zeros in A +integer(kind=8),intent(in):: sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. +real(kind=8),intent(in):: A(nnz) +integer(kind=8),intent(in):: jA(nnz), iA(n+1) + +integer(kind=8),intent(out):: ierr +! integer,intent(out):: ierr ! =0 no error, < 0 error (memory, singular, etc.) + +TYPE(DMUMPS_STRUC),pointer:: pmumps_par +TYPE(DMUMPS_STRUC):: mumps_par + +pointer ( pm, mumps_par ) + +allocate(pmumps_par) +pm = loc(pmumps_par) + +call init( mumps_par, sym) +call convert_to_mumps_format(mumps_par, n, A,jA,iA, ierr ) +if (ierr < 0) return ! allocation error + +call factor_matrix(mumps_par, ierr) + +if ( ierr < 0 ) then + ! Error in factorization. + call destroy(mumps_par) + pm_out = 0 +else + pm_out = pm +end if + +return +end subroutine factor_mumps + +!------------------------------------------------- + +subroutine solve_mumps( pm_in, nrhs, lrhs, rhs, x, transpose ) +! Solve A*x = rhs + +use mumps_mod, only: solve + +implicit none + +INCLUDE 'dmumps_struc.h' + +integer(kind=8),intent(in):: pm_in ! mumps pointer +integer(kind=8),intent(in):: nrhs ! # of right-hand-sides +integer(kind=8),intent(in):: lrhs ! # length of right-hand-sides +real(kind=8),intent(in):: rhs(lrhs) ! right-hand-side +real(kind=8),intent(out):: x(lrhs) ! solution +integer(kind=8),intent(in):: transpose ! =1 for transpose + +TYPE(DMUMPS_STRUC):: mumps_par +pointer ( pm, mumps_par ) +pm = pm_in + +call solve(mumps_par, nrhs, rhs, x, (transpose==1) ) + +return +end subroutine solve_mumps + +!------------------------------------------------- + +subroutine destroy_mumps( pm_in ) +! Destroy the instance (deallocate internal data structures) + +!DIR$ ATTRIBUTES DLLEXPORT :: destroy_mumps +!DIR$ ATTRIBUTES ALIAS: 'destroy_mumps_':: destroy_mumps +use mumps_mod, only: destroy + +implicit none +INCLUDE 'dmumps_struc.h' + +integer(kind=8),intent(in):: pm_in ! mumps pointer +TYPE(DMUMPS_STRUC):: mumps_par +pointer ( pm, mumps_par ) +pm = pm_in + +call destroy(mumps_par) + +return +end subroutine destroy_mumps + +!------------------------------------------------- +!------------- COMPLEX ------------ +!------------------------------------------------- + +subroutine factor_mumps_cmplx( n, nnz, sym, A, jA, iA, ierr, pm_out ) + +use mumps_cmplx_mod, only: init, convert_to_mumps_format, factor_matrix, destroy +implicit none + +INCLUDE 'zmumps_struc.h' + +integer(kind=8),intent(out):: pm_out ! mumps pointer +integer(kind=8),intent(in):: n ! # of rows in A +integer(kind=8),intent(in):: nnz ! # of non zeros in A +integer(kind=8),intent(in):: sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. +complex(kind=8),intent(in):: A(nnz) +integer(kind=8),intent(in):: jA(nnz), iA(n+1) + +integer(kind=8),intent(out):: ierr +! integer,intent(out):: ierr ! =0 no error, < 0 error (memory, singular, etc.) + +TYPE(ZMUMPS_STRUC),pointer:: pmumps_par +TYPE(ZMUMPS_STRUC):: mumps_par + +pointer ( pm, mumps_par ) + +allocate(pmumps_par) +pm = loc(pmumps_par) + +call init( mumps_par, sym) +call convert_to_mumps_format(mumps_par, n, A,jA,iA, ierr ) +if (ierr < 0) return ! allocation error + +call factor_matrix(mumps_par, ierr) + +if ( ierr < 0 ) then + ! Error in factorization. + call destroy(mumps_par) + pm_out = 0 +else + pm_out = pm +end if + +return +end subroutine factor_mumps_cmplx + +!------------------------------------------------- + +subroutine solve_mumps_cmplx( pm_in, nrhs, lrhs, rhs, x, transpose ) +! Solve A*x = rhs + +use mumps_cmplx_mod, only: solve + +implicit none + +INCLUDE 'zmumps_struc.h' + +integer(kind=8),intent(in):: pm_in ! mumps pointer +integer(kind=8),intent(in):: nrhs ! # of right-hand-sides +integer(kind=8),intent(in):: lrhs ! # length of right-hand-sides +complex(kind=8),intent(in):: rhs(lrhs) ! right-hand-side +complex(kind=8),intent(out):: x(lrhs) ! solution +integer(kind=8),intent(in):: transpose ! =1 for transpose + +TYPE(ZMUMPS_STRUC):: mumps_par +pointer ( pm, mumps_par ) +pm = pm_in + +call solve(mumps_par, nrhs, rhs, x, (transpose==1) ) + +return +end subroutine solve_mumps_cmplx + +!------------------------------------------------- + +subroutine destroy_mumps_cmplx( pm_in ) +! Destroy the instance (deallocate internal data structures) + +use mumps_cmplx_mod, only: destroy + +implicit none +INCLUDE 'zmumps_struc.h' + +integer(kind=8),intent(in):: pm_in ! mumps pointer +TYPE(ZMUMPS_STRUC):: mumps_par +pointer ( pm, mumps_par ) +pm = pm_in + +call destroy(mumps_par) + +return +end subroutine destroy_mumps_cmplx + +!------------------------------------------------- + diff --git a/pymatsolver/mumps/mumps_p.f90 b/pymatsolver/mumps/mumps_p.f90 new file mode 100755 index 0000000..710e0c8 --- /dev/null +++ b/pymatsolver/mumps/mumps_p.f90 @@ -0,0 +1,266 @@ + + + +module mumps_mod + + save + ! integer,parameter:: ui_out = 10 ! unit for mumps log file + !integer,parameter:: st_unit=11 ! unit for output statistics file + + + INCLUDE 'dmumps_struc.h' + + !logical,private:: mumps_initialized = .false. + + +contains + +!--------------------------------------------------------------- + +subroutine init( mumps_par, sym ) + +implicit none + +TYPE (DMUMPS_STRUC),intent(inout):: mumps_par +integer(kind=8),intent(in):: sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. + + !mumps_par%COMM = MPI_COMM_WORLD + mumps_par%SYM = sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. + mumps_par%PAR = 1 + + mumps_par%JOB = -1 ! initialize + CALL DMUMPS(mumps_par) + + + !if (mumps_par%MYID == 0) then + ! open(unit=ui_out, file='mumps.log', action='write') + !end if + + !mumps_initialized = .true. + +return +end subroutine init + +!--------------------------------------------------------------- + +subroutine convert_to_mumps_format( mumps_par, n, A,jA,iA, ierr ) +! A is actually a transpose. +implicit none + +TYPE (DMUMPS_STRUC),intent(inout):: mumps_par +integer(kind=8),intent(in):: n ! # of rows in A +real(kind=8),intent(in):: A(*) +integer(kind=8),intent(in):: jA(*), iA(n+1) +integer(kind=8),intent(out):: ierr + +integer nonz, i,j, j1,j2, ind, jcol, istat + + +nonz = iA(n+1) - 1 +!n = size(iA) - 1 + +mumps_par%N = n + +if (mumps_par%SYM == 0) then ! not symmetric matrix + mumps_par%NZ = nonz + allocate( mumps_par%IRN(nonz), mumps_par%JCN(nonz), mumps_par%A(nonz), stat=istat) + if ( istat /= 0 ) then + ierr = -13 ! allocation error + return + end if + + + do i = 1, n + mumps_par%JCN( iA(i) : iA(i+1)-1 ) = i + end do ! i + + mumps_par%IRN = jA(1:nonz) + mumps_par%A = A(1:nonz) + +else ! symmetric matrix + ! Keep only the lower half of the matrix (row >= column). + + nonz = nonz/2 + n ! should be +n/2, but I'm using n just in case. + allocate( mumps_par%IRN(nonz), mumps_par%JCN(nonz), mumps_par%A(nonz), stat=istat) + if ( istat /= 0 ) then + ierr = -13 ! allocation error + return + end if + + ind = 0 + do i = 1, n + + j1 = iA(i) + j2 = iA(i+1) - 1 + do j = j1, j2 + jcol = jA(j) + + if (i >= jcol) then + ind = ind + 1 + mumps_par%A(ind) = A(j) + mumps_par%JCN(ind) = jcol + mumps_par%IRN(ind) = i + end if + + end do ! j + end do ! i + + mumps_par%NZ = ind + !if (nonz < ind) call errorstop('nonz < ind') ! debug +end if + +ierr = 0 +return +end subroutine convert_to_mumps_format + +!--------------------------------------------------------------- + +subroutine factor_matrix( mumps_par, ierr ) + +implicit none + +TYPE (DMUMPS_STRUC),intent(inout):: mumps_par +integer(kind=8),intent(out):: ierr + +!mumps_par%icntl(2) = 6 ! output stream for diagnostics +mumps_par%icntl(4) = 0 ! 1 ! amount of output + + +mumps_par%icntl(2) = 0 ! ui_out ! output stream for diagnostic printing +mumps_par%icntl(3) = 0 ! ui_out ! output stream for global information + +!mumps_par%icntl(14) = 40 ! % increase in working space + +!mumps_par%icntl(11) = 0 ! return statistics +! mumps_par%icntl(7) = 5 ! ordering type +!mumps_par%cntl(1) = 0.d0 ! rel. threshold for numerical pivoting + + +mumps_factorization: do + + mumps_par%JOB = 1 ! analysis + CALL DMUMPS(mumps_par) + + mumps_par%JOB = 2 ! factorization + CALL DMUMPS(mumps_par) + + ierr = mumps_par%INFOG(1) + + + if ( mumps_par%INFOG(1) == -9 .or. mumps_par%INFOG(1) == -8 ) then + ! Main internal real/complex workarray S too small. + mumps_par%icntl(14) = mumps_par%icntl(14) + 10 + if ( mumps_par%MYID == 0 ) then + write(*,30) mumps_par%icntl(14) + 30 format(/'MUMPS percentage increase in the estimated working space', & + /'increased to',i4) + end if + + !else if (mumps_par%INFOG(1) == -13) then + ! if ( mumps_par%MYID == 0 ) then + ! write(*,40) + ! 40 format(/'MUMPS memory allocation error.'/) + ! end if + ! stop + + !else if (mumps_par%INFOG(1) == -10) then + ! if ( mumps_par%MYID == 0 ) then + ! write(*,45) + ! 45 format(/'MUMPS ERROR: Numerically singular matrix.'/) + ! end if + ! stop + + !else if (mumps_par%INFOG(1) == -40) then + ! if ( mumps_par%MYID == 0 ) then + ! write(*,46) + ! 46 format(/'MUMPS ERROR: matrix is not positive definite.'/) + ! end if + ! stop + + !else if ( mumps_par%INFOG(1) < 0 ) then + ! if ( mumps_par%MYID == 0 ) then + ! write(*,20) mumps_par%INFOG(1), mumps_par%INFOG(2) + ! 20 format(/'ERROR occured in MUMPS!',/,'INFOG(1), INFOG(2) ', 2i6,/) + ! end if + ! stop + + else + exit mumps_factorization ! factorization successful + end if + +end do mumps_factorization + + +if ( mumps_par%MYID == 0 ) then + !flush(ui_out) + ! Turn off mumps output. + mumps_par%icntl(2) = 0 ! output stream for diagnostic printing + mumps_par%icntl(3) = 0 ! output stream for global information +end if + + +return +end subroutine factor_matrix + +!-------------------------------------------------------------- + +subroutine solve( mumps_par, nrhs, rhs, x, transpose ) +! Solve A*x = rhs + +implicit none + +TYPE (DMUMPS_STRUC),intent(inout):: mumps_par +integer(kind=8),intent(in):: nrhs ! # of right-hand-sides +real(kind=8),intent(in):: rhs(nrhs * mumps_par%N) +real(kind=8),intent(out),target:: x(nrhs * mumps_par%N) ! solution +logical,intent(in):: transpose ! if .true. take the transpose + + x = rhs + + mumps_par%RHS => x + + ! The following is significant only on the host cpu. + mumps_par%NRHS = nrhs ! # of right-hand-sides + mumps_par%LRHS = mumps_par%N ! size of system + + if (transpose) then + mumps_par%icntl(9) = 0 ! for solving A'x = b + else + mumps_par%icntl(9) = 1 ! for solving Ax = b + end if + + + mumps_par%JOB = 3 ! Solve the system. + CALL DMUMPS(mumps_par) + ! At this point mumps_par%RHS (rhs) contains the solution. + + +return +end subroutine solve + +!-------------------------------------------------------------- + +subroutine destroy(mumps_par) + +implicit none +TYPE (DMUMPS_STRUC),intent(inout):: mumps_par + +!if (.not. mumps_initialized) return + +! Destroy the instance (deallocate internal data structures) +mumps_par%JOB = -2 +CALL DMUMPS(mumps_par) + +if (associated(mumps_par%A)) then + deallocate(mumps_par%IRN, mumps_par%JCN, mumps_par%A) +end if + +!if (mumps_par%MYID == 0) then +! close(ui_out) +!end if + +return +end subroutine destroy + + +end module mumps_mod diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py new file mode 100644 index 0000000..a617615 --- /dev/null +++ b/tests/test_Mumps.py @@ -0,0 +1,85 @@ +import unittest +import numpy as np +import scipy.sparse as sp +import pymatsolver + +TOL = 1e-11 + + +class TestMumps(unittest.TestCase): + + if pymatsolver.AvailableSolvers['Mumps']: + + def setUp(self): + 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] + a = np.r_[ + 3.0, -3.0, 2.0, 1.0, 3.0, 2.0, + 4.0, 2.0, 6.0, -1.0, 4.0, 1.0 + ] + rhs = np.r_[20.0, 24.0, 9.0, 6.0, 13.0] + rhs = np.c_[rhs, 10*rhs, 100*rhs] + 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 + + def test_1to5(self): + rhs = self.rhs + sol = self.sol + Ainv = pymatsolver.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) + + def test_1to5_cmplx(self): + rhs = self.rhs.astype(complex) + sol = self.sol.astype(complex) + self.A = self.A.astype(complex) + Ainv = pymatsolver.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) + + def test_1to5_T(self): + rhs = self.rhs + sol = self.sol + Ainv = pymatsolver.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) + + def test_singular(self): + A = sp.identity(5).tocsr() + A[-1, -1] = 0 + self.assertRaises(Exception, pymatsolver.Mumps, A) + + def test_multiFactorsInMem(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)] + + for Ainv in solvers: + self.assertLess( + 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 + ) + +if __name__ == '__main__': + unittest.main() From 614d54f705ff05cada0ffa5e3aed386e1c51057f Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 7 Jan 2017 20:42:19 -0800 Subject: [PATCH 2/6] travis and failure for mumps --- .travis.yml | 11 +- pymatsolver/mumps/__init__.py | 194 +++++++++++++++++----------------- 2 files changed, 108 insertions(+), 97 deletions(-) diff --git a/.travis.yml b/.travis.yml index b11c483..35c949b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -8,6 +8,15 @@ matrix: allow_failures: - python: 3.6 sudo: false +addons: + apt: + packages: + - gcc + - gfortran + - libopenmpi-dev + - libmumps-seq-dev + - libblas-dev + - liblapack-dev before_install: - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh @@ -21,8 +30,8 @@ install: - pip install nose-cov python-coveralls pyMKL - git clone https://github.com/simpeg/discretize.git; cd discretize; python setup.py install; cd ..; -- python setup.py install - python setup.py build_ext --inplace +- make mumps script: - nosetests --with-cov --cov pymatsolver --cov-config .coveragerc -v -s after_success: diff --git a/pymatsolver/mumps/__init__.py b/pymatsolver/mumps/__init__.py index a554a0e..926d3fa 100644 --- a/pymatsolver/mumps/__init__.py +++ b/pymatsolver/mumps/__init__.py @@ -1,9 +1,12 @@ from __future__ import unicode_literals from __future__ import print_function from __future__ import division -# from __future__ import absolute_import -import MumpsInterface +try: + import MumpsInterface as _MUMPSINT + _mumpsExists = True +except ImportError: + _mumpsExists = False import gc from pymatsolver.solvers import Base @@ -60,112 +63,111 @@ } -# if _mumpsExists: +if _mumpsExists: -class _Pointer(object): - """Gets an int and a destroy call that gets called on garbage collection. + 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. + 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. + 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 __init__(self, pointerINT, destroyCall): - self.INT = pointerINT - self.destroyCall = destroyCall - - def __del__(self): - self.destroyCall(self.INT) - - -class Mumps(Base): - """ - - documentation:: + You can always clean the Solver, and it will explicitly call the destroy method. + """ + def __init__(self, pointerINT, destroyCall): + self.INT = pointerINT + self.destroyCall = destroyCall - http://mumps.enseeiht.fr/doc/userguide_4.10.0.pdf + def __del__(self): + self.destroyCall(self.INT) - """ + class Mumps(Base): + """ - transpose = False - symmetric = False + documentation:: - @property - def T(self): - newMS = self.__class__( - self.A, - symmetric=self.symmetric, - fromPointer=self.pointer - ) - newMS.transpose = not self.transpose - return newMS + http://mumps.enseeiht.fr/doc/userguide_4.10.0.pdf - def __init__(self, A, symmetric=False, fromPointer=None): - self.A = A.tocsc() - self.symmetric = symmetric + """ - if fromPointer is None: + transpose = False + symmetric = False + + @property + def T(self): + newMS = self.__class__( + self.A, + symmetric=self.symmetric, + fromPointer=self.pointer + ) + newMS.transpose = not self.transpose + return newMS + + def __init__(self, A, symmetric=False, fromPointer=None): + self.A = A.tocsc() + self.symmetric = symmetric + + if fromPointer is None: + self.factor() + elif isinstance(fromPointer, _Pointer): + self.pointer = fromPointer + else: + raise Exception('Unknown pointer for construction.') + + @property + def isfactored(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'] + + means factor, solve, destroy + """ + if self.A.dtype == float: + return {'F': _MUMPSINT.factor_mumps, + 'S': _MUMPSINT.solve_mumps, + 'D': _MUMPSINT.destroy_mumps}[ftype] + elif self.A.dtype == complex: + return {'F': _MUMPSINT.factor_mumps_cmplx, + 'S': _MUMPSINT.solve_mumps_cmplx, + 'D': _MUMPSINT.destroy_mumps_cmplx}[ftype] + + def factor(self): + if self.isfactored: + return + + sym = 1 if self.symmetric else 0 + ierr, p = self._funhandle('F')( + sym, + self.A.data, + self.A.indices+1, + self.A.indptr+1 + ) + if ierr < 0: + raise Exception("Mumps Exception [{}] - {}".format(ierr, _mumpsErrors[ierr])) + elif ierr > 0: + print("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) + + self.pointer = _Pointer(p, self._funhandle('D')) + + def _solveM(self, rhs): self.factor() - elif isinstance(fromPointer, _Pointer): - self.pointer = fromPointer - else: - raise Exception('Unknown pointer for construction.') + 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 - @property - def isfactored(self): - return getattr(self, 'pointer', None) is not None + _solve1 = _solveM - def _funhandle(self, ftype): - """ - switches the function handle between real and complex. + def clean(self): + self._funhandle('D')(self.pointer.INT) + del self.pointer + gc.collect() - ftype in ['F','S','D'] - means factor, solve, destroy - """ - if self.A.dtype == float: - return {'F': MumpsInterface.factor_mumps, - 'S': MumpsInterface.solve_mumps, - 'D': MumpsInterface.destroy_mumps}[ftype] - elif self.A.dtype == complex: - return {'F': MumpsInterface.factor_mumps_cmplx, - 'S': MumpsInterface.solve_mumps_cmplx, - 'D': MumpsInterface.destroy_mumps_cmplx}[ftype] - - def factor(self): - if self.isfactored: - return - - sym = 1 if self.symmetric else 0 - ierr, p = self._funhandle('F')( - sym, - self.A.data, - self.A.indices+1, - self.A.indptr+1 - ) - if ierr < 0: - raise Exception("Mumps Exception [{}] - {}".format(ierr, _mumpsErrors[ierr])) - elif ierr > 0: - print("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) - - self.pointer = _Pointer(p, self._funhandle('D')) - - 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 - - _solve1 = _solveM - - def clean(self): - self._funhandle('D')(self.pointer.INT) - del self.pointer - gc.collect() - - -# del _mumpsExists +del _mumpsExists From 69f8dce8d4d402bcfc00340d18e8da327e89fe06 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 7 Jan 2017 20:56:53 -0800 Subject: [PATCH 3/6] Relative imports for the Mumps --- pymatsolver/mumps/Makefile | 1 + pymatsolver/mumps/__init__.py | 188 ++++++++++++++++------------------ tests/test_Mumps.py | 8 +- 3 files changed, 95 insertions(+), 102 deletions(-) diff --git a/pymatsolver/mumps/Makefile b/pymatsolver/mumps/Makefile index f6a809b..e6501e4 100644 --- a/pymatsolver/mumps/Makefile +++ b/pymatsolver/mumps/Makefile @@ -40,6 +40,7 @@ build: clean: rm -f *.o *.mod *.so + rm -rf *.dSYM test: cd ..;nosetests diff --git a/pymatsolver/mumps/__init__.py b/pymatsolver/mumps/__init__.py index 926d3fa..6826e25 100644 --- a/pymatsolver/mumps/__init__.py +++ b/pymatsolver/mumps/__init__.py @@ -1,14 +1,11 @@ from __future__ import unicode_literals from __future__ import print_function from __future__ import division +from __future__ import absolute_import -try: - import MumpsInterface as _MUMPSINT - _mumpsExists = True -except ImportError: - _mumpsExists = False import gc from pymatsolver.solvers import Base +from . import MumpsInterface as _MUMPSINT _mumpsErrors = { -1: "An error occurred on processor INFO(2).", @@ -63,111 +60,106 @@ } -if _mumpsExists: +class _Pointer(object): + """Gets an int and a destroy call that gets called on garbage collection. - 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. - 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. - 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 __init__(self, pointerINT, destroyCall): + self.INT = pointerINT + self.destroyCall = destroyCall - You can always clean the Solver, and it will explicitly call the destroy method. - """ - def __init__(self, pointerINT, destroyCall): - self.INT = pointerINT - self.destroyCall = destroyCall + def __del__(self): + self.destroyCall(self.INT) - def __del__(self): - self.destroyCall(self.INT) +class Mumps(Base): + """ - class Mumps(Base): - """ + documentation:: - documentation:: + http://mumps.enseeiht.fr/doc/userguide_4.10.0.pdf - http://mumps.enseeiht.fr/doc/userguide_4.10.0.pdf + """ - """ + transpose = False + symmetric = False + + @property + def T(self): + newMS = self.__class__( + self.A, + symmetric=self.symmetric, + fromPointer=self.pointer + ) + newMS.transpose = not self.transpose + return newMS - transpose = False - symmetric = False - - @property - def T(self): - newMS = self.__class__( - self.A, - symmetric=self.symmetric, - fromPointer=self.pointer - ) - newMS.transpose = not self.transpose - return newMS - - def __init__(self, A, symmetric=False, fromPointer=None): - self.A = A.tocsc() - self.symmetric = symmetric - - if fromPointer is None: - self.factor() - elif isinstance(fromPointer, _Pointer): - self.pointer = fromPointer - else: - raise Exception('Unknown pointer for construction.') - - @property - def isfactored(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'] - - means factor, solve, destroy - """ - if self.A.dtype == float: - return {'F': _MUMPSINT.factor_mumps, - 'S': _MUMPSINT.solve_mumps, - 'D': _MUMPSINT.destroy_mumps}[ftype] - elif self.A.dtype == complex: - return {'F': _MUMPSINT.factor_mumps_cmplx, - 'S': _MUMPSINT.solve_mumps_cmplx, - 'D': _MUMPSINT.destroy_mumps_cmplx}[ftype] - - def factor(self): - if self.isfactored: - return - - sym = 1 if self.symmetric else 0 - ierr, p = self._funhandle('F')( - sym, - self.A.data, - self.A.indices+1, - self.A.indptr+1 - ) - if ierr < 0: - raise Exception("Mumps Exception [{}] - {}".format(ierr, _mumpsErrors[ierr])) - elif ierr > 0: - print("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) - - self.pointer = _Pointer(p, self._funhandle('D')) - - def _solveM(self, rhs): + def __init__(self, A, symmetric=False, fromPointer=None): + self.A = A.tocsc() + self.symmetric = symmetric + + if fromPointer is None: 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 + elif isinstance(fromPointer, _Pointer): + self.pointer = fromPointer + else: + raise Exception('Unknown pointer for construction.') - _solve1 = _solveM + @property + def isfactored(self): + return getattr(self, 'pointer', None) is not None - def clean(self): - self._funhandle('D')(self.pointer.INT) - del self.pointer - gc.collect() + def _funhandle(self, ftype): + """ + switches the function handle between real and complex. + ftype in ['F','S','D'] -del _mumpsExists + means factor, solve, destroy + """ + if self.A.dtype == float: + return {'F': _MUMPSINT.factor_mumps, + 'S': _MUMPSINT.solve_mumps, + 'D': _MUMPSINT.destroy_mumps}[ftype] + elif self.A.dtype == complex: + return {'F': _MUMPSINT.factor_mumps_cmplx, + 'S': _MUMPSINT.solve_mumps_cmplx, + 'D': _MUMPSINT.destroy_mumps_cmplx}[ftype] + + def factor(self): + if self.isfactored: + return + + sym = 1 if self.symmetric else 0 + ierr, p = self._funhandle('F')( + sym, + self.A.data, + self.A.indices+1, + self.A.indptr+1 + ) + if ierr < 0: + raise Exception("Mumps Exception [{}] - {}".format(ierr, _mumpsErrors[ierr])) + elif ierr > 0: + print("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) + + self.pointer = _Pointer(p, self._funhandle('D')) + + 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 + + _solve1 = _solveM + + def clean(self): + self._funhandle('D')(self.pointer.INT) + del self.pointer + gc.collect() diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index a617615..4eacc88 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -59,10 +59,10 @@ def test_1to5_T(self): ) self.assertLess(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 + # self.assertRaises(Exception, pymatsolver.Mumps, A) def test_multiFactorsInMem(self): n = 100 From 0932d7432a61497219c8c98659eb85ef8d220480 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 7 Jan 2017 21:06:30 -0800 Subject: [PATCH 4/6] bringing back the readme for mumps install --- README.rst | 53 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/README.rst b/README.rst index 1968e4d..d611982 100644 --- a/README.rst +++ b/README.rst @@ -37,6 +37,59 @@ All solvers work with :code:`scipy.sparse` matricies, and a single or multiple r * L/U Triangular Solves * Wrapping of SciPy matrix solvers (direct and indirect) * Pardiso solvers now that MKL comes with conda! +* Mumps solver with nice error messages + + +Installing Mumps +================ + +I have not been able to get the pip install to work because of multiple dependencies on fortran libraries. +However, the linux and mac installs are relatively easy. Note that you must have mumps pre-installed, +currently I have only got this working for the sequential version, so when you are installing, +you will need to point to that one. You can also look at the `.travis.yml` file for how to get it working on TravisCI. + +Linux +----- + +From a clean install on Ubuntu: + +.. code-block:: bash + + apt-get update + apt-get -y install gcc gfortran git libopenmpi-dev libmumps-seq-dev libblas-dev liblapack-dev + + # Install all the python you need! + wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh -O miniconda.sh; + chmod +x miniconda.sh + ./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 + + git clone https://github.com/rowanc1/pymatsolver.git + cd pymatsolver + make mumps + +Mac +--- + +This assumes that you have Brew and some python installed (numpy, scipy): + +.. code-block:: bash + + brew install mumps --with-scotch5 --without-mpi + + git clone https://github.com/rowanc1/pymatsolver.git + cd pymatsolver + make mumps_mac + +If you have problems you may have to go into the Makefile and update the pointers to Lib and Include for the various libraries. + +This command is helpful for finding dependencies. You should also take note of have happens when brew installs mumps. + +.. code-block:: bash + + mpicc --showme Code: From dbf16993bd772d5b78e0c47c9311d200114316b3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 7 Jan 2017 21:12:32 -0800 Subject: [PATCH 5/6] Add two shortcuts for backwards compatibility. --- pymatsolver/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pymatsolver/__init__.py b/pymatsolver/__init__.py index db3bc41..603ef2b 100644 --- a/pymatsolver/__init__.py +++ b/pymatsolver/__init__.py @@ -27,6 +27,7 @@ try: from pymatsolver.direct import Pardiso AvailableSolvers['Pardiso'] = True + ParadisoSolver = Pardiso # backwards compatibility except ImportError: SolverHelp['Pardiso'] = """Pardiso is not working @@ -37,6 +38,7 @@ try: from pymatsolver.mumps import Mumps AvailableSolvers['Mumps'] = True + MumpsSolver = Mumps # backwards compatibility except ImportError: SolverHelp['Mumps'] = """Mumps is not working. From 12cbab0fcae1fc2aa0c2c20fec173a09a383a4cb Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 7 Jan 2017 21:12:54 -0800 Subject: [PATCH 6/6] =?UTF-8?q?Bump=20version:=200.0.3=20=E2=86=92=200.1.0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- docs/conf.py | 4 ++-- pymatsolver/__init__.py | 2 +- setup.py | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index e7121ea..de3a7e3 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.0.3 +current_version = 0.1.0 files = setup.py pymatsolver/__init__.py docs/conf.py diff --git a/docs/conf.py b/docs/conf.py index 67d417a..d602818 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -60,9 +60,9 @@ # built documents. # # The short X.Y version. -version = u'0.0.3' +version = u'0.1.0' # The full version, including alpha/beta/rc tags. -release = u'0.0.3' +release = u'0.1.0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/pymatsolver/__init__.py b/pymatsolver/__init__.py index 603ef2b..d612936 100644 --- a/pymatsolver/__init__.py +++ b/pymatsolver/__init__.py @@ -55,7 +55,7 @@ mpicc --showme """ -__version__ = '0.0.3' +__version__ = '0.1.0' __author__ = 'Rowan Cockett' __license__ = 'MIT' __copyright__ = 'Copyright 2017 Rowan Cockett' diff --git a/setup.py b/setup.py index 074f4b3..75e42ad 100644 --- a/setup.py +++ b/setup.py @@ -29,12 +29,12 @@ setup( name="pymatsolver", - version="0.0.3", + version="0.1.0", packages=find_packages(), install_requires=[ 'numpy>=1.7', 'scipy>=0.13', - 'pyMKL>=0.0.3' + 'pyMKL>=0.1.0' ], author="Rowan Cockett", author_email="rowanc1@gmail.com",