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.

Reply via email to