coin-or / Ipopt

COIN-OR Interior Point Optimizer IPOPT
https://coin-or.github.io/Ipopt
Other
1.4k stars 279 forks source link

interface to MKL PARDISO #216

Closed svigerske closed 4 years ago

svigerske commented 5 years ago

Issue created by migration from Trac.

Original creator: akalinki

Original creation time: 2013-08-20 07:59:12

Assignee: ipopt-team

Version: 3.11

CC: jonnyz007 @tkelman

Hello, I am engineer in Intel MKL, responsible for developing Sparse component of Intel MKL. We are interested in support interface of MKL PARDISO by ipopt project but don't know how we can implement interface to our functionality in your code. Is there any way to collaborate with you in such activity? Thanks, Alexander Kalinkin

svigerske commented 5 years ago

Comment by @svigerske created at 2013-08-23 17:02:07

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 MKLDOMAIN* 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@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@math.hu-berlin.de] Sent: Friday, August 23, 2013 1:01 AM To: jonathan.currie@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_pardisowsmp". 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@berkeley.edu] Sent: Thursday, 22 August 2013 2:30 p.m. To: jonathan.currie@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@math.hu-berlin.de] Sent: Wednesday, 21 August 2013 9:41 p.m. To: Tony Kelman Cc: jonathan.currie@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@list.coin-or.org [mailto:ipopt-bounces@list.coin-or.org] On Behalf Of Tony Kelman Sent: Wednesday, 21 August 2013 4:52 p.m. To: Kalinkin, Alexander A Cc: ipopt@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@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@berkeley.edu] Sent: Wednesday, August 21, 2013 10:52 AM To: Kalinkin, Alexander A Cc: ipopt@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

svigerske commented 5 years ago

Attachment IpPardisoSolverInterface.cpp by @svigerske created at 2013-08-23 17:03:36

patched Ipopt PARDISO interface by J. Currie to work with MKL PARDISO

svigerske commented 5 years ago

Comment by @svigerske created at 2013-08-23 17:13:49

