[CUSOLVER] Avoid the conversion to CSR format for reordering routines (#1786)

- index is now an optional argument for reordering routines;
- we don't need to convert CSC matrices into CSR matrices anymore because reordering routines are based on the pattern of A + A';
- `zfd` now returns a column permutation for `SparseMatrixCSC`. This change makes sense because the CUSOLVER routine `csrzfd` computes a row permutation for sparse matrices in CSR format.

Author: amontoison

lib/cusolver/sparse.jl

@@ -231,60 +231,52 @@ for (fname, elty, relty) in ((:cusolverSpScsreigsHost, :ComplexF32, :Float32),
 end
 
 #csrsymrcm
-function symrcm(A::SparseMatrixCSC, index::Char)
+function symrcm(A::SparseMatrixCSC, index::Char='O')
     n, m = size(A)
     (m ≠ m) && throw(DimensionMismatch("A must be square, but has dimensions ($n,$m)!"))
     descA = CuMatrixDescriptor('G', 'L', 'N', index)
     nnzA = nnz(A)
-    Mat = similar(A)
-    transpose!(Mat, A)
-    colsA = convert(Vector{Cint}, Mat.rowval)
-    rowsA = convert(Vector{Cint}, Mat.colptr)
+    colsA = convert(Vector{Cint}, A.rowval)
+    rowsA = convert(Vector{Cint}, A.colptr)
     p = zeros(Cint, n)
     cusolverSpXcsrsymrcmHost(sparse_handle(), n, nnzA, descA, rowsA, colsA, p)
     return p
 end
 
 #csrsymmdq
-function symmdq(A::SparseMatrixCSC, index::Char)
+function symmdq(A::SparseMatrixCSC, index::Char='O')
     n, m = size(A)
     (m ≠ m) && throw(DimensionMismatch("A must be square, but has dimensions ($n,$m)!"))
     descA = CuMatrixDescriptor('G', 'L', 'N', index)
     nnzA = nnz(A)
-    Mat = similar(A)
-    transpose!(Mat, A)
-    colsA = convert(Vector{Cint}, Mat.rowval)
-    rowsA = convert(Vector{Cint}, Mat.colptr)
+    colsA = convert(Vector{Cint}, A.rowval)
+    rowsA = convert(Vector{Cint}, A.colptr)
     p = zeros(Cint, n)
     cusolverSpXcsrsymmdqHost(sparse_handle(), n, nnzA, descA, rowsA, colsA, p)
     return p
 end
 
 #csrsymamd
-function symamd(A::SparseMatrixCSC, index::Char)
+function symamd(A::SparseMatrixCSC, index::Char='O')
     n, m = size(A)
     (m ≠ m) && throw(DimensionMismatch("A must be square, but has dimensions ($n,$m)!"))
     descA = CuMatrixDescriptor('G', 'L', 'N', index)
     nnzA = nnz(A)
-    Mat = similar(A)
-    transpose!(Mat, A)
-    colsA = convert(Vector{Cint}, Mat.rowval)
-    rowsA = convert(Vector{Cint}, Mat.colptr)
+    colsA = convert(Vector{Cint}, A.rowval)
+    rowsA = convert(Vector{Cint}, A.colptr)
     p = zeros(Cint, n)
     cusolverSpXcsrsymamdHost(sparse_handle(), n, nnzA, descA, rowsA, colsA, p)
     return p
 end
 
 #csrmetisnd
-function metisnd(A::SparseMatrixCSC, index::Char)
+function metisnd(A::SparseMatrixCSC, index::Char='O')
     n, m = size(A)
     (m ≠ m) && throw(DimensionMismatch("A must be square, but has dimensions ($n,$m)!"))
     descA = CuMatrixDescriptor('G', 'L', 'N', index)
     nnzA = nnz(A)
-    Mat = similar(A)
-    transpose!(Mat, A)
-    colsA = convert(Vector{Cint}, Mat.rowval)
-    rowsA = convert(Vector{Cint}, Mat.colptr)
+    colsA = convert(Vector{Cint}, A.rowval)
+    rowsA = convert(Vector{Cint}, A.colptr)
     p = zeros(Cint, n)
     cusolverSpXcsrmetisndHost(sparse_handle(), n, nnzA, descA, rowsA, colsA, C_NULL, p)
     return p
@@ -296,16 +288,14 @@ for (fname, elty) in ((:cusolverSpScsrzfdHost, :Float32),
                       (:cusolverSpCcsrzfdHost, :ComplexF32),
                       (:cusolverSpZcsrzfdHost, :ComplexF64))
     @eval begin
-        function zfd(A::SparseMatrixCSC{$elty}, index::Char)
+        function zfd(A::SparseMatrixCSC{$elty}, index::Char='O')
             n, m = size(A)
             (m ≠ m) && throw(DimensionMismatch("A must be square, but has dimensions ($n,$m)!"))
             descA = CuMatrixDescriptor('G', 'L', 'N', index)
             nnzA = nnz(A)
-            Mat = similar(A)
-            transpose!(Mat, A)
-            colsA = convert(Vector{Cint}, Mat.rowval)
-            rowsA = convert(Vector{Cint}, Mat.colptr)
-            valsA = Mat.nzval
+            colsA = convert(Vector{Cint}, A.rowval)
+            rowsA = convert(Vector{Cint}, A.colptr)
+            valsA = A.nzval
             p = zeros(Cint, n)
             numnz = Ref{Cint}(0)
             $fname(sparse_handle(), n, nnzA, descA, valsA, rowsA, colsA, p, numnz)

test/cusolver/sparse.jl

@@ -13,31 +13,31 @@ A = rand(n, n)
 A = SparseMatrixCSC(A)
 
 @testset "symrcm" begin
-    p = symrcm(A, 'O')
+    p = symrcm(A)
     p .+= 1
     @test minimum(p) == 1
     @test maximum(p) == n
     @test isperm(p)
 end
 
 @testset "symmdq" begin
-    p = symmdq(A, 'O')
+    p = symmdq(A)
     p .+= 1
     @test minimum(p) == 1
     @test maximum(p) == n
     @test isperm(p)
 end
 
 @testset "symamd" begin
-    p = symamd(A, 'O')
+    p = symamd(A)
     p .+= 1
     @test minimum(p) == 1
     @test maximum(p) == n
     @test isperm(p)
 end
 
 @testset "metisnd" begin
-    p = metisnd(A, 'O')
+    p = metisnd(A)
     p .+= 1
     @test minimum(p) == 1
     @test maximum(p) == n
@@ -51,12 +51,12 @@ end
             A[i,i] = 0
         end
         A = SparseMatrixCSC(A)
-        p = zfd(A, 'O')
+        p = zfd(A)
         p .+= 1
         @test minimum(p) == 1
         @test maximum(p) == n
         @test isperm(p)
-        @test 0 ∉ diag(A[p,:])
+        @test 0 ∉ diag(A[:,p])
     end
 
     @testset "csrlsvlu!" begin