Hi,
I am trying to write some sparse matrix code. I do the following:
costruct some (very sparse) matrices, pack them up into csc (that has
been solved in an earlier post) and call UMFpack.
For handling the matrices, I came up with the following code:
(defclass sparsematrix ()
((rows :initarg :rows)
(columns :initarg :columns)
(elements :initform (make-hash-table :test 'equal))))
(defun make-sm (rows columns)
"Create a sparse matrix of size rows x columns."
(assert (and (plusp rows) (plusp columns)))
(make-instance 'sparsematrix :rows rows :columns columns))
(defun sm-ref (sm row column)
"Retrieve element (row,column) from sm. An error is signaled
for indices out of range."
(with-slots (rows columns elements) sm
(assert (and (<= 0 row) (< row rows)))
(assert (and (<= 0 column) (< column columns)))
(print "sm-ref called")
(multiple-value-bind (value found) (gethash (cons row column) elements)
(if found
value
0))))
(defun (setf sm-ref) (value sm row column)
"Set element (row,column) in sm. An error is signaled
for indices out of range."
;; !!! should remove element for (zerop value)
(with-slots (rows columns elements) sm
(assert (and (<= 0 row) (< row rows)))
(assert (and (<= 0 column) (< column columns)))
(print "(setf sm-ref) called")
(setf (gethash (cons row column) elements) value)))
;; example
(defparameter *sm* (make-sm 5 5))
(sm-ref *sm* 2 1)
(setf (sm-ref *sm* 2 1) 3)
(incf (sm-ref *sm* 2 1))
Now when I run the example:
;;;; (defparameter *sm* (make-sm 5 5)) ...
;;;; (sm-ref *sm* 2 1) ...
"sm-ref called"
;;;; (setf (sm-ref *sm* 2 1) 3) ...
"(setf sm-ref) called"
;;;; (incf (sm-ref *sm* 2 1)) ...
"sm-ref called"
"(setf sm-ref) called"
Two questions:
1. I understand that incf calls sm-ref to retrieve the value, and
(setf sm-ref) to set it. Yet if the value is nonzero (ie it is in the
hash table), the "place" is retrieved by the first call already.
Would it be possible to have just one call to gethash in this case?
2. How can I define "multf", a function that multiplies the value in
place, ie (multf a b) would set a to (* a b)?
Thanks,
Tamas
On Jun 29, 10:04 am, Tamas Papp <······@gmail.com> wrote:
[...]
> 2. How can I define "multf", a function that multiplies the value in
> place, ie (multf a b) would set a to (* a b)?
I see Pascal already took a stab at answering this, but I'll expand on
it a little.
First, you can't write a function to do stuff like this as a general
rule. There's no way to pass things by reference in Common Lisp, which
is what allows you to use functions for this purpose in other
languages. Instead, you have to define a macro. SETF and friends
(INCF, PUSH, et c.) are all macros to begin with.
Unfortunately, the na�ve approach to this macro won't work. If you try
(defmacro multf (place factor)
`(setf ,place (* ,place ,factor)))
it will appear to work in some cases:
* (let ((x 4))
(multf x 2)
x)
=> 8
Don't be fooled. That's just the language being sneaky and trying to
lull you into a false sense of security. This can screw up in all
sorts of hopefully entertaining ways:
* (let ((a (make-array 3 :initial-contents '(1 2 3))
(b 1))
(multf (aref a (incf b)) 2)
=> #(1 6 3)
Macroexpanding the MULTF form, you get
(setf (aref a (incf b)) (* (aref a (incf b)) 2))
which probably isn't what you wanted to do at all. This is a common
bug in Lisp macros---something you only wanted to evaluate once ends
up being evaluated more than once. In most cases it's possible to code
around in a straightforward way[1], but not with SETF-based forms. The
usual solution is to establish temporary bindings (in a LET) for all
the things you don't want to evaluate more than once, but that doesn't
work at all for SETF.
These issues are so tricky to deal with that, in fact, Common Lisp
provides it's own special set of functions to help you out. Pascal
mentioned GET-SETF-EXPANSION, and that's the one that will let you
write MULTF. Unfortunately, it's still sort of tricky to use GET-SETF-
EXPANSION right, since it returns five values and you need to figure
out what to do with them all. For now, I'd probably recommend just
living without MULTF, but if you're feeling ambitious, there's a good
discussion of these issues in /On Lisp/[2].
[1] /Practical Common Lisp/ has a chapter on defining macros that
covers easier to deal with cases beautifully. See
http://www.gigamonkeys.com/book/macros-defining-your-own.html
[2] /On Lisp/ is a substantially more advanced book than /PCL/ is, but
if you want to learn lots and lots about macros, it's the thing to
read. The chapter in question is available online at
http://www.bookshelf.jp/texi/onlisp/onlisp_13.html
One thing to be careful of is that the author is the book is a little
old, and the author says GET-SETF-METHOD where we'd say GET-SETF-
EXPANSION. I assume the name was changed during the ANSI
standardization process.
Pillsy <·········@gmail.com> wrote:
> These issues are so tricky to deal with that, in fact, Common Lisp
> provides it's own special set of functions to help you out. Pascal
> mentioned GET-SETF-EXPANSION, and that's the one that will let you
> write MULTF. Unfortunately, it's still sort of tricky to use GET-SETF-
> EXPANSION right, since it returns five values and you need to figure
> out what to do with them all.
What's wrong with the non-tricky (define-modify-macro multf (value) *)?
--
Juho Snellman
On Jun 29, 12:38 pm, Juho Snellman <······@iki.fi> wrote:
> Pillsy <·········@gmail.com> wrote:
> > These issues are so tricky to deal with that, in fact, Common Lisp
> > provides it's own special set of functions to help you out. Pascal
> > mentioned GET-SETF-EXPANSION, and that's the one that will let you
> > write MULTF. Unfortunately, it's still sort of tricky to use GET-SETF-
> > EXPANSION right, since it returns five values and you need to figure
> > out what to do with them all.
> What's wrong with the non-tricky (define-modify-macro multf (value) *)?
Nothing would be wrong with it, if I hadn't been too dumb to remember
what it does. :-/
Cheers,
Pillsy
Tamas Papp wrote:
> I am trying to write some sparse matrix code. I do the following:
> costruct some (very sparse) matrices, pack them up into csc (that has
> been solved in an earlier post) and call UMFpack.
>
> For handling the matrices, I came up with the following code:
>
> (defclass sparsematrix ()
> ((rows :initarg :rows)
> (columns :initarg :columns)
> (elements :initform (make-hash-table :test 'equal))))
I'm curious why you would want to store the elements of the sparse
matrix in a hash table. The usual applications in which you need to add
elements dynamically can be handled by extensible arrays (the choice of
the Sparse BLAS reference implementation) or lists, and it's pretty rare
that you would want to _remove_ elements (usually direct factorizations
don't bother removing numerically canceled elements).
I can provide some references on typical sparse matrix data structures
if you're interested. There are some decent online references, e.g.,
Saad's "Iterative Methods for Sparse Linear Systems" has an overview:
http://www-users.cs.umn.edu/~saad/books.html
Best,
mfh
Mark Hoemmen <············@gmail.com> writes:
> Tamas Papp wrote:
>> I am trying to write some sparse matrix code. I do the following:
>> costruct some (very sparse) matrices, pack them up into csc (that has
>> been solved in an earlier post) and call UMFpack.
>>
>> For handling the matrices, I came up with the following code:
>>
>> (defclass sparsematrix () ((rows :initarg :rows)
>> (columns :initarg :columns)
>> (elements :initform (make-hash-table :test 'equal))))
Pascal, Juho and Pillsy: thanks for the followup.
> I'm curious why you would want to store the elements of the sparse
> matrix in a hash table. The usual applications in which you need to
> add elements dynamically can be handled by extensible arrays (the
> choice of the Sparse BLAS reference implementation) or lists, and it's
> pretty rare that you would want to _remove_ elements (usually direct
> factorizations don't bother removing numerically canceled elements).
Even if it is unusual, is there any problem with it? I am using hash
tables not because I want to remove elements but because finding and
adding elements is quick (I could use balanced trees and such things,
but probably I would end up reimplementing something like hash tables
;-)
Wouldn't plain lists, for example, impose a large overhead when I want
to find an existing element?
> I can provide some references on typical sparse matrix data structures
> if you're interested. There are some decent online references, e.g.,
> Saad's "Iterative Methods for Sparse Linear Systems" has an overview:
>
> http://www-users.cs.umn.edu/~saad/books.html
Thanks, I will look into it.
Tamas
Tamas Papp wrote:
> Mark Hoemmen <············@gmail.com> writes:
>> I'm curious why you would want to store the elements of the sparse
>> matrix in a hash table. The usual applications in which you need to
>> add elements dynamically can be handled by extensible arrays (the
>> choice of the Sparse BLAS reference implementation) or lists, and it's
>> pretty rare that you would want to _remove_ elements (usually direct
>> factorizations don't bother removing numerically canceled elements).
>
> Even if it is unusual, is there any problem with it? I am using hash
> tables not because I want to remove elements but because finding and
> adding elements is quick (I could use balanced trees and such things,
> but probably I would end up reimplementing something like hash tables
> ;-)
There's no problem with it, of course ;-) It will probably affect
performance for typical tasks involving iterating over all the elements
of the matrix. It also makes it harder to give your sparse matrix to
some C or Fortran program, which will likely use CSC, CSR, or coordinate
storage as an input format. But if you don't plan on calling some
external routine, then the format doesn't matter so much. Hash tables
would be especially useful if you do a lot of tests of "is there an
element at (i,j)?", and the (i,j) queries aren't correlated much at all.
Here's an example: say you're writing a finite element code. Here's
how typical FEM codes build the stiffness matrix:
* Start with an array of all the elements
* For each element:
Construct its element matrix
Add it to the array of element matrices
* Assemble the sparse matrix from the element matrices
Each element matrix is a small clique in the graph of the sparse matrix.
Some codes chose to leave the sparse matrix in unassembled form as an
array of cliques; certain direct factorization techniques ("frontal")
evolved around this data structure. However, it was (and is) generally
considered faster if one first assembles the sparse matrix into some
CSC-like or CSR-like format. This is what modern factorization codes,
as well as iterative solvers, tend to do.
It might be interesting to see whether the assembly process could be
speeded up by using some sorted data structure to store the cliques.
>> I can provide some references on typical sparse matrix data structures
>> if you're interested. There are some decent online references, e.g.,
>> Saad's "Iterative Methods for Sparse Linear Systems" has an overview:
>>
>> http://www-users.cs.umn.edu/~saad/books.html
>
> Thanks, I will look into it.
You may also want to read the following:
T. Davis, "Direct Methods for Sparse Linear Systems," SIAM, 2006.
Tim Davis demonstrates a concise yet decently fast factorization code here:
http://www.cise.ufl.edu/research/sparse/CSparse/
Duff, Ersiman, and Reid, "Direct Methods for Sparse Matrices," Oxford
University Press, 1986.
Let me know if you want some more resources -- our group might have some
tutorials, diagrams of typical data structures, and the like.
Best,
mfh
Mark Hoemmen <············@gmail.com> writes:
> You may also want to read the following:
>
> T. Davis, "Direct Methods for Sparse Linear Systems," SIAM, 2006.
>
> Tim Davis demonstrates a concise yet decently fast factorization code here:
> http://www.cise.ufl.edu/research/sparse/CSparse/
>
> Duff, Ersiman, and Reid, "Direct Methods for Sparse Matrices," Oxford
> University Press, 1986.
>
> Let me know if you want some more resources -- our group might have
> some tutorials, diagrams of typical data structures, and the like.
Hi Mark,
I decided to implement a rudimentary sparse matrix solver in Lisp,
partly for fun and partly to figure out how I could apply partial
evaluation to speed it up.
Can you please recommend some online tutorials or papers to get me
started? I don't need the fastest or most advanced method, just
something "good enough". The matrices are similar to the ones that
result from solving 2 or 3 dimensional PDEs by finite differences
(they are used in the Kushner-Dupuis method of solving variational
equations). I am asking for something online because currently I
don't have access to a good library (but I can get articles through my
university's library proxy).
Thanks,
Tamas
Tamas Papp wrote:
> I decided to implement a rudimentary sparse matrix solver in Lisp,
> partly for fun and partly to figure out how I could apply partial
> evaluation to speed it up.
Ooo, that would be very interesting ;-) I'm curious to see how it turns
out!
> Can you please recommend some online tutorials or papers to get me
> started? I don't need the fastest or most advanced method, just
> something "good enough". The matrices are similar to the ones that
> result from solving 2 or 3 dimensional PDEs by finite differences
> (they are used in the Kushner-Dupuis method of solving variational
> equations). I am asking for something online because currently I
> don't have access to a good library (but I can get articles through my
> university's library proxy).
CSparse is a good bet:
http://www.cise.ufl.edu/research/sparse/CSparse/
It's C, but very concise and written as a tutorial code (but still
efficient).
mfh
Mark Hoemmen wrote:
>> Can you please recommend some online tutorials or papers to get me
>> started? I don't need the fastest or most advanced method, just
>> something "good enough". The matrices are similar to the ones that
>> result from solving 2 or 3 dimensional PDEs by finite differences
>> (they are used in the Kushner-Dupuis method of solving variational
>> equations). I am asking for something online because currently I
>> don't have access to a good library (but I can get articles through my
>> university's library proxy).
Oh, oops, you wanted tutorials / papers and not libraries, sorry about
that ;-)
Let me see, there are some old papers that talk about elimination trees:
@article{ demmel99supernodal,
author = "James W. Demmel and Stanley C. Eisenstat and John R.
Gilbert and Xiaoye S. Li and Joseph W. H. Liu",
title = "A Supernodal Approach to Sparse Partial Pivoting",
journal = "SIAM Journal on Matrix Analysis and Applications",
volume = "20",
number = "3",
pages = "720--755",
year = "1999",
url = "citeseer.ist.psu.edu/demmel95supernodal.html"
}
and a lot of the bibliography in Citeseer for the above paper, including
e.g.,
@article{ gilbert94predicting,
author = "John R. Gilbert",
title = "Predicting Structure in Sparse Matrix Computations",
journal = "SIAM Journal on Matrix Analysis and Applications",
volume = "15",
number = "1",
pages = "62--79",
year = "1994",
url = "citeseer.ist.psu.edu/gilbert94predicting.html" }
J.W.H. LIU. The role of elimination trees in sparse factorization. SIAM
Journal of Matrix Analysis and Applications 11 (1990), 134--172.
and other papers that talk about "predicting structure" or "elimination
trees." I've read those three above.
A lot of these are SIAM papers and thus probably available via SIAM's
website through your library's proxy.
Best,
mfh
Tamas Papp <······@gmail.com> writes:
> Hi Mark,
>
> I decided to implement a rudimentary sparse matrix solver in Lisp,
> partly for fun and partly to figure out how I could apply partial
> evaluation to speed it up.
>
> Can you please recommend some online tutorials or papers to get me
> started? I don't need the fastest or most advanced method, just
> something "good enough". The matrices are similar to the ones that
> result from solving 2 or 3 dimensional PDEs by finite differences
> (they are used in the Kushner-Dupuis method of solving variational
> equations). I am asking for something online because currently I
> don't have access to a good library (but I can get articles through my
> university's library proxy).
>
> Thanks,
>
> Tamas
Maybe you could also look at Femlisp's (www.femlisp.org) sparse matrix
approach, which also uses hash-tables (but the entries are assumed to be
full matrix blocks). There is also code for standard sparse storage
schemes.
femlisp/src/matlisp/compressed.lisp
femlisp/src/algebra/sparse.lisp
femlisp/src/algebra/sparse-lu.lisp
I'm not completely happy with it. But I don't have anything better at the
moment.
Nicolas
Tamas Papp wrote:
> I am trying to write some sparse matrix code. I do the following:
> costruct some (very sparse) matrices, pack them up into csc (that has
> been solved in an earlier post) and call UMFpack.
>
> For handling the matrices, I came up with the following code:
>
> (defclass sparsematrix ()
> ((rows :initarg :rows)
> (columns :initarg :columns)
> (elements :initform (make-hash-table :test 'equal))))
I'm curious why you would want to store the elements of the sparse
matrix in a hash table. The usual applications in which you need to add
elements dynamically can be handled by extensible arrays (the choice of
the Sparse BLAS reference implementation) or lists, and it's pretty rare
that you would want to _remove_ elements (usually direct factorizations
don't bother removing numerically canceled elements).
I can provide some references on typical sparse matrix data structures
if you're interested. There are some decent online references, e.g.,
Saad's "Iterative Methods for Sparse Linear Systems" has an overview:
http://www-users.cs.umn.edu/~saad/books.html
Best,
mfh
Tamas Papp <······@gmail.com> writes:
> (multiple-value-bind (value found) (gethash (cons row column) elements)
> (if found
> value
> 0))))
(gethash (cons row column) elements 0)
> Now when I run the example:
>
> ;;;; (defparameter *sm* (make-sm 5 5)) ...
> ;;;; (sm-ref *sm* 2 1) ...
>
> "sm-ref called"
> ;;;; (setf (sm-ref *sm* 2 1) 3) ...
>
> "(setf sm-ref) called"
> ;;;; (incf (sm-ref *sm* 2 1)) ...
>
> "sm-ref called"
> "(setf sm-ref) called"
>
> Two questions:
>
> 1. I understand that incf calls sm-ref to retrieve the value, and
> (setf sm-ref) to set it. Yet if the value is nonzero (ie it is in the
> hash table), the "place" is retrieved by the first call already.
> Would it be possible to have just one call to gethash in this case?
It is not possible. Well, it's implementation dependant, and any gain
you may expect even more so, if the implementation doesn't already
optimize the double gethash call.
> 2. How can I define "multf", a function that multiplies the value in
> place, ie (multf a b) would set a to (* a b)?
Defining a multf macro that would use the setf expander for the place.
See GET-SETF-EXPANSION and DEFINE-SETF-EXPANDER
You could also define a specific setf expander for sm-ref so incf
doesn't call gethash twice, but this is actually out of your hands:
the setf expander is implementation specific, and while it's
theorically possible for an implementation to expose the internal
objects as setf expander values that would definitely avoid computing
the hash slot twice, it can more easily, safely, and generally provide
the same benefit by just caching one hash slot. That is, calling
twice gethash with the same arguments will compute the hash, and the
hash slot only once.
For example, here is clisp's expansion for gethash:
C/USER[234]> (get-setf-expansion '(gethash key table))
(#:G10908 #:G10909) ;
(KEY TABLE) ;
(#:G10910) ;
(SYSTEM::PUTHASH #:G10908 #:G10909 #:G10910) ; ; writer form
(GETHASH #:G10908 #:G10909) ; reader form
It computes gethash only once, and use system::puthash to write.
This is the function actually used when you write (setf (gethash key table) value).
An implementation could return instead:
C/USER[234]> (get-setf-expansion '(gethash key table))
(#:LOCATION-123) ;
((SYSTEM::HASH-SLOT-LOCATION KEY TABLE)) ;
(#:G10910) ;
(SYSTEM::LOCATION-STORE #:LOCATION-123 #:G10910) ; ; writer form
(SYSTEM::LOCATION-LOAD #:LOCATION-123) ; reader form
but it would involve more complexities for the implementation (eg to
keep the location consistent accross garbage collection, etc), than
simply caching the last gethash call.
--
__Pascal Bourguignon__ http://www.informatimago.com/
NOTE: The most fundamental particles in this product are held
together by a "gluing" force about which little is currently known
and whose adhesive power can therefore not be permanently
guaranteed.
Pascal Bourguignon <···@informatimago.com> writes:
> Tamas Papp <······@gmail.com> writes:
>
> (gethash (cons row column) elements 0)
good point, thanks
>
> It is not possible. Well, it's implementation dependant, and any gain
> you may expect even more so, if the implementation doesn't already
> optimize the double gethash call.
Thanks, if you say that the implementations already use caching, then
the probelm is solved.
Another question: I tried to write a print-object method like this (I
want dots for zero elements):
(defmethod print-object ((obj sparsematrix) stream)
(print-unreadable-object (obj stream :type t)
(with-slots (rows columns elements) obj
(dotimes (row rows)
(format stream "~%")
(dotimes (column columns)
(let ((elt (sm-ref obj row column)))
(format stream "~8,,,@A "
(if (zerop elt)
"."
elt))))))))
But I have a problem with numbers that are too long:
(defparameter *sm* (make-sm 4 4))
(setf (sm-ref *sm* 0 0) 4.1)
(setf (sm-ref *sm* 3 1) 1234567890)
*sm*
#<SPARSEMATRIX
4.1 . . .
. . . .
. . . .
. 1234567890 . . >
Any suggestions to arrange the values in the matrix properly (padded
on the left)?
Thanks,
Tamas
Tamas Papp <······@gmail.com> writes:
> Another question: I tried to write a print-object method like this (I
> want dots for zero elements):
>
> (defmethod print-object ((obj sparsematrix) stream)
> (print-unreadable-object (obj stream :type t)
> (with-slots (rows columns elements) obj
> (dotimes (row rows)
> (format stream "~%")
> (dotimes (column columns)
> (let ((elt (sm-ref obj row column)))
> (format stream "~8,,,@A "
> (if (zerop elt)
> "."
> elt))))))))
>
> But I have a problem with numbers that are too long:
>
> (defparameter *sm* (make-sm 4 4))
> (setf (sm-ref *sm* 0 0) 4.1)
> (setf (sm-ref *sm* 3 1) 1234567890)
> *sm*
>
> #<SPARSEMATRIX
> 4.1 . . .
> . . . .
> . . . .
> . 1234567890 . . >
>
> Any suggestions to arrange the values in the matrix properly (padded
> on the left)?
Yes. Do two passes, the first one to compute the width of each columns
(or all cells).
(defmethod print-object ((obj sparsematrix) stream)
(print-unreadable-object (obj stream :type t)
(with-slots (rows columns elements) obj
(let ((widths (make-array columns :element-type 'integer
:initial-element 0)))
(dotimes (column columns)
(dotimes (row rows)
(setf (aref widths column)
(max (length (format nil "~A"
(gethash (cons row column) elements 0)))
(aref widths column)))))
(dotimes (row rows)
(format stream "~%")
(dotimes (column columns)
(let ((elt (sm-ref obj row column)))
(format stream "~V,,,@A "
(aref widths column)
(if (nth-value 1 (gethash (cons row column) elements))
; to print dots only on empty cells.
elt
".")))))))))
But getting it to print nicely portably is more complex. For example
clisp will insert spaces before the first line:
C/USER[289]> *sm*
#<SPARSEMATRIX
4.1 . . .
. . . .
. . 0.0 .
. 1234567890 . . >
--
__Pascal Bourguignon__ http://www.informatimago.com/
NOTE: The most fundamental particles in this product are held
together by a "gluing" force about which little is currently known
and whose adhesive power can therefore not be permanently
guaranteed.