Replying to [comment:1 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 think this ticket would be a good place to continue the discussion. It's probably too technical and mail-intensive for the mailing list.

When I find some time, I can look into extending the build system and incorporating the patch from Jonathan.

svigerske commented 5 years ago

Comment by @svigerske created at 2013-08-23 17:14:51

Changing assignee from ipopt-team to @svigerske.

svigerske commented 5 years ago

Comment by @svigerske created at 2013-08-23 17:14:51

Changing status from new to assigned.

svigerske commented 5 years ago

Comment by @tkelman created at 2013-08-31 22:45:03

I took a crack at the autotools part of this. I renamed the define Jonathan added to HAVE_PARDISO_MKL, so it would end up next to the other Pardiso-related defines in config.h.in (which is created in alphabetical order by autoheader).

The additions to configure.ac are pretty short. I check for MKL based on whether the function mkl_get_version_string is available in any of the BLAS_LIBS, LAPACK_LIBS, or Pardiso libs (if any are specified by the user). If MKL is available, then we should always try to use Pardiso (unless we think there are still folks using old enough versions of MKL that don't have Pardiso available, in which case we'd have to require anyone who wants to use MKL-Pardiso to specify --with-pardiso).

The HAVE_PARDISO_MKL define should only get set when we've determined we aren't using Basel Pardiso, which would have pardiso_ipopt_newinterface available.

I don't have access to MKL to test this anywhere except Linux - can anyone test on Mac and/or Windows via MinGW/MSYS? You'll need to run BuildTools/run_autotools to propagate the changes from configure.ac to configure and config.h.in.

svigerske commented 5 years ago

Attachment mkl_pardiso.patch by @tkelman created at 2013-09-01 14:25:32

Add new define to config_default.h, and catch a typo about the linear solver loader

svigerske commented 5 years ago

Comment by @svigerske created at 2013-09-01 16:35:26

I applied this patch to trunk with 6d0b2e67ce.

On my Linux, with a fairly old MKL, it works like a charm.

I tried MKL PARDISO on one of the scalabe example and set

linear_solver pardiso
pardiso_msglvl 1

This gave me

 < Parallel Direct Factorization with #processors: >         2

in the statistics output of PARDISO.

So it's maybe confusing that configure says that it is not the parallel version.

However, enabling HAVE_PARDISO_PARALLEL would let Ipopt raise an error if a user does not set OMP_NUM_THREADS. I see that with the patch, one lets MKL PARDISO choose the number of threads as it likes, which seems to be the number of processors on the machine.

Would it make sense to have the interface behave the same with any PARDISO library? That is, always set IPARM_4b93956e64 to the value of OMP_NUM_THREADS, if set (yes, that seems redundant, at least in case of MKL PARDISO), and leave it unset if the variable has not been set, but also not reporting an error then. It seems odd to me that one is forced to set the OMP_NUM_THREADS variable for Basel PARDISO.

svigerske commented 5 years ago

Comment by @svigerske created at 2013-09-01 17:13:22

Replying to [comment:4 kelman]:

I don't have access to MKL to test this anywhere except Linux - can anyone test on Mac and/or Windows via MinGW/MSYS?

I tried

I don't have an Intel MKL for Mac so far.

svigerske commented 5 years ago

Comment by @tkelman created at 2013-09-01 17:17:16

I guess the configure message about HAVE_PARDISO_PARALLEL is a bit misleading, since MKL Pardiso doesn't have pardiso_exist_parallel() defined but it does multi-threading with its own system.

My understanding is MKL lets you use either OMP_NUM_THREADS or MKL_NUM_THREADS, or uses all processors if neither environment variable is set.

The MKL documentation says: "iparm(3) Reserved. Set to zero." We can test setting it to the value of OMP_NUM_THREADS and see if anything bad happens?

Maybe a warning instead of an error if OMP_NUM_THREADS isn't set when using Basel Pardiso? The code is setting num_procs = 1 by default anyway.

svigerske commented 5 years ago

Comment by @svigerske created at 2013-09-01 17:56:02

OK, just not setting it seems fine to me.

cb02e2ac87 disables the check for pardiso_exist_parallel() if PARDISO is from MKL and changes the error to a warning.

svigerske commented 5 years ago

Attachment parallel_pardiso_check.patch by @tkelman created at 2013-09-01 19:03:31

svigerske commented 5 years ago

Comment by @tkelman created at 2013-09-01 19:03:50

Replying to [comment:8 stefan]:

OK, just not setting it seems fine to me.

cb02e2ac87 disables the check for pardiso_exist_parallel() if PARDISO is from MKL and changes the error to a warning.

I think the logic in configure from cb02e2ac87 is a tiny bit off... have_mkl = yes means we're at least using MKL for Blas, but might still be using MKL Blas with Basel Pardiso. If it's possible that any non-parallel version of Basel Pardiso is still floating around, then we would still want to check for pardiso_exist_parallel() in that case. Although a false positive for HAVE_PARDISO_PARALLEL looks pretty harmless so I doubt it matters.

svigerske commented 5 years ago

Comment by @svigerske created at 2013-09-01 19:38:57

Well, in the previous version, we were not recognizing a multithreaded PARDISO if it was the PARDISO MKL.

Also now, if one is using MKL and has a PARDISO without pardiso_ipopt_newinterface(), the system assumes that PARDISO comes from MKL, even though it could still be the Basel MKL.

I'm fine with assuming that one either has a recent Basel PARDISO (new interface, parallelized) or an MKL PARDISO (or none at all).

svigerske commented 5 years ago

Comment by @tkelman created at 2013-09-02 03:01:58

Since MKL has its own threading control system, HAVE_PARDISO_PARALLEL wasn't necessary for it. I think your fix for that part is fine, my concern was MKL Blas triggering a false positive for HAVE_PARDISO_PARALLEL using an old non-parallel Basel Pardiso.

Good point that you could still have old-interface Basel Pardiso with MKL Blas, the pardiso_ipopt_newinterface() is not a perfect proxy for whether or not we're using MKL Pardiso. Though that combination may have possible namespace shadowing of the libraries, it might be possible that MKL Pardiso could get called there instead of Basel Pardiso when the interfaces are the same.

I don't know how many combinations of old Basel Pardiso versions are around, in terms of old/new interface and parallel/not. There will always be old versions of Ipopt still available, if any of these old Pardiso versions no longer works properly with the most up-to-date version of Ipopt.

svigerske commented 5 years ago

Comment by akalinki created at 2013-09-17 03:10:47

Hi Colleagues, Sorry for delay in answering but could i currently help you somehow? Thanks, Alexander Kalinkin

svigerske commented 5 years ago

Comment by @svigerske created at 2013-09-17 08:18:54

I think all that would be needed is some more testing and maybe tuning of PARDISO parameters.

The changes in the !Ipopt/Pardiso interface from Jonathan and Tony have been committed, they are also in the latest release. With Ipopt/trunk, Ipopt's configure should also recognize a MKL/Pardiso if MKL is specified (or automatically found) for Blas. (I disabled this for the stable branch, since it would have changed the default linear solver for people who build Ipopt with MKL but without HSL.)

svigerske commented 5 years ago

Comment by akalinki created at 2013-09-17 09:39:12

Hi Stefan, Good news for us! Are there any ways to download last release to check how its work with MKL functionality? Thanks, Alexander Kalinkin

svigerske commented 5 years ago

Comment by @svigerske created at 2013-09-17 09:45:47

The easiest would be to try it with the current trunk. You can get it via subversion:

  svn co https://projects.coin-or.org/svn/Ipopt/trunk Ipopt-trunk

If configure finds your MKL libs, then Pardiso should be used automatically.

If you prefer the latest release, then download http://www.coin-or.org/download/source/Ipopt/Ipopt-3.11.4.tgz For this one, you would have to specify the MKL libs for the --with-pardiso configure flag and add -DHAVE_PARDISO_MKL -DHAVE_PARDISO_PARALLEL to the compiler flags, e.g., I'm running configure via

configure \
--with-blas="-L${MKLDIR} -Wl,--start-group -lmkl_def -lmkl_lapack -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -lippcoreem64t -Wl,--end-group -liomp5 -lpthread" \
--with-pardiso="-L${MKLDIR} -Wl,--start-group -lmkl_def -lmkl_lapack -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -lippcoreem64t -Wl,--end-group -liomp5 -lpthread" \
ADD_CXXFLAGS="-DHAVE_PARDISO_MKL -DHAVE_PARDISO_PARALLEL"

(there may be easier ways)

svigerske commented 5 years ago

Comment by akalinki created at 2013-09-17 10:07:01

Thanks again. Am i right that documentation still doesn't update? For example I found following strings:

''The Parallel Sparse Direct Solver (PARDISO) (see http://www.pardiso-project.org). Note: The Pardiso version in Intel's MKL library does not yet support the features necessary for Ipopt.''

svigerske commented 5 years ago

Comment by @svigerske created at 2013-09-17 10:37:54

Yes, I forgot to update the docu.

I now updated the documentation source regarding Ipopt/trunk (87dbd9fe2a), but maybe I should keep the previous version online while we still test if things are really working fine.

svigerske commented 5 years ago

Comment by akalinki created at 2013-10-18 02:46:08

Hi Stefan, I've provided your linkline to our customer and diff of revisions but they are interested about official release as one archive with updated version of documentation. Can you provide some estimation date if it is possible? Thanks, Alex

svigerske commented 5 years ago

Comment by @svigerske created at 2013-10-19 18:44:34

The archive for the official 3.11.4 is there.

I've now (308fa7f485) merged the changes in the docu to the 3.11 branch and changed them according to what is written in comment:15.

I also updated the HTML docu accordingly, in particular the expert installation options.

An updated pdf docu will be available with the next release, which could come within the next weeks (I have one open point on my list).

svigerske commented 5 years ago

Comment by @svigerske created at 2013-11-24 16:19:14

With 3866cc73ec I committed a few parameter changes for MKL Pardiso. Before that, I had some NLPs where PARDISO exited with error -4 (Zero pivot, numerical fact. or iterative refinement problem.), which then made also Ipopt fail on that instance (in most cases).

Changing IPARM_7902a27040 to 3 (which isn't really documented in the MKL Pardiso manual) as it is set for Basel Pardiso seems to help a bit. Also setting IPARM_f97428a6b8 (pivot permutation) to 12 helps a bit.

Switching IPARM_6afa9fa918 (matrix ordering algorithm) from its default 1 (metis) to 0 (minimum degree ordering) also helped on some cases, but I left this one at 1. For Basel Pardiso, Ipopt uses the undocumented option 5, which is not available with MKL Pardiso.

Enabling iterative refinement (IPARM_6a14a728fe seems to help most, so I enabled it, too. For Basel Pardiso, there was a comment that one should decide whether it should be enabled. As I don't have the latter, I kept it at 0 there.

svigerske commented 5 years ago

Comment by @svigerske created at 2013-11-26 10:32:03

With 195a4b65e7, I reset IPARM_6a14a728fe back to 0 when using MKL Pardiso (so it's now the same again as with Basel Pardiso). While setting it to 1 seemed to help to avoid some numerical difficulties on some instances, it also seems to decrease performance on another testset quite much, i.e., in a way that it became even worse than using Ipopt with Mumps. Further, on that testset, it also increased the number of instances where Ipopt failed to converge (from 12 to 28 out of 1314).

svigerske commented 5 years ago

Comment by @svigerske created at 2013-12-18 14:35:35

I should note also here that c1719e17a4 sets IPARM_6a14a728fe to 1 when using MKL Pardiso. A closer look has shown that this is numerically more robust. This will also be the setting in the upcoming maintenance release (3.11.7).

svigerske commented 4 years ago

I don't remember whether there was still something to do here. Some more testing and tuning would always be useful.