Skip to content

Commit

Permalink
Added utilities for sparse and symmetric matrices
Browse files Browse the repository at this point in the history 
Added Blas functions for multiplication by symmetric matrices.  The
Blas.symv function is used in the arpack.jl code that will soon be
committed.

Combined mminfo and mmread functions in sparse.jl.  Also added the
fkeep!, dropzeros! and droptol! functions based on corresponding
functions in CSparse documented in section 2.7 of Tim Davis's book
"Direct Methods for Sparse Linear Systems".
  • Loading branch information
@dmbates
dmbates committed Nov 9, 2012
1 parent 752fc21 commit 3ab2610
Show file tree
Hide file tree
Showing 2 changed files with 116 additions and 64 deletions.
64 changes: 63 additions & 1 deletion base/blas.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,11 @@ export copy!,
sbmv!,
sbmv,
gemm!,
gemm
gemm,
symm!,
symm,
symv!,
symv

# SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
for (fname, elty) in ((:dcopy_,:Float64), (:scopy_,:Float32),
Expand Down Expand Up @@ -321,6 +325,64 @@ for (fname, elty) in ((:dgemv_,:Float64), (:sgemv_,:Float32),
end
end

# (SY) symmetric matrix-matrix and matrix-vector multiplication

# SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
# .. Scalar Arguments ..
# DOUBLE PRECISION ALPHA,BETA
# INTEGER LDA,LDB,LDC,M,N
# CHARACTER SIDE,UPLO
# .. Array Arguments ..
# DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)

# SUBROUTINE DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
# .. Scalar Arguments ..
# DOUBLE PRECISION ALPHA,BETA
# INTEGER INCX,INCY,LDA,N
# CHARACTER UPLO
# .. Array Arguments ..
# DOUBLE PRECISION A(LDA,*),X(*),Y(*)

for (vfname, mfname, elty) in
((:dsymv_,:dsymm_,:Float64),
(:ssymv_,:ssymm_,:Float32),
(:zsymv_,:zsymm_,:Complex128),
(:csymv_,:csymm_,:Complex64))
@eval begin
function symv!(uplo, alpha::($elty), A::StridedMatrix{$elty}, X::StridedVector{$elty},
beta::($elty), Y::StridedVector{$elty})
m, n = size(A)
if m != n error("symm!: matrix A is $m by $n but must be square") end
if m != length(X) || m != length(Y) error("symm!: dimension mismatch") end
ccall(dlsym(Base.libblas, $(string(vfname))), Void,
(Ptr{Uint8}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{Int32},
Ptr{$elty}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{Int32}),
&uplo, &n, &alpha, A, &stride(A,2), X, &stride(X,1), &beta, Y, &stride(Y,1))
Y
end
function symv(uplo, alpha::($elty), A::StridedMatrix{$elty}, X::StridedVector{$elty})
symv!(uplo, alpha, A, X, zero($elty), similar(X))
end
function symm!(side, uplo, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty},
beta::($elty), C::StridedMatrix{$elty})
side = uppercase(convert(Char, side))
m, n = size(C)
k, j = size(A)
if k != j error("symm!: matrix A is $k by $j but must be square") end
if j != (side == 'L' ? m : n) || size(B,2) != n error("symm!: Dimension mismatch") end
ccall(dlsym(Base.libblas, $(string(mfname))), Void,
(Ptr{Uint8}, Ptr{Uint8}, Ptr{Int32}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{Int32},
Ptr{$elty}, Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{Int32}),
&side, &uplo, &m, &n, &alpha, A, &stride(A,2), B, &stride(B,2),
&beta, C, &stride(C,2))
C
end
function symm(side, uplo, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty})
symm!(side, uplo, alpha, A, B, zero($elty), similar(B))
end
end
end

end # module

# Use BLAS copy for small arrays where it is faster than memcpy, and for strided copying
Expand Down
116 changes: 53 additions & 63 deletions extras/sparse.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1067,60 +1067,17 @@ function assign(S::SparseAccumulator, v, i::Integer)
return S
end

function mminfo(filename::ASCIIString)
# function [rows, cols, entries, rep, field, symmetry] = mminfo(filename)
#
function mmread(filename::ASCIIString, infoonly::Bool)
# Reads the contents of the Matrix Market file 'filename'
# and extracts size and storage information.
#
# In the case of coordinate matrices, entries refers to the
# number of coordinate entries stored in the file. The number
# of non-zero entries in the final matrix cannot be determined
# until the data is read (and symmetrized, if necessary).
#
# In the case of array matrices, entries is the product
# rows*cols, regardless of whether symmetry was used to
# store the matrix efficiently.

