Closed
Description
I encounter a problem using LinearSolve to solve a sparse linear system with PardisoJL():
The test case, test_PardisoJL.jl:
using Pkg
Pkg.activate(; temp=true)
Pkg.add("LinearSolve")
Pkg.add("Pardiso")
using LinearAlgebra
using LinearSolve
using Random
using SparseArrays
import Pardiso
include("A.jl")
A = makeA()
b = ones(size(A, 1))
linprob = LinearProblem(A, b)
x1 = LinearSolve.solve(linprob, PardisoJL())
x2 = A\b
@show norm(x1 - x2)
the result is clearly wrong. Matrix A is a 60x60 SparseMatrixCSC constructed in A.jl:
function makeA()
n = 60
colptr = [1, 4, 7, 11, 15, 17, 22, 26, 30, 34, 38, 40, 46, 50, 54, 58, 62, 64, 70, 74, 78, 82, 86, 88, 94, 98, 102, 106, 110, 112, 118, 122, 126, 130, 134, 136, 142, 146, 150, 154, 158, 160, 166, 170, 174, 178, 182, 184, 190, 194, 198, 202, 206, 208, 214, 218, 222, 224, 226, 228, 232]
rowval = [1, 3, 4, 1, 2, 4, 2, 4, 9, 10, 3, 5, 11, 12, 1, 3, 2, 4, 6, 11, 12, 2, 7, 9, 10, 2, 7, 8, 10, 8, 10, 15, 16, 9, 11, 17, 18, 7, 9, 2, 8, 10, 12, 17, 18, 8, 13, 15, 16, 8, 13, 14, 16, 14, 16, 21, 22, 15, 17, 23, 24, 13, 15, 8, 14, 16, 18, 23, 24, 14, 19, 21, 22, 14, 19, 20, 22, 20, 22, 27, 28, 21, 23, 29, 30, 19, 21, 14, 20, 22, 24, 29, 30, 20, 25, 27, 28, 20, 25, 26, 28, 26, 28, 33, 34, 27, 29, 35, 36, 25, 27, 20, 26, 28, 30, 35, 36, 26, 31, 33, 34, 26, 31, 32, 34, 32, 34, 39, 40, 33, 35, 41, 42, 31, 33, 26, 32, 34, 36, 41, 42, 32, 37, 39, 40, 32, 37, 38, 40, 38, 40, 45, 46, 39, 41, 47, 48, 37, 39, 32, 38, 40, 42, 47, 48, 38, 43, 45, 46, 38, 43, 44, 46, 44, 46, 51, 52, 45, 47, 53, 54, 43, 45, 38, 44, 46, 48, 53, 54, 44, 49, 51, 52, 44, 49, 50, 52, 50, 52, 57, 58, 51, 53, 59, 60, 49, 51, 44, 50, 52, 54, 59, 60, 50, 55, 57, 58, 50, 55, 56, 58, 56, 58, 57, 59, 55, 57, 50, 56, 58, 60]
nzval = [-0.64, 1.0, -1.0, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -0.03510101010101016, -0.975, -1.0806825309567203, 1.0, -0.95, -0.025, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0, -0.025, -0.95, -0.3564, -0.64, 1.0, -1.0, 13.792569659442691, 0.8606811145510832, -13.792569659442691, 1.0, 0.03475000000000006, 1.0, -1.0806825309567203, 1.0, 2.370597639417811, -2.3705976394178108, 10.698449178570607, -11.083604432603583, -0.2770901108150896, 1.0]
A = SparseMatrixCSC(n, n, colptr, rowval, nzval)
return(A)
end
test_pardiso.jl shows that using Pardiso directly works fine:
using Pkg
Pkg.activate(; temp=true)
Pkg.add("MKL")
Pkg.add("Pardiso")
using LinearAlgebra
using MKL
using Random
using SparseArrays
using Pardiso
include("A.jl")
A = makeA()
b = ones(size(A, 1))
ps = MKLPardisoSolver()
x1 = similar(b)
Pardiso.solve!(ps, x1, A, b)
x2 = A\b
@show norm(x1 - x2)
Environment :
- Output of
using Pkg; Pkg.status()
Status `/tmp/jl_MTVh52/Project.toml`
[7ed4a6bd] LinearSolve v2.29.1
[46dd5b70] Pardiso v0.5.7- Output of
using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
Status `/tmp/jl_MTVh52/Manifest.toml`
[47edcb42] ADTypes v1.0.0
[7d9f7c33] Accessors v0.1.36
[79e6a3ab] Adapt v4.0.4
[4fba245c] ArrayInterface v7.10.0
[4c555306] ArrayLayouts v1.9.2
[62783981] BitTwiddlingConvenienceFunctions v0.1.5
[2a0fbf3d] CPUSummary v0.2.4
[d360d2e6] ChainRulesCore v1.23.0
[fb6a15b2] CloseOpenIntervals v0.1.12
[38540f10] CommonSolve v0.2.4
[34da2185] Compat v4.14.0
[a33af91c] CompositionsBase v0.1.2
[2569d6c7] ConcreteStructs v0.2.3
[187b0558] ConstructionBase v1.5.5
[adafc99b] CpuId v0.3.1
[9a962f9c] DataAPI v1.16.0
[e2d170a0] DataValueInterfaces v1.0.0
[ffbed154] DocStringExtensions v0.9.3
[4e289a0a] EnumX v1.0.4
[e2ba6199] ExprTools v0.1.10
[29a986be] FastLapackInterface v2.0.2
[1a297f60] FillArrays v1.10.2
[069b7b12] FunctionWrappers v1.1.3
[77dc65aa] FunctionWrappersWrappers v0.1.3
[46192b85] GPUArraysCore v0.1.6
[3e5b6fbb] HostCPUFeatures v0.1.16
[615f187c] IfElse v0.1.1
[3587e190] InverseFunctions v0.1.14
[82899510] IteratorInterfaceExtensions v1.0.0
[692b3bcd] JLLWrappers v1.5.0
[ef3ab10e] KLU v0.6.0
[ba0b0d4f] Krylov v0.9.5
[10f19ff3] LayoutPointers v0.1.15
[5078a376] LazyArrays v1.10.0
[7ed4a6bd] LinearSolve v2.29.1
[bdcacae8] LoopVectorization v0.12.169
[1914dd2f] MacroTools v0.5.13
[d125e4d3] ManualMemory v0.1.8
[a3b82374] MatrixFactorizations v2.2.0
[6fe1bfb0] OffsetArrays v1.14.0
[bac558e1] OrderedCollections v1.6.3
[46dd5b70] Pardiso v0.5.7
[f517fe37] Polyester v0.7.13
[1d0040c9] PolyesterWeave v0.2.1
[aea7be01] PrecompileTools v1.2.1
[21216c6a] Preferences v1.4.3
[3cdcf5f2] RecipesBase v1.3.4
[731186ca] RecursiveArrayTools v3.13.0
[f2c3362d] RecursiveFactorization v0.2.23
[189a3867] Reexport v1.2.2
[ae029012] Requires v1.3.0
[7e49a35a] RuntimeGeneratedFunctions v0.5.13
[94e857df] SIMDTypes v0.1.0
[476501e8] SLEEFPirates v0.6.42
[0bca4576] SciMLBase v2.35.0
[c0aeaf25] SciMLOperators v0.3.8
[53ae85a6] SciMLStructures v1.1.0
[efcf1570] Setfield v1.1.1
[e56a9233] Sparspak v0.3.9
[aedffcd0] Static v0.8.10
[0d7ed370] StaticArrayInterface v1.5.0
[1e83bf80] StaticArraysCore v1.4.2
[7792a7ef] StrideArraysCore v0.5.6
[2efcf032] SymbolicIndexingInterface v0.3.16
[3783bdb8] TableTraits v1.0.1
[bd369af6] Tables v1.11.1
[8290d209] ThreadingUtilities v0.5.2
[d5829a12] TriangularSolve v0.2.0
[3a884ed6] UnPack v1.0.2
[3d5dd08c] VectorizationBase v0.21.66
[1d5cc7b8] IntelOpenMP_jll v2024.1.0+0
[856f044c] MKL_jll v2024.1.0+0
[1317d2d5] oneTBB_jll v2021.12.0+0
[0dad84c5] ArgTools v1.1.1
[56f22d72] Artifacts
[2a0f44e3] Base64
[ade2ca70] Dates
[8ba89e20] Distributed
[f43a241f] Downloads v1.6.0
[7b1f6079] FileWatching
[9fa8497b] Future
[b77e0a4c] InteractiveUtils
[4af54fe1] LazyArtifacts
[b27032c2] LibCURL v0.6.4
[76f85450] LibGit2
[8f399da3] Libdl
[37e2e46d] LinearAlgebra
[56ddb016] Logging
[d6f4376e] Markdown
[ca575930] NetworkOptions v1.2.0
[44cfe95a] Pkg v1.10.0
[de0858da] Printf
[3fa0cd96] REPL
[9a3f8284] Random
[ea8e919c] SHA v0.7.0
[9e88b42a] Serialization
[6462fe0b] Sockets
[2f01184e] SparseArrays v1.10.0
[10745b16] Statistics v1.10.0
[4607b0f0] SuiteSparse
[fa267f1f] TOML v1.0.3
[a4e569a6] Tar v1.10.0
[8dfed614] Test
[cf7118a7] UUIDs
[4ec0a83e] Unicode
[e66e0078] CompilerSupportLibraries_jll v1.1.0+0
[deac9b47] LibCURL_jll v8.4.0+0
[e37daf67] LibGit2_jll v1.6.4+0
[29816b5a] LibSSH2_jll v1.11.0+1
[c8ffd9c3] MbedTLS_jll v2.28.2+1
[14a3606d] MozillaCACerts_jll v2023.1.10
[4536629a] OpenBLAS_jll v0.3.23+4
[bea87d4a] SuiteSparse_jll v7.2.1+1
[83775a58] Zlib_jll v1.2.13+1
[8e850b90] libblastrampoline_jll v5.8.0+1
[8e850ede] nghttp2_jll v1.52.0+1
[3f19e933] p7zip_jll v17.4.0+2- Output of
versioninfo()
Julia Version 1.10.2
Commit bd47eca2c8a (2024-03-01 10:14 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 20 × 13th Gen Intel(R) Core(TM) i7-1370P
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, goldmont)
Threads: 1 default, 0 interactive, 1 GC (on 20 virtual cores)