@ -1,78 +1,97 @@
# SPDX-License-Identifier: MIT
module HYPRE
module HYPRE
module LibHYPRE
using MPI : MPI
include ( " ../lib/LibHYPRE.jl " )
using PartitionedArrays : IndexRange , MPIData , PSparseMatrix , PVector , PartitionedArrays ,
SequentialData , map_parts
using SparseArrays : SparseArrays , SparseMatrixCSC , nnz , nonzeros , nzrange , rowvals
using SparseMatricesCSR : SparseMatrixCSR , colvals , getrowptr
# Add manual methods for some ::Function signatures where the library wants function
export HYPREMatrix , HYPREVector
# pointers. Instead of creating function pointers to the Julia wrappers we can just look
# up the pointer in the library and pass that.
# TODO: Maybe this can be done automatically as post-process pass in Clang.jl
import Libdl : dlsym
function HYPRE_PCGSetPrecond ( solver , precond :: Function , precond_setup :: Function , precond_solver )
# Clang.jl auto-generated bindings and some manual methods
precond_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond ) )
include ( " LibHYPRE.jl " )
precond_setup_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond_setup ) )
using . LibHYPRE
return HYPRE_PCGSetPrecond ( solver , precond_ptr , precond_setup_ptr , precond_solver )
using . LibHYPRE : @check
end
function HYPRE_GMRESSetPrecond ( solver , precond :: Function , precond_setup :: Function , precond_solver )
# Internal namespace to hide utility functions
precond_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond ) )
include ( " Internals.jl " )
precond_setup_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond_setup ) )
return HYPRE_GMRESSetPrecond ( solver , precond_ptr , precond_setup_ptr , precond_solver )
end
function HYPRE_FlexGMRESSetPrecond ( solver , precond :: Function , precond_setup :: Function , precond_solver )
precond_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond ) )
precond_setup_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond_setup ) )
return HYPRE_FlexGMRESSetPrecond ( solver , precond_ptr , precond_setup_ptr , precond_solver )
end
function HYPRE_LGMRESSetPrecond ( solver , precond :: Function , precond_setup :: Function , precond_solver )
precond_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond ) )
precond_setup_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond_setup ) )
return HYPRE_LGMRESSetPrecond ( solver , precond_ptr , precond_setup_ptr , precond_solver )
end
function HYPRE_COGMRESSetPrecond ( solver , precond :: Function , precond_setup :: Function , precond_solver )
precond_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond ) )
precond_setup_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond_setup ) )
return HYPRE_COGMRESSetPrecond ( solver , precond_ptr , precond_setup_ptr , precond_solver )
end
function HYPRE_BiCGSTABSetPrecond ( solver , precond :: Function , precond_setup :: Function , precond_solver )
precond_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond ) )
precond_setup_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond_setup ) )
return HYPRE_BiCGSTABSetPrecond ( solver , precond_ptr , precond_setup_ptr , precond_solver )
end
function HYPRE_CGNRSetPrecond ( solver , precond :: Function , precondT :: Function , precond_setup :: Function , precond_solver )
precond_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond ) )
precondT_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precondT ) )
precond_setup_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond_setup ) )
return HYPRE_CGNRSetPrecond ( solver , precond_ptr , precondT_ptr , precond_setup_ptr , precond_solver )
end
function HYPRE_LOBPCGSetPrecond ( solver , precond :: Function , precond_setup :: Function , precond_solver )
precond_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond ) )
precond_setup_ptr = dlsym ( HYPRE_jll . libHYPRE_handle , Symbol ( precond_setup ) )
return HYPRE_LOBPCGSetPrecond ( solver , precond_ptr , precond_setup_ptr , precond_solver )
end
###############################
# HYPREMatrix and HYPREVector #
###############################
# Export everything with HYPRE_ prefix
mutable struct HYPREMatrix # <: AbstractMatrix{HYPRE_Complex}
for name in names ( @__MODULE__ ; all = true )
IJMatrix :: HYPRE_IJMatrix
if startswith ( string ( name ) , " HYPRE_ " )
ParCSRMatrix :: HYPRE_ParCSRMatrix
@eval export $ name
HYPREMatrix ( ) = new ( C_NULL , C_NULL )
end
end
mutable struct HYPREVector # <: AbstractVector{HYPRE_Complex}
IJVector :: HYPRE_IJVector
ParVector :: HYPRE_ParVector
HYPREVector ( ) = new ( C_NULL , C_NULL )
end
end
# Create a new IJMatrix, set the object type, prepare for setting values
function Internals . init_matrix ( comm :: MPI . Comm , ilower , iupper )
# Create the IJ matrix
A = HYPREMatrix ( )
IJMatrixRef = Ref { HYPRE_IJMatrix } ( C_NULL )
@check HYPRE_IJMatrixCreate ( comm , ilower , iupper , ilower , iupper , IJMatrixRef )
A . IJMatrix = IJMatrixRef [ ]
# Attach a finalizer
finalizer ( x -> HYPRE_IJMatrixDestroy ( x . IJMatrix ) , A )
# Set storage type
@check HYPRE_IJMatrixSetObjectType ( A . IJMatrix , HYPRE_PARCSR )
# Initialize to make ready for setting values
@check HYPRE_IJMatrixInitialize ( A . IJMatrix )
return A
end
end
using SparseArrays : SparseArrays , SparseMatrixCSC , nnz , rowvals , nonzeros , nzrange
# Finalize the matrix and fetch the assembled matrix
using SparseMatricesCSR : SparseMatrixCSR , colvals
# This should be called after setting all the values
using MPI : MPI
function Internals . assemble_matrix ( A :: HYPREMatrix )
using . LibHYPRE
# Finalize after setting all values
@check HYPRE_IJMatrixAssemble ( A . IJMatrix )
# Fetch the assembled CSR matrix
ParCSRMatrixRef = Ref { Ptr { Cvoid } } ( C_NULL )
@check HYPRE_IJMatrixGetObject ( A . IJMatrix , ParCSRMatrixRef )
A . ParCSRMatrix = convert ( Ptr { HYPRE_ParCSRMatrix } , ParCSRMatrixRef [ ] )
return A
end
export HYPREMatrix , HYPREVector
function Internals . init_vector ( comm :: MPI . Comm , ilower , iupper )
# Create the IJ vector
b = HYPREVector ( )
b_ref = Ref { HYPRE_IJVector } ( C_NULL )
@check HYPRE_IJVectorCreate ( comm , ilower , iupper , b_ref )
b . IJVector = b_ref [ ]
# Attach a finalizer
finalizer ( x -> HYPRE_IJVectorDestroy ( x . IJVector ) , b ) # Set storage type
# Set storage type
@check HYPRE_IJVectorSetObjectType ( b . IJVector , HYPRE_PARCSR )
# Initialize to make ready for setting values
@check HYPRE_IJVectorInitialize ( b . IJVector )
return b
end
module Internals
function Internals . assemble_vector ( b :: HYPREVector )
function check_n_rows end
# Finalize after setting all values
function to_hypre_data end
@check HYPRE_IJVectorAssemble ( b . IJVector )
# Fetch the assembled vector
par_b_ref = Ref { Ptr { Cvoid } } ( C_NULL )
@check HYPRE_IJVectorGetObject ( b . IJVector , par_b_ref )
b . ParVector = convert ( Ptr { HYPRE_ParVector } , par_b_ref [ ] )
return b
end
end
######################################
# SparseMatrixCS(C|R) -> HYPREMatrix #
######################################
function Internals . check_n_rows ( A , ilower , iupper )
function Internals . check_n_rows ( A , ilower , iupper )
if size ( A , 1 ) != ( iupper - ilower + 1 )
if size ( A , 1 ) != ( iupper - ilower + 1 )
throw ( ArgumentError ( " number of rows in matrix does not match global start/end rows ilower and iupper " ) )
throw ( ArgumentError ( " number of rows in matrix does not match global start/end rows ilower and iupper " ) )
@ -85,7 +104,7 @@ function Internals.to_hypre_data(A::SparseMatrixCSC, ilower, iupper)
A_rows = rowvals ( A )
A_rows = rowvals ( A )
A_vals = nonzeros ( A )
A_vals = nonzeros ( A )
# Initialize data as HYPRE expec ts
# Initialize the data buffers HYPRE wan ts
nrows = HYPRE_Int ( iupper - ilower + 1 ) # Total number of rows
nrows = HYPRE_Int ( iupper - ilower + 1 ) # Total number of rows
ncols = zeros ( HYPRE_Int , nrows ) # Number of colums for each row
ncols = zeros ( HYPRE_Int , nrows ) # Number of colums for each row
rows = collect ( HYPRE_BigInt , ilower : iupper ) # The row indices
rows = collect ( HYPRE_BigInt , ilower : iupper ) # The row indices
@ -123,7 +142,7 @@ function Internals.to_hypre_data(A::SparseMatrixCSR, ilower, iupper)
A_cols = colvals ( A )
A_cols = colvals ( A )
A_vals = nonzeros ( A )
A_vals = nonzeros ( A )
# Initialize data as HYPRE expec ts
# Initialize the data buffers HYPRE wan ts
nrows = HYPRE_Int ( iupper - ilower + 1 ) # Total number of rows
nrows = HYPRE_Int ( iupper - ilower + 1 ) # Total number of rows
ncols = Vector { HYPRE_Int } ( undef , nrows ) # Number of colums for each row
ncols = Vector { HYPRE_Int } ( undef , nrows ) # Number of colums for each row
rows = collect ( HYPRE_BigInt , ilower : iupper ) # The row indices
rows = collect ( HYPRE_BigInt , ilower : iupper ) # The row indices
@ -147,42 +166,17 @@ function Internals.to_hypre_data(A::SparseMatrixCSR, ilower, iupper)
return nrows , ncols , rows , cols , values
return nrows , ncols , rows , cols , values
end
end
mutable struct HYPREMatrix # <: AbstractMatrix{HYPRE_Complex}
IJMatrix :: HYPRE_IJMatrix
ParCSRMatrix :: HYPRE_ParCSRMatrix
HYPREMatrix ( ) = new ( C_NULL , C_NULL )
end
function HYPREMatrix ( B :: Union { SparseMatrixCSC , SparseMatrixCSR } , ilower , iupper , comm :: MPI . Comm = MPI . COMM_WORLD )
function HYPREMatrix ( B :: Union { SparseMatrixCSC , SparseMatrixCSR } , ilower , iupper , comm :: MPI . Comm = MPI . COMM_WORLD )
# Compute indices/values in the format SetValues expect
A = Internals . init_matrix ( comm , ilower , iupper )
nrows , ncols , rows , cols , values = Internals . to_hypre_data ( B , ilower , iupper )
nrows , ncols , rows , cols , values = Internals . to_hypre_data ( B , ilower , iupper )
# Create the IJ matrix
@check HYPRE_IJMatrixSetValues ( A . IJMatrix , nrows , ncols , rows , cols , values )
A = HYPREMatrix ( )
Internals . assemble_matrix ( A )
IJMatrixRef = Ref { HYPRE_IJMatrix } ( C_NULL )
HYPRE_IJMatrixCreate ( comm , ilower , iupper , ilower , iupper , IJMatrixRef )
A . IJMatrix = IJMatrixRef [ ]
# Attach a finalizer
finalizer ( x -> HYPRE_IJMatrixDestroy ( x . IJMatrix ) , A )
# Set storage type
HYPRE_IJMatrixSetObjectType ( A . IJMatrix , HYPRE_PARCSR )
# Initialize to make ready for setting values
HYPRE_IJMatrixInitialize ( A . IJMatrix )
# Set all the values
HYPRE_IJMatrixSetValues ( A . IJMatrix , nrows , ncols , rows , cols , values )
# Finalize
HYPRE_IJMatrixAssemble ( A . IJMatrix )
# Fetch the assembled CSR matrix
ParCSRMatrixRef = Ref { Ptr { Cvoid } } ( C_NULL )
HYPRE_IJMatrixGetObject ( A . IJMatrix , ParCSRMatrixRef )
A . ParCSRMatrix = convert ( Ptr { HYPRE_ParCSRMatrix } , ParCSRMatrixRef [ ] )
return A
return A
end
end
mutable struct HYPREVector # <: AbstractVector{HYPRE_Complex}
#########################
IJVector :: HYPRE_IJVector
# Vector -> HYPREVector #
ParVector :: HYPRE_ParVector
#########################
HYPREVector ( ) = new ( C_NULL , C_NULL )
end
function Internals . to_hypre_data ( x :: Vector , ilower , iupper )
function Internals . to_hypre_data ( x :: Vector , ilower , iupper )
Internals . check_n_rows ( x , ilower , iupper )
Internals . check_n_rows ( x , ilower , iupper )
@ -190,25 +184,188 @@ function Internals.to_hypre_data(x::Vector, ilower, iupper)
values = convert ( Vector { HYPRE_Complex } , x )
values = convert ( Vector { HYPRE_Complex } , x )
return HYPRE_Int ( length ( indices ) ) , indices , values
return HYPRE_Int ( length ( indices ) ) , indices , values
end
end
# TODO: Internals.to_hypre_data(x::SparseVector, ilower, iupper) (?)
function HYPREVector ( x :: Vector , ilower , iupper , comm = MPI . COMM_WORLD )
function HYPREVector ( x :: Vector , ilower , iupper , comm = MPI . COMM_WORLD )
b = Internals . init_vector ( comm , ilower , iupper )
nvalues , indices , values = Internals . to_hypre_data ( x , ilower , iupper )
nvalues , indices , values = Internals . to_hypre_data ( x , ilower , iupper )
b = HYPREVector ( )
@check HYPRE_IJVectorSetValues ( b . IJVector , nvalues , indices , values )
b_ref = Ref { HYPRE_IJVector } ( C_NULL )
Internals . assemble_vector ( b )
HYPRE_IJVectorCreate ( comm , ilower , iupper , b_ref )
return b
b . IJVector = b_ref [ ]
end
finalizer ( x -> HYPRE_IJVectorDestroy ( x . IJVector ) , b ) # Set storage type
HYPRE_IJVectorSetObjectType ( b . IJVector , HYPRE_PARCSR )
# Initialize to make ready for setting values
##################################################
HYPRE_IJVectorInitialize ( b . IJVector )
# PartitionedArrays.PSparseMatrix -> HYPREMatrix #
# Set the values
##################################################
HYPRE_IJVectorSetValues ( b . IJVector , nvalues , indices , values )
# TODO: This has some duplicated code with to_hypre_data(::SparseMatrixCSC, ilower, iupper)
function Internals . to_hypre_data ( A :: SparseMatrixCSC , r :: IndexRange , c :: IndexRange )
@assert r . oid_to_lid isa UnitRange && r . oid_to_lid . start == 1
ilower = r . lid_to_gid [ r . oid_to_lid . start ]
iupper = r . lid_to_gid [ r . oid_to_lid . stop ]
a_rows = rowvals ( A )
a_vals = nonzeros ( A )
# Initialize the data buffers HYPRE wants
nrows = HYPRE_Int ( iupper - ilower + 1 ) # Total number of rows
ncols = zeros ( HYPRE_Int , nrows ) # Number of colums for each row
rows = collect ( HYPRE_BigInt , ilower : iupper ) # The row indices
# cols = Vector{HYPRE_BigInt}(undef, nnz) # The column indices
# values = Vector{HYPRE_Complex}(undef, nnz) # The values
# First pass to count nnz per row (note that the fact that columns are permuted
# doesn't matter for this pass)
a_rows = rowvals ( A )
a_vals = nonzeros ( A )
@inbounds for j in 1 : size ( A , 2 )
for i in nzrange ( A , j )
row = a_rows [ i ]
row > r . oid_to_lid . stop && continue # Skip ghost rows
# grow = r.lid_to_gid[lrow]
ncols [ row ] += 1
end
end
# Initialize remaining buffers now that nnz is known
nnz = sum ( ncols )
cols = Vector { HYPRE_BigInt } ( undef , nnz )
values = Vector { HYPRE_Complex } ( undef , nnz )
# Keep track of the last index used for every row
lastinds = zeros ( Int , nrows )
cumsum! ( ( @view lastinds [ 2 : end ] ) , ( @view ncols [ 1 : end - 1 ] ) )
# Second pass to populate the output -- here we need to take care of the permutation
# of columns. TODO: Problem that they are not sorted?
@inbounds for j in 1 : size ( A , 2 )
for i in nzrange ( A , j )
row = a_rows [ i ]
row > r . oid_to_lid . stop && continue # Skip ghost rows
k = lastinds [ row ] += 1
val = a_vals [ i ]
cols [ k ] = c . lid_to_gid [ j ]
values [ k ] = val
end
end
return nrows , ncols , rows , cols , values
end
# TODO: Possibly this can be optimized if it is possible to pass overlong vectors to HYPRE.
# At least values should be possible to directly share, but cols needs to translated
# to global ids.
function Internals . to_hypre_data ( A :: SparseMatrixCSR , r :: IndexRange , c :: IndexRange )
@assert r . oid_to_lid isa UnitRange && r . oid_to_lid . start == 1
ilower = r . lid_to_gid [ r . oid_to_lid . start ]
iupper = r . lid_to_gid [ r . oid_to_lid . stop ]
a_cols = colvals ( A )
a_vals = nonzeros ( A )
nnz = getrowptr ( A ) [ r . oid_to_lid . stop + 1 ] - 1
# Initialize the data buffers HYPRE wants
nrows = HYPRE_Int ( iupper - ilower + 1 ) # Total number of rows
ncols = zeros ( HYPRE_Int , nrows ) # Number of colums for each row
rows = collect ( HYPRE_BigInt , ilower : iupper ) # The row indices
cols = Vector { HYPRE_BigInt } ( undef , nnz ) # The column indices
values = Vector { HYPRE_Complex } ( undef , nnz ) # The values
# Loop over the (owned) rows and collect all values
k = 0
@inbounds for i in r . oid_to_lid
nzr = nzrange ( A , i )
ncols [ i ] = length ( nzr )
for j in nzr
k += 1
col = a_cols [ j ]
val = a_vals [ j ]
cols [ k ] = c . lid_to_gid [ col ]
values [ k ] = val
end
end
@assert nnz == k
return nrows , ncols , rows , cols , values
end
function Internals . get_comm ( A :: Union { PSparseMatrix { <: Any , <: M } , PVector { <: Any , <: M } } ) where M <: MPIData
return A . rows . partition . comm
end
Internals . get_comm ( _ :: Union { PSparseMatrix , PVector } ) = MPI . COMM_WORLD
function Internals . get_proc_rows ( A :: Union { PSparseMatrix { <: Any , <: M } , PVector { <: Any , <: M } } ) where M <: MPIData
r = A . rows . partition . part
ilower :: HYPRE_BigInt = r . lid_to_gid [ r . oid_to_lid [ 1 ] ]
iupper :: HYPRE_BigInt = r . lid_to_gid [ r . oid_to_lid [ end ] ]
return ilower , iupper
end
function Internals . get_proc_rows ( A :: Union { PSparseMatrix { <: Any , <: S } , PVector { <: Any , <: S } } ) where S <: SequentialData
ilower :: HYPRE_BigInt = typemax ( HYPRE_BigInt )
iupper :: HYPRE_BigInt = typemin ( HYPRE_BigInt )
for r in A . rows . partition . parts
ilower = min ( r . lid_to_gid [ r . oid_to_lid [ 1 ] ] , ilower )
iupper = max ( r . lid_to_gid [ r . oid_to_lid [ end ] ] , iupper )
end
return ilower , iupper
end
function HYPREMatrix ( B :: PSparseMatrix )
# Use the same communicator as the matrix
comm = Internals . get_comm ( B )
# Fetch rows owned by this process
ilower , iupper = Internals . get_proc_rows ( B )
# Create the IJ matrix
A = Internals . init_matrix ( comm , ilower , iupper )
# Set all the values
map_parts ( B . values , B . rows . partition , B . cols . partition ) do Bv , Br , Bc
nrows , ncols , rows , cols , values = Internals . to_hypre_data ( Bv , Br , Bc )
@check HYPRE_IJMatrixSetValues ( A . IJMatrix , nrows , ncols , rows , cols , values )
return nothing
end
# Finalize
# Finalize
HYPRE_IJVectorAssemble ( b . IJVector )
Internals . assemble_matrix ( A )
# Fetch the assembled object
return A
par_b_ref = Ref { Ptr { Cvoid } } ( C_NULL )
end
HYPRE_IJVectorGetObject ( b . IJVector , par_b_ref )
b . ParVector = convert ( Ptr { HYPRE_ParVector } , par_b_ref [ ] )
############################################
# PartitionedArrays.PVector -> HYPREVector #
############################################
#
function HYPREVector ( v :: PVector )
# Use the same communicator as the matrix
comm = Internals . get_comm ( v )
# Fetch rows owned by this process
ilower , iupper = Internals . get_proc_rows ( v )
# Create the IJ vector
b = Internals . init_vector ( comm , ilower , iupper )
# Set all the values
map_parts ( v . values , v . owned_values , v . rows . partition ) do vv , vo , vr
ilower_part = vr . lid_to_gid [ vr . oid_to_lid . start ]
iupper_part = vr . lid_to_gid [ vr . oid_to_lid . stop ]
# Option 1: Set all values
nvalues = HYPRE_Int ( iupper_part - ilower_part + 1 )
indices = collect ( HYPRE_BigInt , ilower_part : iupper_part )
# TODO: Could probably just pass the full vector even if it is too long
# values = convert(Vector{HYPRE_Complex}, vv)
values = collect ( HYPRE_Complex , vo )
# # Option 2: Set only non-zeros
# indices = HYPRE_BigInt[]
# values = HYPRE_Complex[]
# for (i, vi) in zip(ilower_part:iupper_part, vo)
# if !iszero(vi)
# push!(indices, i)
# push!(values, vi)
# end
# end
# nvalues = length(indices)
@check HYPRE_IJVectorSetValues ( b . IJVector , nvalues , indices , values )
return nothing
end
# Finalize
Internals . assemble_vector ( b )
return b
return b
end
end