Skip to content
Merged
Show file tree
Hide file tree
Changes from all 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
2 changes: 1 addition & 1 deletion .buildkite/pipeline-julia.yml
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
asteps:
steps:
- label: "CPUs -- StaticArrays.jl"
plugins:
- JuliaCI/julia#v1:
Expand Down
8 changes: 4 additions & 4 deletions .buildkite/pipeline.yml
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ steps:
- label: "Nvidia GPUs -- CUDA.jl"
plugins:
- JuliaCI/julia#v1:
version: "1.12"
version: "1.13"
agents:
queue: "cuda"
command: |
Expand All @@ -22,7 +22,7 @@ steps:
- label: "AMD GPUs -- AMDGPU.jl"
plugins:
- JuliaCI/julia#v1:
version: "1.12"
version: "1.13"
agents:
queue: "rocm"
rocmgpu: "*"
Expand All @@ -43,7 +43,7 @@ steps:
- label: "Intel GPUs -- oneAPI.jl"
plugins:
- JuliaCI/julia#v1:
version: "1.12"
version: "1.13"
agents:
queue: "oneapi"
command: |
Expand All @@ -59,7 +59,7 @@ steps:
- label: "Apple M1 GPUs -- Metal.jl"
plugins:
- JuliaCI/julia#v1:
version: "1.12"
version: "1.13"
agents:
queue: "metal"
command: |
Expand Down
4 changes: 2 additions & 2 deletions .github/workflows/release-libkrylov.yml
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ jobs:

- uses: julia-actions/setup-julia@v2
with:
version: '1.12'
version: '1.13'

