I wrote this up more cleanly here http://matthewrocklin.com/blog/work/2012/10/29/Matrix-Computations/
In it I describe a bit more background as well such as why BLAS is important and how NumPy can be improved. On Sun, Oct 28, 2012 at 10:25 AM, Matthew Rocklin <mrock...@gmail.com>wrote: > I pressed 'send' too early. > > This is the simplest non-trivial problem with which to use blas > computations. The system is able to deal with more complex computations and > predicates. Ideally we will be able to use rules/pattern matching to select > specialized routines that match the predicates. I apologize about the short > function names 'gemm', 'trsm', these are from BLAS. > > The goal of the project is to allow users to skip the middle step of > building the computation. The first and last (math and use) are > comprehensible to mathematical users, the second (gemm, trsm, etc...) is > not. I think that we should be able to automate the second step. > > The generated Fortran code is below. It mostly just calls BLAS. Note that > this code is about as optimized as you will find. BLAS should easily beat > numpy and most hand-written C (unless the hand-written C is very very > optimized). I'll post timings in a bit. > > subroutine f(alpha, A, B, beta, C, x, n) > > real*8, intent(in) :: A(n, n) > real*8, intent(in) :: B(n, n) > real*8, intent(inout) :: C(n, n) > real*8, intent(in) :: alpha > real*8, intent(in) :: beta > integer, intent(in) :: n > real*8, intent(inout) :: x(n, 1) > > call dgemm('N', 'N', n, n, n, alpha, A, n, B, n, beta, C, n) > call dtrsv('L', 'N', 'N', n, C, n, x, 1) > > RETURN > END > > We now have three intermediate representations > > 1. MatrixExprs (MatrixSymbol, MatMul, Inverse, ...) > 2. MatrixComputations (gemm, trsv, symm, cholesky) > 3. Fortran > > We lack a compiler from (1) to (2). I hope that rules, pattern matching, > and assumptions can perform this transformation. Once we have this (and > once bugs in MatrixComputations are cleaned up) I think that we'll have a > pretty awesome dense linear algebra library in SymPy. > > > On Sun, Oct 28, 2012 at 8:14 AM, Matthew Rocklin <mrock...@gmail.com>wrote: > >> Update: >> >> I have basic code generation completed. I found that the Routine class >> didn't meet my needs to I built something new. There is some duplication. >> An example: >> >> # Mathematical Problem >> A = MatrixSymbol('A', n, n) >> B = MatrixSymbol('B', n, n) >> C = MatrixSymbol('C', n, n) >> x = MatrixSymbol('x', n, 1) >> context = (Q.lower_triangular(A) & >> Q.lower_triangular(B) & >> Q.lower_triangular(C)) >> expr = (alpha*A*B + beta*C).I*x >> >> # Build a computation object from this expression >> # This should be replaced by a clever pattern matching/rules system >> gemm = GEMM(alpha, A, B, beta, C) # produces alpha*A*B + beta*C and >> stores result in C >> trsv = TRSV(alpha*A*B + beta*C, x) # computes (alpha*A*B+beta*C).I*x >> and stores result in x >> comp = MatrixRoutine((gemm, trsv), (alpha, A, B, beta, C, x), >> (expr,)) # a computation to do both gemm and trsv >> >> # build callable object >> f = comp.build(str, context) # Write fortran code (below) and compile >> with f2py >> >> # Use f >> import numpy as np >> nalpha, nbeta = 2.0, 3.0 >> nA,nB,nC = [np.asarray([[1,0],[3,4]], order='F', dtype='float64') >> for i in range(3)] >> nx = np.asarray((2.3, 4.5), order='F', dtype='float64') >> >> f(nalpha, nA, nB, nbeta, nC, nx) >> >> >> >> >> >> >> On Thu, Oct 25, 2012 at 1:09 PM, Matthew Rocklin <mrock...@gmail.com>wrote: >> >>> For those interested in following along I'm working off of my blas >>> branch. Here is a link to the comparison master->blas >>> >>> https://github.com/mrocklin/sympy/compare/blas >>> >>> In particular I recommend looking at matrices/expressions/blas.py. This >>> branch also includes my work on matrix assumptions (separate PR). >>> >>> The goal is to create Fortran code from MatrixExprs. The work left to do >>> can be decomposed into the following tasks >>> >>> 1. Write a good Computation class (the current structure isn't yet ideal) >>> 2. Write down enough of BLAS to get going (close!) >>> 3. Write down rules that translate MatrixExprs into BLAS computations >>> (this is easier if we have good pattern matching) >>> 4. Write down strategies to handle non-confluent rules (there are lots >>> of possible solutions) >>> 5. Use f2py to reconnect the Fortran code to numpy interface >>> 6. Compare SymPy MatExpr times to NumPy times. I expect we'll be able to >>> generate code that is several times faster than NumPy while still >>> maintaining a nice interface. >>> >>> >>> On Thu, Oct 25, 2012 at 12:58 PM, Øyvind Jensen <jensen.oyv...@gmail.com >>> > wrote: >>> >>>> >>>> >> If you don't mind, maybe we should CC the mailing list? >>>> > >>>> > Go for it. I've been flooding the list recently which is why I didn't >>>> do >>>> > this earlier. >>>> > >>>> > I've recently been working on a Theano-like graph for SymPy (a >>>> Computation >>>> > object). I'm may be reinventing something. I'm working off of this >>>> branch >>>> > https://github.com/mrocklin/sympy/tree/blas . What is there now will >>>> soon >>>> > change however. I'm still trying to figure out the best way to do >>>> this. >>>> > >>>> > Question about codegen. I have to compose routines that have multiple >>>> > out-params and inplace operations. Is SymPy codegen equipped to >>>> handle this? >>>> > Does it have tools to reason about these sorts of problems? Maybe >>>> this will >>>> > become clearer to me when I do the mat-vec multiply example. >>>> >>>> Multiple output has at least some support. In test_codegen.py you 'll >>>> find test_multiple_results_f(). What kind of inplace operations are you >>>> referring to? >>>> >>>> I took a quick glance at your blas branch. I think that your >>>> Computation class is somewhat comparable with codegen.Routine, except that >>>> you are doing it properly. The Routine class aims to describe how a given >>>> mathematical expression can be programmed as a function or subroutine. But >>>> with the Computation class, you aim to describe also the relationship >>>> between several computations, isn't that right? It might be a good idea to >>>> remove or at least shrink the Routine class in favor of your stuff. >>>> >>>> If you want codegen to support MatrixExpr objects, you would need to >>>> implement them in the code printers, and I guess you will also need to hack >>>> the constructor of Routine. >>>> >>>> >>>> Cheers, >>>> Øyvind >>>> >>>> -- >>>> You received this message because you are subscribed to the Google >>>> Groups "sympy" group. >>>> To view this discussion on the web visit >>>> https://groups.google.com/d/msg/sympy/-/9PvoswQHp6kJ. >>>> >>>> To post to this group, send email to sympy@googlegroups.com. >>>> To unsubscribe from this group, send email to >>>> sympy+unsubscr...@googlegroups.com. >>>> For more options, visit this group at >>>> http://groups.google.com/group/sympy?hl=en. >>>> >>> >>> >> > -- You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to sympy@googlegroups.com. To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.