From c8fb75cd347068363d2bfb708071869c79f1368a Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Thu, 14 Nov 2019 11:29:12 +0100 Subject: [PATCH 1/5] update dependencies and compatibility to SparseMatricesCSR v0.4.0 package --- Manifest.toml | 48 ++++++++++++++++++++++++------------------------ Project.toml | 2 +- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 60e834f..9d0fc49 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -16,10 +16,10 @@ uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee" version = "0.8.10" [[BinaryProvider]] -deps = ["Libdl", "Logging", "SHA"] -git-tree-sha1 = "c7361ce8a2129f20b0e05a89f7070820cfed6648" +deps = ["Libdl", "SHA"] +git-tree-sha1 = "5b08ed6036d9d3f0ee6369410b830f8873d4024c" uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" -version = "0.5.6" +version = "0.5.8" [[Calculus]] deps = ["Compat"] @@ -52,10 +52,10 @@ uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" version = "2.2.0" [[DataStructures]] -deps = ["InteractiveUtils", "OrderedCollections", "Random", "Serialization", "Test"] -git-tree-sha1 = "ca971f03e146cf144a9e2f2ce59674f5bf0e8038" +deps = ["InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "1fe8fad5fc84686dcbc674aa255bc867a64f8132" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.15.0" +version = "0.17.5" [[Dates]] deps = ["Printf"] @@ -67,9 +67,9 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" [[DiffEqDiffTools]] deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] -git-tree-sha1 = "21b855cb29ec4594f9651e0e9bdc0cdcfdcd52c1" +git-tree-sha1 = "81edfb3a8b55154772bb6080b5db40868e1778ed" uuid = "01453d9d-ee7c-5054-8395-0335cb756afa" -version = "1.3.0" +version = "1.4.0" [[DiffResults]] deps = ["Compat", "StaticArrays"] @@ -78,10 +78,10 @@ uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" version = "0.0.4" [[DiffRules]] -deps = ["Random", "Test"] -git-tree-sha1 = "dc0869fb2f5b23466b32ea799bd82c76480167f7" +deps = ["NaNMath", "Random", "SpecialFunctions"] +git-tree-sha1 = "f734b5f6bc9c909027ef99f6d91d5d9e4b111eed" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "0.0.10" +version = "0.1.0" [[Distances]] deps = ["LinearAlgebra", "Statistics"] @@ -100,10 +100,10 @@ uuid = "442a2c76-b920-505d-bb47-c5924d526838" version = "0.3.3" [[ForwardDiff]] -deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"] -git-tree-sha1 = "4c4d727f1b7e0092134fabfab6396b8945c1ea5b" +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"] +git-tree-sha1 = "4407e7b76999eca2646abdb68203bd4302476168" uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "0.10.3" +version = "0.10.6" [[Gridap]] deps = ["Combinatorics", "FastGaussQuadrature", "JSON", "LineSearches", "LinearAlgebra", "NLsolve", "QuadGK", "Reexport", "SparseArrays", "StaticArrays", "TensorPolynomialBases", "TensorValues", "Test", "UnstructuredGrids", "WriteVTK"] @@ -161,9 +161,9 @@ version = "7.5.0" [[NLsolve]] deps = ["DiffEqDiffTools", "Distances", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"] -git-tree-sha1 = "c9578878f56f425a2160f5b436c7f79a154d154c" +git-tree-sha1 = "be25a8486c7be4b075165315e201fa5d199073f1" uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" -version = "4.1.0" +version = "4.2.0" [[NaNMath]] deps = ["Compat"] @@ -185,9 +185,9 @@ version = "0.12.0" [[Parsers]] deps = ["Dates", "Test"] -git-tree-sha1 = "ef0af6c8601db18c282d092ccbd2f01f3f0cd70b" +git-tree-sha1 = "c56ecb484f286639f161e712b8311f5ab77e8d32" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "0.3.7" +version = "0.3.8" [[Pkg]] deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] @@ -195,9 +195,9 @@ uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [[Polynomials]] deps = ["LinearAlgebra", "RecipesBase"] -git-tree-sha1 = "f7c0c07e82798aef542d60a6e6e85e39f4590750" +git-tree-sha1 = "ae71c2329790af97b7682b11241b3609e4d48626" uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" -version = "0.5.3" +version = "0.6.0" [[Printf]] deps = ["Unicode"] @@ -205,9 +205,9 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" [[QuadGK]] deps = ["DataStructures", "LinearAlgebra"] -git-tree-sha1 = "438630b843c210b375b2a246329200c113acc61b" +git-tree-sha1 = "1af46bf083b9630a5b27d4fd94f496c5fca642a8" uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -version = "2.1.0" +version = "2.1.1" [[REPL]] deps = ["InteractiveUtils", "Markdown", "Sockets"] @@ -253,9 +253,9 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[SparseMatricesCSR]] deps = ["LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "c32c12d65a48919e6c34582ac0be99e4a5faa7a6" +git-tree-sha1 = "c714a2ad5ef5472b28523fa89e77323173eaa135" uuid = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" -version = "0.3.0" +version = "0.4.0" [[SpecialFunctions]] deps = ["BinDeps", "BinaryProvider", "Libdl"] diff --git a/Project.toml b/Project.toml index 828b180..44f0080 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,7 @@ SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" [compat] Gridap = "0.5.0" -SparseMatricesCSR = "0.3.0" +SparseMatricesCSR = "0.4.0" julia = "1.1" [extras] From 043f2d6e590823930cb91a36f5fd9cafbce0a0ab Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Thu, 14 Nov 2019 11:30:05 +0100 Subject: [PATCH 2/5] Add parameters to control Pardiso 0/1-based indexing --- src/PardisoParameters.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/PardisoParameters.jl b/src/PardisoParameters.jl index 2624842..b3250d3 100644 --- a/src/PardisoParameters.jl +++ b/src/PardisoParameters.jl @@ -101,6 +101,13 @@ const PARDISO_SOLVE_LINEAR_SYSTEM = 0 # linear system. const PARDISO_SOLVE_CONJUGATE_TRANSPOSED_SYSTEM = 1 # conjugate transposed system const PARDISO_SOLVE_TRANSPOSED_SYSTEM = 2 # transposed system +############################################################### +# Parameter values for IPARM[35] +# IPARM_ONE_OR_ZERO_BASED_INDEXING +############################################################### +const PARDISO_ONE_BASED_INDEXING = 0 # One-based indexing +const PARDISO_ZERO_BASED_INDEXING = 1 # Zero-based indexing + ############################################################### # Iparm parameters # From 4e6e4d1196b9650ee3f6ad8827a13874a81ef7de Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Thu, 14 Nov 2019 11:32:15 +0100 Subject: [PATCH 3/5] update LinearSolver implementation to support SparseMatricesCSR v.0.4.0 (0/1-based sparse matrices indexing) --- src/LinearSolver.jl | 152 +++++++++++++++++++++++-------------------- test/LinearSolver.jl | 80 ++++++++++++----------- 2 files changed, 126 insertions(+), 106 deletions(-) diff --git a/src/LinearSolver.jl b/src/LinearSolver.jl index 7501983..be36d7f 100644 --- a/src/LinearSolver.jl +++ b/src/LinearSolver.jl @@ -54,9 +54,8 @@ Gridap SymbolicSetup implementation for Intel Pardiso MKL solver. """ mutable struct PardisoSymbolicSetup{T,Ti} <: SymbolicSetup phase :: Int - mat :: SparseMatrixCSC{T,Ti} + mat :: AbstractSparseMatrix{T,Ti} solver :: PardisoSolver - transpose :: Bool end """ @@ -65,9 +64,8 @@ Gridap NumericalSetup implementation for Intel Pardiso MKL solver. """ mutable struct PardisoNumericalSetup{T,Ti} <: NumericalSetup phase :: Int - mat :: SparseMatrixCSC{T,Ti} + mat :: AbstractSparseMatrix{T,Ti} solver :: PardisoSolver - transpose :: Bool end """ @@ -106,91 +104,95 @@ PardisoSolver(mtype, iparm, msglvl) = """ - function PardisoSymbolicSetup(phase::Integer, mat::AbstractSparseMatrix{T,Ti}, solver::PardisoSolver) where {T,Ti} + function build_PardisoSymbolicSetup(phase::Integer, mat::AbstractSparseMatrix{T,Ti}, solver::PardisoSolver) where {T,Ti} PardisoSymbolicSetup constructor overloading with default values. Returns a PardisoSymbolicSetup from a given AbstractSparseMatrix. """ -function PardisoSymbolicSetup(phase::Integer, +function build_PardisoSymbolicSetup(phase::Integer, mat::AbstractSparseMatrix{T,Ti}, solver::PardisoSolver) where {T,Ti} - PardisoSymbolicSetup(phase,SparseMatrixCSC{T,Ti}(mat),solver,true) + PardisoSymbolicSetup(phase,SparseMatrixCSC{T,Ti}(mat),solver) end """ - function PardisoNumericalSetup(phase::Integer, mat::AbstractSparseMatrix{T,Ti}, solver::PardisoSolver) where {T,Ti} + function build_PardisoNumericalSetup(phase::Integer, mat::AbstractSparseMatrix{T,Ti}, solver::PardisoSolver) where {T,Ti} PardisoNumericalSetup constructor overloading with default values. Returns a PardisoNumericalSetup from a given AbstractSparseMatrix. """ -function PardisoNumericalSetup(phase::Integer, +function build_PardisoNumericalSetup(phase::Integer, mat::AbstractSparseMatrix{T,Ti}, solver::PardisoSolver) where {T,Ti} - PardisoNumericalSetup(phase,SparseMatrixCSC{T,Ti}(mat),solver,true) + PardisoNumericalSetup(phase,SparseMatrixCSC{T,Ti}(mat),solver) end """ - function PardisoSymbolicSetup(phase::Integer, mat::SparseMatrixCSC{T,Ti}, solver::PardisoSolver) where {T,Ti} + function build_PardisoSymbolicSetup(phase::Integer, mat::SparseMatrixCSC{T,Ti}, solver::PardisoSolver) where {T,Ti} PardisoSymbolicSetup constructor overloading. Returns a PardisoSymbolicSetup from a given SparseMatrixCSC. """ -function PardisoSymbolicSetup(phase::Integer, +function build_PardisoSymbolicSetup(phase::Integer, mat::SparseMatrixCSC{T,Ti}, solver::PardisoSolver) where {T,Ti} - PardisoSymbolicSetup(phase,mat,solver,true) + PardisoSymbolicSetup(phase,mat,solver) end """ - function PardisoNumericalSetup(phase::Integer, mat::SparseMatrixCSC{T,Ti}, solver::PardisoSolver) where {T,Ti} + function build_PardisoNumericalSetup(phase::Integer, mat::SparseMatrixCSC{T,Ti}, solver::PardisoSolver) where {T,Ti} PardisoNumericalSetup constructor overloading. Returns a PardisoNumericalSetup from a given SparseMatrixCSC. """ -function PardisoNumericalSetup(phase::Integer, +function build_PardisoNumericalSetup(phase::Integer, mat::SparseMatrixCSC{T,Ti}, solver::PardisoSolver) where {T,Ti} - PardisoNumericalSetup(phase,mat,solver,true) + PardisoNumericalSetup(phase,mat,solver) end """ - function PardisoSymbolicSetup(phase::Integer, mat::SparseMatrixCSR{T,Ti}, solver::PardisoSolver) where {T,Ti} + function build_PardisoSymbolicSetup(phase::Integer, mat::SparseMatrixCSR{T,Ti}, solver::PardisoSolver) where {T,Ti} PardisoSymbolicSetup constructor overloading. Returns a PardisoSymbolicSetup from a given SparseMatrixCSR. """ -function PardisoSymbolicSetup(phase::Integer, - mat::SparseMatrixCSR{T,Ti}, - solver::PardisoSolver) where {T,Ti} - PardisoSymbolicSetup(phase,mat.transpose,solver,false) +function build_PardisoSymbolicSetup(phase::Integer, + mat::SparseMatrixCSR{Bi}, + solver::PardisoSolver) where {Bi} + solver.iparm[IPARM_ONE_OR_ZERO_BASED_INDEXING] = Bi == 0 ? PARDISO_ZERO_BASED_INDEXING : PARDISO_ONE_BASED_INDEXING + PardisoSymbolicSetup(phase,mat,solver) end """ - function PardisoNumericalSetup(phase::Integer, mat::SparseMatrixCSR{T,Ti}, solver::PardisoSolver) where {T,Ti} + function build_PardisoNumericalSetup(phase::Integer, mat::SparseMatrixCSR{T,Ti}, solver::PardisoSolver) where {T,Ti} PardisoNumericalSetup constructor overloading. Returns a PardisoNumericalSetup from a given SparseMatrixCSR. """ -function PardisoNumericalSetup(phase::Integer, - mat::SparseMatrixCSR{T,Ti}, - solver::PardisoSolver) where {T,Ti} - PardisoNumericalSetup(phase,mat.transpose,solver,false) +function build_PardisoNumericalSetup(phase::Integer, + mat::SparseMatrixCSR{Bi}, + solver::PardisoSolver) where {Bi} + solver.iparm[IPARM_ONE_OR_ZERO_BASED_INDEXING] = Bi == 0 ? PARDISO_ZERO_BASED_INDEXING : PARDISO_ONE_BASED_INDEXING + PardisoNumericalSetup(phase,mat,solver) end """ - function PardisoSymbolicSetup(phase::Integer, mat::SymSparseMatrixCSR{T,Ti}, solver::PardisoSolver) where {T,Ti} + function build_PardisoSymbolicSetup(phase::Integer, mat::SymSparseMatrixCSR{T,Ti}, solver::PardisoSolver) where {T,Ti} PardisoSymbolicSetup constructor overloading. Returns a PardisoSymbolicSetup from a given SymSparseMatrixCSR. """ -function PardisoSymbolicSetup(phase::Integer, - mat::SymSparseMatrixCSR{T,Ti}, - solver::PardisoSolver) where {T,Ti} - PardisoSymbolicSetup(phase,mat.lowertrian,solver,false) +function build_PardisoSymbolicSetup(phase::Integer, + mat::SymSparseMatrixCSR{Bi}, + solver::PardisoSolver) where {Bi} + solver.iparm[IPARM_ONE_OR_ZERO_BASED_INDEXING] = Bi == 0 ? PARDISO_ZERO_BASED_INDEXING : PARDISO_ONE_BASED_INDEXING + PardisoSymbolicSetup(phase,mat.uppertrian,solver) end """ - function PardisoNumericalSetup(phase::Integer, mat::SymSparseMatrixCSR{T,Ti}, solver::PardisoSolver) where {T,Ti} + function build_PardisoNumericalSetup(phase::Integer, mat::SymSparseMatrixCSR{T,Ti}, solver::PardisoSolver) where {T,Ti} PardisoNumericalSetup constructor overloading. Returns a PardisoNumericalSetup from a given SparseMatrixCSR. """ -function PardisoNumericalSetup(phase::Integer, - mat::SymSparseMatrixCSR{T,Ti}, - solver::PardisoSolver) where {T,Ti} - PardisoNumericalSetup(phase,mat.lowertrian,solver,false) +function build_PardisoNumericalSetup(phase::Integer, + mat::SymSparseMatrixCSR{Bi}, + solver::PardisoSolver) where {Bi} + solver.iparm[IPARM_ONE_OR_ZERO_BASED_INDEXING] = Bi == 0 ? PARDISO_ZERO_BASED_INDEXING : PARDISO_ONE_BASED_INDEXING + PardisoNumericalSetup(phase,mat.uppertrian,solver) end """ @@ -202,17 +204,18 @@ function symbolic_setup( ps::PardisoSolver{Ti}, mat::M) where {T<:Float64,Ti<:Int32,M<:AbstractSparseMatrix{T,Ti}} - pss = PardisoSymbolicSetup(GridapPardiso.PHASE_ANALYSIS, mat, ps) + pss = build_PardisoSymbolicSetup(GridapPardiso.PHASE_ANALYSIS, mat, ps) + m,n = size(pss.mat) err = pardiso!( pss.solver.pt, # Handle to internal data structure. The entries must be set to zero prior to the first call to pardiso maxfct, # Maximum number of factors with identical sparsity structure that must be kept in memory at the same time mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. pss.solver.mtype, # Defines the matrix type, which influences the pivoting method pss.phase, # Controls the execution of the solver (11 == Analysis) - pss.mat.n, # Number of equations in the sparse linear systems of equations - pss.mat.nzval, # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja - pss.mat.colptr, # Pointers to columns in CSR format - pss.mat.rowval, # Column indices of the CSR sparse matrix + n, # Number of equations in the sparse linear systems of equations + nonzeros(pss.mat), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja + getptr(pss.mat), # Pointers to columns in CSR format + getindices(pss.mat), # Column indices of the CSR sparse matrix Vector{Ti}(), # Permutation vector nrhs, # Number of right-hand sides that need to be solved for pss.solver.iparm, # This array is used to pass various parameters to Intel MKL PARDISO @@ -232,17 +235,18 @@ function symbolic_setup( ps::PardisoSolver{Ti}, mat::M) where {T<:Float64,Ti<:Int64,M<:AbstractSparseMatrix{T,Ti}} - pss = PardisoSymbolicSetup(GridapPardiso.PHASE_ANALYSIS, mat, ps) + pss = build_PardisoSymbolicSetup(GridapPardiso.PHASE_ANALYSIS, mat, ps) + m,n = size(pss.mat) err = pardiso_64!( pss.solver.pt, # Handle to internal data structure. The entries must be set to zero prior to the first call to pardiso maxfct, # Maximum number of factors with identical sparsity structure that must be kept in memory at the same time mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. pss.solver.mtype, # Defines the matrix type, which influences the pivoting method pss.phase, # Controls the execution of the solver (11 == Analysis) - pss.mat.n, # Number of equations in the sparse linear systems of equations - pss.mat.nzval, # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja - pss.mat.colptr, # Pointers to columns in CSR format - pss.mat.rowval, # Column indices of the CSR sparse matrix + n, # Number of equations in the sparse linear systems of equations + nonzeros(pss.mat), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja + getptr(pss.mat), # Pointers to columns in CSR format + getindices(pss.mat), # Column indices of the CSR sparse matrix Vector{Ti}(), # Permutation vector nrhs, # Number of right-hand sides that need to be solved for pss.solver.iparm, # This array is used to pass various parameters to Intel MKL PARDISO @@ -262,13 +266,14 @@ function symbolic_setup_finalize( pss::PardisoSymbolicSetup{T,Ti}) where {T,Ti<:Int32} pss.phase = GridapPardiso.PHASE_RELEASE_INTERNAL_MEMORY + m,n = size(pss.mat) err = pardiso!( pss.solver.pt, # Handle to internal data structure. The entries must be set to zero prior to the first call to pardiso maxfct, # Maximum number of factors with identical sparsity structure that must be kept in memory at the same time mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. pss.solver.mtype, # Defines the matrix type, which influences the pivoting method pss.phase, # Controls the execution of the solver (11 == Analysis) - pss.mat.n, # Number of equations in the sparse linear systems of equations + n, # Number of equations in the sparse linear systems of equations Vector{T}(), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja Vector{Ti}(), # Pointers to columns in CSR format Vector{Ti}(), # Column indices of the CSR sparse matrix @@ -290,13 +295,14 @@ function symbolic_setup_finalize( pss::PardisoSymbolicSetup{T,Ti}) where {T,Ti<:Int64} pss.phase = GridapPardiso.PHASE_RELEASE_INTERNAL_MEMORY + m,n = size(pss.mat) err = pardiso_64!( pss.solver.pt, # Handle to internal data structure. The entries must be set to zero prior to the first call to pardiso maxfct, # Maximum number of factors with identical sparsity structure that must be kept in memory at the same time mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. pss.solver.mtype, # Defines the matrix type, which influences the pivoting method pss.phase, # Controls the execution of the solver (11 == Analysis) - pss.mat.n, # Number of equations in the sparse linear systems of equations + n, # Number of equations in the sparse linear systems of equations Vector{T}(), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja Vector{Ti}(), # Pointers to columns in CSR format Vector{Ti}(), # Column indices of the CSR sparse matrix @@ -318,15 +324,17 @@ function numerical_setup!( pns::PardisoNumericalSetup{T,Ti}, mat::M) where {T<:Float64,Ti<:Int32,M<:AbstractSparseMatrix{T,Ti}} + m,n = size(pns.mat) + err = pardiso!( pns.solver.pt, # Handle to internal data structure. The entries must be set to zero prior to the first call to pardiso maxfct, # Maximum number of factors with identical sparsity structure that must be kept in memory at the same time mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. pns.solver.mtype, # Defines the matrix type, which influences the pivoting method pns.phase, # Controls the execution of the solver (11 == Analysis) - pns.mat.n, # Number of equations in the sparse linear systems of equations - pns.mat.nzval, # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja - pns.mat.colptr, # Pointers to columns in CSR format - pns.mat.rowval, # Column indices of the CSR sparse matrix + n, # Number of equations in the sparse linear systems of equations + nonzeros(pns.mat), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja + getptr(pns.mat), # Pointers to columns in CSR format + getindices(pns.mat), # Column indices of the CSR sparse matrix Vector{Ti}(), # Permutation vector nrhs, # Number of right-hand sides that need to be solved for pns.solver.iparm, # This array is used to pass various parameters to Intel MKL PARDISO @@ -346,15 +354,17 @@ function numerical_setup!( pns::PardisoNumericalSetup{T,Ti}, mat::M) where {T<:Float64,Ti<:Int64,M<:AbstractSparseMatrix{T,Ti}} + m,n = size(pns.mat) + err = pardiso_64!( pns.solver.pt, # Handle to internal data structure. The entries must be set to zero prior to the first call to pardiso maxfct, # Maximum number of factors with identical sparsity structure that must be kept in memory at the same time mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. pns.solver.mtype, # Defines the matrix type, which influences the pivoting method pns.phase, # Controls the execution of the solver (11 == Analysis) - pns.mat.n, # Number of equations in the sparse linear systems of equations - pns.mat.nzval, # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja - pns.mat.colptr, # Pointers to columns in CSR format - pns.mat.rowval, # Column indices of the CSR sparse matrix + n, # Number of equations in the sparse linear systems of equations + nonzeros(pns.mat), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja + getptr(pns.mat), # Pointers to columns in CSR format + getindices(pns.mat), # Column indices of the CSR sparse matrix Vector{Ti}(), # Permutation vector nrhs, # Number of right-hand sides that need to be solved for pns.solver.iparm, # This array is used to pass various parameters to Intel MKL PARDISO @@ -374,13 +384,14 @@ function numerical_setup_finalize!( pns::PardisoNumericalSetup{T,Ti}) where {T, Ti<:Int32} pns.phase = GridapPardiso.PHASE_RELEASE_INTERNAL_MEMORY + m,n = size(pns.mat) err = pardiso!( pns.solver.pt, # Handle to internal data structure. The entries must be set to zero prior to the first call to pardiso maxfct, # Maximum number of factors with identical sparsity structure that must be kept in memory at the same time mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. pns.solver.mtype, # Defines the matrix type, which influences the pivoting method pns.phase, # Controls the execution of the solver (11 == Analysis) - pns.mat.n, # Number of equations in the sparse linear systems of equations + n, # Number of equations in the sparse linear systems of equations Vector{T}(), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja Vector{Ti}(), # Pointers to columns in CSR format Vector{Ti}(), # Column indices of the CSR sparse matrix @@ -402,13 +413,14 @@ function numerical_setup_finalize!( pns::PardisoNumericalSetup{T,Ti}) where {T, Ti<:Int64} pns.phase = GridapPardiso.PHASE_RELEASE_INTERNAL_MEMORY + m,n = size(pns.mat) err = pardiso_64!( pns.solver.pt, # Handle to internal data structure. The entries must be set to zero prior to the first call to pardiso maxfct, # Maximum number of factors with identical sparsity structure that must be kept in memory at the same time mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. pns.solver.mtype, # Defines the matrix type, which influences the pivoting method pns.phase, # Controls the execution of the solver (11 == Analysis) - pns.mat.n, # Number of equations in the sparse linear systems of equations + n, # Number of equations in the sparse linear systems of equations Vector{T}(), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja Vector{Ti}(), # Pointers to columns in CSR format Vector{Ti}(), # Column indices of the CSR sparse matrix @@ -429,7 +441,7 @@ Create the PardisoSymbolicSetup object and use Intel Pardiso MKL to perform the function numerical_setup( pss::PardisoSymbolicSetup{T,Ti}, mat::M) where {T<:Float64,Ti<:Integer,M<:AbstractSparseMatrix{T,Ti}} - pns = PardisoNumericalSetup(GridapPardiso.PHASE_NUMERICAL_FACTORIZATION, mat, pss.solver) + pns = build_PardisoNumericalSetup(GridapPardiso.PHASE_NUMERICAL_FACTORIZATION, mat, pss.solver) return numerical_setup!(pns, mat) end @@ -444,10 +456,11 @@ function solve!( b::AbstractVector{T}) where {T<:Float64,Ti<:Int32} phase = GridapPardiso.PHASE_SOLVE_ITERATIVE_REFINEMENT + m,n = size(ns.mat) # Here we assume that the users defines iparm for CSR matrix iparmcopy = copy(ns.solver.iparm) - if ns.transpose + if hascolmajororder(ns.mat) if ns.solver.iparm[IPARM_TRANSPOSED_OR_CONJUGATED_TRANSPOSED] == PARDISO_SOLVE_LINEAR_SYSTEM iparmcopy[IPARM_TRANSPOSED_OR_CONJUGATED_TRANSPOSED] = PARDISO_SOLVE_TRANSPOSED_SYSTEM elseif ns.solver.iparm[IPARM_TRANSPOSED_OR_CONJUGATED_TRANSPOSED] == PARDISO_SOLVE_TRANSPOSED_SYSTEM @@ -465,10 +478,10 @@ function solve!( mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. ns.solver.mtype, # Defines the matrix type, which influences the pivoting method phase, # Controls the execution of the solver (11 == Analysis) - ns.mat.n, # Number of equations in the sparse linear systems of equations - ns.mat.nzval, # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja - ns.mat.colptr, # Pointers to columns in CSR format - ns.mat.rowval, # Column indices of the CSR sparse matrix + n, # Number of equations in the sparse linear systems of equations + nonzeros(ns.mat), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja + getptr(ns.mat), # Pointers to columns in CSR format + getindices(ns.mat), # Column indices of the CSR sparse matrix Vector{Ti}(), # Permutation vector nrhs, # Number of right-hand sides that need to be solved for iparmcopy, # This array is used to pass various parameters to Intel MKL PARDISO @@ -490,10 +503,11 @@ function solve!( b::AbstractVector{T}) where {T<:Float64,Ti<:Int64} phase = GridapPardiso.PHASE_SOLVE_ITERATIVE_REFINEMENT + m,n = size(ns.mat) # Here we assume that the users defines iparm for CSR matrix iparmcopy = copy(ns.solver.iparm) - if ns.transpose + if hascolmajororder(ns.mat) if ns.solver.iparm[IPARM_TRANSPOSED_OR_CONJUGATED_TRANSPOSED] == PARDISO_SOLVE_LINEAR_SYSTEM iparmcopy[IPARM_TRANSPOSED_OR_CONJUGATED_TRANSPOSED] = PARDISO_SOLVE_TRANSPOSED_SYSTEM elseif ns.solver.iparm[IPARM_TRANSPOSED_OR_CONJUGATED_TRANSPOSED] == PARDISO_SOLVE_TRANSPOSED_SYSTEM @@ -511,10 +525,10 @@ function solve!( mnum, # Actual matrix for the solution phase. The value must be: 1 <= mnum <= maxfct. ns.solver.mtype, # Defines the matrix type, which influences the pivoting method phase, # Controls the execution of the solver (11 == Analysis) - ns.mat.n, # Number of equations in the sparse linear systems of equations - ns.mat.nzval, # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja - ns.mat.colptr, # Pointers to columns in CSR format - ns.mat.rowval, # Column indices of the CSR sparse matrix + n, # Number of equations in the sparse linear systems of equations + nonzeros(ns.mat), # Contains the non-zero elements of the coefficient matrix A corresponding to the indices in ja + getptr(ns.mat), # Pointers to columns in CSR format + getindices(ns.mat), # Column indices of the CSR sparse matrix Vector{Ti}(), # Permutation vector nrhs, # Number of right-hand sides that need to be solved for iparmcopy, # This array is used to pass various parameters to Intel MKL PARDISO diff --git a/test/LinearSolver.jl b/test/LinearSolver.jl index 93147b4..409a6d1 100644 --- a/test/LinearSolver.jl +++ b/test/LinearSolver.jl @@ -95,29 +95,14 @@ rows = 5 cols = 5 # pardiso! -I = Vector{Int32}(); J = Vector{Int32}(); V = Vector{Float64}() -for (ik, jk, vk) in zip(I_,J_,V_) - push_coo!(SparseMatrixCSR,I,J,V,ik,jk,vk) -end -finalize_coo!(SparseMatrixCSR,I,J,V,rows, cols) -A = sparsecsr(I,J,V,rows, cols) -b = ones(size(A)[2]) -x = similar(b) -ps = PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE) -ss = symbolic_setup(ps, A) -ns = numerical_setup(ss, A) -solve!(x, ns, b) -@test maximum(abs.(A*x-b)) < tol -test_linear_solver(ps, A, b, x) -if Int == Int64 - # pardiso_64! - I = Vector{Int64}(); J = Vector{Int64}(); V = Vector{Float64}() +for Bi in (0,1) + I = Vector{Int32}(); J = Vector{Int32}(); V = Vector{Float64}() for (ik, jk, vk) in zip(I_,J_,V_) push_coo!(SparseMatrixCSR,I,J,V,ik,jk,vk) end - finalize_coo!(SparseMatrixCSR,I,J,V,rows, cols) - A = sparsecsr(I,J,V,rows, cols) + finalize_coo!(SparseMatrixCSR,I,J,V,rows,cols) + A = sparsecsr(SparseMatrixCSR{Bi},I,J,V,rows,cols) b = ones(size(A)[2]) x = similar(b) ps = PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE) @@ -126,6 +111,24 @@ if Int == Int64 solve!(x, ns, b) @test maximum(abs.(A*x-b)) < tol test_linear_solver(ps, A, b, x) + + if Int == Int64 + # pardiso_64! + I = Vector{Int64}(); J = Vector{Int64}(); V = Vector{Float64}() + for (ik, jk, vk) in zip(I_,J_,V_) + push_coo!(SparseMatrixCSR,I,J,V,ik,jk,vk) + end + finalize_coo!(SparseMatrixCSR,I,J,V,rows,cols) + A = sparsecsr(SparseMatrixCSR{Bi},I,J,V,rows,cols) + b = ones(size(A)[2]) + x = similar(b) + ps = PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE) + ss = symbolic_setup(ps, A) + ns = numerical_setup(ss, A) + solve!(x, ns, b) + @test maximum(abs.(A*x-b)) < tol + test_linear_solver(ps, A, b, x) + end end ##################################################### @@ -159,29 +162,14 @@ rows = 8 cols = 8 # pardiso! -I = Vector{Int32}(); J = Vector{Int32}(); V = Vector{Float64}() -for (ik, jk, vk) in zip(I_,J_,V_) - push_coo!(SymSparseMatrixCSR,I,J,V,ik,jk,vk) -end -finalize_coo!(SymSparseMatrixCSR,I,J,V,rows,cols) -A = symsparsecsr(I,J,V,rows,cols) -b = ones(size(A)[2]) -x = similar(b) -ps = PardisoSolver(GridapPardiso.MTYPE_REAL_SYMMETRIC_INDEFINITE, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE) -ss = symbolic_setup(ps, A) -ns = numerical_setup(ss, A) -solve!(x, ns, b) -@test maximum(abs.(A*x-b)) < tol -test_linear_solver(ps, A, b, x) -if Int == Int64 - # pardiso_64! - I = Vector{Int64}(); J = Vector{Int64}(); V = Vector{Float64}() +for Bi in (0,1) + I = Vector{Int32}(); J = Vector{Int32}(); V = Vector{Float64}() for (ik, jk, vk) in zip(I_,J_,V_) push_coo!(SymSparseMatrixCSR,I,J,V,ik,jk,vk) end finalize_coo!(SymSparseMatrixCSR,I,J,V,rows,cols) - A = symsparsecsr(I,J,V,rows,cols) + A = symsparsecsr(SymSparseMatrixCSR{Bi},I,J,V,rows,cols) b = ones(size(A)[2]) x = similar(b) ps = PardisoSolver(GridapPardiso.MTYPE_REAL_SYMMETRIC_INDEFINITE, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE) @@ -190,6 +178,24 @@ if Int == Int64 solve!(x, ns, b) @test maximum(abs.(A*x-b)) < tol test_linear_solver(ps, A, b, x) + + if Int == Int64 + # pardiso_64! + I = Vector{Int64}(); J = Vector{Int64}(); V = Vector{Float64}() + for (ik, jk, vk) in zip(I_,J_,V_) + push_coo!(SymSparseMatrixCSR,I,J,V,ik,jk,vk) + end + finalize_coo!(SymSparseMatrixCSR,I,J,V,rows,cols) + A = symsparsecsr(SymSparseMatrixCSR{Bi},I,J,V,rows,cols) + b = ones(size(A)[2]) + x = similar(b) + ps = PardisoSolver(GridapPardiso.MTYPE_REAL_SYMMETRIC_INDEFINITE, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE) + ss = symbolic_setup(ps, A) + ns = numerical_setup(ss, A) + solve!(x, ns, b) + @test maximum(abs.(A*x-b)) < tol + test_linear_solver(ps, A, b, x) + end end end From f55c946552a8f93b26f15dc00d7667e23690c936 Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Thu, 14 Nov 2019 12:37:52 +0100 Subject: [PATCH 4/5] update Gridapt compat to latest version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 44f0080..d0fb9a2 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" [compat] -Gridap = "0.5.0" +Gridap = "0.5.2" SparseMatricesCSR = "0.4.0" julia = "1.1" From 0c37bb09999cfe5145fb1430237e2ea8975952a5 Mon Sep 17 00:00:00 2001 From: victorsndvg Date: Thu, 14 Nov 2019 12:40:31 +0100 Subject: [PATCH 5/5] Prepare version 0.3.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d0fb9a2..8298fb6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GridapPardiso" uuid = "34aa2546-dee6-11e9-014e-739fa02ec06f" authors = ["Francesc Verdugo ", "Víctor Sande "] -version = "0.3.0" +version = "0.3.1" [deps] Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"