JuliaInv / MUMPSjInv.jl

MUMPS interface for Julia
18 stars 15 forks source link

Updated MUMPS.jl to use MUMPS_5.0.0, metis-5.1.0, and added sparse rhs capability #6

Closed Pbellive closed 8 years ago

Pbellive commented 9 years ago

Hi Lars, I've found your MUMPS package very useful. Thanks a lot for putting it together! I forked the JuliaSparse/MUMPS.jl repo and made some changes. I've submitted a pull request in case you're interested in merging them into your code. My main changes were:

1) Added source code for MUMPS_5.0.0 and metis-5.1.0 and updated the Makefile to build MUMPS_5.0.0 and metis-5.1.0 as static libraries then build the julia interface as a shared library.

2) Changed the build procedure as follows: User must add their BLAS library location in src/MUMPS_5.0.0/Makefile.inc instead of in the main Makefile. This allows MUMPS example programs to build correctly.

3) Changed build settings to compile MUMPS with 64 bit integers. Since MUMPS is coded with 32 bit integers this needs to be done (as far as I know) by compiling with 64 bit as the default integer type (e.g. by using the compiler flag -fdefault-integer-8 in gfortran). I also built metis with 64 bit integers. In order to have the Julia interface Fortran code work properly with the 64bit integer MUMPS, I changed integer(kind=8) declarations to default integers and compiled the interface code with -fdefault-integer-8.

4) Added sparse right hand side capabilities without changing MUMPS.jl api. I wrote new Fortran interface routines that use MUMPS' sparse rhs functionality (icntl(20)=1). The right hand side(s) are input to MUMPS as a CSC sparse matrix and MUMPS decides automatically whether to exploit sparsity in the solution procedure. The user calls solveMUMPS or applyMUMPS as normal and the package dispatches based on right hand side type (Array or SparseMatrixCSC). For extremely sparse right hand sides, this resulted in significant speedup (>factor of 2) compared to dense rhs. I've included some benchmarking results in the file sparse-rhs-benchmark-results.txt. I also added tests of the sparse rhs functionality in src/tests.

Let me know if you have any questions about the changes I made.

All the best, Patrick

lruthotto commented 9 years ago

Thanks, Patrick! This looks very cool. I would like to test it before merging the pull request. When trying to install it on my Mac, I somehow get a weird warning (in the step where I try to build MUMPS).

lruthot$ make
/Applications/Xcode.app/Contents/Developer/usr/bin/make ARITH=d mumps_lib
(cd libseq; /Applications/Xcode.app/Contents/Developer/usr/bin/make)
make[2]: Nothing to be done for `all'.
if [ "./PORD/lib/" != "" ] ; then \
      cd ./PORD/lib/; \
      /Applications/Xcode.app/Contents/Developer/usr/bin/make CC="gcc" CFLAGS="-O3 -fPIC -DINTSIZE64" AR="ar vr " RANLIB="echo" OUTC="-o" LIBEXT=.a; \
    fi;
make[2]: `libpord.a' is up to date.
if [ "./PORD/lib/" != "" ] ; then \
      cp ./PORD/lib//libpord.a lib/libpord.a; \
    fi;
cp: lib/libpord.a: No such file or directory
make[1]: *** [lib/libpord.a] Error 1
make: *** [d] Error 2

I probably just need to adjust an option here or there, but couldn't find it. Do you have an idea?

Thanks. Lars

Pbellive commented 8 years ago

Hi Lars, I'm sorry, the build failure was due to my sloppy use of git. I forgot that git doesn't track empty directories so the MUMPS_5.0.0/lib folder didn't get pushed to github from my local repo. The MUMPS makefile seems to require that directory to exist before build time. I've added a .gitignore file to the folder so it's now tracked. I also noticed another error in the metis build settings. I've corrected those two errors and pushed them to my fork of your repo. I'm now able to successfully build and test a fresh clone of my github MUMPS repo on kubuntu. As long as you're compiling with gcc and gfortran you should now be able to build successfully without changing anything but the path to your MUMPS.jl/src directory in the main makefile and the path to your blas libraries in MUMPS_5.0.0/Makefile.inc. Let me know if you have any other issues!

Patrick

P.S. I didn't mention before, I'm a Ph.D. student with Eldad Haber at UBC. Eldad mentioned that you'll be in Vancouver soon. We could chat about this then as well.

lruthotto commented 8 years ago

Thanks for the hint. I was now able to compile MUMPS on my system. However, I am running into Segmentation Faults in all test files.

               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "help()" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.4.0-dev+5576 (2015-06-24 15:42 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit 92f8c17* (16 days old master)
|__/                   |  x86_64-apple-darwin13.4.0

julia> include("testDivGrad.jl")
Test for real SPD matrix: one rhs

signal (11): Segmentation fault: 11
libmetis__CreateCoarseGraphNoMask at /tmp/MUMPS.jl/lib/MUMPS (unknown line)
Segmentation fault: 11

Any ideas? Of course we can also try to sort things out next week in Vancouver.

Pbellive commented 8 years ago

Hi Lars, I'm not positive but my guess is that the metis integer type is set incorrectly in your build. Please insure that line 33 of .../MUMPS.jl/src/metis-5.1.0/include/metis.h is

define IDXTYPEWIDTH 64

This sets metis to be built with 64 bit integers. I can reproduce your segmentation fault on my laptop when compiling with

define IDXTYPEWIDTH 32

IDXTYPEWIDTH is set correctly on my fork now but for some reason when I initially pushed my local changes to github and sent the PR the incorrect version of metis.h was sent to github. That's the second error I mentioned in my last message. It'll teach me to submit a pull request before testing a fresh clone of my own github repo rather than just doing the tests on my local branch.

Did you download a fresh clone of my fork before your most recent round of tests or did you just add the MUMPS_5.0.0/lib directory manually on your machine? If you just added the directory manually then the segmentation fault you experienced would be expected. I'm sorry if that wasn't clear from my last message. If the IDXTYPEWIDTH setting is not your problem then I'm not sure what could be causing the segfault.

Patrick

lruthotto commented 8 years ago

Thanks, Patrick. With the new version everything seems to work. Before I merge this I want to do some more extensive testing and compare some timings. Will merge this, in the next couple of days.

looking forward to meeting you in Vancouver. Lars