PreallocatedLinearOperator docs

.travis.yml

@@ -28,7 +28,7 @@ branches:
 jobs:
   include:
     - stage: Documentation
-      julia: 1.0
+      julia: 1.1
       os: linux
       script:
          - julia --project=docs -e 'using Pkg; Pkg.instantiate(); Pkg.add(PackageSpec(path=pwd()))'

docs/Project.toml

@@ -4,4 +4,4 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
 
 [compat]
-Documenter = "~0.20"
+Documenter = "~0.21"

docs/make.jl

@@ -6,7 +6,7 @@ makedocs(
   doctest = true,
   strict = true,
   assets = ["assets/style.css"],
-  format = :html,
+  format = Documenter.HTML(),
   sitename = "LinearOperators.jl",
   pages = Any["Home" => "index.md",
               "Tutorial" => "tutorial.md",

docs/src/index.md

@@ -17,24 +17,25 @@ Pkg.add("LinearOperators")
 
 ## Operators Available
 
-Operator               | Description
------------------------|------------
-`LinearOperator`       | Base class. Useful to define operators from functions
-`opEye`                | Identity operator
-`opOnes`               | All ones operator
-`opZeros`              | All zeros operator
-`opDiagonal`           | Square (equivalent to `diagm()`) or rectangular diagonal operator
-`opInverse`            | Equivalent to `\`
-`opCholesky`           | More efficient than `opInverse` for symmetric positive definite matrices
-`opLDL`                | Similar to `opCholesky`, for general sparse symmetric matrices
-`opHouseholder`        | Apply a Householder transformation `I-2hh'`
-`opHermitian`          | Represent a symmetric/hermitian operator based on the diagonal and strict lower triangle
-`opRestriction`        | Represent a selection of "rows" when composed on the left with an existing operator
-`opExtension`          | Represent a selection of "columns" when composed on the right with an existing operator
-`LBFGSOperator`        | Limited-memory BFGS approximation in operator form (damped or not)
-`InverseLBFGSOperator` | Inverse of a limited-memory BFGS approximation in operator form (damped or not)
-`LSR1Operator`         | Limited-memory SR1 approximation in operator form
-`kron`                 | Kronecker tensor product in linear operator form
+Operator                     | Description
+-----------------------------|------------
+`LinearOperator`             | Base class. Useful to define operators from functions
+`PreallocatedLinearOperator` | Define operators with preallocation for efficient use of memory
+`opEye`                      | Identity operator
+`opOnes`                     | All ones operator
+`opZeros`                    | All zeros operator
+`opDiagonal`                 | Square (equivalent to `diagm()`) or rectangular diagonal operator
+`opInverse`                  | Equivalent to `\`
+`opCholesky`                 | More efficient than `opInverse` for symmetric positive definite matrices
+`opLDL`                      | Similar to `opCholesky`, for general sparse symmetric matrices
+`opHouseholder`              | Apply a Householder transformation `I-2hh'`
+`opHermitian`                | Represent a symmetric/hermitian operator based on the diagonal and strict lower triangle
+`opRestriction`              | Represent a selection of "rows" when composed on the left with an existing operator
+`opExtension`                | Represent a selection of "columns" when composed on the right with an existing operator
+`LBFGSOperator`              | Limited-memory BFGS approximation in operator form (damped or not)
+`InverseLBFGSOperator`       | Inverse of a limited-memory BFGS approximation in operator form (damped or not)
+`LSR1Operator`               | Limited-memory SR1 approximation in operator form
+`kron`                       | Kronecker tensor product in linear operator form
 
 ## Utility Functions
 

docs/src/reference.md

@@ -2,6 +2,7 @@
 
 ```@docs
 LinearOperator
+PreallocatedLinearOperator
 opEye
 opOnes
 opZeros

docs/src/tutorial.md

@@ -82,48 +82,10 @@ A = [im 1.0; 0.0 1.0]
 op = LinearOperator(Float64, 2, 2, false, false,
                      v -> A * v, u -> transpose(A) * u, w -> A' * w)
 Matrix(op) # Tries to create Float64 matrix with contents of A
-
 # output
-
-ERROR: InexactError: Float64(Float64, 0.0 + 1.0im)
+ERROR: InexactError: Float64(0.0 + 1.0im)
 [...]
 ```
-
-### Preallocation
-
-**Preallocation is being updated. Watch the repo for updates.**
-
-Currently, operators created from matrices allocate internal memory for storage of
-matrix vector products.
-```@example ex1
-A = rand(5, 3)
-op = LinearOperator(A)
-v = rand(3)
-Av = op * v
-w = rand(3)
-al = @allocated op * w
-println("al = $al")
-Aw = op * w
-Aw === Av
-```
-it is possible to recreate this functionality by preallocating the required elements and
-using appropriate inplace functions.  For instance:
-```@example ex1
-function Aprod!(Av, v)
-  n = length(v)
-  Av[1] = 4v[1] - v[2]
-  for i = 2:n-1
-    Av[i] = -v[i-1] + 4*v[i] - v[i+1]
-  end
-  Av[n] = -v[n-1] + 4*v[n]
-  return Av
-end
-Av = Array{Float64}(undef, 100)
-op = LinearOperator(100, 100, true, true, v -> Aprod!(Av, v))
-w = op * ones(100)
-Av === w
-```
-
 ## Limited memory BFGS and SR1
 
 Two other useful operators are the Limited-Memory BFGS in forward and inverse form.
@@ -191,3 +153,47 @@ opA[:,1] * ones(1)
 ```@example ex1
 opA[1,1] * ones(1)
 ```
+
+## Preallocated Operators
+
+Operators created from matrices are very practical, however, it is often useful to reuse
+the memory used by the operator. For that use, we can use
+`PreallocatedLinearOperator(A)` to create an operator that reuses the memory.
+
+```@example ex2
+using LinearOperators # hide
+m, n = 50, 30
+A = rand(m, n) + im * rand(m, n)
+op1 = PreallocatedLinearOperator(A)
+op2 = LinearOperator(A)
+v = rand(n)
+al = @allocated op1 * v
+println("Allocation of PreallocatedLinearOperator product = $al")
+v = rand(n)
+al = @allocated op2 * v
+println("Allocation of LinearOperator product = $al")
+```
+Notice the memory overwrite:
+```@example ex2
+Av = op1 * v
+w = rand(n)
+Aw = op1 * w
+Aw === Av
+```
+which doesn't happen on `LinearOperator`.
+```@example ex2
+Av = op2 * v
+w = rand(n)
+Aw = op2 * w
+Aw === Av
+```
+
+You can also provide the memory to be used.
+```@example ex2
+Mv  = Array{ComplexF64}(undef, m)
+Mtu = Array{ComplexF64}(undef, n)
+Maw = Array{ComplexF64}(undef, n)
+op  = PreallocatedLinearOperator(Mv, Mtu, Maw, A)
+v, u, w = rand(n), rand(m), rand(m)
+(Mv === op * v, Mtu === transpose(op) * u, Maw === adjoint(op) * w)
+```