Matlab mex interface to external direct solver

Hi All,

I am trying to use an external direct solver in Matlab. The matrices I am dealing with are sparse and unsymmetric, and the size can be large (several hundred thousands). Further, the problem is dynamic and nonlinear and so the matrix formation, assembly and solution occur multiple times.

Matlab's internal direct solver uses UMFPACK and seems to require very large amounts of memory (In this case, it requires more than 24GB of memory even for a matrix of size 60,000. So for memory usage efficiency and solution times to be reasonable, I am trying to use a direct external solver (MUMPS) with Matlab.

I checked the Matlab interface that comes with the MUMPS distribution and it works fine. However, it only works in sequential mode. Just to get a feel for it I wrote my own simple Matlab mex interface in Fortran 90 for sequential MUMPS and it works fine in comparison with Matlab's in-built solver. The problem starts when I try to compile a parallel version of MUMPS to work with Matlab. The parallel version of MUMPS works perfectly well when I call MUMPS from a FORTRAN program. But every time I interface it with Matlab, Matlab crashes during the matrix factorization stage.

I have tried compiling the parallel version of MUMPS with the following libraries that are required:

1. BLAS
2. LAPACK (version used): 3.4.2
3. SCALAPACK (version used): 2.02 (BLACS is built-in to this)
4. MPICH2 (version used): 1.2.1

There can possibly be a conflict of the libraries being used by Matlab and MUMPS. I am not sure how to address this issue, if it is indeed the case. Please find my simple mex code and the crash dump log below.

Thanks in advance for your help and your time. Any suggestion is appreciated.

Best Regards,

Anand

************************

Mex Code in Fortran 90: ( I am new to Fortran 90, so my code can be flawed and redundant in places).

************************

#include  "fintrf.h"

!     Gateway routine
SUBROUTINE MEXFUNCTION(nlhs, plhs, nrhs, prhs)

#include "mpif.h"
#include  "dmumps_struc.h"
type (dmumps_struc) id

integer*4 :: ierr

!     mexFunction arguments:
mwPointer :: plhs(*), prhs(*)
integer*4  :: nlhs, nrhs

!     Function declarations:
mwPointer  :: mxCreateDoubleMatrix, mxGetPr
integer*4  :: mxIsNumeric
mwSize     :: mxGetM, mxGetN
!mwPointer :: mxCreateSparse

!     Pointers to input/output mxArrays:
real*8, pointer :: A_pr, b_pr, x_pr

!     Array information:
mwSize :: mA, nA
!mwSize :: sizeA
mwSize, external :: mxGetNzmax
mwSize ::nzA

!mwPointer :: mxGetIr, mxGetJc
mwPointer :: IrA, JcA

integer*4 :: mexPrintf
integer*4 :: kwrite
character(len=180) :: line

integer*8 nrows, ncols, ctyp
!mwPointer mxCalloc
integer*8 i

integer*4, target, allocatable :: IRN(:), JCN(:)
real*8, target, allocatable :: A(:), RHS(:)

integer*4 :: mxIsSparse, stat, error
!-----------------------------------------------------------------------
!     Check for proper number of arguments. 
if(nrhs /= 4) then
 call mexErrMsgIdAndTxt ('MATLAB:matsq:nInput', &
                           'Four inputs required.')
elseif(nlhs /= 1) then
 call mexErrMsgIdAndTxt ('MATLAB:matsq:nOutput', &
                           'One output required.')
endif

!     Get the size of the input array.
mA = mxGetM(prhs(1))
nA = mxGetN(prhs(1))
!sizeA = mA*nA

nzA = mxGetNzmax(prhs(1))

!     Check that the array is numeric (not strings).
if(mxIsNumeric(prhs(1)) .eq. 0) then
 call mexErrMsgIdAndTxt ('MATLAB:matsq:NonNumeric', &
                           'Input must be a numeric array.')
endif

nrows=mA
ncols=1
ctype=0
!     Create matrix for the return argument.
plhs(1) = mxCreateDoubleMatrix(nrows,ncols,ctype)

A_pr = mxGetPr(prhs(1))
b_pr = mxGetPr(prhs(2))
x_pr = mxGetPr(plhs(1))

!IrA = mxGetIr(prhs(1))
!JcA = mxGetJc(prhs(1))

IrA = mxGetPr(prhs(3))
JcA = mxGetPr(prhs(4))

