Add documentation on operators (#63)

* Add documentation on operators

docs/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 KrylovPreconditioners = "45d422c2-293f-44ce-8315-2cb988662dec"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 Documenter = "1.0"

docs/make.jl

@@ -1,7 +1,8 @@
-using Documenter, KrylovPreconditioners
+using Documenter, KrylovPreconditioners, SparseArrays
 
 makedocs(
   modules = [KrylovPreconditioners],
+  checkdocs = :exports,
   doctest = true,
   linkcheck = true,
   format = Documenter.HTML(assets = ["assets/style.css"],
@@ -10,6 +11,8 @@ makedocs(
                            collapselevel = 1),
   sitename = "KrylovPreconditioners.jl",
   pages = ["Home" => "index.md",
+           "Krylov operators" => "krylov_operators.md",
+           "Triangular operators" => "triangular_operators.md",
            "Reference" => "reference.md"
           ]
 )

docs/src/krylov_operators.md

@@ -0,0 +1,77 @@
+# [Krylov operators](@id krylov_operators)
+
+```@docs
+KrylovOperator
+update!(::AbstractKrylovOperator, ::Any)
+```
+
+## Nvidia GPUs
+
+Sparse matrices have a specific storage on Nvidia GPUs (`CuSparseMatrixCSC`, `CuSparseMatrixCSR` or `CuSparseMatrixCOO`):
+
+```julia
+using CUDA, CUDA.CUSPARSE
+using SparseArrays
+using KrylovPreconditioners
+
+if CUDA.functional()
+  # CPU Arrays
+  A_cpu = sprand(200, 100, 0.3)
+
+  # GPU Arrays
+  A_csc_gpu = CuSparseMatrixCSC(A_cpu)
+  A_csr_gpu = CuSparseMatrixCSR(A_cpu)
+  A_coo_gpu = CuSparseMatrixCOO(A_cpu)
+
+  # Krylov operators
+  op_csc = KrylovOperator(A_csc_gpu; nrhs=1, transa='N')
+  op_csr = KrylovOperator(A_csr_gpu; nrhs=1, transa='T')
+  op_coo = KrylovOperator(A_coo_gpu; nrhs=5, transa='N')
+end
+```
+
+## AMD GPUs
+
+Sparse matrices have a specific storage on AMD GPUs (`ROCSparseMatrixCSC`, `ROCSparseMatrixCSR` or `ROCSparseMatrixCOO`):
+
+```julia
+using AMDGPU, AMDGPU.rocSPARSE
+using SparseArrays
+using KrylovPreconditioners
+
+if AMDGPU.functional()
+  # CPU Arrays
+  A_cpu = sprand(200, 100, 0.3)
+
+  # GPU Arrays
+  A_csc_gpu = ROCSparseMatrixCSC(A_cpu)
+  A_csr_gpu = ROCSparseMatrixCSR(A_cpu)
+  A_coo_gpu = ROCSparseMatrixCOO(A_cpu)
+
+  # Krylov operators
+  op_csc = KrylovOperator(A_csc_gpu; nrhs=1, transa='N')
+  op_csr = KrylovOperator(A_csr_gpu; nrhs=1, transa='T')
+  op_coo = KrylovOperator(A_coo_gpu; nrhs=5, transa='N')
+end
+```
+
+## Intel GPUs
+
+Sparse matrices have a specific storage on Intel GPUs (`oneSparseMatrixCSR`):
+
+```julia
+using oneAPI, oneAPI.oneMKL
+using SparseArrays
+using KrylovPreconditioners
+
+if oneAPI.functional()
+  # CPU Arrays
+  A_cpu = sprand(Float32, 20, 10, 0.3)
+
+  # GPU Arrays
+  A_csr_gpu = oneSparseMatrixCSR(A_cpu)
+
+  # Krylov operator
+  op_csr = KrylovOperator(A_csr_gpu; nrhs=1, transa='N')
+end
+```

docs/src/reference.md

@@ -5,7 +5,9 @@
 ```@index
 ```
 
-```@autodocs
-Modules = [KrylovPreconditioners]
-Order   = [:function, :type]
+```@docs
+KrylovPreconditioners.update!(::BlockJacobiPreconditioner, ::SparseMatrixCSC)
+KrylovPreconditioners.BlockJacobiPreconditioner
+KrylovPreconditioners.backward_substitution!
+KrylovPreconditioners.forward_substitution!
 ```

docs/src/triangular_operators.md

@@ -0,0 +1,77 @@
+# [Triangular operators](@id triangular_operators)
+
+```@docs
+TriangularOperator
+update!(::AbstractTriangularOperator, ::Any)
+```
+
+## Nvidia GPUs
+
+Sparse matrices have a specific storage on Nvidia GPUs (`CuSparseMatrixCSC`, `CuSparseMatrixCSR` or `CuSparseMatrixCOO`):
+
+```julia
+using CUDA, CUDA.CUSPARSE
+using SparseArrays
+using KrylovPreconditioners
+
+if CUDA.functional()
+  # CPU Arrays
+  A_cpu = sprand(100, 100, 0.3)
+
+  # GPU Arrays
+  A_csc_gpu = CuSparseMatrixCSC(A_cpu)
+  A_csr_gpu = CuSparseMatrixCSR(A_cpu)
+  A_coo_gpu = CuSparseMatrixCOO(A_cpu)
+
+  # Triangular operators
+  op_csc = TriangularOperator(A_csc_gpu; uplo='L', diag='U', nrhs=1, transa='N')
+  op_csr = TriangularOperator(A_csr_gpu; uplo='U', diag='N', nrhs=1, transa='T')
+  op_coo = TriangularOperator(A_coo_gpu; uplo='L', diag='N', nrhs=5, transa='N')
+end
+```
+
+## AMD GPUs
+
+Sparse matrices have a specific storage on AMD GPUs (`ROCSparseMatrixCSC`, `ROCSparseMatrixCSR` or `ROCSparseMatrixCOO`):
+
+```julia
+using AMDGPU, AMDGPU.rocSPARSE
+using SparseArrays
+using KrylovPreconditioners
+
+if AMDGPU.functional()
+  # CPU Arrays
+  A_cpu = sprand(200, 100, 0.3)
+
+  # GPU Arrays
+  A_csc_gpu = ROCSparseMatrixCSC(A_cpu)
+  A_csr_gpu = ROCSparseMatrixCSR(A_cpu)
+  A_coo_gpu = ROCSparseMatrixCOO(A_cpu)
+
+  # Triangular operators
+  op_csc = TriangularOperator(A_csc_gpu; uplo='L', diag='U', nrhs=1, transa='N')
+  op_csr = TriangularOperator(A_csr_gpu; uplo='L', diag='U', nrhs=1, transa='T')
+  op_coo = TriangularOperator(A_coo_gpu; uplo='L', diag='U', nrhs=5, transa='N')
+end
+```
+
+## Intel GPUs
+
+Sparse matrices have a specific storage on Intel GPUs (`oneSparseMatrixCSR`):
+
+```julia
+using oneAPI, oneAPI.oneMKL
+using SparseArrays
+using KrylovPreconditioners
+
+if oneAPI.functional()
+  # CPU Arrays
+  A_cpu = sprand(T, 20, 10, 0.3)
+
+  # GPU Arrays
+  A_csr_gpu = oneSparseMatrixCSR(A_cpu)
+
+  # Triangular operator
+  op_csr = TriangularOperator(A_csr_gpu; uplo='L', diag='U', nrhs=1, transa='N')
+end
+```

src/KrylovPreconditioners.jl

@@ -12,32 +12,9 @@ import LinearAlgebra: ldiv!
 abstract type AbstractKrylovPreconditioner end
 export AbstractKrylovPreconditioner
 
-abstract type AbstractKrylovOperator{T} end
-export AbstractKrylovOperator
-
-abstract type AbstractTriangularOperator{T} end
-export AbstractTriangularOperator
-
-update!(p::AbstractKrylovPreconditioner, A::SparseMatrixCSC) = error("update!() for $(typeof(p)) is not implemented.")
-update!(p::AbstractKrylovPreconditioner, A) = error("update!() for $(typeof(p)) is not implemented.")
-update!(p::AbstractKrylovOperator, A::SparseMatrixCSC) = error("update!() for $(typeof(p)) is not implemented.")
-update!(p::AbstractKrylovOperator, A) = error("update!() for $(typeof(p)) is not implemented.")
-
-export update!, get_timer, reset_timer!
-
-function get_timer(p::AbstractKrylovPreconditioner)
-    return p.timer_update
-end
-
-function reset_timer!(p::AbstractKrylovPreconditioner)
-    p.timer_update = 0.0
-end
-
-function KrylovOperator end
-export KrylovOperator
-
-function TriangularOperator end
-export TriangularOperator
+# Operators
+include("krylov_operators.jl")
+include("triangular_operators.jl")
 
 # Preconditioners
 include("ic0.jl")
@@ -52,4 +29,16 @@ export scaling_csr!
 # Ordering
 # include(ordering.jl)
 
+update!(op::AbstractKrylovPreconditioner, A) = error("update!() for $(typeof(op)) is not implemented.")
+
+export update!, get_timer, reset_timer!
+
+function get_timer(p::AbstractKrylovPreconditioner)
+    return p.timer_update
+end
+
+function reset_timer!(p::AbstractKrylovPreconditioner)
+    p.timer_update = 0.0
+end
+
 end # module KrylovPreconditioners

src/krylov_operators.jl

@@ -0,0 +1,35 @@
+export AbstractKrylovOperator, KrylovOperator
+
+abstract type AbstractKrylovOperator{T} end
+
+"""
+    update!(op::AbstractKrylovOperator, A)
+
+Update the sparse matrix `A` associated with the given `AbstractKrylovOperator` without the need to reallocate buffers 
+or repeat the structural analysis phase for detecting parallelism for sparse matrix-vector or matrix-matrix products.
+`A` and the operator `op` must have the same sparsity pattern, enabling efficient reuse of existing resources.
+
+#### Input arguments
+
+* `op`: The Krylov operator to update;
+* `A`: The new sparse matrix to associate with the operator.
+"""
+update!(op::AbstractKrylovOperator, A) = error("update!() for $(typeof(op)) is not implemented.")
+
+"""
+    KrylovOperator(A; nrhs::Int=1, transa::Char='N')
+
+Create a Krylov operator to accelerate sparse matrix-vector or matrix-matrix products on GPU architectures.
+The operator is compatible with sparse matrices stored on NVIDIA, AMD, and Intel GPUs.
+
+#### Input arguments
+
+* `A`: The sparse matrix on the GPU that serves as the operator for matrix-vector or matrix-matrix products;
+* `nrhs`: Specifies the number of columns for the right-hand sides. Defaults to `1` for standard matrix-vector products;
+* `transa`: Determines how the matrix `A` is applied during the products; `'N'` for no transposition, `'T'` for transpose, and `'C'` for conjugate transpose.
+
+#### Output argument
+
+* `op`: An instance of `AbstractKrylovOperator` representing the Krylov operator for the specified sparse matrix and parameters.
+"""
+function KrylovOperator end

src/triangular_operators.jl

@@ -0,0 +1,37 @@
+export AbstractTriangularOperator, TriangularOperator
+
+abstract type AbstractTriangularOperator{T} end
+
+"""
+    update!(op::AbstractTriangularOperator, A)
+
+Update the sparse matrix `A` associated with the given `AbstractTriangularOperator` without the need to reallocate buffers 
+or repeat the structural analysis phase for detecting parallelism for sparse triangular solves.
+`A` and the operator `op` must have the same sparsity pattern, enabling efficient reuse of existing resources.
+
+#### Input arguments
+
+* `op`: The triangular operator to update;
+* `A`: The new sparse matrix to associate with the operator.
+"""
+update!(op::AbstractTriangularOperator, A) = error("update!() for $(typeof(op)) is not implemented.")
+
+"""
+    TriangularOperator(A, uplo::Char, diag::Char; nrhs::Int=1, transa::Char='N')
+
+Create a triangular operator for efficient solution of sparse triangular systems on GPU architectures. 
+Supports sparse matrices stored on NVIDIA, AMD, and Intel GPUs.
+
+#### Input arguments
+
+* `A`: A sparse matrix on the GPU representing the triangular system to be solved;
+* `uplo`: Specifies whether the triangular matrix `A` is upper triangular (`'U'`) or lower triangular (`'L'`);
+* `diag`: Indicates whether the diagonal is unit (`'U'`) or non-unit (`'N'`);
+* `nrhs`: Specifies the number of columns for the right-hand side(s). Defaults to 1, corresponding to solving triangular systems with a single vector as the right-hand side;
+* `transa`: Determines how the matrix `A` is applied during the triangle solves; `'N'` for no transposition, `'T'` for transpose, and `'C'` for conjugate transpose.
+
+#### Output argument
+
+* `op`: An instance of `AbstractTriangularOperator` representing the triangular operator for the specified sparse matrix and parameters.
+"""
+function TriangularOperator end