Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adaptations for Panua Pardiso #75

Merged
merged 17 commits into from
Mar 1, 2024
Merged
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,15 +4,21 @@ repo = "https://github.com/JuliaSparse/Pardiso.jl.git"
version = "0.5.5"

[deps]
Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MKL_jll = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[compat]
Compat = "4.10.0"
j-fu marked this conversation as resolved.
Show resolved Hide resolved
MKL_jll = "2021, 2022, 2023, 2024"
Libdl= "1.6"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Formatting

LinearAlgebra = "1.6"
SparseArrays = "1.6"
julia = "1.6"


j-fu marked this conversation as resolved.
Show resolved Hide resolved
[extras]
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Expand Down
56 changes: 25 additions & 31 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -31,35 +31,29 @@ If you instead use a self installed MKL, follow these instructions:
executing something like `source /opt/intel/oneapi/setvars.sh intel64` or
running `"C:\Program Files (x86)\IntelSWTools\compilers_and_libraries\windows\mkl\bin\mklvars.bat" intel64`
* Run `Pkg.build("Pardiso", verbose=true)`
* Run `Pardiso.show_build_log()` to see the build log for additional
information.
* Note that the `MKLROOT` environment variable must be set, and `LD_LIBRARY_PATH` must contain `$MKLROOT/lib`
whenever using the library.

### PARDISO from [panua.ch](https://pardiso-project.org) ("ProjectPardiso")

* Unzip the download file `panua-pardiso-yyyymmdd-os.zip` to some folder and set the environment variable `JULIA_PARDISO`
to the `lib` subdirectory of this folder.
* (Project-Pardiso 6.0): Put the PARDISO library `libpardisoVVV-WIN-X86-64.dll`, `libpardisoVVV-GNUXXX-X86-64.so` or
`libpardisoVVV-MACOS-X86-64.dylib` in a folder somewhere and set the environment variable `JULIA_PARDISO` to that folder.
For example, create an entry `ENV["JULIA_PARDISO"] = "/Users/Someone/Pardiso"` in the
`.julia/config/startup.jl` file and download the Pardiso library to that folder.
* Eventually, run `Pardiso.show_build_log()` to see the build log for additional information.
* Note that the `MKLROOT` environment variable must be set, and `LD_LIBRARY_PATH` must contain `$MKLROOT/lib` whenever using the library.

### PARDISO from [panua.ch](https://panua.ch) ("PanuaPardiso", formerly "ProjectPardiso")

* Unzip the download file `panua-pardiso-yyyymmdd-os.zip` to some folder and set the environment variable `JULIA_PARDISO` to the `lib` subdirectory of this folder. For example, create an entry `ENV["JULIA_PARDISO"] = "/Users/Someone/panua-pardiso-yyyymmdd-os/lib"` in `.julia/config/startup.jl`. If you have a valid license for the predecessor from pardiso-project.org, put the PARDISO library to a subdirectory denoted by `ENV["JULIA_PARDISO"]`,
evenutally rename it to `libpardiso.so`.
* Perform the platform specific steps below
* Run `Pkg.build("Pardiso")`
* Run `Pardiso.show_build_log()` to see the build log for additional information.
* Run `Pkg.build("Pardiso", verbose=true)`
* Eventually, run `Pardiso.show_build_log()` to see the build log for additional information.

Note: Weird errors and problems with MKL Pardiso has been observed when ProjectPardiso is enabled
(likely because some library that is needed by ProjectPardiso is problematic with MKL).
If you want to use MKL Pardiso it is better ot just disable ProjectPardiso by not setting
the environment variable `JULIA_PARDISO` (and rerunning `build Pardiso`).
Note: In the past, weird errors and problems with MKL Pardiso had been observed when PanuaPardiso is enabled
(likely because some library that is needed by PanauaPardiso was problematic with MKL).
In that case of you want to use MKL Pardiso it is better ot just disable PanuaPardiso by not setting
the environment variable `JULIA_PARDISO` (and rerunning `Pkg.build("Pardiso")`).

##### Linux / macOS specific

