From 8ac7f9a642f719d1c5037b7922883c36385427d7 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 10 Sep 2026 20:42:55 -0600 Subject: [PATCH 1/3] [CI] Upgrade the version of Julia and JuliaC --- .buildkite/pipeline-julia.yml | 2 +- .buildkite/pipeline.yml | 8 ++++---- .github/workflows/release-libkrylov.yml | 4 ++-- .github/workflows/test-libkrylov.yml | 4 ++-- docs/src/interfaces/building.md | 22 ++++------------------ interfaces/README.md | 18 +++++------------- 6 files changed, 18 insertions(+), 40 deletions(-) diff --git a/.buildkite/pipeline-julia.yml b/.buildkite/pipeline-julia.yml index 4ee7039cd..e6cae4a51 100644 --- a/.buildkite/pipeline-julia.yml +++ b/.buildkite/pipeline-julia.yml @@ -1,4 +1,4 @@ -asteps: +steps: - label: "CPUs -- StaticArrays.jl" plugins: - JuliaCI/julia#v1: diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 35edc2351..971a1c225 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -2,7 +2,7 @@ steps: - label: "Nvidia GPUs -- CUDA.jl" plugins: - JuliaCI/julia#v1: - version: "1.12" + version: "1.13" agents: queue: "cuda" command: | @@ -22,7 +22,7 @@ steps: - label: "AMD GPUs -- AMDGPU.jl" plugins: - JuliaCI/julia#v1: - version: "1.12" + version: "1.13" agents: queue: "rocm" rocmgpu: "*" @@ -43,7 +43,7 @@ steps: - label: "Intel GPUs -- oneAPI.jl" plugins: - JuliaCI/julia#v1: - version: "1.12" + version: "1.13" agents: queue: "oneapi" command: | @@ -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: | diff --git a/.github/workflows/release-libkrylov.yml b/.github/workflows/release-libkrylov.yml index 90d899c31..be72f6fc3 100644 --- a/.github/workflows/release-libkrylov.yml +++ b/.github/workflows/release-libkrylov.yml @@ -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) @@ -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 diff --git a/.github/workflows/test-libkrylov.yml b/.github/workflows/test-libkrylov.yml index 5804cb6f5..1588af730 100644 --- a/.github/workflows/test-libkrylov.yml +++ b/.github/workflows/test-libkrylov.yml @@ -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) @@ -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 diff --git a/docs/src/interfaces/building.md b/docs/src/interfaces/building.md index 07aeefcfb..bd8066d50 100644 --- a/docs/src/interfaces/building.md +++ b/docs/src/interfaces/building.md @@ -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. @@ -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) @@ -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: ``` diff --git a/interfaces/README.md b/interfaces/README.md index a4d12f49f..a5f41062f 100644 --- a/interfaces/README.md +++ b/interfaces/README.md @@ -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). @@ -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) @@ -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: @@ -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/). From f805b68891dcc36143af39e230d81ad212114a6d Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 10 Sep 2026 21:59:09 -0600 Subject: [PATCH 2/3] Fix a bug in the C interface --- interfaces/scripts/generate_stores.jl | 11 ++++++----- interfaces/src/c_stores.jl | 11 ++++++----- interfaces/test/test_libkrylov.jl | 9 +++++---- 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/interfaces/scripts/generate_stores.jl b/interfaces/scripts/generate_stores.jl index 5da1df321..b34a0318c 100644 --- a/interfaces/scripts/generate_stores.jl +++ b/interfaces/scripts/generate_stores.jl @@ -395,14 +395,15 @@ 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...) @@ -410,13 +411,13 @@ function _typed_solve_ls_m_radius!(ws::Krylov.KrylovWorkspace{T, FC, S}, fptr_A, 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...) diff --git a/interfaces/src/c_stores.jl b/interfaces/src/c_stores.jl index 0a7beb5b5..6268d1b16 100644 --- a/interfaces/src/c_stores.jl +++ b/interfaces/src/c_stores.jl @@ -442,14 +442,15 @@ 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...) @@ -457,13 +458,13 @@ function _typed_solve_ls_m_radius!(ws::Krylov.KrylovWorkspace{T, FC, S}, fptr_A, 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...) diff --git a/interfaces/test/test_libkrylov.jl b/interfaces/test/test_libkrylov.jl index 5b7ed5dab..f16f4350b 100644 --- a/interfaces/test/test_libkrylov.jl +++ b/interfaces/test/test_libkrylov.jl @@ -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 @@ -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 From 6c17e63ab78f563828b7ffd20999bcf1c2572419 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 10 Sep 2026 21:59:32 -0600 Subject: [PATCH 3/3] Fix a few typos in the documentation --- docs/src/preconditioners.md | 25 ++++++++++++++++--------- src/block_krylov_utils.jl | 2 +- src/car.jl | 2 +- src/cgls.jl | 2 +- src/cgls_lanczos_shift.jl | 22 ++++++---------------- src/cgne.jl | 2 +- src/crls.jl | 2 +- src/crmr.jl | 2 +- src/workspace_accessors.jl | 2 +- 9 files changed, 29 insertions(+), 32 deletions(-) diff --git a/docs/src/preconditioners.md b/docs/src/preconditioners.md index a77d8d553..672280b68 100644 --- a/docs/src/preconditioners.md +++ b/docs/src/preconditioners.md @@ -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. @@ -49,7 +49,7 @@ 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$ | @@ -57,11 +57,11 @@ However, there is no need to specify $L$ and one may specify $P_c = LL^H$ or its | 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 | |:---------------------:|:------------------------------------:|:-------------------------------------------:| @@ -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 | |:--------------------:|:----------------------------------------------------:|:----------------------------------------------------:| @@ -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 @@ -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 @@ -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 diff --git a/src/block_krylov_utils.jl b/src/block_krylov_utils.jl index 45f5a9605..93fd0059d 100644 --- a/src/block_krylov_utils.jl +++ b/src/block_krylov_utils.jl @@ -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`. # """ diff --git a/src/car.jl b/src/car.jl index 5d0e8970d..77ce0d8ce 100644 --- a/src/car.jl +++ b/src/car.jl @@ -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 # diff --git a/src/cgls.jl b/src/cgls.jl index e6021d071..17414f42f 100644 --- a/src/cgls.jl +++ b/src/cgls.jl @@ -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; diff --git a/src/cgls_lanczos_shift.jl b/src/cgls_lanczos_shift.jl index 13168b060..b7a4c7e45 100644 --- a/src/cgls_lanczos_shift.jl +++ b/src/cgls_lanczos_shift.jl @@ -13,8 +13,8 @@ export cgls_lanczos_shift, cgls_lanczos_shift! """ (x, stats) = cgls_lanczos_shift(A, b::AbstractVector{FC}, shifts::AbstractVector{T}; - M=I, λ::T=zero(T), atol::T=√eps(T), rtol::T=√eps(T), - radius::T=zero(T), itmax::Int=0, verbose::Int=0, + atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0, + timemax::Float64=Inf, verbose::Int=0, history::Bool=false, callback=workspace->false, iostream::IO=kstdout) @@ -50,8 +50,6 @@ For an in-place variant that reuses memory across solves, see [`cgls_lanczos_shi #### Keyword arguments -* `M`: linear operator that models a Hermitian positive-definite matrix of size `n` used for preconditioning; -* `ldiv`: define whether the preconditioner uses `ldiv!` or `mul!`; * `atol`: absolute stopping tolerance based on the residual norm; * `rtol`: relative stopping tolerance based on the residual norm; * `itmax`: the maximum number of iterations. If `itmax=0`, the default number of iterations is set to `m+n`; @@ -88,9 +86,7 @@ def_args_cgls_lanczos_shift = (:(A ), :(b::AbstractVector{FC} ), :(shifts::AbstractVector{T})) -def_kwargs_cgls_lanczos_shift = (:(; M = I ), - :(; ldiv::Bool = false ), - :(; atol::T = √eps(T) ), +def_kwargs_cgls_lanczos_shift = (:(; atol::T = √eps(T) ), :(; rtol::T = √eps(T) ), :(; itmax::Int = 0 ), :(; timemax::Float64 = Inf ), @@ -102,7 +98,7 @@ def_kwargs_cgls_lanczos_shift = (:(; M = I ), def_kwargs_cgls_lanczos_shift = extract_parameters.(def_kwargs_cgls_lanczos_shift) args_cgls_lanczos_shift = (:A, :b, :shifts) -kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) +kwargs_cgls_lanczos_shift = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) @eval begin function cgls_lanczos_shift!(workspace :: CglsLanczosShiftWorkspace{T,FC,Sm,Sn}, $(def_args_cgls_lanczos_shift...); $(def_kwargs_cgls_lanczos_shift...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, Sm <: AbstractVector{FC}, Sn <: AbstractVector{FC}} @@ -119,10 +115,6 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose nshifts == workspace.nshifts || error("workspace.nshifts = $(workspace.nshifts) is inconsistent with length(shifts) = $nshifts") (verbose > 0) && @printf(iostream, "CGLS-LANCZOS-SHIFT: system of %d equations in %d variables with %d shifts\n", m, n, nshifts) - # Tests M = Iₙ - MisI = (M === I) - !MisI && error("Preconditioner `M` is not supported.") - # Check type consistency eltype(A) == FC || @warn "eltype(A) ≠ $FC. This could lead to errors or additional allocations in operator-vector products." ktypeof(b) == Sm || error("ktypeof(b) must be equal to $Sm") @@ -131,7 +123,6 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose Aᴴ = A' # Set up workspace. - allocate_if(!MisI, workspace, :v, Sn, workspace.Mv) # The length of v is n v, u_prev, u, u_next = workspace.Mv, workspace.u_prev, workspace.u, workspace.u_next x, p, σ, δhat = workspace.x, workspace.p, workspace.σ, workspace.δhat ω, γ, rNorms, converged = workspace.ω, workspace.γ, workspace.rNorms, workspace.converged @@ -148,7 +139,7 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose kcopy!(m, u, b) # u ← b kfill!(u_prev, zero(FC)) kmul!(v, Aᴴ, u) # v₁ ← Aᴴ * b - β = knorm_elliptic(n, v, v) # β₁ = v₁ᵀ M v₁ + β = knorm_elliptic(n, v, v) # β₁ = ‖v₁‖ kfill!(rNorms, β) if history for i = 1 : nshifts @@ -211,13 +202,12 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose kaxpy!(m, -δ, u, u_next) # uₖ₊₁ = u_nextₖ - δₖuₖ - βₖuₖ₋₁ kaxpy!(m, -β, u_prev, u_next) kmul!(v, Aᴴ, u_next) # vₖ₊₁ = Aᴴuₖ₊₁ - β = knorm_elliptic(n, v, v) # βₖ₊₁ = vₖ₊₁ᵀ M vₖ₊₁ + β = knorm_elliptic(n, v, v) # βₖ₊₁ = ‖vₖ₊₁‖ kdiv!(n, v, β) # vₖ₊₁ = vₖ₊₁ / βₖ₊₁ kdiv!(m, u_next, β) # uₖ₊₁ = uₖ₊₁ / βₖ₊₁ kcopy!(m, u_prev, u) # u_prev ← u kcopy!(m, u, u_next) # u ← u_next - MisI || (ρ = kdotr(n, v, v)) for i = 1 : nshifts δhat[i] = δ + ρ * shifts[i] γ[i] = inv(δhat[i] - ω[i] / γ[i]) diff --git a/src/cgne.jl b/src/cgne.jl index 053015642..3eb1ea1a2 100644 --- a/src/cgne.jl +++ b/src/cgne.jl @@ -74,7 +74,7 @@ For an in-place variant that reuses memory across solves, see [`cgne!`](@ref). #### Keyword arguments -* `N`: linear operator that models a Hermitian positive-definite matrix of size `n` used for preconditioning; +* `N`: linear operator that models a Hermitian positive-definite matrix of size `m` used for preconditioning; * `ldiv`: define whether the preconditioner uses `ldiv!` or `mul!`; * `λ`: regularization parameter; * `atol`: absolute stopping tolerance based on the residual norm; diff --git a/src/crls.jl b/src/crls.jl index 75c8ce742..50c5bfba1 100644 --- a/src/crls.jl +++ b/src/crls.jl @@ -59,7 +59,7 @@ For an in-place variant that reuses memory across solves, see [`crls!`](@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; diff --git a/src/crmr.jl b/src/crmr.jl index c8e290193..6d281a7b1 100644 --- a/src/crmr.jl +++ b/src/crmr.jl @@ -72,7 +72,7 @@ For an in-place variant that reuses memory across solves, see [`crmr!`](@ref). #### Keyword arguments -* `N`: linear operator that models a Hermitian positive-definite matrix of size `n` used for preconditioning; +* `N`: linear operator that models a Hermitian positive-definite matrix of size `m` used for preconditioning; * `ldiv`: define whether the preconditioner uses `ldiv!` or `mul!`; * `λ`: regularization parameter; * `atol`: absolute stopping tolerance based on the residual norm; diff --git a/src/workspace_accessors.jl b/src/workspace_accessors.jl index 900aef7b2..a8f8685c4 100644 --- a/src/workspace_accessors.jl +++ b/src/workspace_accessors.jl @@ -69,7 +69,7 @@ Return the number of iterations performed by the Krylov method associated with ` The number of iterations alone is not a reliable basis for comparing different Krylov methods, since the work performed in each iteration can vary significantly. -For a fairer performance comparison, use the total number of operator-vector products with `A` and `A'` (see [Aprod_count](@ref) and [Atprod_count](@ref)). +For a fairer performance comparison, use the total number of operator-vector products with `A` and `A'` (see [`Aprod_count`](@ref) and [`Atprod_count`](@ref)). """ function iteration_count end