Use CUBLAS.geam for elementwise + and - on CuMatrix (JuliaGPU#1775)

lib/cublas/linalg.jl

@@ -591,3 +591,15 @@ end
         error("only supports BLAS type, got $T")
     end
 end
+
+op_wrappers = ((identity, T -> 'N', identity),
+               (T -> :(Transpose{T, <:$T}), T -> 'T', A -> :(parent($A))),
+               (T -> :(Adjoint{T, <:$T}), T -> T <: Real ? 'T' : 'C', A -> :(parent($A))))
+
+for op in (:(+), :(-))
+    for (wrapa, transa, unwrapa) in op_wrappers, (wrapb, transb, unwrapb) in op_wrappers
+        TypeA = wrapa(:(CuMatrix{T}))
+        TypeB = wrapb(:(CuMatrix{T}))
+        @eval Base.$op(A::$TypeA, B::$TypeB) where {T <: CublasFloat} = CUBLAS.geam($transa(T), $transb(T), one(T), $(unwrapa(:A)), $(op)(one(T)), $(unwrapb(:B)))
+    end
+end

test/cublas.jl

@@ -1198,6 +1198,27 @@ end
             h_C = Array(d_C)
             @test D ≈ h_C
         end
+        @testset "CuMatrix -- A ± B -- $elty" begin
+            for opa in (identity, transpose, adjoint)
+                for opb in (identity, transpose, adjoint)
+                    n = 10
+                    m = 20
+                    geam_A = opa == identity ? rand(elty, n, m) : rand(elty, m, n)
+                    geam_B = opb == identity ? rand(elty, n, m) : rand(elty, m, n)
+
+                    geam_dA = CuMatrix{elty}(geam_A)
+                    geam_dB = CuMatrix{elty}(geam_B)
+
+                    geam_C = opa(geam_A) + opb(geam_B)
+                    geam_dC = opa(geam_dA) + opb(geam_dB)
+                    @test geam_C ≈ collect(geam_dC)
+
+                    geam_C = opa(geam_A) - opb(geam_B)
+                    geam_dC = opa(geam_dA) - opb(geam_dB)
+                    @test geam_C ≈ collect(geam_dC)
+                end
+            end
+        end
         A = rand(elty,m,k)
         d_A = CuArray(A)
         @testset "syrkx!" begin