Improved applyMUMPS type stability

[Pbellive]

Hi Lars, I've changed applyMUMPS to make it more type stable. It now has two methods one, for Float64 matrices and the other for Complex128 matrices. Each method branches on the rhs type, calling different fortran functions for sparse and dense rhs. That branching doesn't seem to introduce any type instability. There is one remaining bit of type instability in the function that only seems to occur when the factor is stored on a different worker than solve is called from. And even then I think it's a minor instability that probably shouldn't be fixed. After factoring a real SPD matrix and generating a random right hand side, running pre>@code_warntype</pre give the output:

Variables:
  factor::MUMPS.MUMPSfactorizationReal
  rhs::Array{Float64,1}
  ##dims#6650::Tuple{Int64,Int64}

Body:
  begin  # /home/patrick/Julia/v0.4/MUMPS.jl/src/MUMPSfuncs.jl, line 93:
      return (MUMPS.applyMUMPS)(factor::MUMPS.MUMPSfactorizationReal,rhs::Array{Float64,1},(Base.fill!)((Base.Array)(Base.Float64,0,0)::Array{Float64,2},(Base.box)(Float64,(Base.sitofp)(Float64,0)::Any)::Float64)::Array{Float64,2},0)::Array{Float64,2}
  end::Array{Float64,2}

My changes slightly reduced the flexibility of the API. In my current version, for complex matrices, the rhs must be of type Array{Complex128}. Also, the syntax x = applyMUMPS(factor,rhs,[]) is no longer supported. Instead use applyMUMPS(factor,rhs) or applyMUMPS(factor,rhs,zeros(eltype(rhs),0,0)).

[lruthotto]

Hey Patrick, I tried something based on your new version. Basically I introduced a function applyMUMPS! that overwrites a provided vector or array x. You can take a look in : https://github.com/lruthotto/MUMPS.jl which is a fork of your branch. let me know what you think.

[Pbellive]

Hi Lars, I merged your changes into my MUMPS.jl repo and tested them. I really like the new way you've structured the code. It's very clear and eliminates any branching based on type. I also added a string constant at the top of MUMPSfuncs.jl where the path to the MUMPS library is written so that the path only has to be changed in one place in MUMPSfuncs.jl. I looked a little further into calling MUMPS directly from Julia with no glue code. I think it's impossible with julia 0.4 due to limitations in how ccall handles composite types. I asked a question about it on the julia-users google group in case I'm just missing something in my understanding of ccall but I'm not optimistic.

[Pbellive]

Hi Lars, I merged all your changes into my MUMPS repo and updated it to use MUMPS 5.0.1 and Roman's new Makefile. I also added type assertions for the integer input arguments (sym, tr, and ooc). Without the type assertions, calling with the old syntax solveMUMPS(factorization,rhs,[]) yields a cryptic error message from convert rather than a no method error. I updated the benchmark codes so that they no longer use that syntax. If you accept the pull request as is, it should give the latest MUMPS library, the latest MUMPSfuncs.jl and Roman's simplified makefile.

[lruthotto]

Hi Patrick, thanks a lot. This looks really nice. However, when building it (I did not change anything in the makefile), I get the following error. Can you help:

lruthot$ make
(make -C /Users/lruthot/Projects/MUMPS.jl/src/metis-5.1.0 config )
rm -rf build/Darwin-x86_64
mkdir -p build/Darwin-x86_64
cd build/Darwin-x86_64 && cmake /Users/lruthot/Projects/MUMPS.jl/src/metis-5.1.0 -DCMAKE_VERBOSE_MAKEFILE=1 -DGKLIB_PATH=/Users/lruthot/Projects/MUMPS.jl/src/metis-5.1.0/GKlib
-- The C compiler identification is AppleClang 6.0.0.6000057
-- The CXX compiler identification is AppleClang 6.0.0.6000057
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Looking for execinfo.h
-- Looking for execinfo.h - found
-- Looking for getline
-- Looking for getline - found
-- checking for thread-local storage - found
-- Configuring done
-- Generating done
-- Build files have been written to: /Users/lruthot/Projects/MUMPS.jl/src/metis-5.1.0/build/Darwin-x86_64
(make -C /Users/lruthot/Projects/MUMPS.jl/src/metis-5.1.0 )
/bin/sh: -c: line 1: syntax error: unexpected end of file
make[1]: *** [all] Error 2
make: *** [all] Error 2

[Pbellive]

Hi Lars, I just tested a fresh clone of my repo and it builds correctly for me on Linux. I think I know why it's not working for you on Mac though. It looks like Roman edited the Metis Makefile to a form that his and my machines are able to parse but yours is not. Try deleting lines 67 and 68 of src/metis-5.1.0/Makefile (which are commented out) and building again. After you delete those comments, lines 63-68 of the metis Makefile should read

all clean install:
    @if [ ! -f $(BUILDDIR)/Makefile ]; then \
        more BUILD.txt; \
    else \
        make -C $(BUILDDIR) $@; \
    fi

If that works let me know and I'll push that change in the metis makefile to my focus-shift/MUMPS.jl repo.

[lruthotto]

Thanks Patrick, the code compiles now, after following your instructions and also changing the path in LMETISDIR. For Mac Users the libmetis is in a different subfolders. We should probably be able to check this in the makefile.

I also changed the way the pre-allocated output x is handled in applyMUMPS. My change is motivated by the fact that we often want to manipulate tr but not specify x. This can now be done be setting x=zeros(0). (BTW: keyword arguments are no solution here as then we could not do any remotecalls on applyMUMPS anymore)

You can merge the code from my branch!

[Pbellive]

Hi Lars, my focus-shift/MUMPS.jl repo has now incorporated your change in applyMUMPS and deleted the commented lines from the metis Makefile. I also edited the root Makefile so that it will choose the correct metis build directory based on the operating system. I tested this and it works on Ubuntu. Would you mind testing it on OS X? I believe the user should now be able to build without manually setting the metis build directory as long as he or she is on Linux or OS X with an x86_64 processor.

[lruthotto]

Thanks Patrick, I just tested it and everything runs fine for me! I like that the makefile became even simpler. Great job!