bmxam / petscwrap.jl Goto Github PK
View Code? Open in Web Editor NEWParallel Julia wrapper for PETSc
License: MIT License
Parallel Julia wrapper for PETSc
License: MIT License
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!
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)
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!
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
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.