reintegrate old changes, drop numeric resolves

examples/highlevel.jl

@@ -9,7 +9,7 @@ using MatrixMarket
 using SparseBase
 using LinearAlgebra
 MPI.Init()
-nprow, npcol, nrhs = 1, 1, 4
+nprow, npcol, nrhs = 2, 2, 4
 root = 0
 comm = MPI.COMM_WORLD
 grid = Grid{Int32}(nprow, npcol, comm)
@@ -22,6 +22,11 @@ isroot = iam == root
 coo, b, xtrue = SuperLUDIST.mmread_and_generatesolution(
     Float64, Int32, nrhs, joinpath(@__DIR__, "add32.mtx"), grid; root
 )
+
+# A second rhs and xtrue for testing different sized b.
+b2 = [b;; b]
+x2 = [xtrue;; xtrue]
+
 csr = isroot ? convert(SparseBase.CSRStore, coo) : nothing
 chunksizes = isroot ? distribute_evenly(size(csr, 1), nprow * npcol) : nothing
 
@@ -37,19 +42,25 @@ A = Communication.scatterstore!(
 
 b_local = b[A.first_row : A.first_row + localsize(A, 1) - 1, :] # shrink b
 xtrue_local = xtrue[A.first_row : A.first_row + localsize(A, 1) - 1, :] # shrink xtrue
+F = lu!(A);
 
-options = SuperLUDIST.Options()
-stat = SuperLUDIST.LUStat{Int32}()
-b1 = Matrix{Float64}(undef, localsize(A, 2), 0)
-_, F = pgssvx!(A, b1; options, stat);
+b_local = F \ b_local
 
-b_local, F = pgssvx!(F, b_local);
-GC.gc()
 if !(iam == root) || (nprow * npcol == 1)
     SuperLUDIST.inf_norm_error_dist(b_local, xtrue_local, grid)
+    SuperLUDIST.PStatPrint(F) # printing may be messy.
 end
-SuperLUDIST.PStatPrint(options, stat, grid)
+
+b2_local = b2[A.first_row : A.first_row + localsize(A, 1) - 1, :] # shrink b2
+xtrue2_local = x2[A.first_row : A.first_row + localsize(A, 1) - 1, :] # shrink xtrue2
+
+b2_local = F \ b2_local
+
+if !(iam == root) || (nprow * npcol == 1)
+    SuperLUDIST.inf_norm_error_dist(b2_local, xtrue2_local, grid)
+    SuperLUDIST.PStatPrint(F) #printing may be messy.
+end
+
 
 # @show iam b_local xtrue_local
 MPI.Finalize()
-

examples/pdrive.jl

@@ -9,7 +9,7 @@ using MatrixMarket
 using SparseBase
 using LinearAlgebra
 MPI.Init()
-nprow, npcol, nrhs = 1, 1, 2
+nprow, npcol, nrhs = 1, 1, 3
 root = 0
 comm = MPI.COMM_WORLD
 grid = Grid{Int32}(nprow, npcol, comm)
@@ -44,17 +44,15 @@ xtrue_local = xtrue[A.first_row : A.first_row + localsize(A, 1) - 1, :] # shrink
 # for each user.
 
 # creating options and stat is optional, they will be created if not provided.
-options = SuperLUDIST.Options()
-stat = SuperLUDIST.LUStat{Int32}()
+
 # b1 = Matrix{Float64}(undef, localsize(A, 2), 2)
-b1 = rand(localsize(A, 1))
-_, F = pgssvx!(A, b1; options, stat);
-@show F.options
+b1 = rand(localsize(A, 1), 0)
+_, F = pgssvx!(A, b1);
 b_local, F = pgssvx!(F, b_local);
 if !(iam == root) || (nprow * npcol == 1)
     SuperLUDIST.inf_norm_error_dist(b_local, xtrue_local, grid)
 end
-SuperLUDIST.PStatPrint(options, stat, grid)
+SuperLUDIST.PStatPrint(F)
 
 # @show iam b_local xtrue_local
 MPI.Finalize()

examples/pdrive_ABglobal.jl

@@ -1,10 +1,10 @@
 using MPI
 using SuperLUDIST
-using SuperLUDIST: Grid, Options, LUStat, ScalePermStruct,
+using SuperLUDIST: Grid, Options, LUStat, ScalePerm,
     ReplicatedSuperMatrix, pgssvx!
 using SuperLUDIST.Common
 using MatrixMarket
-nprow, npcol, nrhs = (2, 2, 1)
+nprow, npcol, nrhs = (1, 1, 1)
 root = 0
 MPI.Init()
 comm = MPI.COMM_WORLD
@@ -13,10 +13,10 @@ iam = grid.iam
 
 # This function handles broadcasting internally!
 A = MatrixMarket.mmread(
-    ReplicatedSuperMatrix{Float64, Int}, 
+    SuperLUDIST.ReplicatedSuperMatrix{Float64, Int}, 
     joinpath(@__DIR__, "add32.mtx"),
     grid
-)
+) ;
 
 # on single nodes this will help prevent oversubscription of threads.
 SuperLUDIST.superlu_set_num_threads(Int64, 2)

