WIP: continuing higher-level wrappers

src/MAGMA.jl

@@ -2,17 +2,32 @@ module MAGMA
 using CUDAdrv, CUDAapi, CUDAnative, CuArrays
 using CEnum
 
+# some useful functions from LinearAlgebra for some intermediate processes
 using LinearAlgebra: triu, tril, dot, checksquare, chkstride1
 
+# some useful const definitions or utility functions
+export magmaTypeTuple, magmaTypeDict, magmaTypeList
+export chkfinite
+
 # export wrappers in svds
 export magma_gels!, magma_gesvd!, magma_gesdd!, magmaInit, magmaFinalize, magma_gebrd!
 
 # export wrappers in linearsystemsolver
-export magma_gesv!, magma_getri!, magma_getrs!, magma_getrf!
+export magma_gesv!, magma_getri!, magma_getrs!, magma_getrf!, magma_gerbt!, magma_gesv_rbt!, magma_geqrsv!
+export magma_posv!, magma_hesv!, magma_sysv!
+
+# export wrappers in factorization
+export magma_geqrf!, magma_geqlf!, magma_gelqf!
+
+# export wrappers in eigenvalues
+export magma_geev!, magma_gehrd!
 
 # export some wrappers in clang auto-generation
 export magma_init, magma_finalize
 
+#* some extra equipments
+require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
+
 include("clang/libmagma_common.jl")
 include("clang/libmagma_api.jl")
 
@@ -34,11 +49,26 @@ macro magmafunc_nb(function_name)
 	return Expr(:quote, Symbol("magma_get_", function_name, "_nb"))
 end
 
+macro magmafunc_generic(function_name, str)
+	return Expr(:quote, Symbol(function_name, str))
+end
+
 const magmaTypeTuple= (Float32, Float64, ComplexF32, ComplexF64)
 const magmaTypeList = ["Float32", "Float64", "ComplexF32", "ComplexF64"]
+const magmaTypeList_real = ["Float32", "Float64"]
+const magmaTypeList_complex = ["ComplexF32", "ComplexF64"]
 const magmaTypeDict = Dict(Float32=>"s", Float64=>"d", ComplexF32=>"c", ComplexF64=>"z",
 "Float32"=>"s", "Float64"=>"d", "ComplexF32"=>"c", "ComplexF64"=>"z",)
 
+function chkfinite(A::AbstractMatrix)
+    for a in A
+        if !isfinite(a)
+            throw(ArgumentError("matrix contains Infs or NaNs"))
+        end
+    end
+    return true
+end
+
 # include the files of subroutines
 include("dense/dense.jl")
 

src/dense/dense.jl

@@ -15,3 +15,5 @@ subsetrows(X::AbstractMatrix, Y::AbstractArray, k) = Y[1:k, :]
 
 include("svd.jl")
 include("linearsystemsolver.jl")
+include("factorization.jl")
+include("eigenvalues.jl")

src/dense/eigenvalues.jl

