[Ipopt-tickets] [Ipopt] #216: interface to MKL PARDISO
Ipopt
coin-trac at coin-or.org
Fri Aug 23 13:02:07 EDT 2013
#216: interface to MKL PARDISO
-----------------------+------------------------
Reporter: akalinki | Owner: ipopt-team
Type: defect | Status: new
Priority: normal | Component: Ipopt
Version: 3.11 | Severity: normal
Resolution: | Keywords:
-----------------------+------------------------
Comment (by stefan):
In the following some E-Mail correspondence from the last days:
Alexander,
Fantastic, thanks for following up in so much detail through so many back-
and-forth emails. One thing to note first, I can log in to the Ipopt Trac
page again, so perhaps we might want to post the contents of this
discussion up there (or on the Ipopt mailing list) so other users of Ipopt
can reference this information in the future.
I can send you my test problems (and my interface code you need to compile
in addition to Ipopt to run them) to look at, they are nonlinear
optimization problems from MPC on HVAC systems, with a few different model
variations - Pardiso gets called at every iteration of Ipopt (possibly
multiple times on iterations where perturbation/regularization are
required to get the right inertia for descent properties) on the symmetric
indefinite KKT matrix (refer to main Ipopt implementation paper
http://dx.doi.org/10.1007/s10107-004-0559-y). If the Ipopt interface to
Pardiso has a dump-matrix option we could save the exact matrices being
sent to Pardiso at each iteration for testing.
Thanks for the info on the MKL_DOMAIN_* variables, I didn't know about
those. Makes sense that Pardiso just uses sequential BLAS/LAPACK for small
dense sub-blocks.
I was thinking something similar regarding checking for mkl_get_version to
identify if an Ipopt user is building Ipopt linked to MKL. I see a
function mkl_get_version_string in the libraries, but which library is
mkl_get_version located in? Stefan, Jonathan, and I were discussing how we
could incorporate this check into the build system, it may end up being
straightforward.
It seems like with Jonathan's simple modifications to the Pardiso options,
we're mostly there. We should look through them in a little more detail to
be sure we set appropriate defaults for use with Ipopt, and possibly
advertise to the Ipopt mailing list to find a few more testers.
-Tony
-----Original Message----- From: Kalinkin, Alexander A
Sent: Thursday, August 22, 2013 8:39 PM
To: Stefan Vigerske ; jonathan.currie at aut.ac.nz
Cc: 'Tony Kelman' ; David I Wilson ; Andreas Waechter
Subject: RE: [Ipopt] Ticket 216 - interface to MKL PARDISO
Hi,
I’ve tried to combine all questions/answers below:
Jonathan: “- Pardiso 4.1.2 is slightly faster on large problems (4e6
elements in the Jacobian) than MKL Pardiso, although MKL Pardiso does not
seem to be using all cores currently (4.1.2 is). Difference of 0.3s in 9s
runs. Numerical results are very similar.”
[akalinki] Is it comparison between serial or parallel version of
implementations? If parallel can I ask number of threads?
Jonathan: “As far as documentation goes, it does appear as though the
missing functionality (inertia, number of negative pivots) is now
implemented in MKL Pardiso. Not sure on the matching stuff although I am
finding similar terms between the paper:
"weighted symmetric matchings, and the static factorization with
Supernodal-Bunch–Kaufman pivoting"
and the MKL Pardiso documentation:
"Intel MKL PARDISO can use a maximum weighted matching algorithm to
permute large elements close the diagonal", "the solver uses a complete
supernode pivoting approach", "Apply 1x1 and 2x2 Bunch-Kaufman pivoting
during the factorization process"
I guess only Intel can tell us what functionality has been implemented
from later versions of Pardiso?”
[akalink] I’ve investigated paper mentioned in this topic (http://list
.coin-or.org/pipermail/ipopt/2009-June/001583.html) and didn’t find any
significant difference between MKL PARDISO implementation of matching.
It’s clear that algorithm explained only on high level, but MKL similar
implementation of matching.
Tony: “but this is consistent with the results I got last year on the same
problems comparing to non-MKL Pardiso (HSL and WSMP solvers had noticeably
better performance).”
[akalinki] We doesn’t expect performance worse than this 2 solvers… Can I
asked about matrix and way of using pardiso to understand lower
performance of MKL PARDISO implementation?
Stefan: “The best may be if Alexander just tries how the current
Ipopt/PARDISO interface works with the MKL PARDISO.”
[akalinki] will check it.
Tony: “ 1. Regarding threading, MKL Pardiso presumably uses MKL BLAS
and/or LAPACK subroutines, right? How is nested parallelism handled
between the Pardiso routines and any lower-level BLAS/LAPACK subroutines
it calls? There's only one MKL_NUM_THREADS environment variable right, it
doesn't get any finer-grained to control BLAS/LAPACK number of threads
independently of Pardiso number of threads? My experience with Ipopt has
shown it's better to allocate all threads to the linear solver rather than
BLAS/LAPACK, since the dense submatrices that occur during sparse
factorizations are often too small for multithreading to be worth the
synchronization overhead or reducing the number of outer-level threads for
Pardiso to use. Though the maximum size of problems I've looked at in
Ipopt is only a few tens of thousands of variables, so I'm not sure if my
conclusions hold in all cases.”
[akalinki] In MKL PARDISO functionality sequential LAPACK and BLAS kernels
calls from PARDISO omp region that allow us to achieve performance on
small size of LAPACK and BLAS sizes. This approach is a bit different with
mumps for example which use shared parallelism from LAPACK and Blas
functionality. You can use different threading approach for different part
by MKL using environment variable by following line (for example):
MKL_DOMAIN_ALL 2: MKL_DOMAIN_BLAS 4, MKL_DOMAIN_LAPACK 4
but such it doesn’t the point for blas and lapack kernels in pardiso
functionality – they call sequentially.
Tony: “2. Is there a programmatic way for Ipopt's build system to know
whether or not the version of Pardiso that it is linking to and calling
comes from MKL or not? For example, currently the Ipopt build system
checks whether calling
pardiso_exist_parallel() or pardiso_ipopt_newinterface() work to determine
the appropriate defines to use. If we need to make source-level changes to
Ipopt or its build system to get MKL Pardiso to work, it would be ideal if
the build system could automatically determine whether or not those
changes need to be applied (via defines or some other method).”
[akalinki] Maybe I didn’t understood your question – from my point of view
mkl_get_version() can be used here. This function return version of
current MKL library. Also I want to underline that we didn’t make change
in our PARDISO interface from its implementation a years ago because such
changes could affect functionality of our customer.
Thanks,
Alexander Kalinkin
-----Original Message-----
From: Stefan Vigerske [mailto:stefan at math.hu-berlin.de]
Sent: Friday, August 23, 2013 1:01 AM
To: jonathan.currie at aut.ac.nz
Cc: 'Tony Kelman'; Kalinkin, Alexander A; David I Wilson; Andreas Waechter
Subject: Re: [Ipopt] Ticket 216 - interface to MKL PARDISO
Hi,
so I guess what you guys say is that the MKL PARDISO seems usable so far
and we should integrate Jonathan's patch and define MKL_PARDISO if
configure figures out that PARDISO is from MKL ?
We could also automatically activate the PARDISO interface, if MKL
Blas/Lapack is used.
Stefan
On 08/22/2013 12:05 PM, Jonathan Currie wrote:
> Quick update comparing MKL Pardiso vs Pardiso 4.1.2 (official):
>
> - Pardiso 4.1.2 is slightly faster on large problems (4e6 elements in
the Jacobian) than MKL Pardiso, although MKL Pardiso does not seem to be
using all cores currently (4.1.2 is). Difference of 0.3s in 9s runs.
Numerical results are very similar.
> - Both are 2-3 faster than MA57 on large problems (I guess that’s
> parallelism for you)
>
> I haven't been able to compile with "PARDISO_MATCHING_PREPROCESS"
defined as I can't find "smat_reordering_pardiso_wsmp_". Anyone know where
that lives or its function?
>
> Also, I can't reproduce the failing behaviour of MKL Pardiso (or 4.1.2)
with IPOPT in MATLAB that I saw in BLOM. I suspect I had something wrong
in compilation of BLOM?
>
>
> As far as documentation goes, it does appear as though the missing
functionality (inertia, number of negative pivots) is now implemented in
MKL Pardiso. Not sure on the matching stuff although I am finding similar
terms between the paper:
> "weighted symmetric matchings, and the static factorization with
Supernodal-Bunch–Kaufman pivoting"
>
> and the MKL Pardiso documentation:
> "Intel MKL PARDISO can use a maximum weighted matching algorithm to
> permute large elements close the diagonal", "the solver uses a
> complete supernode pivoting approach", "Apply 1x1 and 2x2 Bunch-Kaufman
pivoting during the factorization process"
>
> I guess only Intel can tell us what functionality has been implemented
from later versions of Pardiso?
>
> Jonathan
>
> -----Original Message-----
> From: Tony Kelman [mailto:kelman at berkeley.edu]
> Sent: Thursday, 22 August 2013 2:30 p.m.
> To: jonathan.currie at aut.ac.nz; 'Stefan Vigerske'
> Cc: 'Kalinkin, Alexander A'
> Subject: Re: [Ipopt] Ticket 216 - interface to MKL PARDISO
>
> I'm actually seeing different results than Jonathan on Linux. Enabling
his initial set of options with -DMKL_PARDISO (though his modification
will need some revision to avoid changing behavior with the non-MKL
Pardiso) things actually appeared to work and converge reasonably well. On
the set of optimization problems I'm looking at there is not a convincing
performance difference between Mumps and MKL Pardiso, but this is
consistent with the results I got last year on the same problems comparing
to non-MKL Pardiso (HSL and WSMP solvers had noticeably better
performance).
>
> We might just need some different test problems that are numerically
challenging for non-MKL Pardiso, or problems where Pardiso has previously
shown much better performance than other available solvers, to see how MKL
Pardiso handles them.
>
> -Tony
>
> -----Original Message-----
> From: Jonathan Currie
> Sent: Wednesday, August 21, 2013 4:43 AM
> To: 'Stefan Vigerske' ; 'Tony Kelman'
> Cc: 'Kalinkin, Alexander A' ; 'Andreas Waechter' ; David I Wilson
> Subject: RE: [Ipopt] Ticket 216 - interface to MKL PARDISO
>
> Hmmm initial testing on small problems with PARDISO was OK, however
larger problems using Tony's BLOM models basically went straight to
restoration (from about iteration 2) then failure. Not quite sure what is
going on.
>
> I've attached my quick and dirty modification to the Ipopt PARDISO
interface, just define "MKL_PARDISO" to swop the settings over, and
obviously "HAVE_PARDISO" and "HAVE_PARDISO_OLDINTERFACE" to use it. I
followed the Intel recommendations on settings, but I know a lot more time
should be spent here (e.g. double parameters not supported, thus many
IPOPT settings there are redundant). Also not sure whether I should be
using "PARDISO_MATCHING_PREPROCESS" or not?
>
> Anyway if you want to have a go I've uploaded mexw32 and mexw64 ipopt
builds, with MUMPS and PARDISO (MKL v11 R5) statically linked, and MA57
dynamically linked to MATLAB's version. MA57 is chosen by default, so as
per normal, change the linear_solver option. Mex files available below:
> http://www.i2c2.aut.ac.nz/Downloads/Temp/IPOPT_PARDISO.zip
>
> As with my builds you will need VC++ 2012 to run them. I don't have
access to the non-MKL version of PARDISO so I can't really compare
performance.
>
> Jonathan
>
> p.s. Tony Sparse++ was a mission to get going! Had to hack quite a
number of the source files before it would compile, and good old VC++ 2012
doesn't support the new C++11 features you were using so they had to come
out too!
> Got BLOM compiled after an hour so... Works fine with mumps, virtually
no blom problems would work with pardiso.
>
>
> -----Original Message-----
> From: Stefan Vigerske [mailto:stefan at math.hu-berlin.de]
> Sent: Wednesday, 21 August 2013 9:41 p.m.
> To: Tony Kelman
> Cc: jonathan.currie at aut.ac.nz; 'Kalinkin, Alexander A'; Andreas
> Waechter
> Subject: Re: [Ipopt] Ticket 216 - interface to MKL PARDISO
>
> Hi,
>
> Toni, Jonathan: Thanks for picking this up.
> I don't have much to add.
>
> According to Andreas, the current PARDISO interface in Ipopt was mainly
written by Olaf Schenk, so he should be able to comment on it, too.
>
> In the past, the inertia were the main problem with the MKL PARDISO.
> The best may be if Alexander just tries how the current Ipopt/PARDISO
interface works with the MKL PARDISO.
>
> Best,
> Stefan
>
> On 08/21/2013 11:05 AM, Tony Kelman wrote:
>> Jonathan,
>>
>> That didn't take long, and figured you'd also be interested in this.
>> Though your build environment with Visual Studio is unlike the other
>> 3 or 4 common environments that use autotools-based builds (Linux,
>> OSX, MSYS/MinGW-w64, and MSYS/MinGW-w32, all potentially multiplied
>> by GNU vs Intel compiler choice). I only have MKL on Linux and I'm
>> pretty sure my non-MKL Pardiso license is expired.
>>
>> Besides the difference in managing number of threads (where it makes
>> sense that MKL would handle this differently), any of the options
>> that have currently set values in Ipopt that are either non-default
>> or not implemented in MKL Pardiso should be very carefully looked at.
>> Not just by someone familiar with MKL Pardiso, but also by someone
>> who has more experience than I do with using the non-MKL Pardiso with
>> Ipopt and how these Pardiso options influence Ipopt's performance and
>> reliability in solving optimization problems.
>>
>> Regarding patches, if you're on Windows with TortoiseSVN installed
>> and you're working with a copy of Ipopt checked out from SVN, then
>> it's quite easy. Make the source changes, then just right-click on
>> the top-level Ipopt directory and go to TortoiseSVN->Create patch...
>> We can of course just do diffs manually and compare to see what
>> modifications you make, but with patches the files you send around
>> are smaller and can be somewhat easier to read.
>>
>>
>> Questions for Alexander:
>>
>> 1. Regarding threading, MKL Pardiso presumably uses MKL BLAS and/or
>> LAPACK subroutines, right? How is nested parallelism handled between
>> the Pardiso routines and any lower-level BLAS/LAPACK subroutines it
calls?
>> There's only one MKL_NUM_THREADS environment variable right, it
>> doesn't get any finer-grained to control BLAS/LAPACK number of
>> threads independently of Pardiso number of threads? My experience
>> with Ipopt has shown it's better to allocate all threads to the
>> linear solver rather than BLAS/LAPACK, since the dense submatrices
>> that occur during sparse factorizations are often too small for
>> multithreading to be worth the synchronization overhead or reducing
>> the number of outer-level threads for Pardiso to use. Though the
>> maximum size of problems I've looked at in Ipopt is only a few tens
>> of thousands of variables, so I'm not sure if my conclusions hold in
all cases.
>>
>> 2. Is there a programmatic way for Ipopt's build system to know
>> whether or not the version of Pardiso that it is linking to and
>> calling comes from MKL or not? For example, currently the Ipopt build
>> system checks whether calling pardiso_exist_parallel() or
>> pardiso_ipopt_newinterface() work to determine the appropriate
>> defines to use. If we need to make source-level changes to Ipopt or
>> its build system to get MKL Pardiso to work, it would be ideal if the
>> build system could automatically determine whether or not those
>> changes need to be applied (via defines or some other method).
>>
>> Thanks,
>> Tony
>>
>>
>> -----Original Message----- From: Jonathan Currie
>> Sent: Wednesday, August 21, 2013 1:00 AM
>> To: 'Tony Kelman' ; 'Kalinkin, Alexander A'
>> Subject: RE: [Ipopt] Ticket 216 - interface to MKL PARDISO
>>
>> Quick update, got it working on the HS71 example (no idea re inertia
>> stuff for now...). Few things different:
>>
>> - Most of the options setup in IpPardisoSolverInterface between lines
>> 386 and 417 are wrong or different from the Intel specified ones. Can
>> basically comment them all except IPARM_[5] = 1.
>> - All options are listed here:
>> http://software.intel.com/sites/products/documentation/doclib/mkl_sa/
>> 1
>> 1/mklman/index.htm
>>
>> - You will need to define HAVE_PARDISO and HAVE_PARDISO_OLDINTERFACE.
>> I don't think HAVE_PARDISO_PARALLEL has any effect on the solver as
>> the threading is controlled via MKL_NUM_THREADS. Best to check though.
>>
>> It will definitely pay to go through the available Intel MKL PARDISO
>> options carefully with someone who knows the pros and cons of them all.
>> I will upload a couple of mex builds shortly you can have a play with
>> and a copy of some slightly modified source (if all working in Matlab
>> OK). Still no expert on this patch system!
>>
>> Jonathan
>>
>> -----Original Message-----
>> From: ipopt-bounces at list.coin-or.org
>> [mailto:ipopt-bounces at list.coin-or.org] On Behalf Of Tony Kelman
>> Sent: Wednesday, 21 August 2013 4:52 p.m.
>> To: Kalinkin, Alexander A
>> Cc: ipopt at list.coin-or.org
>> Subject: Re: [Ipopt] Ticket 216 - interface to MKL PARDISO
>>
>> From the discussion in 2011 it appears inertia functionality should
>> be fine, the matching-based preprocessing is a separate issue.
>>
>> You can see exactly how Ipopt is currently calling Pardiso here, to
>> check if anything looks like it would not work with the MKL version
>> at a source
>> level:
>> https://projects.coin-or.org/Ipopt/browser/releases/3.11.3/Ipopt/src/
>> A lgorithm/LinearSolvers/IpPardisoSolverInterface.cpp
>>
>> and how Ipopt's configure handles detection of Pardiso and building
>> Ipopt's interface to it here:
>> https://projects.coin-or.org/Ipopt/browser/releases/3.11.3/Ipopt/conf
>> i
>> gure.ac#L278
>>
>>
>> Although curiously, the line in IpPardisoSolverInterface.hpp that
>> says #define PARDISO_MATCHING_PREPROCESS has been commented out since
>> it was added in SVN revision 1791, so it would seem all the code
>> under an #ifdef PARDISO_MATCHING_PREPROCESS is not currently being
>> used? I wonder if that might explain the not-fantastic performance of
>> the official Pardiso the last time I tried to run it on one of my
problems.
>>
>> If I have an extra minute I might check what version of MKL I have
>> installed and see if anything obvious breaks when trying to build/run
>> Ipopt with MKL Pardiso currently.
>>
>> -Tony
>>
>>
>> -----Original Message-----
>> From: Kalinkin, Alexander A
>> Sent: Tuesday, August 20, 2013 9:22 PM
>> To: Tony Kelman
>> Cc: ipopt at list.coin-or.org
>> Subject: RE: Ticket 216 - interface to MKL PARDISO
>>
>> Hi Tony,
>> I see. I will check current implementation of inertia functionality
>> to provide deeper answer on this question. Nevertheless, is inertia
>> is only one barrier or there is another requests?
>> Thanks,
>> Alexander
>>
>> -----Original Message-----
>> From: Tony Kelman [mailto:kelman at berkeley.edu]
>> Sent: Wednesday, August 21, 2013 10:52 AM
>> To: Kalinkin, Alexander A
>> Cc: ipopt at list.coin-or.org
>> Subject: Ticket 216 - interface to MKL PARDISO
>>
>> Alexander,
>>
>> I would’ve responded to this on the ticket page you opened
>> (https://projects.coin-or.org/Ipopt/ticket/216), but it appears there
>> was an update to Trac today and I can’t seem to log in.
>>
>> You’ve probably done the searching as well regarding past discussions
>> on this topic, but I’ll include links here anyway to a few that I was
>> able to find for reference's sake:
>>
>> https://projects.coin-or.org/Ipopt/ticket/88
>> http://list.coin-or.org/pipermail/ipopt/2009-June/001583.html
>> http://list.coin-or.org/pipermail/ipopt/2008-January/000990.html
>> http://list.coin-or.org/pipermail/ipopt/2011-September/002567.html
>>
>> In the 2008-2009 time frame of the first three links, the biggest
>> issue was that MKL's PARDISO implementation did not calculate and/or
>> output the inertia of the factorized matrix. The fourth link from
>> 2011 indicates the inertia was added in a relatively recent version of
MKL.
>> Olaf Schenk provided a reference there to a paper "Matching-based
>> Preprocessing Algorithms to the Solution of Saddle-Point Problems in
>> Large-Scale Nonconvex Interior-Point Optimization" that describes the
>> remaining feature present in the primary implementation of PARDISO
>> that is not yet available (or was not at the time) in the MKL version.
>> Can you comment on whether the pivoting algorithm described in that
>> paper is or could be implemented in Intel MKL's version of PARDISO?
>>
>> Glad you posted, this could be a quite helpful alternative to provide
>> users of Ipopt.
>>
>> -Tony
--
Ticket URL: <https://projects.coin-or.org/Ipopt/ticket/216#comment:1>
Ipopt <http://projects.coin-or.org/Ipopt>
Interior-point optimizer for nonlinear programs.
More information about the Ipopt-tickets
mailing list