Merge pull request #2 from bmxam/dev

Dev

Project.toml

@@ -1,14 +1,15 @@
 name = "SlepcWrap"
 uuid = "c3679e3b-785e-4ccc-b734-b7685cbb935e"
 authors = ["bmxam <maxime.bouyges@gmail.com>"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 PetscWrap = "5be22e1c-01b5-4697-96eb-ef9ccdc854b8"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 MPI = "0.16.1"
-PetscWrap = "0.1.2"
-julia = "1.5"
+PetscWrap = "0.1.3"
+julia = "1.5"

README.md

@@ -21,9 +21,11 @@ pkg> add SlepcWrap
 Any contribution(s) and/or remark(s) are welcome! If you need a function that is not wrapped yet but you don't think you are capable of contributing, post an issue with a minimum working example.
 
 ## SLEPc compat.
-This version of PetscWrap.jl has been tested with slepc-3.13.1 Complex numbers are not supported yet.
+This version of PetscWrap.jl has been tested with slepc-3.13.1. Complex numbers are supported.
 
 ## How to use it
+SLEPc methods wrappers share the same name as their C equivalent : for instance `EPSCreate` or `EPSGetEigenvalue`. Furthermore, an optional "higher level" API, referred to as "fancy", is exposed : for instance `create_eps` or `get_eig`). Note that this second way of manipulating SLEPc will evolve according the package's author needs while the first one will try to follow SLEPc official API.
+
 You will find examples of use by building the documentation: `julia SlepcWrap.jl/docs/make.jl`. Here is one of the examples:
 ### Helmholtz equation
 In this example, we use the SLEPc to find the eigenvalues of the following Helmholtz equation:

docs/make.jl

@@ -8,8 +8,9 @@ using Literate
 example_src = joinpath(@__DIR__,"..","example")
 example_dir = joinpath(@__DIR__,"src","example")
 Sys.rm(example_dir; recursive=true, force=true)