@@ -0,0 +1,141 @@
+const function_list_eigenvalues = ("geev", "gehrd", 
+                                   "heevx", "heevd", "heevdx", "heevr", "hetrd",
+                                   "hegvx", "hegvd", "hegvdx", "hegvr")
+macro char_to_vec(op)
+    return op == 'V' ? MagmaVec : MagmaNoVec
+end
+macro type_to_rtype(type)
+    return type == Float64 || type == ComplexF64 ? Float64 : Float32
+end
+
+for type in magmaTypeList_real
+    # create the symbols for element types
+    elty = Symbol(type)
+    relty= Symbol(@type_to_rtype(type))
+    # generate the symbol variables for our wrappers
+    for func_name in function_list_eigenvalues
+        @eval $(Symbol(func_name)) = (Symbol(magmaTypeDict[$type], $func_name))
+    end
+
+    @eval begin
+
+        function magma_geev!(jobvl::AbstractChar, jobvr::AbstractChar, A::AbstractMatrix{$elty})
+            magma_init()
+            chkstride1(A)
+            n = checksquare(A)
+            chkfinite(A) # balancing routines don't support NaNs and Infs
+            lvecs = jobvl == 'V'
+            rvecs = jobvr == 'V'
+            VL    = similar(A, $elty, (n, lvecs ? n : 0))
+            VR    = similar(A, $elty, (n, rvecs ? n : 0))
+            WR    = similar(A, $elty, n)
+            WI    = similar(A, $elty, n)
+            work  = Vector{$elty}(undef, 1)
+            lwork = Cint(-1)
+            info  = Ref{Cint}()
+            lda   = max(1, stride(A, 2))
+
+            jobvl = @char_to_vec(jobvl)
+            jobvr = @char_to_vec(jobvr)
+            func  = eval(@magmafunc($geev))
+            for i = 1:2  # first call returns lwork as work[1]
+                func(jobvl, jobvr, n, A, lda, WR, WI, VL, n,
+                        VR, n, work, lwork, info)
+                if i == 1
+                    lwork = ceil(Int, real(work[1]))
+                    resize!(work, lwork)
+                end
+            end
+            magma_finalize()
+            (WR, WI, VL, VR)
+        end
+
+    end
+
+end
+
+for type in magmaTypeList_complex
+    # create the symbols for element types
+    elty = Symbol(type)
+    relty= Symbol(@type_to_rtype(type))
+    # generate the symbol variables for our wrappers
+    for func_name in function_list_eigenvalues
+        @eval $(Symbol(func_name)) = (Symbol(magmaTypeDict[$type], $func_name))
+    end
+
+    @eval begin
+
+        function magma_geev!(jobvl::AbstractChar, jobvr::AbstractChar, A::AbstractMatrix{$elty})
+            magma_init()
+            chkstride1(A)
+            n = checksquare(A)
+            chkfinite(A) # balancing routines don't support NaNs and Infs
+            lvecs = jobvl == 'V'
+            rvecs = jobvr == 'V'
+            VL    = similar(A, $elty, (n, lvecs ? n : 0))
+            VR    = similar(A, $elty, (n, rvecs ? n : 0))
+            W     = similar(A, $elty, n)
+            rwork = similar(A, $relty, 2n)
+            work  = Vector{$elty}(undef, 1)
+            lwork = Cint(-1)
+            info  = Ref{Cint}()
+            lda   = max(1, stride(A, 2))
+
+            jobvl = @char_to_vec(jobvl)
+            jobvr = @char_to_vec(jobvr)
+            func  = eval(@magmafunc($geev))
+            for i = 1:2  # first call returns lwork as work[1]
+                func(jobvl, jobvr, n, A, lda, W, VL, n, VR, n,
+                    work, lwork, rwork, info)
+                if i == 1
+                    lwork = ceil(Int, real(work[1]))
+                    resize!(work, lwork)
+                end
+            end
+            magma_finalize()
+            (W, VL, VR)
+        end
+
+    end
+
+end
+
+for type in magmaTypeList
+    # create the symbols for element types
+    elty = Symbol(type)
+    relty= Symbol(@type_to_rtype(type))
+    # generate the symbol variables for our wrappers
+    for func_name in function_list_eigenvalues
+        @eval $(Symbol(func_name)) = (Symbol(magmaTypeDict[$type], $func_name))
+    end
+
+    @eval begin
+
+        function magma_gehrd!(ilo::Integer, ihi::Integer, A::AbstractMatrix{$elty})
+            magma_init()
+            chkstride1(A)
+            n = checksquare(A)
+            chkfinite(A) # balancing routines don't support NaNs and Infs
+            tau = similar(A, $elty, max(0,n - 1))
+            work = Vector{$elty}(undef, 1)
+            lwork = Cint(-1)
+            info = Ref{Cint}()
+            func  = eval(@magmafunc($gehrd))    
+            nb    = eval(@magmafunc_nb($gehrd))(n)
+            T     = CuArray{$elty}(undef, nb*n)
+            for i = 1:2  # first call returns lwork as work[1]
+                func(n, ilo, ihi, A,
+                    max(1, stride(A, 2)), tau, work, lwork, T, 
+                    info)
+                if i == 1
+                    lwork = ceil(Int, real(work[1]))
+                    resize!(work, lwork)
+                end
+            end
+            magma_finalize()
+            A, tau
+        end
+
+    end
+
+end

src/dense/factorization.jl

