Use opEye() as a default preconditioner

REQUIRE

@@ -1,2 +1,2 @@
 julia 0.7 2.0
-LinearOperators 0.5.1
+LinearOperators 0.5.2

src/cg.jl

@@ -19,7 +19,7 @@ A preconditioner M may be provided in the form of a linear operator and is
 assumed to be symmetric and positive definite.
 """
 function cg(A :: AbstractLinearOperator, b :: AbstractVector{T};
-            M :: AbstractLinearOperator=opEye(size(A,1)), atol :: Float64=1.0e-8,
+            M :: AbstractLinearOperator=opEye(), atol :: Float64=1.0e-8,
             rtol :: Float64=1.0e-6, itmax :: Int=0, radius :: Float64=0.0,
             verbose :: Bool=false) where T <: Number
 

src/cgls.jl

@@ -38,7 +38,7 @@ It is formally equivalent to LSQR, though can be slightly less accurate,
 but simpler to implement.
 """
 function cgls(A :: AbstractLinearOperator, b :: AbstractVector{T};
-              M :: AbstractLinearOperator=opEye(size(b,1)), λ :: Float64=0.0,
+              M :: AbstractLinearOperator=opEye(), λ :: Float64=0.0,
               atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, radius :: Float64=0.0,
               itmax :: Int=0, verbose :: Bool=false) where T <: Number
 

src/cgne.jl

@@ -54,7 +54,7 @@ but simpler to implement. Only the x-part of the solution is returned.
 A preconditioner M may be provided in the form of a linear operator.
 """
 function cgne(A :: AbstractLinearOperator, b :: AbstractVector{T};
-              M :: AbstractLinearOperator=opEye(size(A,1)), λ :: Float64=0.0,
+              M :: AbstractLinearOperator=opEye(), λ :: Float64=0.0,
               atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, itmax :: Int=0,
               verbose :: Bool=false) where T <: Number
 

src/cgs.jl

@@ -33,7 +33,7 @@ TFQMR and BICGSTAB were developed to remedy this difficulty.»
 This implementation allows a right preconditioner M.
 """
 function cgs(A :: AbstractLinearOperator, b :: AbstractVector{T};
-             M :: AbstractLinearOperator=opEye(size(A,1)),
+             M :: AbstractLinearOperator=opEye(),
              atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6,
              itmax :: Int=0, verbose :: Bool=false) where {T <: Number}
 

src/cr.jl

@@ -15,7 +15,7 @@ assumed to be symmetric and positive definite.
 In a linesearch context, 'linesearch' must be set to 'true'.
 """
 function cr(A :: AbstractLinearOperator, b :: AbstractVector{T};
-            M :: AbstractLinearOperator=opEye(size(A,1)), atol :: Float64=1.0e-8,
+            M :: AbstractLinearOperator=opEye(), atol :: Float64=1.0e-8,
             rtol :: Float64=1.0e-6, γ :: Float64=1.0e-6, itmax :: Int=0,
             radius :: Float64=0.0, verbose :: Bool=false, linesearch :: Bool=false) where T <: Number
 
@@ -29,7 +29,7 @@ function cr(A :: AbstractLinearOperator, b :: AbstractVector{T};
   # Initial state.
   x = zeros(T, n) # initial estimation x = 0
   xNorm = 0.0
-  r = M * b # initial residual r = M * (b - Ax) = M * b
+  r = copy(M * b) # initial residual r = M * (b - Ax) = M * b
   Ar = A * r
   ρ = @kdot(n, r, Ar)
   ρ == 0.0 && return (x, SimpleStats(true, false, [0.0], [], "x = 0 is a zero-residual solution"))

src/crls.jl

@@ -37,7 +37,7 @@ It is formally equivalent to LSMR, though can be slightly less accurate,
 but simpler to implement.
 """
 function crls(A :: AbstractLinearOperator, b :: AbstractVector{T};
-              M :: AbstractLinearOperator=opEye(size(b,1)), λ :: Float64=0.0,
+              M :: AbstractLinearOperator=opEye(), λ :: Float64=0.0,
               atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, radius :: Float64=0.0,
               itmax :: Int=0, verbose :: Bool=false) where T <: Number
 

src/crmr.jl