!---------------------------
!MUMPS routine for solution
!---------------------------
call mpi_init(ierr)

! Define a communicator for the package
id%comm = mpi_comm_world

! Unsymmetric code
id%SYM=0
id%PAR=1

!Initialize an instance of the MUMPS package
id%JOB=-1
call dmumps(id)

!id%ICNTL(7) = 5     ! Sequential analysis (METIS ordering for value=5)
!id%ICNTL(28) = 2    ! Parallel analysis
!id%ICNTL(29) = 2    !ParMETIS ordering
!id%ICNTL(6) = 1
!id%ICNTL(8) = 7
!id%ICNTL(14) = 80

if (id%MYID == 0) then

!Matrix parameters
id%N = mA               ! Matrix Size
id%NZ = nzA            ! Number of non-zeros

!------------------------------------------------------------------------
! Allocate memory
!------------------------------------------------------------------------
allocate(IRN(nzA),stat=error)
if (stat.ne.0) then
call mexErrMsgIdAndTxt ('error: couldn nott allocate memory for array IRN')
endif

allocate(JCN(nzA),stat=error)
if (stat.ne.0) then
call mexErrMsgIdAndTxt ('error: could not allocate memory for array JCN')
endif

allocate(A(nzA), stat=error)
if (stat.ne.0) then
call mexErrMsgIdAndTxt ('error: could not allocate memory for array A')
endif

allocate(RHS(mA), stat=error)
if (stat.ne.0) then
call mexErrMsgIdAndTxt ('error: could not allocate memory for array RHS')
endif

! Get data in Fortran format
call mxCopyPtrToInteger4(IrA, IRN, nzA)
call mxCopyPtrToInteger4(JcA, JCN, nzA)
call mxCopyPtrToReal8(A_pr, A, nzA)
call mxCopyPtrToReal8(b_pr, RHS, mA)