# -----------------------------------------------------------------------
# Install JuliaC.jl as a Julia app (provides the juliac CLI)
Expand All @@ -104,7 +104,7 @@ jobs:
julia --startup-file=no -e "
import Pkg
Pkg.Registry.add(\"General\")
Pkg.Apps.add(url=\"https://github.com/JuliaLang/JuliaC.jl\", rev=\"v0.3.8\")
Pkg.Apps.add(url=\"https://github.com/JuliaLang/JuliaC.jl\", rev=\"v0.3.10\")
"

- name: Add juliac to PATH
Expand Down
4 changes: 2 additions & 2 deletions .github/workflows/test-libkrylov.yml
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ jobs:

- uses: julia-actions/setup-julia@v2
with:
version: '1.12'
version: '1.13'

# -----------------------------------------------------------------------
# Install gfortran on macOS (not pre-installed on GitHub runners)
Expand All @@ -105,7 +105,7 @@ jobs:
julia --startup-file=no -e "
import Pkg
Pkg.Registry.add(\"General\")
Pkg.Apps.add(url=\"https://github.com/JuliaLang/JuliaC.jl\", rev=\"v0.3.8\")
Pkg.Apps.add(url=\"https://github.com/JuliaLang/JuliaC.jl\", rev=\"v0.3.10\")
"

- name: Add juliac to PATH
Expand Down
22 changes: 4 additions & 18 deletions docs/src/interfaces/building.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@ This page describes how to build the library from source and how to compile and

| Tool | Version |
|------|---------|
| Julia | ≥ 1.12 |
| JuliaC.jl | ≥ 0.3.8 |
| Julia | ≥ 1.13 |
| JuliaC.jl | ≥ 0.3.10 |
| C / Fortran compiler | gcc / clang, gfortran |

[JuliaC.jl](https://github.com/JuliaLang/JuliaC.jl) wraps Julia's `juliac` compiler and adds `--bundle`. This produces a self-contained library that embeds the Julia runtime, so no separate Julia installation is required at run time.
Expand All @@ -19,7 +19,7 @@ All commands run from the root of the Krylov.jl repository.

```bash
# Install JuliaC.jl once (it installs juliac into ~/.julia/bin)
julia -e 'import Pkg; Pkg.Apps.add(url="https://github.com/JuliaLang/JuliaC.jl", rev="v0.3.8")'
julia -e 'import Pkg; Pkg.Apps.add(url="https://github.com/JuliaLang/JuliaC.jl", rev="v0.3.10")'
export PATH="$HOME/.julia/bin:$PATH" # add to ~/.bashrc to make it permanent

# Build the bundle (library + embedded Julia runtime)
Expand All @@ -33,25 +33,11 @@ juliac \

# Generate the headers and copy them next to the library
julia --startup-file=no --project=. interfaces/scripts/generate_header.jl
mkdir -p interfaces/build/include
cp interfaces/include/krylov.h interfaces/build/include/
cp interfaces/include/krylov.f90 interfaces/build/include/

# Copy the SuiteSparse libraries into the bundle (see the note below)
JLIB="$(julia --startup-file=no -e 'print(joinpath(Sys.BINDIR, "..", "lib", "julia"))')"
for name in amd btf camd ccolamd cholmod colamd klu ldl rbio spqr suitesparseconfig umfpack; do
cp -a "$JLIB"/lib"$name".* interfaces/build/lib/julia/
done
```

!!! warning "SuiteSparse must be copied manually"
`juliac --bundle` only copies the libraries it can trace statically. The
SuiteSparse stack (`libbtf`, `libcholmod`, `libumfpack`, ...) is loaded
dynamically by Julia at startup (`SparseArrays` is a dependency of Krylov),
so `juliac` does not see it. Without the copy step above, the bundle runs
fine on a machine that has Julia installed. On a clean machine it fails with
`could not load library "libbtf.so.2"`. On Windows the libraries live in
`Sys.BINDIR` (the `bin/` folder) instead of `lib/julia`.

The `--bundle` flag produces a relocatable directory:

```
Expand Down
25 changes: 16 additions & 9 deletions docs/src/preconditioners.md
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ There exist three variants of preconditioning:
where $P_{\ell}$ and $P_r$ are square and nonsingular.

The left preconditioning preserves the error $x_k - x^{\star}$ whereas the right preconditioning keeps invariant the residual $b - A x_k$.
Two-sided preconditioning is the only variant that allows to preserve the hermicity of a linear system.
Two-sided preconditioning is the only variant that allows to preserve the Hermiticity of a linear system.

!!! note
Because det$(P^{-1}A - \lambda I)$ = det$(A - \lambda P)$ det$(P^{-1})$ = det$(AP^{-1} - \lambda I)$, the eigenvalues of $P^{-1}A$ and $AP^{-1}$ are identical. If $P = LL^{H}$, $L^{-1}AL^{-H}$ also has the same spectrum.
Expand Down Expand Up @@ -49,19 +49,19 @@ A Krylov method dedicated to non-Hermitian linear systems allows the three varia
Methods concerned: [`SYMMLQ`](@ref symmlq), [`CG`](@ref cg), [`CG-LANCZOS`](@ref cg_lanczos), [`CG-LANCZOS-SHIFT`](@ref cg_lanczos_shift), [`CR`](@ref cr), [`CAR`](@ref car), [`MINRES`](@ref minres), [`BLOCK-MINRES`](@ref block_minres), [`MINRES-QLP`](@ref minres_qlp) and [`MINARES`](@ref minares).

When $A$ is Hermitian, we can only use centered preconditioning $L^{-1}AL^{-H}y = L^{-1}b$ with $x = L^{-H}y$.
Centered preconditioning is a special case of two-sided preconditioning with $P_{\ell} = L = P_r^H$ that maintains hermicity.
Centered preconditioning is a special case of two-sided preconditioning with $P_{\ell} = L = P_r^H$ that maintains Hermiticity.
However, there is no need to specify $L$ and one may specify $P_c = LL^H$ or its inverse directly.

| Preconditioners | $P_c^{-1}$ | $P_c$ |
|:---------------:|:-------------------------:|:--------------------:|
| Arguments | `M` with `ldiv=false` | `M` with `ldiv=true` |

!!! warning
The preconditioner `M` must be hermitian and positive definite.
The preconditioner `M` must be Hermitian and positive definite.

### Linear least-squares problems

Methods concerned: [`CGLS`](@ref cgls), [`CGLS-LANCZOS-SHIFT`](@ref cgls_lanczos_shift), [`CRLS`](@ref crls), [`LSLQ`](@ref lslq), [`LSQR`](@ref lsqr) and [`LSMR`](@ref lsmr).
Methods concerned: [`CGLS`](@ref cgls), [`CRLS`](@ref crls), [`LSLQ`](@ref lslq), [`LSQR`](@ref lsqr) and [`LSMR`](@ref lsmr).

| Formulation | Without preconditioning | With preconditioning |
|:---------------------:|:------------------------------------:|:-------------------------------------------:|
Expand All @@ -82,11 +82,14 @@ Methods concerned: [`CGLS`](@ref cgls), [`CGLS-LANCZOS-SHIFT`](@ref cgls_lanczos
| Arguments | `M` with `ldiv=false` | `M` with `ldiv=true` | `N` with `ldiv=false` | `N` with `ldiv=true` |

!!! warning
The preconditioners `M` and `N` must be hermitian and positive definite.
The preconditioners `M` and `N` must be Hermitian and positive definite.

!!! note
[`CGLS`](@ref cgls) and [`CRLS`](@ref crls) only accept `M`; they do not support `N`.

### Linear least-norm problems

Methods concerned: [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`LNLQ`](@ref lnlq), [`CRAIG`](@ref craig) and [`CRAIGMR`](@ref craigmr).
Methods concerned: [`LNLQ`](@ref lnlq), [`CRAIG`](@ref craig) and [`CRAIGMR`](@ref craigmr).

| Formulation | Without preconditioning | With preconditioning |
|:--------------------:|:----------------------------------------------------:|:----------------------------------------------------:|
Expand All @@ -107,7 +110,11 @@ Methods concerned: [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`LNLQ`](@ref lnlq)
| Arguments | `M` with `ldiv=false` | `M` with `ldiv=true` | `N` with `ldiv=false` | `N` with `ldiv=true` |

!!! warning
The preconditioners `M` and `N` must be hermitian and positive definite.
The preconditioners `M` and `N` must be Hermitian and positive definite.

!!! note
[`CGNE`](@ref cgne) and [`CRMR`](@ref crmr) only accept `N`, and it does not play the role of $F^{-1}$ above.
It preconditions the normal equations of the second kind $(AA^H + \lambda^2 I)y = b$, so it is $m \times m$ and leaves the minimized norm unchanged.

### Saddle-point and symmetric quasi-definite systems

Expand All @@ -120,7 +127,7 @@ Methods concerned: [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`LNLQ`](@ref lnlq)
| Arguments | `M` with `ldiv=false` | `M` with `ldiv=true` | `N` with `ldiv=false` | `N` with `ldiv=true` |

!!! warning
The preconditioners `M` and `N` must be hermitian and positive definite.
The preconditioners `M` and `N` must be Hermitian and positive definite.

### Generalized saddle-point and unsymmetric partitioned systems

Expand All @@ -133,7 +140,7 @@ Methods concerned: [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`LNLQ`](@ref lnlq)
| Arguments | `C` and `E` with `ldiv=false` | `C` and `E` with `ldiv=true` | `D` and `F` with `ldiv=false` | `D` and `F` with `ldiv=true` |

!!! note
Our implementations of [`BiLQ`](@ref bilq), [`QMR`](@ref qmr), [`BiLQR`](@ref bilqr), [`USYMLQ`](@ref usymlq), [`USYMQR`](@ref usymqr), [`USYMLQR`](@ref usymlqr) and [`TriLQR`](@ref trilqr) don't support preconditioning.
Our implementations of [`BiLQ`](@ref bilq), [`QMR`](@ref qmr), [`BiLQR`](@ref bilqr), [`USYMLQ`](@ref usymlq), [`USYMQR`](@ref usymqr), [`USYMLQR`](@ref usymlqr), [`TriLQR`](@ref trilqr) and [`CGLS-LANCZOS-SHIFT`](@ref cgls_lanczos_shift) don't support preconditioning.

## Packages that provide preconditioners

Expand Down
18 changes: 5 additions & 13 deletions interfaces/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,8 @@ Pre-built, self-contained bundles for Linux, macOS and Windows are attached to e

| Tool | Version |
|------|---------|
| Julia | ≥ 1.12 |
| [JuliaC.jl](https://github.com/JuliaLang/JuliaC.jl) | ≥ 0.3.8 |
| Julia | ≥ 1.13 |
| [JuliaC.jl](https://github.com/JuliaLang/JuliaC.jl) | ≥ 0.3.10 |
| C / Fortran compiler | gcc / clang, gfortran |

[JuliaC.jl](https://github.com/JuliaLang/JuliaC.jl) wraps Julia's `juliac` compiler and adds `--bundle`, which produces a self-contained library that embeds the Julia runtime (no separate Julia installation required at run time).
Expand All @@ -22,7 +22,7 @@ All commands run from the **root of the Krylov.jl repository**.

```bash
# Install JuliaC.jl once (installs juliac into ~/.julia/bin)
julia -e 'import Pkg; Pkg.Apps.add(url="https://github.com/JuliaLang/JuliaC.jl", rev="v0.3.8")'
julia -e 'import Pkg; Pkg.Apps.add(url="https://github.com/JuliaLang/JuliaC.jl", rev="v0.3.10")'
export PATH="$HOME/.julia/bin:$PATH"

# Build the bundle (library + embedded Julia runtime)
Expand All @@ -36,17 +36,9 @@ juliac \

# Generate the headers and copy them next to the library
julia --startup-file=no --project=. interfaces/scripts/generate_header.jl
mkdir -p interfaces/build/include
cp interfaces/include/krylov.h interfaces/build/include/
cp interfaces/include/krylov.f90 interfaces/build/include/

# Copy the SuiteSparse libraries into the bundle.
# juliac --bundle does not trace libraries that Julia dlopen's at startup
# (SparseArrays -> SuiteSparse_jll). Without this the bundle runs only on a
# machine that already has Julia. On Windows the libs are in Sys.BINDIR (bin/).
JLIB="$(julia --startup-file=no -e 'print(joinpath(Sys.BINDIR, "..", "lib", "julia"))')"
for name in amd btf camd ccolamd cholmod colamd klu ldl rbio spqr suitesparseconfig umfpack; do
cp -a "$JLIB"/lib"$name".* interfaces/build/lib/julia/
done
```

The `--bundle` flag produces a relocatable directory:
Expand All @@ -61,7 +53,7 @@ interfaces/build/
└── krylov.f90
```

> **Windows:** use `--output-lib interfaces/build/bin/libkrylov.dll`; the bundle lands in `build/bin/`, and the SuiteSparse libraries are taken from `Sys.BINDIR` (`bin/`).
> **Windows:** use `--output-lib interfaces/build/bin/libkrylov.dll`; the bundle lands in `build/bin/`.
> **macOS:** replace `.so` with `.dylib`.

Compiling and linking a C or Fortran program against the bundle is documented in the [building guide](https://jso.dev/Krylov.jl/dev/interfaces/building/).
Expand Down
11 changes: 6 additions & 5 deletions interfaces/scripts/generate_stores.jl
Original file line number Diff line number Diff line change
Expand Up @@ -395,28 +395,29 @@ function _typed_solve_ls_mn!(ws::Krylov.KrylovWorkspace{T, FC, S}, fptr_A, fptr_
Cint(0)
end

# λ + radius + M(n) — CGLS, CRLS (CG/CR on the normal equations; single
# preconditioner on the n-space, passed via matvec_M).
# λ + radius + M(m) — CGLS, CRLS (CG/CR on the normal equations; single
# preconditioner on the m-space, passed via matvec_M: M is applied to the
# residual r = b - Ax and to q = A*p, which both have length m).
function _typed_solve_ls_m_radius!(ws::Krylov.KrylovWorkspace{T, FC, S}, fptr_A, fptr_At, fptr_M, fptr_N, b_ptr, c_ptr, userdata, opts) where {T, FC, S}
A = COperator{FC}(ws.m, ws.n, fptr_A, fptr_At, userdata)
b = unsafe_wrap(Vector{FC}, Ptr{FC}(b_ptr), ws.m)
kw = _opts_kw(opts, T); λ = T(opts.lambda); rad = T(opts.radius)
if fptr_M != C_NULL
M = CPreconditioner{FC}(ws.n, fptr_M, userdata)
M = CPreconditioner{FC}(ws.m, fptr_M, userdata)
Krylov.krylov_solve!(ws, A, b; M=M, λ=λ, radius=rad, kw...)
else
Krylov.krylov_solve!(ws, A, b; λ=λ, radius=rad, kw...)
end
Cint(0)
end

# λ + N(n) — CGNE, CRMR (single preconditioner on the n-space, via matvec_N).
# λ + N(m) — CGNE, CRMR (single preconditioner on the m-space, via matvec_N).
function _typed_solve_ls_n!(ws::Krylov.KrylovWorkspace{T, FC, S}, fptr_A, fptr_At, fptr_M, fptr_N, b_ptr, c_ptr, userdata, opts) where {T, FC, S}
A = COperator{FC}(ws.m, ws.n, fptr_A, fptr_At, userdata)
b = unsafe_wrap(Vector{FC}, Ptr{FC}(b_ptr), ws.m)
kw = _opts_kw(opts, T); λ = T(opts.lambda)
if fptr_N != C_NULL
N = CPreconditioner{FC}(ws.n, fptr_N, userdata)
N = CPreconditioner{FC}(ws.m, fptr_N, userdata)
Krylov.krylov_solve!(ws, A, b; N=N, λ=λ, kw...)
else
Krylov.krylov_solve!(ws, A, b; λ=λ, kw...)
Expand Down
11 changes: 6 additions & 5 deletions interfaces/src/c_stores.jl
Original file line number Diff line number Diff line change
Expand Up @@ -442,28 +442,29 @@ function _typed_solve_ls_mn!(ws::Krylov.KrylovWorkspace{T, FC, S}, fptr_A, fptr_
Cint(0)
end

# λ + radius + M(n) — CGLS, CRLS (CG/CR on the normal equations; single
# preconditioner on the n-space, passed via matvec_M).
# λ + radius + M(m) — CGLS, CRLS (CG/CR on the normal equations; single
# preconditioner on the m-space, passed via matvec_M: M is applied to the
# residual r = b - Ax and to q = A*p, which both have length m).
function _typed_solve_ls_m_radius!(ws::Krylov.KrylovWorkspace{T, FC, S}, fptr_A, fptr_At, fptr_M, fptr_N, b_ptr, c_ptr, userdata, opts) where {T, FC, S}
A = COperator{FC}(ws.m, ws.n, fptr_A, fptr_At, userdata)
b = unsafe_wrap(Vector{FC}, Ptr{FC}(b_ptr), ws.m)
kw = _opts_kw(opts, T); λ = T(opts.lambda); rad = T(opts.radius)
if fptr_M != C_NULL
M = CPreconditioner{FC}(ws.n, fptr_M, userdata)
M = CPreconditioner{FC}(ws.m, fptr_M, userdata)
Krylov.krylov_solve!(ws, A, b; M=M, λ=λ, radius=rad, kw...)
else
Krylov.krylov_solve!(ws, A, b; λ=λ, radius=rad, kw...)
end
Cint(0)
end

# λ + N(n) — CGNE, CRMR (single preconditioner on the n-space, via matvec_N).
# λ + N(m) — CGNE, CRMR (single preconditioner on the m-space, via matvec_N).
function _typed_solve_ls_n!(ws::Krylov.KrylovWorkspace{T, FC, S}, fptr_A, fptr_At, fptr_M, fptr_N, b_ptr, c_ptr, userdata, opts) where {T, FC, S}
A = COperator{FC}(ws.m, ws.n, fptr_A, fptr_At, userdata)
b = unsafe_wrap(Vector{FC}, Ptr{FC}(b_ptr), ws.m)
kw = _opts_kw(opts, T); λ = T(opts.lambda)
if fptr_N != C_NULL
N = CPreconditioner{FC}(ws.n, fptr_N, userdata)
N = CPreconditioner{FC}(ws.m, fptr_N, userdata)
Krylov.krylov_solve!(ws, A, b; N=N, λ=λ, kw...)
else
Krylov.krylov_solve!(ws, A, b; λ=λ, kw...)
Expand Down
9 changes: 5 additions & 4 deletions interfaces/test/test_libkrylov.jl
Original file line number Diff line number Diff line change
Expand Up @@ -560,8 +560,9 @@ function test_ls_preconditioner()
c_workspace_free(ws)
end

# CGLS: single preconditioner on the n-space (normal equations), via matvec_M.
_Minv_f64[] = fill(2.0, n)
# CGLS: single preconditioner on the m-space — M is applied to the residual
# r = b - Ax and to q = A*p, both of length m — via matvec_M.
_Minv_f64[] = fill(2.0, m)
ws = c_workspace_create(solver_enum("cgls"), m, n, KRYLOV_FLOAT64)
try
@test solve_with(ws, cb_A, b; cb_At=cb_At, cb_M=CB_M_F64) == 0
Expand All @@ -570,8 +571,8 @@ function test_ls_preconditioner()
c_workspace_free(ws)
end

# CGNE: single preconditioner on the n-space, via matvec_N.
_Ninv_f64[] = fill(2.0, n)
# CGNE: single preconditioner on the m-space, via matvec_N.
_Ninv_f64[] = fill(2.0, m)
ws = c_workspace_create(solver_enum("cgne"), m, n, KRYLOV_FLOAT64)
try
@test solve_with(ws, cb_A, b; cb_At=cb_At, cb_N=CB_N_F64) == 0
Expand Down
2 changes: 1 addition & 1 deletion src/block_krylov_utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -329,7 +329,7 @@ kunmqr!(side :: Char, trans :: Char, A :: AbstractMatrix{T}, tau :: AbstractVect
# β, τ = larfg!(α, x)
#
# Generate an elementary Householder reflector `H = I - τ vvᴴ` such that
# `Hᴴ * [α; x] = [β; 0]` (LAPACK convention; `H` is not hermitian in the complex
# `Hᴴ * [α; x] = [β; 0]` (LAPACK convention; `H` is not Hermitian in the complex
# case), where `v = [1; y]` and `β` is the (signed) Euclidean norm of `[α; x]`.
# On output `x` is overwritten by the tail `y` of `v`.
# """
Expand Down
2 changes: 1 addition & 1 deletion src/car.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# An implementation of CAR for the solution of hermitian positive definite linear systems.
# An implementation of CAR for the solution of Hermitian positive definite linear systems.
#
# This method is described in
#
Expand Down
2 changes: 1 addition & 1 deletion src/cgls.jl
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ For an in-place variant that reuses memory across solves, see [`cgls!`](@ref).

#### Keyword arguments

* `M`: linear operator that models a Hermitian positive-definite matrix of size `n` used for preconditioning;
* `M`: linear operator that models a Hermitian positive-definite matrix of size `m` used for preconditioning;
* `ldiv`: define whether the preconditioner uses `ldiv!` or `mul!`;
* `radius`: add the trust-region constraint ‖x‖ ≤ `radius` if `radius > 0`. Useful to compute a step in a trust-region method for optimization;
* `λ`: regularization parameter;
Expand Down
Loading
Loading