lib/common.jl

@@ -33,7 +33,7 @@ struct superlu_scope_t
     Iam::Cint
 end
 
-struct gridinfo_t{I}
+mutable struct gridinfo_t{I}
     comm::MPI_Comm
     rscp::superlu_scope_t
     cscp::superlu_scope_t
@@ -42,7 +42,7 @@ struct gridinfo_t{I}
     npcol::I
 end
 
-struct gridinfo3d_t{I}
+mutable struct gridinfo3d_t{I}
     comm::MPI_Comm
     rscp::superlu_scope_t
     cscp::superlu_scope_t
@@ -54,4 +54,4 @@ struct gridinfo3d_t{I}
     npdep::I
     rankorder::Cint
 end
-end
+end

src/drivers.jl

@@ -12,13 +12,21 @@ function pgssvx!(
     A::ReplicatedSuperMatrix{Tv, Ti},
     b::VecOrMat{Tv};
     options = Options(),
-    perm = ScalePermStruct{Tv, Ti}(size(A)...),
-    LU = LUStruct{Tv, Ti}(size(A, 2), A.grid),
+    perm = ScalePerm{Tv, Ti}(size(A)...),
+    LU = LUFactors{Tv, Ti}(size(A, 2), A.grid),
     stat = LUStat{Ti}(),
     berr = Vector{Tv}(undef, size(b, 2))
 ) where {Tv, Ti}