* Make sure that the version of `gfortran` corresponding to the pardiso library is installed.
* Make sure OpenMP is installed.
* Install a (fast) installation of a BLAS and LAPACK (this should preferably be single threaded since PARDISO handles threading itself), using for example [OpenBLAS](https://github.com/xianyi/OpenBLAS/wiki/Precompiled-installation-packages)

`gfortran` and OpenMP usually come with recent version of gcc/gfortran. On Linux, ProjectPardiso
`gfortran` and OpenMP usually come with recent version of gcc/gfortran. On Linux, Panua Pardiso
looks for libraries `libgfortran.so` and `libgomp.so` . They may be named differently on your system.
In this situation you may try to create links to them with names known to
`Pardiso.jl` (bash; pathnames serve as examples here):
Expand All @@ -76,7 +70,7 @@ This section will explain how to solve equations using `Pardiso.jl` with the def

## Creating the PardisoSolver

A `PardisoSolver` is created with `PardisoSolver()` for solving with ProjectPardiso or `MKLPardisoSolver()` for solving with MKL PARDISO. This object will hold the settings of the solver and will be passed into the solve functions. In the following sections an instance of a `PardisoSolver` or an `MKLPardisoSolver()` will be referred to as `ps`.
A `PardisoSolver` is created with `PardisoSolver()` for solving with PanuaPardiso or `MKLPardisoSolver()` for solving with MKL PARDISO. This object will hold the settings of the solver and will be passed into the solve functions. In the following sections an instance of a `PardisoSolver` or an `MKLPardisoSolver()` will be referred to as `ps`.

### Solving

Expand Down Expand Up @@ -115,7 +109,7 @@ julia> X
-1.17295 8.47922
```

### Schur Complement (ProjectPardiso only)
### Schur Complement (PanuaPardiso only)

Given a partitioned matrix `M = [A B; C D]`, the Schur complement of `A` in `M` is `S = D-CA⁻¹B`.
This can be found with the function `schur_complement` with the following signatures:
Expand Down Expand Up @@ -162,7 +156,7 @@ At present there seems to be an instability in the Schur complement computation

### Setting the number of threads

The number of threads to use is set in different ways for MKL PARDISO and ProjectPardiso.
The number of threads to use is set in different ways for MKL PARDISO and PanuaPardiso.

#### MKL PARDISO

Expand All @@ -171,7 +165,7 @@ set_nprocs!(ps, i) # Sets the number of threads to use
get_nprocs(ps) # Gets the number of threads being used
```

#### ProjectPardiso
#### PanuaPardiso

The number of threads are set at the creation of the `PardisoSolver` by looking for the environment variable `OMP_NUM_THREADS`. This can be done in Julia with for example `ENV["OMP_NUM_THREADS"] = 2`. **Note:** `OMP_NUM_THREADS` must be set *before* `Pardiso` is loaded and can not be changed during runtime.

Expand All @@ -181,7 +175,7 @@ The number of threads used by a `PardisoSolver` can be retrieved with `get_nproc

This section discusses some more advanced usage of `Pardiso.jl`.

For terminology in this section please refer to the [ProjectPardiso manual](http://www.pardiso-project.org/manual/manual.pdf) and the [MKL PARDISO section](https://software.intel.com/en-us/node/470282).
For terminology in this section please refer to the [PanuaPardiso manual](http://panua.ch/manual/manual.pdf) and the [oneMKL PARDISO section](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-0/onemkl-pardiso-parallel-direct-sparse-solver-iface.html).

After using functionality in this section, calls should no longer be made to the `solve` functions but instead directly to the function

Expand Down Expand Up @@ -213,8 +207,8 @@ The matrix type can be explicitly set with `set_matrixtype!(ps, key)` where the

The matrix type for a solver can be retrieved with `get_matrixtype(ps)`.

### Setting the solver (ProjectPardiso only)
ProjectPardiso supports direct and iterative solvers. The solver is set with `set_solver!(ps, key)` where the key has the following meaning:
### Setting the solver (PanuaPardiso only)
PanuatPardiso supports direct and iterative solvers. The solver is set with `set_solver!(ps, key)` where the key has the following meaning:

| enum | integer | Solver |
|--------------------|---------|----------------------------------|
Expand Down Expand Up @@ -242,15 +236,15 @@ Depending on the phase calls to `solve` (and `pardiso` which is mentioned later)
| RELEASE_ALL | -1 | Release all internal memory for all matrices |

### Setting `IPARM` and `DPARM` explicitly
Advanced users likely want to explicitly set and retrieve the `IPARM` and `DPARM` (ProjectPardiso only) parameters.
Advanced users likely want to explicitly set and retrieve the `IPARM` and `DPARM` (PanuaPardiso only) parameters.
This can be done with the getters and setters:

```jl
get_iparm(ps, i) # Gets IPARM[i]
get_iparms(ps) # Gets IPARM
set_iparm!(ps, i, v) # Sets IPARM[i] = v

# ProjectPardiso only
# PanuaPardiso only
get_dparm(ps, i) # Gets DPARM[i]
get_dparms(ps) # Gets DPARM
set_dparm!(ps, i, v) # Sets DPARM[i] = v
Expand All @@ -269,7 +263,7 @@ It is possible for Pardiso to print out timings and statistics when solving. Thi

### Matrix and vector checkers

ProjectPardiso comes with a few matrix and vector checkers to check the consistency and integrity of the input data. These can be called with the functions:
PanuaPardiso comes with a few matrix and vector checkers to check the consistency and integrity of the input data. These can be called with the functions:

```jl
printstats(ps, A, B)
Expand All @@ -294,7 +288,7 @@ get_perm(ps)
set_perm!(ps, perm) # Perm is a Vector{Int}
```

### Schur Complement (ProjectPardiso only)
### Schur Complement (PanuaPardiso only)

The `pardiso(ps,...)` syntax can be used to compute the Schur compelement (as described below). The answer can be retrieved with `pardisogetschur(ps)`.

Expand Down
15 changes: 4 additions & 11 deletions src/Pardiso.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
__precompile__()

module Pardiso
using Compat

if !isfile(joinpath(@__DIR__, "..", "deps", "deps.jl"))
error("""please run Pkg.build("Pardiso") before loading the package""")
Expand Down Expand Up @@ -51,6 +52,7 @@
export get_matrix
export schur_complement, pardisogetschur
export fix_iparm!
@compat public mkl_is_available, panua_is_available

struct PardisoException <: Exception
info::String
Expand Down Expand Up @@ -112,6 +114,8 @@
const pardiso_get_schur_f = Ref{Ptr}()
const PARDISO_LOADED = Ref(false)

panua_is_available() = PARDISO_LOADED[]

Check warning on line 117 in src/Pardiso.jl

View check run for this annotation

Codecov / codecov/patch

src/Pardiso.jl#L117

Added line #L117 was not covered by tests

function __init__()
global MKL_LOAD_FAILED
if LOCAL_MKL_FOUND
Expand Down Expand Up @@ -275,18 +279,7 @@
pardisoinit(ps)
fix_iparm!(ps, T)
try
if isa(ps,PardisoSolver)
# Pardiso7 does not report an error when A is not positive definite,
# but nevertheless does not touch X.
h=hash(X)
end
pardiso(ps, X, get_matrix(ps, A, T), B)
if isa(ps,PardisoSolver)
if hash(X)==h
# Pardiso7 remedy.
throw(PardisoPosDefException("Matrix probably not positive definite"))
end
end
catch e
set_phase!(ps, RELEASE_ALL)
pardiso(ps, X, A, B)
Expand Down
4 changes: 2 additions & 2 deletions src/project_pardiso.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@
end

function PardisoSolver()
if !PARDISO_LOADED[]
error("pardiso library was not loaded")
if !panua_is_available()
error("Panua pardiso library was not loaded")

Check warning on line 16 in src/project_pardiso.jl

View check run for this annotation

Codecov / codecov/patch

src/project_pardiso.jl#L15-L16

Added lines #L15 - L16 were not covered by tests
end

pt = zeros(Int, 64)
Expand Down
Loading