Merge pull request #171 from haampie/fix-maxiter

Make maxiter the size of the matrix

Author: haampie

src/bicgstabl.jl

@@ -28,7 +28,7 @@ bicgstabl_iterator(A, b, l; kwargs...) = bicgstabl_iterator!(zerox(A, b), A, b,
 
 function bicgstabl_iterator!(x, A, b, l::Int = 2;
     Pl = Identity(),
-    max_mv_products = min(30, size(A, 1)),
+    max_mv_products = size(A, 2),
     initial_zero = false,
     tol = sqrt(eps(real(eltype(b))))
 )
@@ -161,7 +161,7 @@ bicgstabl(A, b, l = 2; kwargs...) = bicgstabl!(zerox(A, b), A, b, l; initial_zer
 
 ## Keywords
 
-- `max_mv_products::Int = min(30, size(A, 1))`: maximum number of matrix vector products.
+- `max_mv_products::Int = size(A, 2)`: maximum number of matrix vector products.
 For BiCGStab(l) this is a less dubious term than "number of iterations";
 - `Pl = Identity()`: left preconditioner of the method;
 - `tol::Real = sqrt(eps(real(eltype(b))))`: tolerance for stopping condition `|r_k| / |r_0| ≤ tol`. 
@@ -182,7 +182,7 @@ For BiCGStab(l) this is a less dubious term than "number of iterations";
 """
 function bicgstabl!(x, A, b, l = 2;
     tol = sqrt(eps(real(eltype(b)))),
-    max_mv_products::Int = min(20, size(A, 1)),
+    max_mv_products::Int = size(A, 2),
     log::Bool = false,
     verbose::Bool = false,
     Pl = Identity(),

src/cg.jl

@@ -96,7 +96,7 @@ end
 
 function cg_iterator!(x, A, b, Pl = Identity();
     tol = sqrt(eps(real(eltype(b)))),
-    maxiter = min(20, length(b)),
+    maxiter::Int = size(A, 2),
     initially_zero::Bool = false
 )
     u = zeros(x)
@@ -155,7 +155,7 @@ cg(A, b; kwargs...) = cg!(zerox(A, b), A, b; initially_zero = true, kwargs...)
 - `Pl = Identity()`: left preconditioner of the method. Should be symmetric, 
   positive-definite like `A`.
 - `tol::Real = sqrt(eps(real(eltype(b))))`: tolerance for stopping condition `|r_k| / |r_0| ≤ tol`;
-- `maxiter::Integer = size(A,2)`: maximum number of iterations;
+- `maxiter::Int = size(A,2)`: maximum number of iterations;
 - `verbose::Bool = false`: print method information;
 - `log::Bool = false`: keep track of the residual norm in each iteration;
 
@@ -177,7 +177,7 @@ cg(A, b; kwargs...) = cg!(zerox(A, b), A, b; initially_zero = true, kwargs...)
 """
 function cg!(x, A, b;
     tol = sqrt(eps(real(eltype(b)))),
-    maxiter::Integer = min(20, size(A, 1)),
+    maxiter::Int = size(A, 2),
     log::Bool = false,
     verbose::Bool = false,
     Pl = Identity(),

src/chebyshev.jl

@@ -58,7 +58,7 @@ chebyshev_iterable(A, b, λmin::Real, λmax::Real; kwargs...) =
     chebyshev_iterable!(zerox(A, b), A, b, λmin, λmax; kwargs...)
 
 function chebyshev_iterable!(x, A, b, λmin::Real, λmax::Real;
-    tol = sqrt(eps(real(eltype(b)))), maxiter = size(A, 1), Pl = Identity(), initially_zero = false)
+    tol = sqrt(eps(real(eltype(b)))), maxiter = size(A, 2), Pl = Identity(), initially_zero = false)
 
     λ_avg = (λmax + λmin) / 2
     λ_diff = (λmax - λmin) / 2
@@ -117,7 +117,7 @@ Solve Ax = b for symmetric, definite matrices A using Chebyshev iteration.
   matrix-vector product can be saved when computing the initial 
   residual vector;
 - `tol`: tolerance for stopping condition `|r_k| / |r_0| ≤ tol`.
-- `maxiter::Int`: maximum number of inner iterations of GMRES;
+- `maxiter::Int = size(A, 2)`: maximum number of inner iterations of GMRES;
 - `Pl = Identity()`: left preconditioner;
 - `log::Bool = false`: keep track of the residual norm in each iteration;
 - `verbose::Bool = false`: print convergence information during the iterations.
@@ -136,7 +136,7 @@ Solve Ax = b for symmetric, definite matrices A using Chebyshev iteration.
 function chebyshev!(x, A, b, λmin::Real, λmax::Real;
     Pl = Identity(),
     tol::Real=sqrt(eps(real(eltype(b)))),
-    maxiter::Int=size(A, 1),
+    maxiter::Int=size(A, 2),
     log::Bool=false,
     verbose::Bool=false,
     initially_zero::Bool=false

src/gmres.jl

@@ -104,8 +104,8 @@ function gmres_iterable!(x, A, b;
     Pl = Identity(),
     Pr = Identity(),
     tol = sqrt(eps(real(eltype(b)))),
-    restart::Int = min(20, length(b)),
-    maxiter::Int = restart,
+    restart::Int = min(20, size(A, 2)),
+    maxiter::Int = size(A, 2),
     initially_zero::Bool = false
 )
     T = eltype(x)
@@ -152,8 +152,8 @@ Solves the problem ``Ax = b`` with restarted GMRES.
   matrix-vector product can be saved when computing the initial 
   residual vector;
 - `tol`: relative tolerance;
-- `restart::Int`: restarts GMRES after specified number of iterations;
-- `maxiter::Int`: maximum number of inner iterations of GMRES;
+- `restart::Int = min(20, size(A, 2))`: restarts GMRES after specified number of iterations;
+- `maxiter::Int = size(A, 2)`: maximum number of inner iterations of GMRES;
 - `Pl`: left preconditioner;
 - `Pr`: right preconditioner;
 - `log::Bool`: keep track of the residual norm in each iteration;
@@ -174,8 +174,8 @@ function gmres!(x, A, b;
   Pl = Identity(),
   Pr = Identity(),
   tol = sqrt(eps(real(eltype(b)))),
-  restart::Int = min(20, length(b)),
-  maxiter::Int = restart,
+  restart::Int = min(20, size(A, 2)),
+  maxiter::Int = size(A, 2),
   log::Bool = false,
   initially_zero::Bool = false,
   verbose::Bool = false

src/idrs.jl

@@ -23,7 +23,7 @@ shadow space.
 
 - `s::Integer = 8`: dimension of the shadow space;
 - `tol`: relative tolerance;
-- `maxiter::Int`: maximum number of iterations;
+- `maxiter::Int = size(A, 2)`: maximum number of iterations;
 - `log::Bool`: keep track of the residual norm in each iteration;
 - `verbose::Bool`: print convergence information during the iterations.
 
@@ -39,7 +39,7 @@ shadow space.
 - `history`: convergence history.
 """
 function idrs!(x, A, b;
-    s = 8, tol=sqrt(eps(real(eltype(b)))), maxiter=length(x)^2,
+    s = 8, tol=sqrt(eps(real(eltype(b)))), maxiter=size(A, 2),
     log::Bool=false, kwargs...
     )
     history = ConvergenceHistory(partial=!log)

src/lsmr.jl

@@ -39,7 +39,7 @@ especially if ``A`` is ill-conditioned.
   Maximum precision can be obtained by setting
 - `atol` = `btol` = `conlim` = zero, but the number of iterations
   may then be excessive.
-- `maxiter::Integer = min(20,length(b))`: maximum number of iterations.
+- `maxiter::Int = maximum(size(A))`: maximum number of iterations.
 - `log::Bool`: keep track of the residual norm in each iteration;
 - `verbose::Bool`: print convergence information during the iterations.
 
@@ -65,7 +65,7 @@ especially if ``A`` is ill-conditioned.
 - `:resnom` => `::Vector`: residual norm at each iteration.
 """
 function lsmr!(x, A, b;
-    maxiter::Integer = max(size(A,1), size(A,2)),
+    maxiter::Int = maximum(size(A)),
     log::Bool=false, kwargs...
     )
     history = ConvergenceHistory(partial=!log)
@@ -87,7 +87,7 @@ end
 
 function lsmr_method!(log::ConvergenceHistory, x, A, b, v, h, hbar;
     atol::Number = 1e-6, btol::Number = 1e-6, conlim::Number = 1e8,
-    maxiter::Integer = max(size(A,1), size(A,2)), λ::Number = 0,
+    maxiter::Int = maximum(size(A)), λ::Number = 0,
     verbose::Bool=false
     )
     verbose && @printf("=== lsmr ===\n%4s\t%7s\t\t%7s\t\t%7s\n","iter","anorm","cnorm","rnorm")

src/lsqr.jl

@@ -37,7 +37,7 @@ especially if A is ill-conditioned.
   Maximum precision can be obtained by setting
   `atol` = `btol` = `conlim` = zero, but the number of iterations
   may then be excessive.
-- `maxiter::Integer = min(20,length(b))`: maximum number of iterations.
+- `maxiter::Int = maximum(size(A))`: maximum number of iterations.
 - `verbose::Bool = false`: print method information.
 - `log::Bool = false`: output an extra element of type `ConvergenceHistory`
   containing extra information of the method execution.
@@ -64,7 +64,7 @@ especially if A is ill-conditioned.
 - `:resnom` => `::Vector`: residual norm at each iteration.
 """
 function lsqr!(x, A, b;
-    maxiter::Int=max(size(A,1), size(A,2)), log::Bool=false,
+    maxiter::Int=maximum(size(A)), log::Bool=false,
     kwargs...
     )
     T = Adivtype(A, b)
@@ -90,7 +90,7 @@ end
 function lsqr_method!(log::ConvergenceHistory, x, A, b;
     damp=0, atol=sqrt(eps(real(Adivtype(A,b)))), btol=sqrt(eps(real(Adivtype(A,b)))),
     conlim=real(one(Adivtype(A,b)))/sqrt(eps(real(Adivtype(A,b)))),
-    maxiter::Int=max(size(A,1), size(A,2)), verbose::Bool=false,
+    maxiter::Int=maximum(size(A)), verbose::Bool=false,
     )
     verbose && @printf("=== lsqr ===\n%4s\t%7s\t\t%7s\t\t%7s\t\t%7s\n","iter","resnorm","anorm","cnorm","rnorm")
     # Sanity-checking

src/minres.jl

@@ -42,7 +42,7 @@ function minres_iterable!(x, A, b;
     initially_zero::Bool = false, 
     skew_hermitian::Bool = false, 
     tol = sqrt(eps(real(eltype(b)))), 
-    maxiter = size(A, 1)
+    maxiter = size(A, 2)
 )
     T = eltype(x)
     HessenbergT = skew_hermitian ? T : real(T)
@@ -173,7 +173,7 @@ Solve Ax = b for (skew-)Hermitian matrices A using MINRES.
   residual vector;
 - `skew_hermitian::Bool = false`: if `true` assumes that `A` is skew-symmetric or skew-Hermitian;
 - `tol`: tolerance for stopping condition `|r_k| / |r_0| ≤ tol`. Note that the residual is computed only approximately;
-- `maxiter::Int`: maximum number of inner iterations of GMRES;
+- `maxiter::Int = size(A, 2)`: maximum number of iterations;
 - `Pl`: left preconditioner;
 - `Pr`: right preconditioner;
 - `log::Bool = false`: keep track of the residual norm in each iteration;
@@ -195,7 +195,7 @@ function minres!(x, A, b;
     verbose::Bool = false,
     log::Bool = false,
     tol = sqrt(eps(real(eltype(b)))),
-    maxiter::Int = min(30, size(A, 1)),
+    maxiter::Int = size(A, 2),
     initially_zero::Bool = false
 )
     history = ConvergenceHistory(partial = !log)

test/idrs.jl

@@ -28,7 +28,7 @@ srand(1234567)
 end
 
 @testset "SparseMatrixCSC{$T}" for T in (Float64, Complex128)
-    A = sprand(T, n, n, 0.5) + I
+    A = sprand(T, n, n, 0.5) + n * I
     b = rand(T, n)
     tol = √eps(real(T))
 

test/lsmr.jl

@@ -45,6 +45,12 @@ end
 
 eltype(A::DampenedMatrix) = promote_type(eltype(A.A), eltype(A.diagonal))
 
+function size(A::DampenedMatrix)
+    m, n = size(A.A)
+    l = length(A.diagonal)
+    (m + l, n)
+end
+
 function size(A::DampenedMatrix, dim::Integer)
     m, n = size(A.A)
     l = length(A.diagonal)