mmfile = open(filename,"r")
tokens = split(chomp(readline(mmfile)), ' ')
if length(tokens) != 5 error("Not enough words on header line") end
if tokens[1] != "%%MatrixMarket" error("Not a valid MatrixMarket header.") end
(head1, rep, field, symm) = map(lowercase, tokens[2:5])
if head1 != "matrix"
error("This seems to be a MatrixMarket $head1 file, not a MatrixMarket matrix file")
end
# Read through comments, ignoring them
ll = readline(mmfile)
while length(ll) > 0 && ll[1] == '%'
ll = readline(mmfile)
end
dd = int(split(ll, ' ')) # Read dimensions
rows = dd[1]
cols = dd[2]
entries = (rep == "coordinate" ? dd[3] : rows * cols)
return rows, cols, entries, rep, field, symm
end

function mmread(filename::ASCIIString)
# function [A] = mmread(filename)
#
# function [A,rows,cols,entries,rep,field,symm] = mmread(filename)
#
# Reads the contents of the Matrix Market file 'filename'
# into the matrix 'A'. 'A' will be either sparse or full,
# into a matrix, which will be either sparse or full,
# depending on the Matrix Market format indicated by
# 'coordinate' (coordinate sparse storage), or
# 'array' (dense array storage). The data will be duplicated
# as appropriate if symmetry is indicated in the header.
#
# Optionally, size information about the matrix can be
# obtained by using the return values rows, cols, and
# entries, where entries is the number of nonzero entries
# in the final matrix. Type information can also be retrieved
# using the optional return values rep (representation), field,
# and symm (symmetry).
# as appropriate if symmetry is indicated in the header. (Not yet
# implemented).
#
# If infoonly is true information on the size and structure is
# returned.
mmfile = open(filename,"r")
tokens = split(chomp(readline(mmfile)))
if length(tokens) != 5 error("Not enough words on header line") end
Expand All @@ -1131,29 +1088,62 @@ function mmread(filename::ASCIIString)
end
if field != "real" error("non-float fields not yet allowed") end

ll = readline(mmfile) # Read through comments, ignoring them
while length(ll) > 0 && ll[1] == '%'
ll = readline(mmfile)
end
# Read size information
dd = int(split(ll, ' '))
rows = dd[1]
cols = dd[2]
ll = readline(mmfile) # Read through comments, ignoring them
while length(ll) > 0 && ll[1] == '%' ll = readline(mmfile) end
dd = int(split(ll)) # Read dimensions
rows = dd[1]
cols = dd[2]
entries = rep == "coordinate" ? dd[3] : rows * cols
if infoonly return rows, cols, entries, rep, field, symm end
if rep == "coordinate"
entries = dd[3]
rr = Array(Int32, entries)
cc = Array(Int32, entries)
xx = Array(Float64, entries)
for i in 1:entries
flds = split(chomp(readline(mmfile)))
flds = split(readline(mmfile))
rr[i] = int32(flds[1])
cc[i] = int32(flds[2])
xx[i] = float64(flds[3])
end
return sparse(rr, cc, xx, rows, cols)
elseif rep == "array"
aa = Array(Float64, 0)
for ll in EachLine(mmfile) push(aa, float64(ll)) end
return reshape(aa, (rows, cols))
end
reshape([float64(readline(mmfile)) for i in 1:entries], (rows,cols))
end

mmread(filename::ASCIIString) = mmread(filename, false)

## expand a colptr or rowptr into a full index vector
function expandptr{T<:Integer}(V::Vector{T})

This comment has been minimized.

Sign in to view

Copy link
@tkelman

tkelman Feb 7, 2016

Contributor

@dmbates did you ever use this anywhere? I can think of a few places that might be doing this operation, but they've never called this function in base apparently
if V[1] != 1 error("expandptr: first index must be one") end
res = similar(V, (int64(V[end]-1),))
for i in 1:(length(V)-1), j in V[i]:(V[i+1] - 1) res[j] = i end
res
end

# Based on the function cs_fkeep from the CSparse library
function fkeep!{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, f, other)
nzorig = nnz(A)
nz = 1
for j = 1:A.n
p = A.colptr[j] # record current position
A.colptr[j] = nz # set new position
while p < A.colptr[j+1]
if f(A.rowval[p], j, A.nzval[p], other)
A.nzval[nz] = A.nzval[p]
A.rowval[nz] = A.rowval[p]
nz += 1
end
p += 1
end
end
nz -= 1
A.colptr[A.n + 1] = nz
if nz < nzorig
grow(A.nzval, nz - nzorig)
grow(A.rowval, nz - nzorig)
end
A
end

dropzeros!{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}) = fkeep!(A, (i,j,x,other)->x!=zero(Tv), 0)
droptol!{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, tol::Tv) = fkeep!(A, (i,j,x,other)->abs(x)>other, tol)

0 comments on commit 3ab2610

Please sign in to comment.