-return pgssvx!(SuperLUFactorization(A, options, nothing, perm, LU, stat, berr, b), b)
-
+    return pgssvx!(SuperLUFactorization(A, options, nothing, perm, LU, stat, berr, b), b)
+end
+function pgssvx!(
+    A::ReplicatedSuperMatrix{Tv, Ti};
+    options = Options(),
+    perm = ScalePerm{Tv, Ti}(size(A)...),
+    LU = LUFactors{Tv, Ti}(size(A, 2), A.grid),
+    stat = LUStat{Ti}()
+) where {Tv, Ti}
+    return pgssvx!(SuperLUFactorization(A, options, nothing, perm, LU, stat, berr, b), b)
 end
 
 """
@@ -35,36 +43,47 @@ function pgssvx!(
     A::DistributedSuperMatrix{Tv, Ti},
     b::VecOrMat{Tv};
     options = Options(),
-    Solve = SOLVE{Tv, Ti}(options),
-    perm = ScalePermStruct{Tv, Ti}(size(A)...),
-    LU = LUStruct{Tv, Ti}(size(A, 2), A.grid),
+    Solve = SolveData{Tv, Ti}(options),
+    perm = ScalePerm{Tv, Ti}(size(A)...),
+    LU = LUFactors{Tv, Ti}(size(A, 2), A.grid),
     stat = LUStat{Ti}(),
     berr = Vector{Tv}(undef, size(b, 2))
 ) where {Tv, Ti}
     return pgssvx!(SuperLUFactorization(A, options, Solve, perm, LU, stat, berr, b), b)
 end
+function pgssvx!(
+    A::DistributedSuperMatrix{Tv, Ti};
+    options = Options(),
+    Solve = SolveData{Tv, Ti}(options),
+    perm = ScalePerm{Tv, Ti}(size(A)...),
+    LU = LUFactors{Tv, Ti}(size(A, 2), A.grid),
+    stat = LUStat{Ti}(),
+) where {Tv, Ti}
+    b = Matrix{Tv}(undef, SparseBase.Communication.localsize(A, 1), 0)
+    pgssvx!(A, b; options, Solve, perm, LU, stat)
+end
 
 function pgssvx!(F::SuperLUFactorization{T, I, <:ReplicatedSuperMatrix{T, I}}, b::VecOrMat{T}) where {T, I}
-    (; mat, options, perm, lu, stat, berr) = F
-    b, _ = pgssvx_ABglobal!(options, mat, perm, b, lu, berr, stat)
+    (; mat, options, perm, factors, stat, berr) = F
+    b, _ = pgssvx_ABglobal!(options, mat, perm, b, factors, berr, stat)
     F.options.Fact = Common.FACTORED
     F.b = b
     return b, F
 end
 function pgssvx!(F::SuperLUFactorization{T, I, <:DistributedSuperMatrix{T, I}}, b::VecOrMat{T}) where {T, I}
-    (; mat, options, solve, perm, lu, stat, berr) = F
+    (; mat, options, solve, perm, factors, stat, berr) = F
     currentnrhs = size(F.b, 2)
     if currentnrhs != size(b, 2)
-        # F = pgstrs_prep!(F)
+        F = pgstrs_prep!(F)
         pgstrs_init!(
             F.solve, 
             reverse(Communication.localsize(F.mat))...,
             size(b, 2), F.mat.first_row - 1, F.perm,
-            F.lu, F.mat.grid
+            F.factors, F.mat.grid
         )
     end
 
-    b, _ = pgssvx_ABdist!(options, mat, perm, b, lu, solve, berr, stat)
+    b, _ = pgssvx_ABdist!(options, mat, perm, b, factors, solve, berr, stat)
     F.options.Fact = Common.FACTORED
     F.b = b
     return b, F
@@ -77,9 +96,9 @@ L = Symbol(:SuperLU_, Symbol(I))
 function pgssvx_ABglobal!(
     options,
     A::ReplicatedSuperMatrix{$T, $I},
-    perm::ScalePermStruct{$T, $I},
+    perm::ScalePerm{$T, $I},
     b::Array{$T},
-    LU::LUStruct{$T, $I},
+    LU::LUFactors{$T, $I},
     berr, stat::LUStat{$I}
 )
     info = Ref{Int32}()
@@ -94,10 +113,10 @@ end
 function pgssvx_ABdist!(
     options,
     A::DistributedSuperMatrix{$T, $I},
-    perm::ScalePermStruct{$T, $I},
+    perm::ScalePerm{$T, $I},
     b::Array{$T},
-    LU::LUStruct{$T, $I},
-    Solve::SOLVE{$T, $I},
+    LU::LUFactors{$T, $I},
+    Solve::SolveData{$T, $I},
     berr, stat::LUStat{$I}
 )
     info = Ref{Int32}()
@@ -118,10 +137,10 @@ function inf_norm_error_dist(x::Array{$T}, xtrue::Array{$T}, grid::Grid{$I})
     )
 end
 function pgstrs_init!(
-    solve::SOLVE{$T, $I},
+    solve::SolveData{$T, $I},
     n, m_local, nrhs, first_row,
-    scaleperm::ScalePermStruct{$T, $I},
-    lu::LUStruct{$T, $I},
+    scaleperm::ScalePerm{$T, $I},
+    lu::LUFactors{$T, $I},
     grid::Grid{$I}
 )
     $L.$(Symbol(:p, prefixsymbol(T), :gstrs_init))(
@@ -136,11 +155,11 @@ function pgstrs_prep!(
 )
     if size(F.b, 2) != 0
         gstrs = unsafe_load(F.solve.gstrs_comm)
-        @show gstrs
         $L.superlu_free_dist(gstrs.B_to_X_SendCnt)
         $L.superlu_free_dist(gstrs.X_to_B_SendCnt)
         $L.superlu_free_dist(gstrs.ptr_to_ibuf)
-        @show gstrs
+    else # there has been no gstrs_comm malloc'd, so do that.
+
     end
     return F
 end

src/highlevel.jl

@@ -1,8 +1,8 @@
 function LinearAlgebra.lu!(
-    A::AbstractSuperMatrix{Tv, Ti};
+    A::AbstractSuperMatrix{Tv, Ti}, 
+    b_local = ones(Tv, Communication.localsize(A, 2), 1);
     kwargs...
 ) where {Tv, Ti}
-    b_local = Matrix{Tv}(undef, Communication.localsize(A, 2), 0)
     return pgssvx!(A, b_local; kwargs...)[2] # F, drop b_local.
 end
 

src/lowlevel.jl

@@ -68,7 +68,7 @@ L = Symbol(String(:SuperLU_) * String(I))
     function Destroy_LU(r::Base.RefValue{$(prefixname(T, :LUstruct_t)){$I}}, n, grid)
         $L.$(prefixname(T, :Destroy_LU))(n, grid, r) # why do I have to grid.grid here?!
     end
-    function LUstructInit(r::Base.RefValue{$(prefixname(T, :LUstruct_t)){$I}}, n, grid)
+    function LUstructInit(r::$(prefixname(T, :LUstruct_t)){$I}, n, grid)
         $L.$(prefixname(T, :LUstructInit))(n, r)
         return finalizer(r) do x
             # !MPI.Finalized() && Destroy_LU(x, n, grid)
@@ -88,4 +88,4 @@ function FillRHS_dist!(b::Array, A::AbstractSuperMatrix, xtrue::Array;
     trans = false, nrhs = size(xtrue, 2), ldx = size(xtrue, 1), ldb = size(b, 1))
     FillRHS_dist!(b, A, xtrue, trans, nrhs, ldx, ldb)
     return b
-end
+end

src/matrixmarket.jl

@@ -71,7 +71,7 @@ A = MatrixMarket.mmread(ReplicatedSuperMatrix{Float64, Int32}, )
 """
 function MatrixMarket.mmread(
   ::Type{<:DistributedSuperMatrix{Tv, Ti}}, filename, grid::Grid{Ti}; 
-  desymmetrize = true, root = 0, partitioner = distribute_evenly
+  desymmetrize = true, root = 0, partitioner = Communication.distribute_evenly
 ) where {Tv <: Union{Float32, Float64, ComplexF64}, Ti <: Union{Int32, Int64}}
   comm = grid.comm
   rank = MPI.Comm_rank(comm)