@@ -53,7 +53,7 @@ but simpler to implement. Only the x-part of the solution is returned.
 A preconditioner M may be provided in the form of a linear operator.
 """
 function crmr(A :: AbstractLinearOperator, b :: AbstractVector{T};
-              M :: AbstractLinearOperator=opEye(size(A,1)), λ :: Float64=0.0,
+              M :: AbstractLinearOperator=opEye(), λ :: Float64=0.0,
               atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6,
               itmax :: Int=0, verbose :: Bool=false) where T <: Number
 

src/diom.jl

@@ -23,7 +23,7 @@ and indefinite systems of linear equations can be handled by this single algorit
 This implementation allows a right preconditioner M.
 """
 function diom(A :: AbstractLinearOperator, b :: AbstractVector{T};
-              M :: AbstractLinearOperator=opEye(size(A,1)),
+              M :: AbstractLinearOperator=opEye(),
               atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6,
               itmax :: Int=0, memory :: Int=20, verbose :: Bool=false) where T <: Number
 

src/dqgmres.jl

@@ -21,7 +21,7 @@ and computes a sequence of approximate solutions with the quasi-minimal residual
 This implementation allows a right preconditioner M.
 """
 function dqgmres(A :: AbstractLinearOperator, b :: AbstractVector{T};
-                 M :: AbstractLinearOperator=opEye(size(A,1)), atol :: Float64=1.0e-8,
+                 M :: AbstractLinearOperator=opEye(), atol :: Float64=1.0e-8,
                  rtol :: Float64=1.0e-6, itmax :: Int=0, memory :: Int=20,
                  verbose :: Bool=false) where T <: Number
 

src/lslq.jl

@@ -93,8 +93,8 @@ The iterations stop as soon as one of the following conditions holds true:
 * R. Estrin, D. Orban and M. A. Saunders, *LSLQ: An Iterative Method for Linear Least-Squares with an Error Minimization Property*, Cahier du GERAD G-2017-xx, GERAD, Montreal, 2017.
 """
 function lslq(A :: AbstractLinearOperator, b :: AbstractVector{T};
-              M :: AbstractLinearOperator=opEye(size(A,1)),
-              N :: AbstractLinearOperator=opEye(size(A,2)),
+              M :: AbstractLinearOperator=opEye(),
+              N :: AbstractLinearOperator=opEye(),
               sqd :: Bool=false, λ :: Float64=0.0, σ :: Float64=0.0,
               atol :: Float64=1.0e-8, btol :: Float64=1.0e-8, etol :: Float64=1.0e-8,
               window :: Int=5, utol :: Float64=1.0e-8, itmax :: Int=0,
@@ -104,6 +104,10 @@ function lslq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   size(b, 1) == m || error("Inconsistent problem size")
   verbose && @printf("LSLQ: system of %d equations in %d variables\n", m, n)
 
+  # Tests M == Iₙ and N == Iₘ
+  MisI = isa(M, opEye)
+  NisI = isa(N, opEye)
+
   # If solving an SQD system, set regularization to 1.
   sqd && (λ = 1.0)
   λ² = λ * λ
@@ -125,7 +129,7 @@ function lslq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   β = β₁
 
   @kscal!(m, 1.0/β₁, u)
-  @kscal!(m, 1.0/β₁, Mu)
+  MisI || @kscal!(m, 1.0/β₁, Mu)
   Nv = copy(A' * u)
   v = N * Nv
   α = sqrt(@kdot(n, v, Nv))  # = α₁
@@ -134,7 +138,7 @@ function lslq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   α == 0.0 && return (x_lq, x_cg, err_lbnds, err_ubnds_lq, err_ubnds_cg,
                       SimpleStats(true, false, [β₁], [0.0], "x = 0 is a minimum least-squares solution"))
   @kscal!(n, 1.0/α, v)
-  @kscal!(n, 1.0/α, Nv)
+  NisI || @kscal!(n, 1.0/α, Nv)
 
   Anorm = α
   Anorm² = α * α
@@ -200,7 +204,7 @@ function lslq(A :: AbstractLinearOperator, b :: AbstractVector{T};
     β = sqrt(@kdot(m, u, Mu))
     if β != 0.0
       @kscal!(m, 1.0/β, u)
-      @kscal!(m, 1.0/β, Mu)
+      MisI || @kscal!(m, 1.0/β, Mu)
 
       # 2. αv = A'u - βv
       @kscal!(n, -β, Nv)
@@ -209,7 +213,7 @@ function lslq(A :: AbstractLinearOperator, b :: AbstractVector{T};
       α = sqrt(@kdot(n, v, Nv))
       if α != 0.0
         @kscal!(n, 1.0/α, v)
-        @kscal!(n, 1.0/α, Nv)
+        NisI || @kscal!(n, 1.0/α, Nv)
       end
 
       # rotate out regularization term if present

src/lsmr.jl

@@ -58,8 +58,8 @@ In this case, `N` can still be specified and indicates the norm
 in which `x` should be measured.
 """
 function lsmr(A :: AbstractLinearOperator, b :: AbstractVector{T};
-              M :: AbstractLinearOperator=opEye(size(A,1)),
-              N :: AbstractLinearOperator=opEye(size(A,2)),
+              M :: AbstractLinearOperator=opEye(),
+              N :: AbstractLinearOperator=opEye(),
               sqd :: Bool=false,
               λ :: Float64=0.0, axtol :: Float64=1.0e-8, btol :: Float64=1.0e-8,
               atol :: Float64=0.0, rtol :: Float64=0.0,
@@ -71,6 +71,10 @@ function lsmr(A :: AbstractLinearOperator, b :: AbstractVector{T};
   size(b, 1) == m || error("Inconsistent problem size")
   verbose && @printf("LSMR: system of %d equations in %d variables\n", m, n)
 
+  # Tests M == Iₙ and N == Iₘ
+  MisI = isa(M, opEye)
+  NisI = isa(N, opEye)
+
   # If solving an SQD system, set regularization to 1.
   sqd && (λ = 1.0)
   ctol = conlim > 0.0 ? 1/conlim : 0.0
@@ -85,7 +89,7 @@ function lsmr(A :: AbstractLinearOperator, b :: AbstractVector{T};
   β = β₁
 
   @kscal!(m, 1.0/β₁, u)
-  @kscal!(m, 1.0/β₁, Mu)
+  MisI || @kscal!(m, 1.0/β₁, Mu)
   Nv = copy(A' * u)
   v = N * Nv
   α = sqrt(@kdot(n, v, Nv))
@@ -130,7 +134,7 @@ function lsmr(A :: AbstractLinearOperator, b :: AbstractVector{T};
   # A'b = 0 so x = 0 is a minimum least-squares solution
   α == 0.0 && return (x, SimpleStats(true, false, [β₁], [0.0], "x = 0 is a minimum least-squares solution"))
   @kscal!(n, 1.0/α, v)
-  @kscal!(n, 1.0/α, Nv)
+  NisI || @kscal!(n, 1.0/α, Nv)
 
   h = copy(v)
   hbar = zeros(T, n)
@@ -151,13 +155,13 @@ function lsmr(A :: AbstractLinearOperator, b :: AbstractVector{T};
 
     # Generate next Golub-Kahan vectors.
     # 1. βu = Av - αu
-    @kscal!(m, -α, Mu,)
+    @kscal!(m, -α, Mu)
     @kaxpy!(m, 1.0, A * v, Mu)
     u = M * Mu
     β = sqrt(@kdot(m, u, Mu))
     if β != 0.0
       @kscal!(m, 1.0/β, u)
-      @kscal!(m, 1.0/β, Mu)
+      MisI || @kscal!(m, 1.0/β, Mu)
 
       # 2. αv = A'u - βv
       @kscal!(n, -β, Nv)
@@ -166,7 +170,7 @@ function lsmr(A :: AbstractLinearOperator, b :: AbstractVector{T};
       α = sqrt(@kdot(n, v, Nv))
       if α != 0.0
         @kscal!(n, 1.0/α, v)
-        @kscal!(n, 1.0/α, Nv)
+        NisI || @kscal!(n, 1.0/α, Nv)
       end
     end
 

src/lsqr.jl

@@ -58,8 +58,8 @@ In this case, `N` can still be specified and indicates the norm
 in which `x` should be measured.
 """
 function lsqr(A :: AbstractLinearOperator, b :: AbstractVector{T};
-              M :: AbstractLinearOperator=opEye(size(A,1)),
-              N :: AbstractLinearOperator=opEye(size(A,2)),
+              M :: AbstractLinearOperator=opEye(),
+              N :: AbstractLinearOperator=opEye(),
               sqd :: Bool=false,
               λ :: Float64=0.0, axtol :: Float64=1.0e-8, btol :: Float64=1.0e-8,
               atol :: Float64=0.0, rtol :: Float64=0.0,
@@ -71,6 +71,10 @@ function lsqr(A :: AbstractLinearOperator, b :: AbstractVector{T};
   size(b, 1) == m || error("Inconsistent problem size")
   verbose && @printf("LSQR: system of %d equations in %d variables\n", m, n)
 
+  # Tests M == Iₙ and N == Iₘ
+  MisI = isa(M, opEye)
+  NisI = isa(N, opEye)
+
   # If solving an SQD system, set regularization to 1.
   sqd && (λ = 1.0)
   λ² = λ * λ
@@ -86,7 +90,7 @@ function lsqr(A :: AbstractLinearOperator, b :: AbstractVector{T};
   β = β₁
 
   @kscal!(m, 1.0/β₁, u)
-  @kscal!(m, 1.0/β₁, Mu)
+  MisI || @kscal!(m, 1.0/β₁, Mu)
   Nv = copy(A' * u)
   v = N * Nv
   α = sqrt(@kdot(n, v, Nv))
@@ -113,7 +117,7 @@ function lsqr(A :: AbstractLinearOperator, b :: AbstractVector{T};
   # A'b = 0 so x = 0 is a minimum least-squares solution
   α == 0.0 && return (x, SimpleStats(true, false, [β₁], [0.0], "x = 0 is a minimum least-squares solution"))
   @kscal!(n, 1.0/α, v)
-  @kscal!(n, 1.0/α, Nv)
+  NisI || @kscal!(n, 1.0/α, Nv)
   w = copy(v)
 
   # Initialize other constants.
@@ -154,7 +158,7 @@ function lsqr(A :: AbstractLinearOperator, b :: AbstractVector{T};
     β = sqrt(@kdot(m, u, Mu))
     if β != 0.0
       @kscal!(m, 1.0/β, u)
-      @kscal!(m, 1.0/β, Mu)
+      MisI || @kscal!(m, 1.0/β, Mu)
       Anorm² = Anorm² + α * α + β * β;  # = ‖B_{k-1}‖²
       λ > 0.0 && (Anorm² += λ²)
 
@@ -165,7 +169,7 @@ function lsqr(A :: AbstractLinearOperator, b :: AbstractVector{T};
       α = sqrt(@kdot(n, v, Nv))
       if α != 0.0
         @kscal!(n, 1.0/α, v)
-        @kscal!(n, 1.0/α, Nv)
+        NisI || @kscal!(n, 1.0/α, Nv)
       end
     end
 

src/minres.jl

@@ -43,7 +43,7 @@ A preconditioner M may be provided in the form of a linear operator and is
 assumed to be symmetric and positive definite.
 """
 function minres(A :: AbstractLinearOperator, b :: AbstractVector{T};
-                M :: AbstractLinearOperator=opEye(size(A,1)), λ :: Float64=0.0,
+                M :: AbstractLinearOperator=opEye(), λ :: Float64=0.0,
                 atol :: Float64=1.0e-12, rtol :: Float64=1.0e-12, etol :: Float64=1.0e-8,
                 window :: Int=5, itmax :: Int=0, conlim :: Float64=1.0e+8,
                 verbose :: Bool=false) where T <: Number

src/symmlq.jl

@@ -25,7 +25,7 @@ A preconditioner M may be provided in the form of a linear operator and is
 assumed to be symmetric and positive definite.
 """
 function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
-                M :: AbstractLinearOperator=opEye(size(A,1)), λ :: Float64=0.0,
+                M :: AbstractLinearOperator=opEye(), λ :: Float64=0.0,
                 λest :: Float64=0.0, atol :: Float64=1.0e-8, rtol :: Float64=1.0e-8,
                 etol :: Float64=1.0e-8, window :: Int=0, itmax :: Int=0,
                 conlim :: Float64=1.0e+8, verbose :: Bool=false) where T <: Number
@@ -35,6 +35,9 @@ function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   size(b, 1) == m || error("Inconsistent problem size")
   verbose && @printf("SYMMLQ: system of size %d\n", n)
 
+  # Test M == Iₘ
+  MisI = isa(M, opEye)
+
   ϵM = eps(T)
   x = zeros(T, n)
   ctol = conlim > 0.0 ? 1 ./ conlim : 0.0;
@@ -47,8 +50,8 @@ function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   β₁ == 0.0 && return (x, x, SimpleStats(true, true, [0.0], [0.0], "x = 0 is a zero-residual solution"))
   β₁ = sqrt(β₁)
   β = β₁
-  vold = @kscal!(m, 1 ./ β, vold)
-  p = @kscal!(m, 1 ./ β, p)
+  @kscal!(m, 1 ./ β, vold)
+  MisI || @kscal!(m, 1 ./ β, p)
 
   v = copy(A * p)
   α = @kdot(m, p, v) + λ
@@ -58,7 +61,7 @@ function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   β < 0.0 && error("Preconditioner is not positive definite")
   β = sqrt(β)
   @kscal!(m, 1 ./ β, v)
-  @kscal!(m, 1 ./ β, p)
+  MisI || @kscal!(m, 1 ./ β, p)
 
   # Start QR factorization
   γbar = α
@@ -143,7 +146,7 @@ function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
     β < 0.0 && error("Preconditioner is not positive definite")
     β = sqrt(β)
     @kscal!(m, 1 ./ β, v_next)
-    @kscal!(m, 1 ./ β, p)
+    MisI || @kscal!(m, 1 ./ β, p)
 
     # Continue A norm estimate
     ANorm² = ANorm² + α * α + oldβ * oldβ + β * β