LinSolveAlgo.jl 4.81 KB
Newer Older
1
module LinSolveAlgo
Michel Juillard's avatar
Michel Juillard committed
2

3
4
5
6
import LinearAlgebra.BlasInt
import LinearAlgebra.BLAS.@blasfunc
import LinearAlgebra.BLAS.libblas
import LinearAlgebra.LAPACK: liblapack, chklapackerror
Michel Juillard's avatar
Michel Juillard committed
7

8
export LinSolveWS, linsolve_core!, linsolve_core_no_lu!, lu!
Michel Juillard's avatar
Michel Juillard committed
9
10

struct LinSolveWS
11
    lu::Matrix{Float64}
Michel Juillard's avatar
Michel Juillard committed
12
    ipiv::Vector{BlasInt}
13
14
15
    function LinSolveWS(n::Int64)
        lu = zeros(Float64,n,n)
        ipiv = zeros(BlasInt,n)
16
        new(lu,ipiv)
Michel Juillard's avatar
Michel Juillard committed
17
18
19
    end
end

20
21
22
23
24
25
26
27
28
function linsolve_core_no_lu!(ws::LinSolveWS,trans::Ref{UInt8},a::StridedMatrix{Float64},b::StridedVecOrMat{Float64})
    mm,nn = size(a)
    m = Ref{BlasInt}(mm)
    n = Ref{BlasInt}(nn)
    nhrs = Ref{BlasInt}(size(b,2))
    lda = Ref{BlasInt}(max(1,stride(a,2)))
    ldb = Ref{BlasInt}(max(1,stride(b,2)))
    info = Ref{BlasInt}(0)

29
    ccall((@blasfunc(dgetrs_), liblapack), Cvoid,
30
31
32
33
34
35
36
37
38
          (Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ptr{Float64},Ref{BlasInt},
           Ptr{BlasInt},Ptr{Float64},Ref{BlasInt},Ref{BlasInt}),
          trans,n,nhrs,ws.lu,lda,ws.ipiv,b,ldb,info)
    if info[] != 0
        println("dgetrs ",info[])
        chklapackerror(info[])
    end
end

Michel Juillard's avatar
Michel Juillard committed
39
40
41
42
43
44
45
46
47
function linsolve_core!(ws::LinSolveWS,trans::Ref{UInt8},a::StridedMatrix{Float64},b::StridedVecOrMat{Float64})
    mm,nn = size(a)
    m = Ref{BlasInt}(mm)
    n = Ref{BlasInt}(nn)
    nhrs = Ref{BlasInt}(size(b,2))
    lda = Ref{BlasInt}(max(1,stride(a,2)))
    ldb = Ref{BlasInt}(max(1,stride(b,2)))
    info = Ref{BlasInt}(0)

48
    lu!(ws.lu,a,ws.ipiv)
49
    ccall((@blasfunc(dgetrs_), liblapack), Cvoid,
Michel Juillard's avatar
Michel Juillard committed
50
51
          (Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ptr{Float64},Ref{BlasInt},
           Ptr{BlasInt},Ptr{Float64},Ref{BlasInt},Ref{BlasInt}),
52
          trans,n,nhrs,ws.lu,lda,ws.ipiv,b,ldb,info)
Michel Juillard's avatar
Michel Juillard committed
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
    if info[] != 0
        println("dgetrs ",info[])
        chklapackerror(info[])
    end
end

function linsolve_core!(ws::LinSolveWS,trans::Ref{UInt8},a::StridedMatrix{Float64},b::StridedVecOrMat{Float64},c::StridedVecOrMat{Float64})
    mm,nn = size(a)
    m = Ref{BlasInt}(mm)
    n = Ref{BlasInt}(nn)
    lda = Ref{BlasInt}(max(1,stride(a,2)))
    ldb = Ref{BlasInt}(max(1,stride(b,2)))
    ldc = Ref{BlasInt}(max(1,stride(c,2)))
    info = Ref{BlasInt}(0)

68
    lu!(ws.lu,a,ws.ipiv)
Michel Juillard's avatar
Michel Juillard committed
69
    nhrs = Ref{BlasInt}(size(b,2))
70
    ccall((@blasfunc(dgetrs_), liblapack), Cvoid,
Michel Juillard's avatar
Michel Juillard committed
71
72
          (Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ptr{Float64},Ref{BlasInt},
           Ptr{BlasInt},Ptr{Float64},Ref{BlasInt},Ref{BlasInt}),
73
          trans,n,nhrs,ws.lu,lda,ws.ipiv,b,ldb,info)
