bmxam/petscwrap.jl

Parallel Julia wrapper for PETSc

petscwrap.jl's Issues

Complex number wrapper and examples

Hi,

Thanks for this awesome Julia wrapper for the PETSc. It would help me a lot.

In my regime, instead of real, I deal with complex numbers. I tried to dig out how to use this Julia wrapper to solve complex problems. Unfortunately, I cannot figure it out.

From https://github.com/bmxam/SlepcWrap.jl/blob/master/example/complex.jl#L8, I understood that PetscScalar should be complex datatype to use complex number problems. The default PETSc_jll.libpetsc would be real one. So, I set libpetsc_path = PETSc_jll.libpetsc_Float64_Complex_Int64_path in build.jl, and it showed that PetscScalar is ComplexF64.

However, when I tried to modify the example linear_system_fancy.jl to validate whether it support the complex numbers. It gave me the wrong answer. For example

n = 11
Δx = 1.0 / (n - 1)

A = create_matrix(comm; nrows_glo = n, ncols_glo = n, autosetup = true)
b = create_vector(comm; nrows_glo = n, autosetup = true)

b_start, b_end = get_range(b)

n_loc = length(b) ## Note that n_loc = b_end - b_start + 1...
b[b_start:b_end] = 1im * 2 * ones(n_loc)

A_start, A_end = get_range(A)
for i = A_start:A_end
    A[i, (i - 1):i] = 1im*[-1.0 1.0] / Δx
end

(b_start == 1) && (b[1] = 0.0)

# Assemble matrice and vector
assemble!(A)
assemble!(b)

# Set up the linear solver
ksp = create_ksp(A; autosetup = true)

# Solve the system
x = solve(ksp, b)

@show x

I would expect that it would give us the same x as I did not multiply 1im. However, I got a different and weird x.

and

M = create_matrix(comm; nrows_glo = 3, ncols_glo = 3, autosetup = true)
M_start, M_end = get_range(M)
I = [1, 1, 1, 2, 3]
J = [1, 3, 1, 3, 2]
V = [1+1im, 2+2im, 3+3im, 4+4im, 5+5im]
k = findall(x -> M_start <= x <= M_end, I) # just a trick to allow this example to run in parallel
set_values!(M, I[k], J[k], V[k], ADD_VALUES)
assemble!(M)
@show M

It cannot show the complex part of M.

Could you give me some direction on how to use PetscWrap.jl to solve complex numbers or provide a complex version in the example?

Many thanks!

Replace MPI_COMM_WORLD by PETSC_COMM_WORLD

Note : PETSC_COMM_WORLD can be accessed this way:

ptr = cglobal((:PETSC_COMM_WORLD, libpetsc), Int32)
value = unsafe_load(ptr)
unsafe_store!(ptr, new_value)

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Attempt to wrap MatCreateComposite

Hi there,

Thanks for all the work creating this cool wrapper! I've been trying to wrap MatCreateComposite, and my attempt is shown below. For the original function in PETSc, see https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateComposite.html#MatCreateComposite

function MatCreateComposite(comm::MPI.Comm, nmat::PetscInt, mats::Vector{PetscMat}, mat::PetscMat)   
    error = ccall((:MatCompositeAddMat, libpetsc), PetscErrorCode, (MPI.MPI_Comm, PetscInt, Ptr{CMat}, Ptr{CMat}),                  
                            comm, nmat, [x.ptr for x in mats], mat.ptr)    
    @assert iszero(error)
end

I've been getting segmentation faults when trying to use this function. In particular, I'm not sure that I'm handling the const Mat *mats argument correctly. I've never used C before, so I'm quite the novice to its syntax. I understand that the "*mats" means this input is passed by reference rather than value, so I should be passing in pointers. However, I'm not sure if [x.ptr for x in mats] is the right way to pass in the pointers for the vector of Mat types.

I've been trying to mirror the set up in https://www.mcs.anl.gov/petsc/petsc-current/src/mat/tutorials/ex9.c.html FYI.

Here is the example script I've been trying to use (based on your linear_system.jl)

using PetscWrap
PetscInitialize()

# Number of mesh points and mesh step
n = 11
Δx = 1. / (n - 1)

# Create a matrix and a vector
A = MatCreate() # place holder to which we write the final matrix
A1 = MatCreate() # component 1 of the composite matrix
A2 = MatCreate() # component 2 of the composite matrix
b = VecCreate()

# Set up
MatSetSizes(A1, PETSC_DECIDE, PETSC_DECIDE, n, n) # not absolutely sure this step is necessary
MatSetSizes(A2, PETSC_DECIDE, PETSC_DECIDE, n, n) # same for this step
VecSetSizes(b, PETSC_DECIDE, n)

# We can then use command-line options to set our matrix/vectors.
MatSetFromOptions(A1)  # not absolutely sure this step is necessary
MatSetFromOptions(A2) # same for this step
VecSetFromOptions(b)

# Finish the set up
MatSetUp(A1)  # not absolutely sure this step is necessary
MatSetUp(A2) # same for this step
VecSetUp(b)

# Build the b vector
b_start, b_end = VecGetOwnershipRange(b)

# Now let's build the right hand side vector. Their are various ways to do this, this is just one.
n_loc = VecGetLocalSize(b) # Note that n_loc = b_end - b_start...
VecSetValues(b, collect(b_start:b_end-1), 2 * ones(n_loc))

# And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices.
A_start, A_end = MatGetOwnershipRange(A1) # A1 is going to hold one diagonal of the differentiation matrix
for i in A_start:A_end-1
    MatSetValues(A1, [i], [i-1], [-1.] / Δx, INSERT_VALUES) # MatSetValues(A, I, J, V, INSERT_VALUES)
end

A_start, A_end = MatGetOwnershipRange(A2) # A2 is going to hold another diagonal
for i in A_start:A_end-1
    MatSetValues(A2, [i], [i], [1.] / Δx, INSERT_VALUES) # MatSetValues(A, I, J, V, INSERT_VALUES)
end

# Set boundary condition (only the proc handling index `0` is acting)
(b_start == 0) && VecSetValue(b, 0, 0.)

# Assemble matrices
VecAssemblyBegin(b)
VecAssemblyEnd(b)
MatCreateComposite(2, [A1, A2], A) # segmentation fault occurs here, but I'm trying to do A = A1 + A2