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 <[email protected]> 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 <[email protected]>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
>> <[email protected]>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 [email protected].
>>> To unsubscribe from this group, send email to
>>> [email protected].
>>> 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 [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