-Literate.markdown(joinpath(example_src, "helmholtz1.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
-Literate.markdown(joinpath(example_src, "helmholtz2.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
+Literate.markdown(joinpath(example_src, "helmholtz.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
+Literate.markdown(joinpath(example_src, "helmholtz_fancy.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
+Literate.markdown(joinpath(example_src, "complex.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
 
 makedocs(;
     modules=[SlepcWrap],
@@ -24,13 +25,17 @@ makedocs(;
     pages=[
         "Home" => "index.md",
         "Examples" => Any[
-            "example/helmholtz1.md",
-            "example/helmholtz2.md",
+            "example/helmholtz.md",
+        ],
+        "Fancy examples" => Any[
+            "example/helmholtz_fancy.md",
+            "example/complex.md",
         ],
         "API Reference" => Any[
             "api/init.md",
             "api/eps.md",
-        ]
+        ],
+        "API fancy" => "api/fancy/fancy.md",
     ],
 )
 

docs/src/api/eps.md

@@ -1,5 +1,5 @@
 # Eigenvalue Problem solver (EPS)
 ```@autodocs
 Modules = [SlepcWrap]
-Pages   = ["eps.jl"]
+Pages   = ["src/eps.jl"]
 ```

docs/src/api/fancy/fancy.md

@@ -0,0 +1,5 @@
+# Eigenvalue Problem solver (EPS)
+```@autodocs
+Modules = [SlepcWrap]
+Pages   = ["fancy/eps.jl"]
+```

docs/src/example/complex.md

@@ -0,0 +1,78 @@
+
+# Complex numbers
+Trivial example to demonstrate the capability of handling complex numbers
+
+```julia
+using PetscWrap
+using SlepcWrap
+```
+
+Start by checking that you're using a PETSc/SLEPc build with complex number
+
+```julia
+(PetscScalar <: Real) && error("You must configure PETSc/SLEPc with complex numbers to run this example")
+```
+
+Initialize SLEPc
+
+```julia
+SlepcInitialize()
+```
+
+Create matrix
+
+```julia
+A = create_matrix(4, 4; auto_setup = true)
+```
+
+Get rows handled by the local processor
+
+```julia
+A_rstart, A_rend = get_range(A)
+```
+
+Create diag matrix with complex numbers
+
+```julia
+J = I = A_rstart:A_rend
+V = im .* I
+set_values!(A, I, J, V)
+assemble!(A)
+```
+
+Now we set up the eigenvalue solver
+
+```julia
+eps = create_eps(A; auto_setup = true)
+```
+
+Then we solve
+
+```julia
+solve!(eps)
+```
+
+Show converged eivenvalues
+
+```julia
+@show get_eigenvalues(eps)
+```
+
+Finally, let's free the memory
+
+```julia
+destroy!(A)
+destroy!(eps)
+```
+
+And call finalize when you're done
+
+```julia
+SlepcFinalize()
+
+```
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+

docs/src/example/helmholtz1.md → docs/src/example/helmholtz.md

@@ -64,33 +64,34 @@ B_rstart, B_rend = MatGetOwnershipRange(B)
 Fill matrix A  with second order derivative central scheme
 
 ```julia
-for i in A_rstart:A_rend
-    if(i == 1)
-        A[1, 1:2] = [-2., 1] / Δx^2
-    elseif (i == n)
-        A[n, n-1:n] = [1., -2.] / Δx^2
+for i in A_rstart:A_rend-1
+    if (i == 0)
+        MatSetValues(A, [0], [0, 1], [-2., 1] / Δx^2, INSERT_VALUES)
+    elseif (i == n-1)
+        MatSetValues(A, [n-1], [n-2, n-1], [1., -2.] / Δx^2, INSERT_VALUES)
     else
-        A[i, i-1:i+1] = [1., -2., 1.] / Δx^2
+        MatSetValues(A, [i], i-1:i+1, [1., -2., 1.] / Δx^2, INSERT_VALUES)
     end
 end
 ```
 
 Fill matrix B with identity matrix
 
 ```julia
-for i in B_rstart:B_rend
-    B[i,i] = -1.
+for i in B_rstart:B_rend-1
+    MatSetValue(B, i, i, -1., INSERT_VALUES)
 end
 ```
 
-Set boundary conditions : u(0) = 0 and u(1) = 0. Only the processor handling the corresponding rows are playing a role here.
+Set boundary conditions : u(0) = 0 and u(1) = 0. Only the processor
+handling the corresponding rows are playing a role here.
 
 ```julia
-(A_rstart == 1) && (A[1, 1:2] = [1. 0.] )
-(B_rstart == 1) && (B[1,   1] = 0.      )
+(A_rstart == 0) && MatSetValues(A, [0], [0,1], [1., 0.], INSERT_VALUES)
+(B_rstart == 0) && MatSetValue(B, 0, 0, 0., INSERT_VALUES)
 
-(A_rend == n) && (A[n, n-1:n] = [0. 1.] )
-(B_rend == n) && (B[n,     n] = 0.      )
+(A_rend == n) && MatSetValues(A, [n-1], [n-2,n-1], [0., 1.], INSERT_VALUES)
+(B_rend == n) && MatSetValue(B, n-1, n-1, 0., INSERT_VALUES)
 ```
 
 Assemble the matrices
@@ -126,7 +127,7 @@ nconv = EPSGetConverged(eps)
 Then we can get/display these eigenvalues (more precisely their square root, i.e ``\simeq \omega``)
 
 ```julia
-for ieig in 1:nconv
+for ieig in 0:nconv - 1
     vpr, vpi = EPSGetEigenvalue(eps, ieig)
     @show √(vpr), √(vpi)
 end
@@ -141,7 +142,7 @@ vecr, veci = MatCreateVecs(A)
 Then loop over the eigen pairs and retrieve eigenvectors
 
 ```julia
-for ieig in 1:nconv
+for ieig in 0:nconv-1
     vpr, vpi, vecpr, vecpi = EPSGetEigenpair(eps, ieig, vecr, veci)
 
     # At this point, you can call VecGetArray to obtain a Julia array (see PetscWrap examples).

docs/src/example/helmholtz2.md → docs/src/example/helmholtz_fancy.md

@@ -1,5 +1,5 @@
 
-# Helmholtz equation
+# Helmholtz equation with fancy names
 In this example, we use the SLEPc to find the eigenvalues of the following Helmholtz equation:
 ``u'' + \omega^2 u = 0`` associated to Dirichlet boundary conditions on the domain ``[0,1]``. Hence
 the theoritical eigenvalues are ``\omega = k \pi`` with ``k \in \mathbb{Z}^*``; and the associated
@@ -64,7 +64,7 @@ Fill matrix A  with second order derivative central scheme
 
 ```julia
 for i in A_rstart:A_rend
-    if(i == 1)
+    if (i == 1)
         A[1, 1:2] = [-2., 1] / Δx^2
     elseif (i == n)
         A[n, n-1:n] = [1., -2.] / Δx^2

example/complex.jl

@@ -0,0 +1,41 @@
+module Complex #hide
+# # Complex numbers
+# Trivial example to demonstrate the capability of handling complex numbers
+using PetscWrap
+using SlepcWrap
+
+# Start by checking that you're using a PETSc/SLEPc build with complex number
+(PetscScalar <: Real) && error("You must configure PETSc/SLEPc with complex numbers to run this example")
+
+# Initialize SLEPc
+SlepcInitialize()
+
+# Create matrix
+A = create_matrix(4, 4; auto_setup = true)
+
+# Get rows handled by the local processor
+A_rstart, A_rend = get_range(A)
+
+# Create diag matrix with complex numbers
+J = I = A_rstart:A_rend
+V = im .* I
+set_values!(A, I, J, V)
+assemble!(A)
+
+# Now we set up the eigenvalue solver
+eps = create_eps(A; auto_setup = true)
+
+# Then we solve
+solve!(eps)
+
+# Show converged eivenvalues
+@show get_eigenvalues(eps)
+
+# Finally, let's free the memory
+destroy!(A)
+destroy!(eps)
+
+# And call finalize when you're done
+SlepcFinalize()
+
+end #hide

example/helmholtz1.jl → example/helmholtz.jl

@@ -47,27 +47,28 @@ A_rstart, A_rend = MatGetOwnershipRange(A)
 B_rstart, B_rend = MatGetOwnershipRange(B)
 
 # Fill matrix A  with second order derivative central scheme
-for i in A_rstart:A_rend
-    if(i == 1)
-        A[1, 1:2] = [-2., 1] / Δx^2
-    elseif (i == n)
-        A[n, n-1:n] = [1., -2.] / Δx^2
+for i in A_rstart:A_rend-1
+    if (i == 0)
+        MatSetValues(A, [0], [0, 1], [-2., 1] / Δx^2, INSERT_VALUES)
+    elseif (i == n-1)
+        MatSetValues(A, [n-1], [n-2, n-1], [1., -2.] / Δx^2, INSERT_VALUES)
     else
-        A[i, i-1:i+1] = [1., -2., 1.] / Δx^2
+        MatSetValues(A, [i], i-1:i+1, [1., -2., 1.] / Δx^2, INSERT_VALUES)
     end
 end
 
 # Fill matrix B with identity matrix
-for i in B_rstart:B_rend
-    B[i,i] = -1.
+for i in B_rstart:B_rend-1
+    MatSetValue(B, i, i, -1., INSERT_VALUES)
 end
 
-# Set boundary conditions : u(0) = 0 and u(1) = 0. Only the processor handling the corresponding rows are playing a role here.
-(A_rstart == 1) && (A[1, 1:2] = [1. 0.] )
-(B_rstart == 1) && (B[1,   1] = 0.      )
+# Set boundary conditions : u(0) = 0 and u(1) = 0. Only the processor
+# handling the corresponding rows are playing a role here.
+(A_rstart == 0) && MatSetValues(A, [0], [0,1], [1., 0.], INSERT_VALUES)
+(B_rstart == 0) && MatSetValue(B, 0, 0, 0., INSERT_VALUES)
 
-(A_rend == n) && (A[n, n-1:n] = [0. 1.] )
-(B_rend == n) && (B[n,     n] = 0.      )
+(A_rend == n) && MatSetValues(A, [n-1], [n-2,n-1], [0., 1.], INSERT_VALUES)
+(B_rend == n) && MatSetValue(B, n-1, n-1, 0., INSERT_VALUES)
 
 # Assemble the matrices
 MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)
@@ -88,7 +89,7 @@ EPSSolve(eps)
 nconv = EPSGetConverged(eps)
 
 # Then we can get/display these eigenvalues (more precisely their square root, i.e ``\simeq \omega``)
-for ieig in 1:nconv
+for ieig in 0:nconv - 1
     vpr, vpi = EPSGetEigenvalue(eps, ieig)
     @show √(vpr), √(vpi)
 end
@@ -97,7 +98,7 @@ end
 vecr, veci = MatCreateVecs(A)
 
 # Then loop over the eigen pairs and retrieve eigenvectors
-for ieig in 1:nconv
+for ieig in 0:nconv-1
     vpr, vpi, vecpr, vecpi = EPSGetEigenpair(eps, ieig, vecr, veci)
 
     ## At this point, you can call VecGetArray to obtain a Julia array (see PetscWrap examples).

example/helmholtz2.jl → example/helmholtz_fancy.jl

@@ -1,5 +1,5 @@
 module Helmholtz #hide
-# # Helmholtz equation
+# # Helmholtz equation with fancy names
 # In this example, we use the SLEPc to find the eigenvalues of the following Helmholtz equation:
 # ``u'' + \omega^2 u = 0`` associated to Dirichlet boundary conditions on the domain ``[0,1]``. Hence
 # the theoritical eigenvalues are ``\omega = k \pi`` with ``k \in \mathbb{Z}^*``; and the associated
@@ -47,7 +47,7 @@ B_rstart, B_rend = get_range(B)
 
 # Fill matrix A  with second order derivative central scheme
 for i in A_rstart:A_rend
-    if(i == 1)
+    if (i == 1)
         A[1, 1:2] = [-2., 1] / Δx^2
     elseif (i == n)
         A[n, n-1:n] = [1., -2.] / Δx^2

src/SlepcWrap.jl

@@ -1,7 +1,3 @@
-"""
-    I have decided to remove all exclamation marks `!` at the end of routines "modifying" their arguments
-    when the name is the same as in SLEPc API. It was too confusing.
-"""
 module SlepcWrap
     using PetscWrap
     using MPI
@@ -30,4 +26,12 @@ module SlepcWrap
             EPSSetUp,
             EPSSolve,
             EPSView
+
+    include("fancy/eps.jl")
+    export  create_eps,
+            get_eig, get_eigs,
+            get_eigenvalue, get_eigenvalues,
+            get_eigenpair,
+            get_tolerances,
+            neigs
 end