Add m and n in all Krylov solvers

src/bicgstab.jl

@@ -101,7 +101,7 @@ function bicgstab!(solver :: BicgstabSolver{T,FC,S}, A, b :: AbstractVector{FC};
                    itmax :: Int=0, verbose :: Int=0, history :: Bool=false,
                    ldiv :: Bool=false, callback = solver -> false) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: DenseVector{FC}}
 
-  n, m = size(A)
+  m, n = size(A)
   m == n || error("System must be square")
   length(b) == m || error("Inconsistent problem size")
   (verbose > 0) && @printf("BICGSTAB: system of size %d\n", n)

src/bilq.jl

@@ -89,7 +89,7 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
                itmax :: Int=0, verbose :: Int=0, history :: Bool=false,
                callback = solver -> false) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: DenseVector{FC}}
 
-  n, m = size(A)
+  m, n = size(A)
   m == n || error("System must be square")
   length(b) == m || error("Inconsistent problem size")
   (verbose > 0) && @printf("BILQ: system of size %d\n", n)

src/bilqr.jl

@@ -94,7 +94,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
                 itmax :: Int=0, verbose :: Int=0, history :: Bool=false,
                 callback = solver -> false) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: DenseVector{FC}}
 
-  n, m = size(A)
+  m, n = size(A)
   m == n || error("Systems must be square")
   length(b) == m || error("Inconsistent problem size")
   length(c) == n || error("Inconsistent problem size")

src/cg.jl

@@ -97,7 +97,7 @@ function cg!(solver :: CgSolver{T,FC,S}, A, b :: AbstractVector{FC};
 
   linesearch && (radius > 0) && error("`linesearch` set to `true` but trust-region radius > 0")
 
-  n, m = size(A)
+  m, n = size(A)
   m == n || error("System must be square")
   length(b) == n || error("Inconsistent problem size")
   (verbose > 0) && @printf("CG: system of %d equations in %d variables\n", n, n)

src/cg_lanczos.jl

@@ -89,7 +89,7 @@ function cg_lanczos!(solver :: CgLanczosSolver{T,FC,S}, A, b :: AbstractVector{F
                      check_curvature :: Bool=false, verbose :: Int=0, history :: Bool=false,
                      ldiv :: Bool=false, callback = solver -> false) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: DenseVector{FC}}
 
-  n, m = size(A)
+  m, n = size(A)
   m == n || error("System must be square")
   length(b) == n || error("Inconsistent problem size")
   (verbose > 0) && @printf("CG Lanczos: system of %d equations in %d variables\n", n, n)

src/cg_lanczos_shift.jl

@@ -77,7 +77,7 @@ function cg_lanczos_shift!(solver :: CgLanczosShiftSolver{T,FC,S}, A, b :: Abstr
                            verbose :: Int=0, history :: Bool=false,
                            ldiv :: Bool=false, callback = solver -> false) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: DenseVector{FC}}
 
-  n, m = size(A)
+  m, n = size(A)
   m == n || error("System must be square")
   length(b) == n || error("Inconsistent problem size")
 

src/cr.jl

@@ -95,7 +95,8 @@ function cr!(solver :: CrSolver{T,FC,S}, A, b :: AbstractVector{FC};
              ldiv :: Bool=false, callback = solver -> false) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: DenseVector{FC}}
 
   linesearch && (radius > 0) && error("'linesearch' set to 'true' but radius > 0")
-  n, m = size(A)
+
+  m, n = size(A)
   m == n || error("System must be square")
   length(b) == n || error("Inconsistent problem size")
   (verbose > 0) && @printf("CR: system of %d equations in %d variables\n", n, n)

src/krylov_processes.jl

