Merge pull request #118 from gridap/add/SparseMatricesCSR

Add/SparseMatricesCSR

Manifest.toml

@@ -245,6 +245,12 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[[SparseMatricesCSR]]
+deps = ["LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "7896f9c00e5413cbaecd758ade9fe67bf50ad32e"
+uuid = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
+version = "0.3.1"
+
 [[SpecialFunctions]]
 deps = ["BinDeps", "BinaryProvider", "Libdl"]
 git-tree-sha1 = "3bdd374b6fd78faf0119b8c5d538788dbf910c6e"

Project.toml

@@ -13,6 +13,7 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TensorPolynomialBases = "f77b161d-e03a-5d68-b365-d5306aabeaa6"
 TensorValues = "31c64edf-cdeb-50e4-845d-ed2334360c20"
@@ -28,6 +29,7 @@ LineSearches = "7.0.1"
 NLsolve = "4.1.0"
 QuadGK = "2.1.0"
 Reexport = "0.2.0"
+SparseMatricesCSR = "0.3.1"
 StaticArrays = "0.10.3"
 TensorPolynomialBases = "0.1.4"
 TensorValues = "0.3.5"

src/FESpaces/Assemblers.jl

@@ -5,11 +5,13 @@ using Gridap
 using Gridap.Helpers
 
 using SparseArrays
+using SparseMatricesCSR
 
 export Assembler
 export SparseMatrixAssembler
 export assemble
 export assemble!
+export sparse_from_coo
 
 """
 Abstract assembly operator
@@ -80,21 +82,28 @@ function assemble!(r,a::Assembler,cv::CellMatrix)
 end
 
 """
-Assembler that produces SparseMatrices from the SparseArrays package
+Assembler that produces SparseMatrices from the SparseArrays and SparseMatricesCSR packages
 """
-struct SparseMatrixAssembler{E} <: Assembler{SparseMatrixCSC{E,Int},Vector{E}}
+struct SparseMatrixAssembler{E,M} <: Assembler{M,Vector{E}}
   testfesp::FESpace
   trialfesp::FESpace
+
+  function SparseMatrixAssembler(
+    ::Type{M},
+    testfesp::FESpace{D,Z,T},
+    trialfesp::FESpace{D,Z,T}) where {M<:AbstractSparseMatrix,D,Z,T}
+    E = eltype(T)
+    new{E,M}(testfesp,trialfesp) 
+  end
 end
 
 function SparseMatrixAssembler(test::FESpace{D,Z,T}, trial::FESpace{D,Z,T}) where {D,Z,T}
-  E = eltype(T)
-  SparseMatrixAssembler{E}(trial,test)
+  SparseMatrixAssembler(SparseMatrixCSC, trial,test)
 end
 
 function assemble(
   this::SparseMatrixAssembler{E},
-  vals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where E
+  vals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E}
 
   n = num_free_dofs(this.testfesp)
   vec = zeros(E,n)
@@ -105,7 +114,7 @@ end
 function assemble!(
   vec::Vector{E},
   this::SparseMatrixAssembler{E},
-  allvals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where E
+  allvals::Vararg{Tuple{<:CellVector,<:CellNumber}}) where {E}
 
   vec .= zero(E)
   rows = celldofids(this.testfesp)
@@ -127,35 +136,35 @@ function _assemble_vector!(vec,vals,rows)
 end
 
 function assemble(
-  this::SparseMatrixAssembler{E},
-  allvals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where E
+  this::SparseMatrixAssembler{E,M},
+  allvals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where {E,M}
 
   I = Int
   aux_row = I[]; aux_col = I[]; aux_val = E[]
 
   _rows_m = celldofids(this.testfesp)
   _cols_m = celldofids(this.trialfesp)
-
   for (vals, cellids_row, cellids_col) in allvals
     _vals = apply_constraints_rows(this.testfesp, vals, cellids_row)
     rows_m = reindex(_rows_m, cellids_row)
     vals_m = apply_constraints_cols(this.trialfesp, _vals, cellids_col)
     cols_m = reindex(_cols_m, cellids_col)
-    _assemble_sparse_matrix_values!(
+    _assemble_sparse_matrix_values!(M,
       aux_row,aux_col,aux_val,vals_m,rows_m,cols_m)
   end
-  sparse(aux_row,aux_col,aux_val)
+  num_rows = num_free_dofs(this.testfesp)
+  num_cols = num_free_dofs(this.trialfesp)
+  finalize_coo!(M,aux_row,aux_col,aux_val,num_rows,num_cols)
+  sparse_from_coo(M,aux_row,aux_col,aux_val)
 end
 
-function _assemble_sparse_matrix_values!(aux_row,aux_col,aux_val,vals,rows,cols)
+function _assemble_sparse_matrix_values!(::Type{M},aux_row,aux_col,aux_val,vals,rows,cols) where {M}
   for (rows_c, cols_c, vals_c) in zip(rows,cols,vals)
      for (j,gidcol) in enumerate(cols_c)
        if gidcol > 0
         for (i,gidrow) in enumerate(rows_c)
           if gidrow > 0
-            push!(aux_row, gidrow)
-            push!(aux_col, gidcol)
-            push!(aux_val, vals_c[i,j])
+            push_coo!(M,aux_row, aux_col, aux_val, gidrow, gidcol, vals_c[i,j])
           end
         end
       end
@@ -164,13 +173,26 @@ function _assemble_sparse_matrix_values!(aux_row,aux_col,aux_val,vals,rows,cols)
 end
 
 function assemble!(
-  mat::SparseMatrixCSC{E},
-  this::SparseMatrixAssembler{E},
-  vals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where E
+  mat::AbstractSparseMatrix{E,Int},
+  this::SparseMatrixAssembler{E,M},
+  vals::Vararg{Tuple{<:CellMatrix,<:CellNumber,<:CellNumber}}) where {E,M}
   # This routine can be optimized a lot taking into a count the sparsity graph of mat
   # For the moment we create an intermediate matrix and then transfer the nz values
   m = assemble(this,vals...)
-  mat.nzval .= m.nzval
+  nonzeros(mat) .= nonzeros(m)
+end
+
+function sparse_from_coo(::Type{<:SparseMatrixCSC}, args...)
+  sparse(args...)
+end
+
+function sparse_from_coo(::Type{<:SparseMatrixCSR}, args...)
+  sparsecsr(args...)
 end
 
+function sparse_from_coo(::Type{<:SymSparseMatrixCSR}, args...)
+  symsparsecsr(args...)
+end
+
+
 end # module

src/FESpaces/FEOperators.jl

@@ -4,6 +4,7 @@ using Gridap
 using Gridap.Helpers
 
 using LinearAlgebra
+using SparseArrays
 
 import Gridap: apply
 export apply!
@@ -187,6 +188,16 @@ function LinearFEOperator(
   LinearFEOperator(testfesp,trialfesp,assem,terms...)
 end
 
+function LinearFEOperator(
+  ::Type{M},
+  testfesp::FESpaceLike,
+  trialfesp::FESpaceLike,
+  terms::Vararg{<:AffineFETerm}) where {M}
+
+  assem = SparseMatrixAssembler(M,testfesp,trialfesp)
+  LinearFEOperator(testfesp,trialfesp,assem,terms...)
+end
+
 function LinearFEOperator(
   biform::Function,
   liform::Function,
@@ -290,6 +301,17 @@ function NonLinearFEOperator(
   NonLinearFEOperator(testfesp,trialfesp,assem,terms)
 end
 
+function NonLinearFEOperator(
+  ::Type{M},
+  testfesp::FESpaceLike,
+  trialfesp::FESpaceLike,
+  terms::Vararg{<:FETerm}) where {M}
+
+  assem = SparseMatrixAssembler(M,testfesp,trialfesp)
+  NonLinearFEOperator(testfesp,trialfesp,assem,terms...)
+end
+
+
 function NonLinearFEOperator(
   res::Function,
   jac::Function,

src/MultiField/MultiAssemblers.jl

@@ -5,13 +5,16 @@ using Gridap.Helpers
 using Gridap.MultiFESpaces: _compute_offsets
 
 using SparseArrays
+using SparseMatricesCSR
 
 export MultiAssembler
 export MultiSparseMatrixAssembler
 
 import Gridap.Assemblers: assemble
 import Gridap.Assemblers: assemble!
 import Gridap.Assemblers: SparseMatrixAssembler
+import Gridap.Assemblers: SparseMatrixAssembler
+using  Gridap.Assemblers: sparse_from_coo
 
 abstract type MultiAssembler{M<:AbstractMatrix,V<:AbstractVector} end
 
@@ -70,34 +73,70 @@ end
 """
 MultiAssembler that produces SparseMatrices from the SparseArrays package
 """
-struct MultiSparseMatrixAssembler{E} <: MultiAssembler{SparseMatrixCSC{E,Int},Vector{E}}
+struct MultiSparseMatrixAssembler{E,M} <: MultiAssembler{M,Vector{E}}
   testfesps::MultiFESpace{E}
   trialfesps::MultiFESpace{E}
+
+  function MultiSparseMatrixAssembler(
+    ::Type{M},
+    testfesp::MultiFESpace{E},
+    trialfesp::MultiFESpace{E}) where {M<:AbstractSparseMatrix,E}
+    new{E,M}(testfesp,trialfesp) 
+  end
+end
+
+function MultiSparseMatrixAssembler(
+  testfesp::MultiFESpace{E},
+  trialfesp::MultiFESpace{E}) where {E}
+  MultiSparseMatrixAssembler(SparseMatrixCSC, testfesps, trialfesp)
 end
 
 function MultiSparseMatrixAssembler(
+  ::Type{M},
   testfesps::Vector{<:FESpaceWithDirichletData},
-  trialfesps::Vector{<:FESpaceWithDirichletData})
+  trialfesps::Vector{<:FESpaceWithDirichletData}) where {M<:AbstractSparseMatrix}
   V = MultiFESpace(testfesps)
   U = MultiFESpace(trialfesps)
-  MultiSparseMatrixAssembler(V,U)
+  MultiSparseMatrixAssembler(M,V,U)
+end
+
+function MultiSparseMatrixAssembler(
+  testfesps::Vector{<:FESpaceWithDirichletData},
+  trialfesps::Vector{<:FESpaceWithDirichletData})
+  MultiSparseMatrixAssembler(SparseMatrixCSC,V,U)
+end
+
+function SparseMatrixAssembler(
+  ::Type{M},
+  testfesps::Vector{<:FESpaceWithDirichletData},
+  trialfesps::Vector{<:FESpaceWithDirichletData}) where {M<:AbstractSparseMatrix}
+  MultiSparseMatrixAssembler(M,testfesps,trialfesps)
 end
 
 function SparseMatrixAssembler(
   testfesps::Vector{<:FESpaceWithDirichletData},
   trialfesps::Vector{<:FESpaceWithDirichletData})
-  MultiSparseMatrixAssembler(testfesps,trialfesps)
+  MultiSparseMatrixAssembler(SparseMatrixCSC,testfesps,trialfesps)
 end
 
 function MultiSparseMatrixAssembler(
-  testfesps::MultiFESpace{E,S}, trialfesps::MultiFESpace{E,S}) where {E,S}
+  ::Type{M},
+  testfesps::MultiFESpace{E,S}, 
+  trialfesps::MultiFESpace{E,S}) where {M<:AbstractSparseMatrix,E,S}
   @notimplementedif S != ConsequtiveMultiFieldStyle
-  MultiSparseMatrixAssembler{E}(testfesps,trialfesps)
+  MultiSparseMatrixAssembler{E,M}(M,testfesps,trialfesps)
+end
+
+function MultiSparseMatrixAssembler(
+  testfesps::MultiFESpace{E,S}, 
+  trialfesps::MultiFESpace{E,S}) where {E,S}
+  @notimplementedif S != ConsequtiveMultiFieldStyle
+  MultiSparseMatrixAssembler(SparseMatrixCSC,testfesps,trialfesps)
 end
 
 function assemble(
   this::MultiSparseMatrixAssembler{E},
-  allvals::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where E
+  allvals::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E}
 
   n = num_free_dofs(this.testfesps)
   vec = zeros(E,n)
@@ -108,7 +147,7 @@ end
 function assemble!(
   vec::Vector{E},
   this::MultiSparseMatrixAssembler{E},
-  allmcv::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where E
+  allmcv::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E}
 
   vec .= zero(E)
   V = this.testfesps
@@ -123,8 +162,8 @@ function assemble!(
 end
 
 function assemble(
-  this::MultiSparseMatrixAssembler{E},
-  allmcm::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}}) where {E}
+  this::MultiSparseMatrixAssembler{E,M},
+  allmcm::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}}) where {E,M}
 
   V = this.testfesps
   U = this.trialfesps
@@ -140,21 +179,24 @@ function assemble(
     mf_vals = apply_constraints_cols(U,mf_vals,cellids_col)
     mf_cols = reindex(_mf_cols, cellids_col)
     i_to_fieldid = mf_vals.fieldids
-    _assemble_mat!(
+    _assemble_mat!(M,
       aux_row,aux_col,aux_val,mf_rows,mf_cols,mf_vals,
       i_to_fieldid,offsets_row,offsets_col)
   end
-  sparse(aux_row, aux_col, aux_val)
+  num_rows = num_free_dofs(this.testfesps)
+  num_cols = num_free_dofs(this.trialfesps)
+  finalize_coo!(M,aux_row,aux_col,aux_val,num_rows,num_cols)
+  sparse_from_coo(M,aux_row, aux_col, aux_val)
 end
 
 function assemble!(
-  mat::SparseMatrixCSC{E},
+  mat::AbstractSparseMatrix{E},
   this::MultiSparseMatrixAssembler{E},
   allvals::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}}) where E
   # This routine can be optimized a lot taking into a count the sparsity graph of mat
   # For the moment we create an intermediate matrix and then transfer the nz values
   m = assemble(this,allvals...)
-  mat.nzval .= m.nzval
+  nonzeros(mat) .= nonzeros(m)
 end
 
 function _assemble_vec!(vec,mf_rows,mf_vals,i_to_fieldid,offsets)
@@ -172,9 +214,9 @@ function _assemble_vec!(vec,mf_rows,mf_vals,i_to_fieldid,offsets)
   end
 end
 
-function _assemble_mat!(
+function _assemble_mat!(::Type{M},
   aux_row,aux_col,aux_val,mf_rows,mf_cols,mf_vals,
-  i_to_fieldid,offsets_row,offsets_col)
+  i_to_fieldid,offsets_row,offsets_col) where {M<:AbstractSparseMatrix}
 
   for (mf_rows_c,mf_cols_c,mf_vals_c) in zip(mf_rows,mf_cols,mf_vals)
     for (i,vals_c) in enumerate(mf_vals_c)
@@ -187,9 +229,8 @@ function _assemble_mat!(
         if gidcol > 0
           for (lidrow,gidrow) in enumerate(rows_c)
             if gidrow > 0
-              push!(aux_row, gidrow+offset_row)
-              push!(aux_col, gidcol+offset_col)
-              push!(aux_val, vals_c[lidrow,lidcol])
+              push_coo!(M,aux_row, aux_col, aux_val, 
+                    gidrow+offset_row, gidcol+offset_col, vals_c[lidrow,lidcol])
             end
           end
         end