!write(line,*)(IRN(i), i = 1,12), nzA
!kwrite=mexPrintf(line//achar(10))
!write(line,*)(JCN(i), i = 1,12)
!kwrite=mexPrintf(line//achar(10))

! Assign data to Mumps structure
id%IRN => IRN
id%JCN => JCN
id%A => A
id%RHS => RHS
endif

!Solution (NOTE: MUMPS replaces the id.RHS with the computed solution)
id%JOB = 6
call dmumps(id)

if (id%MYID == 0) then
write(line,*)(id%INFOG(i), i = 1,2)
kwrite=mexPrintf(line//achar(10))

call mxCopyReal8ToPtr(id%RHS, x_pr, mA)

!Deallocate
deallocate(IRN)
deallocate(JCN)
deallocate(A)
deallocate(RHS)

endif

!Destroy instance
id%JOB = -2
call dmumps(id)

call mpi_finalize(ierr)

return
end

***********************
Matlab Crash Report
***********************
------------------------------------------------------------------------
       Segmentation violation detected
------------------------------------------------------------------------

Configuration:
  Crash Decoding  : Disabled
  Current Visual  : None
  Default Encoding: UTF-8
  GNU C Library   : 2.12 stable
  MATLAB Root     : /opt/matlab/R2012b
  MATLAB Version  : 8.0.0.783 (R2012b)
  Operating System: Linux 2.6.32-279.5.1.el6.x86_64 #1 SMP Tue Aug 14 16:11:42 CDT 2012 x86_64
  Processor ID    : x86 Family 6 Model 44 Stepping 2, GenuineIntel
  Virtual Machine : Java 1.6.0_17-b04 with Sun Microsystems Inc. Java HotSpot(TM) 64-Bit Server VM mixed mode
  Window System   : No active display

Fault Count: 1

Abnormal termination:
Segmentation violation

Register State (from fault):
  RAX = 0000000000000001  RBX = 00007fbac43cd120
  RCX = 0000000000000000  RDX = 0000000100000000
  RSP = 00007fbacec70330  RBP = 54534f4300000001
  RSI = 00007fbac43cd120  RDI = 0000000080000000

     R8 = a29af9d2c43cd128   R9 = 4534f43000000010
    R10 = 453573eac43cd130  R11 = 0000000000000000
    R12 = 00007fbacec70400  R13 = 00007fbac43cd120
    R14 = 00007fba5bdd29f4  R15 = 0000000000000001

    RIP = 00007fba5996b390  EFL = 0000000000010a86

     CS = 0033   FS = 0000   GS = 0000

Stack Trace (from fault):
[  0] 0x00007fbae668c1de          /opt/matlab/R2012b/bin/glnxa64/libmwfl.so+00516574 _ZN2fl4diag15stacktrace_base7captureERKNS0_14thread_contextEm+000158
[  1] 0x00007fbae668d4b2          /opt/matlab/R2012b/bin/glnxa64/libmwfl.so+00521394
[  2] 0x00007fbae668effe          /opt/matlab/R2012b/bin/glnxa64/libmwfl.so+00528382 _ZN2fl4diag13terminate_logEPKcRKNS0_14thread_contextE+000174
[  3] 0x00007fbae597a093         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00557203 _ZN2fl4diag13terminate_logEPKcPK8ucontext+000067
[  4] 0x00007fbae5976b9d         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00543645
[  5] 0x00007fbae5978835         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00550965
[  6] 0x00007fbae5978a55         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00551509
[  7] 0x00007fbae59790fe         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00553214
[  8] 0x00007fbae5979295         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00553621
[  9] 0x0000003aa2c0f500                             /lib64/libpthread.so.0+00062720
[ 10] 0x00007fba5996b390              /opt/matlab/R2012b/bin/glnxa64/mkl.so+22602640 mkl_blas_mc3_idamax+000160
[ 11] 0x00007fba585f42a6              /opt/matlab/R2012b/bin/glnxa64/mkl.so+02192038 mkl_blas_idamax+000038
[ 12] 0x00007fba584df996              /opt/matlab/R2012b/bin/glnxa64/mkl.so+01059222 IDAMAX_+000006
[ 13] 0x00007fba5bc623b8          /home/ase5/Fort_Mex_Par_try/check1.mexa64+00660408
[ 14] 0x00007fba5bc32c6e          /home/ase5/Fort_Mex_Par_try/check1.mexa64+00466030
[ 15] 0x00007fba5bbd0ded          /home/ase5/Fort_Mex_Par_try/check1.mexa64+00065005
[ 16] 0x00007fba5bc50ce2          /home/ase5/Fort_Mex_Par_try/check1.mexa64+00589026
[ 17] 0x00007fba5bc1864b          /home/ase5/Fort_Mex_Par_try/check1.mexa64+00357963
[ 18] 0x00007fba5bc47eed          /home/ase5/Fort_Mex_Par_try/check1.mexa64+00552685
[ 19] 0x00007fba5bbe2c44          /home/ase5/Fort_Mex_Par_try/check1.mexa64+00138308
[ 20] 0x00007fba5bbc7c5a          /home/ase5/Fort_Mex_Par_try/check1.mexa64+00027738 mexfunction_+001182
[ 21] 0x00007fbadd8516ac           /opt/matlab/R2012b/bin/glnxa64/libmex.so+00112300 mexRunMexFile+000108
[ 22] 0x00007fbadd84d4e9           /opt/matlab/R2012b/bin/glnxa64/libmex.so+00095465
[ 23] 0x00007fbadd84e33c           /opt/matlab/R2012b/bin/glnxa64/libmex.so+00099132
[ 24] 0x00007fbae56c7a4b /opt/matlab/R2012b/bin/glnxa64/libmwm_dispatcher.so+00596555 _ZN8Mfh_file11dispatch_fhEiPP11mxArray_tagiS2_+000539
[ 25] 0x00007fbae4f5be56 /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+02264662
[ 26] 0x00007fbae4ee7c1d /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+01788957
[ 27] 0x00007fbae4f1024e /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+01954382
[ 28] 0x00007fbae4f0d0d3 /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+01941715
[ 29] 0x00007fbae4f0ded7 /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+01945303
[ 30] 0x00007fbae4f79760 /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+02385760
[ 31] 0x00007fbae56c7a4b /opt/matlab/R2012b/bin/glnxa64/libmwm_dispatcher.so+00596555 _ZN8Mfh_file11dispatch_fhEiPP11mxArray_tagiS2_+000539
[ 32] 0x00007fbae4f4897b /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+02185595
[ 33] 0x00007fbae4f0621c /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+01913372
[ 34] 0x00007fbae4f0324d /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+01901133
[ 35] 0x00007fbae4f03685 /opt/matlab/R2012b/bin/glnxa64/libmwm_interpreter.so+01902213
[ 36] 0x00007fbadda7c22e      /opt/matlab/R2012b/bin/glnxa64/libmwbridge.so+00143918
[ 37] 0x00007fbadda7c391      /opt/matlab/R2012b/bin/glnxa64/libmwbridge.so+00144273
[ 38] 0x00007fbadda7cf6d      /opt/matlab/R2012b/bin/glnxa64/libmwbridge.so+00147309 _Z8mnParserv+000733
[ 39] 0x00007fbae595f472         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00447602 _ZN11mcrInstance30mnParser_on_interpreter_threadEv+000034
[ 40] 0x00007fbae593db69         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00310121
[ 41] 0x00007fbae593dd48         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00310600
[ 42] 0x00007fbae5f53f73    /opt/matlab/R2012b/bin/glnxa64/libmwservices.so+00999283 _ZN10eventqueue18UserEventQueueImpl5flushEv+000371
[ 43] 0x00007fbae5f54695    /opt/matlab/R2012b/bin/glnxa64/libmwservices.so+01001109 _ZN10eventqueue8ReadPipeEib+000053
[ 44] 0x00007fbae5f53321    /opt/matlab/R2012b/bin/glnxa64/libmwservices.so+00996129 _ZN10eventqueue18UserEventQueueImpl9selectFcnEb+000353
[ 45] 0x00007fbadaaeba65         /opt/matlab/R2012b/bin/glnxa64/libmwuix.so+00518757
[ 46] 0x00007fbae5feda11    /opt/matlab/R2012b/bin/glnxa64/libmwservices.so+01628689 _ZSt8for_eachIN9__gnu_cxx17__normal_iteratorIPN5boost8weak_ptrIN4sysq10ws_ppeHookEEESt6vectorIS6_SaIS6_EEEENS4_8during_FIS6_NS2_10shared_ptrIS5_EEEEET0_T_SH_SG_+000081
[ 47] 0x00007fbae5feeaeb    /opt/matlab/R2012b/bin/glnxa64/libmwservices.so+01633003 _ZN4sysq12ppe_for_eachINS_8during_FIN5boost8weak_ptrINS_10ws_ppeHookEEENS2_10shared_ptrIS4_EEEEEET_RKS9_+000251
[ 48] 0x00007fbae5fec5a2    /opt/matlab/R2012b/bin/glnxa64/libmwservices.so+01623458 _ZN4sysq19ppePollingDuringFcnEb+000114
[ 49] 0x00007fbae5fec969    /opt/matlab/R2012b/bin/glnxa64/libmwservices.so+01624425 _ZN4sysq11ppeMainLoopEiib+000121
[ 50] 0x00007fbae5fecb08    /opt/matlab/R2012b/bin/glnxa64/libmwservices.so+01624840 _ZN4sysq11ppeLoopIfOKEiib+000152
[ 51] 0x00007fbae5fecc63    /opt/matlab/R2012b/bin/glnxa64/libmwservices.so+01625187 _ZN4sysq20processPendingEventsEiib+000147
[ 52] 0x00007fbae593e664         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00312932
[ 53] 0x00007fbae593eb3c         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00314172
[ 54] 0x00007fbae5938592         /opt/matlab/R2012b/bin/glnxa64/libmwmcr.so+00288146
[ 55] 0x0000003aa2c07851                             /lib64/libpthread.so.0+00030801
[ 56] 0x0000003aa24e811d                                   /lib64/libc.so.6+00950557 clone+000109

This error was detected while a MEX-file was running. If the MEX-file
is not an official MathWorks function, please examine its source code
for errors. Please consult the External Interfaces Guide for information
on debugging MEX-files.

Answer 1

[Korbi]

Hi Anand,

For the case that you still have not found an answer for your question: Have you tried the PARDISO package? I still have problems with the compatibility with the Linux version on our cluster to provide more experience, but it could be an interesting alternative for you. It also provides a Matlab interface and seems to perform quite well in benchmark tests, both in memory usage and calculation speed.

I basically have the same problem as you. I implemented a simple CFD-Code in Matlab and I want to run it on a Linux cluster with manifold input parameter variations. Acutally it's nothing else then the permanent solution of SLEs. In my case they have the size of about 200.000 rows/columns. For the present memory usage is not an issue. The performance gain through parallelized solution of sparse SLEs in Matlab is what makes me ponder.