@@ -0,0 +1,157 @@
+# wrappers for MAGMA factorization functions
+
+const function_list_factorization = ("geqrf", "gerqf", "geqlf", "gelqf")
+for type in magmaTypeList
+    # create the symbols for element types
+    elty = Symbol(type)
+    # generate the symbol variables for our wrappers
+    for func_name in function_list_factorization
+        @eval $(Symbol(func_name)) = (Symbol(magmaTypeDict[$type], $func_name))
+    end
+
+    @eval begin
+
+        function magma_geqrf!(A::Matrix{$elty})
+            m, n  = size(A)
+            work  = Vector{$elty}(undef, 1)
+            lwork = Cint(-1)
+            info  = Ref{Cint}()
+            lda   = max(1, m)
+            tau   = similar(A, min(m, n))
+
+            func  = eval(@magmafunc($geqrf))
+            magma_init()
+            for i = 1:2                # 
+                func(m, n, A, lda, tau, work, lwork, info)
+                if i == 1
+                    lwork = ceil(Int, real(work[1]))
+                    resize!(work, lwork)
+                end
+            end
+            magma_finalize()
+            A, tau
+        end
+        function magma_geqrf!(A::CuMatrix{$elty})
+            magma_init()
+            m, n  = size(A)
+            nb    = (eval(@magmafunc_nb($geqrf)))( m, n )
+            info  = Ref{Cint}()
+            lda   = max(1, m)
+            tau   = similar(Matrix(A), min(m, n))
+            T     = similar(A, (2*min( m, n ) + ceil(Int, n/32)*32 )*nb)
+
+            func  = eval(@magmafunc_gpu($geqrf))
+            func(m, n, A, lda, tau, T, info)
+            magma_finalize()
+            A, tau
+        end
+
+        # ! Strangely, it seems that MAGMA does not provide any API for gerqf
+        # function magma_gerqf!(A::Matrix{$elty})
+        #     m, n  = size(A)
+        #     work  = Vector{$elty}(undef, 1)
+        #     lwork = Cint(-1)
+        #     info  = Ref{Cint}()
+        #     lda   = max(1, m)
+        #     tau   = similar(A, min(m, n))
+        #     func  = (eval(@magmafunc_generic($geqrf, "_error")))
+        #     magma_init()
+        #     for i = 1:2                # 
+        #         func(m, n, A, lda, tau, work, lwork, info)
+        #         if i == 1
+        #             lwork = ceil(Int, real(work[1]))
+        #             resize!(work, lwork)
+        #         end
+        #     end
+        #     magma_finalize()
+        #     A, tau
+        # end
+        # function magma_gerqf!(A::CuMatrix{$elty})
+        #     #! MAGMA does not provide direct GPU interface for gerqf
+        #     A = Matrix(A)
+        #     magma_gerqf!(A)
+        # end
+
+        function magma_geqlf!(A::Matrix{$elty})
+            m, n  = size(A)
+            work  = Vector{$elty}(undef, 1)
+            lwork = Cint(-1)
+            info  = Ref{Cint}()
+            lda   = max(1, m)
+            tau   = similar(A, min(m, n))
+
+            func  = eval(@magmafunc($geqlf))
+            magma_init()
+            for i = 1:2                # 
+                func(m, n, A, lda, tau, work, lwork, info)
+                if i == 1
+                    lwork = ceil(Int, real(work[1]))
+                    resize!(work, lwork)
+                end
+            end
+            magma_finalize()
+            A, tau
+        end
+        function magma_geqlf!(A::CuMatrix{$elty})
+            A = Matrix(A)
+            magma_geqlf!(A)
+        end
+
+        function magma_gelqf!(A::Matrix{$elty})
+            m, n  = size(A)
+            work  = Vector{$elty}(undef, 1)
+            lwork = Cint(-1)
+            info  = Ref{Cint}()
+            lda   = max(1, m)
+            tau   = similar(A, min(m, n))
+
+            func  = eval(@magmafunc($gelqf))
+            magma_init()
+            for i = 1:2                # 
+                func(m, n, A, lda, tau, work, lwork, info)
+                if i == 1
+                    lwork = ceil(Int, real(work[1]))
+                    resize!(work, lwork)
+                end
+            end
+            magma_finalize()
+            A, tau
+        end
+        function magma_gelqf!(A::CuMatrix{$elty})
+            m, n  = size(A)
+            work  = Vector{$elty}(undef, 1)
+            lwork = Cint(-1)
+            info  = Ref{Cint}()
+            lda   = max(1, m)
+            tau   = similar(Matrix(A), min(m, n))
+
+            func  = eval(@magmafunc_gpu($gelqf))
+            magma_init()
+            for i = 1:2                # 
+                func(m, n, A, lda, tau, work, lwork, info)
+                if i == 1
+                    lwork = ceil(Int, real(work[1]))
+                    resize!(work, lwork)
+                end
+            end
+            magma_finalize()
+            A, tau
+        end
+
+        function magma_geqrs!(A::CuArray{$elty}, B::CuArray{$elty})
+            magma_init()
+
+            m, n = size(A)
+            nrhs = size(B, 2)
+            lda  = max(1, m)
+            ldb  = max(1, m)
+            tau  = Array{$elty}(undef, n)
+
+
+            magma_finalize()
+            B, info[]
+        end
+
+    end
+
+end