I don't know if it was clear from my previous post, but Alan's two points are already implemented in the build in A\b for square A. First, A is promoted to a type which is stable under division and if that type is one that LAPACK supports then LAPACK is called for the remaining work and if not, the generic build in version does the job with pivoting but without blocking.
2014-02-21 19:34 GMT+01:00 andrew cooke <[email protected]>: > i wish i could *any* write code with no effort. > > > On Friday, 21 February 2014 10:07:35 UTC-3, Alan Edelman wrote: >> >> It would take no effort to write a GE for GF2 and even a GE that >> is mathematically correct for all fields, stability aside. >> >> A few thoughts to mull over: >> >> 1) for GF2 one would pivot on non-zero >> as there is no stability issue. For reals, we pivot to largest to avoid >> round-off. >> >> 2) For GF2, it's not clear to me that blocking would gain speed, >> while for reals, one key value of lapack is getting BLAS3 speeds >> for gaussian elimination through blocking >> >> >> >> >> >> >> >> >> >> >> On Friday, February 21, 2014 6:56:45 AM UTC-5, andrew cooke wrote: >>> >>> that's it exactly. thanks! >>> >>> On Friday, 21 February 2014 08:09:27 UTC-3, Andreas Noack Jensen wrote: >>>> >>>> Bitstypes and algebra are not my home ground, but the lufact that Tim >>>> references is supposed to work generically and I think it does. Have a look >>>> at >>>> >>>> https://gist.github.com/andreasnoackjensen/9132511 >>>> >>>> I am not sure the algebra is right. The important thing for lufact is >>>> that that the type is closed under /. Therefore I use rationals over GF(2). >>>> Maybe it is better just to define / for GF2. >>>> >>>> >>>> 2014-02-21 4:11 GMT+01:00 Stefan Karpinski <[email protected]>: >>>> >>>>> Yes, it's entirely possible that his might work for matrices over >>>>> finite fields. If not, we should see if we can modify it so that it does >>>>> work. >>>>> >>>>> >>>>> On Thu, Feb 20, 2014 at 10:02 PM, Tim Holy <[email protected]> wrote: >>>>> >>>>>> It's not Gaussian elimination (it's better!), but aren't these >>>>>> relevant? >>>>>> https://github.com/JuliaLang/julia/pull/5381 >>>>>> https://github.com/JuliaLang/julia/pull/5430 >>>>>> >>>>>> --Tim >>>>>> >>>>>> On Thursday, February 20, 2014 09:12:51 PM Stefan Karpinski wrote: >>>>>> > We ought to have a generic Gaussian elimination algorithm and >>>>>> making sure >>>>>> > that it can be used for finite fields is a good test of its >>>>>> generality. >>>>>> > This is definitely the kind of thing that can be more reasonably in >>>>>> Julia >>>>>> > than in many other systems. The big difference is that while in this >>>>>> > particular case you may have to write Gaussian elimination, at >>>>>> least the >>>>>> > next person with their own integer matrix type won't. >>>>>> > >>>>>> > On Thu, Feb 20, 2014 at 7:29 PM, andrew cooke <[email protected]> >>>>>> wrote: >>>>>> > > Indeed, but then I need to provide my own implementation. I was >>>>>> hoping >>>>>> > > (navely in retrospect, since everyone else wants fast rather than >>>>>> generic >>>>>> > > code) that I could call generic code for "matrix division" with a >>>>>> new >>>>>> > > numerical type and it would "just work". >>>>>> > > >>>>>> > > I think this is still the best approach - but I need to write >>>>>> that code >>>>>> > > (Gaussian elimination) myself, in a generic way, then define my >>>>>> new >>>>>> > > numeric >>>>>> > > type, then call it. >>>>>> > > >>>>>> > > Which is much more work than I had hoped. >>>>>> > > >>>>>> > > Cheers, >>>>>> > > Andrew >>>>>> > > >>>>>> > > On Thursday, 20 February 2014 21:14:38 UTC-3, Patrick O'Leary >>>>>> wrote: >>>>>> > >> You can redefine the \ operator for your type to not do this. >>>>>> > >> >>>>>> > >> On Thursday, February 20, 2014 6:08:31 PM UTC-6, andrew cooke >>>>>> wrote: >>>>>> > >>> While I still don't see a problem with the theory side of this, >>>>>> there is >>>>>> > >>> a practical problem - Gaussian Elimination in Julia seems to be >>>>>> > >>> delegated >>>>>> > >>> to LAPACK and arguments are promoted to Float: >>>>>> > >>> >>>>>> > >>> julia> [1 0; 0 1] \ [1, 2] >>>>>> > >>> >>>>>> > >>> 2-element Array{Float64,1}: >>>>>> > >>> 1.0 >>>>>> > >>> 2.0 >>>>>> > >>> >>>>>> > >>> Andrew >>>>>> > >>> >>>>>> > >>> On Thursday, 20 February 2014 19:45:28 UTC-3, andrew cooke >>>>>> wrote: >>>>>> > >>>> A broad and a narrow question... >>>>>> > >>>> >>>>>> > >>>> If Julia supports the definition of new integer types can I >>>>>> define a >>>>>> > >>>> new type for a finite field and then use existing linear >>>>>> algebra >>>>>> > >>>> libraries >>>>>> > >>>> to do maths with them? Could I define an integer type for >>>>>> polynomials? >>>>>> > >>>> Is >>>>>> > >>>> this the kind of thing that would work in theory but not in >>>>>> practice? >>>>>> > >>>> Has >>>>>> > >>>> anyone done this? >>>>>> > >>>> >>>>>> > >>>> Specifically, I need to solve a problem modulo 2 (GF(2) - >>>>>> addition and >>>>>> > >>>> subtraction are XOR; multiplication is AND; division is >>>>>> trivial). I >>>>>> > >>>> was >>>>>> > >>>> about to write my own Gaussian Elimination and then remembered >>>>>> a >>>>>> > >>>> comment >>>>>> > >>>> from here saying Julia is the first language where you can >>>>>> define new >>>>>> > >>>> integers... >>>>>> > >>>> >>>>>> > >>>> Am I talking rubbish? I'm not a mathematician, so I may be >>>>>> completely >>>>>> > >>>> muddled anyway. >>>>>> > >>>> >>>>>> > >>>> Thanks, >>>>>> > >>>> Andrew >>>>>> > >>>> >>>>>> > >>>> PS I guess for best speed I should use a Uint for my 0s and 1s? >>>>>> > >>>> Assuming the problem is small enough that it will still fit in >>>>>> cache? >>>>>> >>>>> >>>>> >>>> >>>> >>>> -- >>>> Med venlig hilsen >>>> >>>> Andreas Noack Jensen >>>> >>> -- Med venlig hilsen Andreas Noack Jensen