@@ -19,7 +19,7 @@ export hermitian_lanczos, nonhermitian_lanczos, arnoldi, golub_kahan, saunders_s
 * C. Lanczos, [*An Iteration Method for the Solution of the Eigenvalue Problem of Linear Differential and Integral Operators*](https://doi.org/10.6028/jres.045.026), Journal of Research of the National Bureau of Standards, 45(4), pp. 225--280, 1950.
 """
 function hermitian_lanczos(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOrComplex
-  n, m = size(A)
+  m, n = size(A)
   R = real(FC)
   S = ktypeof(b)
   M = vector_to_matrix(S)
@@ -90,7 +90,7 @@ end
 * C. Lanczos, [*An Iteration Method for the Solution of the Eigenvalue Problem of Linear Differential and Integral Operators*](https://doi.org/10.6028/jres.045.026), Journal of Research of the National Bureau of Standards, 45(4), pp. 225--280, 1950.
 """
 function nonhermitian_lanczos(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k::Int) where FC <: FloatOrComplex
-  n, m = size(A)
+  m, n = size(A)
   Aᴴ = A'
   S = ktypeof(b)
   M = vector_to_matrix(S)
@@ -179,7 +179,7 @@ end
 * W. E. Arnoldi, [*The principle of minimized iterations in the solution of the matrix eigenvalue problem*](https://doi.org/10.1090/qam/42792), Quarterly of Applied Mathematics, 9, pp. 17--29, 1951.
 """
 function arnoldi(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOrComplex
-  n, m = size(A)
+  m, n = size(A)
   S = ktypeof(b)
   M = vector_to_matrix(S)
 
@@ -239,7 +239,7 @@ end
 * G. H. Golub and W. Kahan, [*Calculating the Singular Values and Pseudo-Inverse of a Matrix*](https://doi.org/10.1137/0702016), SIAM Journal on Numerical Analysis, 2(2), pp. 225--224, 1965.
 """
 function golub_kahan(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOrComplex
-  n, m = size(A)
+  m, n = size(A)
   R = real(FC)
   Aᴴ = A'
   S = ktypeof(b)
@@ -259,8 +259,8 @@ function golub_kahan(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOrComple
   rowval[2k+1] = k+1
   colptr[k+2] = 2k+2
 
-  V = M(undef, m, k+1)
-  U = M(undef, n, k+1)
+  V = M(undef, n, k+1)
+  U = M(undef, m, k+1)
   L = SparseMatrixCSC(k+1, k+1, colptr, rowval, nzval)
 
   for i = 1:k
@@ -270,21 +270,21 @@ function golub_kahan(A, b::AbstractVector{FC}, k::Int) where FC <: FloatOrComple
     vᵢ₊₁ = p = view(V,:,i+1)
     if i == 1
       wᵢ = vᵢ
-      βᵢ = @knrm2(n, b)
+      βᵢ = @knrm2(m, b)
       uᵢ .= b ./ βᵢ
       mul!(wᵢ, Aᴴ, uᵢ)
-      αᵢ = @knrm2(m, wᵢ)
+      αᵢ = @knrm2(n, wᵢ)
       L[1,1] = αᵢ
       vᵢ .= wᵢ ./ αᵢ
     end
     mul!(q, A, vᵢ)
     αᵢ = L[i,i] 
-    @kaxpy!(n, -αᵢ, uᵢ, q)
-    βᵢ₊₁ = @knrm2(n, q)
+    @kaxpy!(m, -αᵢ, uᵢ, q)
+    βᵢ₊₁ = @knrm2(m, q)
     uᵢ₊₁ .= q ./ βᵢ₊₁
     mul!(p, Aᴴ, uᵢ₊₁)
-    @kaxpy!(m, -βᵢ₊₁, vᵢ, p)
-    αᵢ₊₁ = @knrm2(m, p)
+    @kaxpy!(n, -βᵢ₊₁, vᵢ, p)
+    αᵢ₊₁ = @knrm2(n, p)
     vᵢ₊₁ .= p ./ αᵢ₊₁
     L[i+1,i]   = βᵢ₊₁
     L[i+1,i+1] = αᵢ₊₁
@@ -314,7 +314,7 @@ end
 * M. A. Saunders, H. D. Simon, and E. L. Yip, [*Two Conjugate-Gradient-Type Methods for Unsymmetric Linear Equations*](https://doi.org/10.1137/0725052), SIAM Journal on Numerical Analysis, 25(4), pp. 927--940, 1988.
 """
 function saunders_simon_yip(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k::Int) where FC <: FloatOrComplex
-  n, m = size(A)
+  m, n = size(A)
   Aᴴ = A'
   S = ktypeof(b)
   M = vector_to_matrix(S)
@@ -337,8 +337,8 @@ function saunders_simon_yip(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k::
     end
   end
 
-  V = M(undef, n, k+1)
-  U = M(undef, m, k+1)
+  V = M(undef, m, k+1)
+  U = M(undef, n, k+1)
   T = SparseMatrixCSC(k+1, k, colptr, rowval, nzval_T)
   Tᴴ = SparseMatrixCSC(k+1, k, colptr, rowval, nzval_Tᴴ)
 
@@ -348,8 +348,8 @@ function saunders_simon_yip(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k::
     vᵢ₊₁ = q = view(V,:,i+1)
     uᵢ₊₁ = p = view(U,:,i+1)
     if i == 1
-      β = @knrm2(n, b)
-      γ = @knrm2(m, c)
+      β = @knrm2(m, b)
+      γ = @knrm2(n, c)
       vᵢ .= b ./ β
       uᵢ .= c ./ γ
     end
@@ -360,16 +360,16 @@ function saunders_simon_yip(A, b::AbstractVector{FC}, c::AbstractVector{FC}, k::
       uᵢ₋₁ = view(U,:,i-1)
       βᵢ = T[i,i-1]
       γᵢ = T[i-1,i]
-      @kaxpy!(n, -γᵢ, vᵢ₋₁, q)
-      @kaxpy!(m, -βᵢ, uᵢ₋₁, p)
+      @kaxpy!(m, -γᵢ, vᵢ₋₁, q)
+      @kaxpy!(n, -βᵢ, uᵢ₋₁, p)
     end
-    αᵢ = @kdot(n, vᵢ, q)
+    αᵢ = @kdot(m, vᵢ, q)
     T[i,i] = αᵢ
     Tᴴ[i,i] = conj(αᵢ)
-    @kaxpy!(n, -     αᵢ , vᵢ, q)
-    @kaxpy!(m, -conj(αᵢ), uᵢ, p)
-    βᵢ₊₁ = @knrm2(n, q)
-    γᵢ₊₁ = @knrm2(m, p)
+    @kaxpy!(m, -     αᵢ , vᵢ, q)
+    @kaxpy!(n, -conj(αᵢ), uᵢ, p)
+    βᵢ₊₁ = @knrm2(m, q)
+    γᵢ₊₁ = @knrm2(n, p)
     vᵢ₊₁ .= q ./ βᵢ₊₁
     uᵢ₊₁ .= p ./ γᵢ₊₁
     T[i+1,i] = βᵢ₊₁
@@ -405,7 +405,7 @@ end
 * A. Montoison and D. Orban, [*GPMR: An Iterative Method for Unsymmetric Partitioned Linear Systems*](https://dx.doi.org/10.13140/RG.2.2.24069.68326), Cahier du GERAD G-2021-62, GERAD, Montréal, 2021.
 """
 function montoison_orban(A, B, b::AbstractVector{FC}, c::AbstractVector{FC}, k::Int) where FC <: FloatOrComplex
-  n, m = size(A)
+  m, n = size(A)
   S = ktypeof(b)
   M = vector_to_matrix(S)
 
@@ -424,8 +424,8 @@ function montoison_orban(A, B, b::AbstractVector{FC}, c::AbstractVector{FC}, k::
     end
   end
 
-  V = M(undef, n, k+1)
-  U = M(undef, m, k+1)
+  V = M(undef, m, k+1)
+  U = M(undef, n, k+1)
   H = SparseMatrixCSC(k+1, k, colptr, rowval, nzval_H)
   F = SparseMatrixCSC(k+1, k, colptr, rowval, nzval_F)
 
@@ -435,8 +435,8 @@ function montoison_orban(A, B, b::AbstractVector{FC}, c::AbstractVector{FC}, k::
     vᵢ₊₁ = q = view(V,:,i+1)
     uᵢ₊₁ = p = view(U,:,i+1)
     if i == 1
-      β = @knrm2(n, b)
-      γ = @knrm2(m, c)
+      β = @knrm2(m, b)
+      γ = @knrm2(n, c)
       vᵢ .= b ./ β
       uᵢ .= c ./ γ
     end
@@ -445,14 +445,14 @@ function montoison_orban(A, B, b::AbstractVector{FC}, c::AbstractVector{FC}, k::
     for j = 1:i
       vⱼ = view(V,:,j)
       uⱼ = view(U,:,j)
-      H[j,i] = @kdot(n, vⱼ, q)
+      H[j,i] = @kdot(m, vⱼ, q)
       @kaxpy!(n, -H[j,i], vⱼ, q)
-      F[j,i] = @kdot(m, uⱼ, p)
+      F[j,i] = @kdot(n, uⱼ, p)
       @kaxpy!(m, -F[j,i], uⱼ, p)
     end
-    H[i+1,i] = @knrm2(n, q)
+    H[i+1,i] = @knrm2(m, q)
     vᵢ₊₁ .= q ./ H[i+1,i]
-    F[i+1,i] = @knrm2(m, p)
+    F[i+1,i] = @knrm2(n, p)
     uᵢ₊₁ .= p ./ F[i+1,i]
   end
   return V, H, U, F