Michel Juillard's avatar
Michel Juillard committed
74
75
76
77
78
    if info[] != 0
        println("dgetrs ",info[])
        chklapackerror(info[])
    end
    nhrs = Ref{BlasInt}(size(c,2))
79
    ccall((@blasfunc(dgetrs_), liblapack), Cvoid,
Michel Juillard's avatar
Michel Juillard committed
80
81
          (Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ptr{Float64},Ref{BlasInt},
           Ptr{BlasInt},Ptr{Float64},Ref{BlasInt},Ref{BlasInt}),
82
          trans,n,nhrs,ws.lu,lda,ws.ipiv,c,ldc,info)
Michel Juillard's avatar
Michel Juillard committed
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
    if info[] != 0
        println("dgetrs ",info[])
        chklapackerror(info[])
    end
end

function linsolve_core!(ws::LinSolveWS,trans::Ref{UInt8},a::StridedMatrix{Float64},b::StridedVecOrMat{Float64},c::StridedVecOrMat{Float64},d::StridedVecOrMat{Float64})
    mm,nn = size(a)
    m = Ref{BlasInt}(mm)
    n = Ref{BlasInt}(nn)
    nhrs = Ref{BlasInt}(size(b,2))
    lda = Ref{BlasInt}(max(1,stride(a,2)))
    ldb = Ref{BlasInt}(max(1,stride(b,2)))
    ldc = Ref{BlasInt}(max(1,stride(c,2)))
    ldd = Ref{BlasInt}(max(1,stride(d,2)))
    info = Ref{BlasInt}(0)

100
    lu!(ws.lu,a,ws.ipiv)
101
    ccall((@blasfunc(dgetrs_), liblapack), Cvoid,
Michel Juillard's avatar
Michel Juillard committed
102
103
          (Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ptr{Float64},Ref{BlasInt},
           Ptr{BlasInt},Ptr{Float64},Ref{BlasInt},Ref{BlasInt}),
104
          trans,n,nhrs,ws.lu,lda,ws.ipiv,b,ldb,info)
Michel Juillard's avatar
Michel Juillard committed
105
106
107
108
    if info[] != 0
        println("dgetrs ",info[])
        chklapackerror(info[])
    end
109
    ccall((@blasfunc(dgetrs_), liblapack), Cvoid,
Michel Juillard's avatar
Michel Juillard committed
110
111
          (Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ptr{Float64},Ref{BlasInt},
           Ptr{BlasInt},Ptr{Float64},Ref{BlasInt},Ref{BlasInt}),
112
          trans,n,nhrs,ws.lu,lda,ws.ipiv,c,ldc,info)
Michel Juillard's avatar
Michel Juillard committed
113
114
115
116
    if info[] != 0
        println("dgetrs ",info[])
        chklapackerror(info[])
    end
117
    ccall((@blasfunc(dgetrs_), liblapack), Cvoid,
Michel Juillard's avatar
Michel Juillard committed
118
119
          (Ref{UInt8},Ref{BlasInt},Ref{BlasInt},Ptr{Float64},Ref{BlasInt},
           Ptr{BlasInt},Ptr{Float64},Ref{BlasInt},Ref{BlasInt}),
120
          trans,n,nhrs,ws.lu,lda,ws.ipiv,d,ldd,info)
Michel Juillard's avatar
Michel Juillard committed
121
122
123
124
125
126
    if info[] != 0
        println("dgetrs ",info[])
        chklapackerror(info[])
    end
end

127
128
function lu!(lu,a,ipiv)
    copy!(lu,a)
Michel Juillard's avatar
Michel Juillard committed
129
130
131
132
133
    mm,nn = size(a)
    m = Ref{BlasInt}(mm)
    n = Ref{BlasInt}(nn)
    lda = Ref{BlasInt}(max(1,stride(a,2)))
    info = Ref{BlasInt}(0)
134
    ccall((@blasfunc(dgetrf_), liblapack), Cvoid,
Michel Juillard's avatar
Michel Juillard committed
135
136
          (Ref{BlasInt},Ref{BlasInt},Ptr{Float64},Ref{BlasInt},
           Ptr{BlasInt},Ref{BlasInt}),
137
          m,n,lu,lda,ipiv,info)
Michel Juillard's avatar
Michel Juillard committed
138
139
140
141
142
143
144
    if info[] != 0
        println("dgetrf ",info[])
        chklapackerror(info[])
    end
end

end