[racket] math/matrix
When I change it to operate on (Vectorof FlVector) instead of (Vectorof
(Vectorof Flonum)), I get this:
cpu time: 996 real time: 995 gc time: 22
1.0000000000009335
cpu time: 15387 real time: 15384 gc time: 13006
1.0000000000009335
cpu time: 1057 real time: 1056 gc time: 85
1.0000000000009335
cpu time: 11514 real time: 11510 gc time: 9097
1.0000000000009335
cpu time: 1079 real time: 1079 gc time: 100
1.0000000000009335
cpu time: 15425 real time: 15426 gc time: 13072
1.0000000000009335
When I use `racket/unsafe/ops` instead of `math/private/unsafe`, I get this:
cpu time: 591 real time: 591 gc time: 17
1.0000000000016622
cpu time: 11514 real time: 11509 gc time: 9195
1.0000000000016622
cpu time: 604 real time: 604 gc time: 31
1.0000000000016622
cpu time: 15739 real time: 15737 gc time: 13358
1.0000000000016622
cpu time: 596 real time: 595 gc time: 24
1.0000000000016622
cpu time: 11498 real time: 11493 gc time: 9154
1.0000000000016622
Racket's floating-point math is fast if the flonums aren't allocated
separately on the heap.
Neil ⊥
On 05/12/2014 02:17 PM, Jens Axel Søgaard wrote:
> Hi Eric,
>
> You were absolute right. The version below cuts the time in half.
> It is mostly cut and paste from existing functions and removing
> non-Flonum cases.
>
> /Jens Axel
>
> #lang typed/racket
> (require math/matrix
> math/array
> math/private/matrix/utils
> math/private/vector/vector-mutate
> math/private/unsafe
> (only-in racket/unsafe/ops unsafe-fl/)
> racket/fixnum
> racket/flonum
> racket/list)
>
> (define-type Pivoting (U 'first 'partial))
>
> (: flonum-matrix-gauss-elim
> (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Listof Index)))
> ((Matrix Flonum) Any -> (Values (Matrix Flonum) (Listof Index)))
> ((Matrix Flonum) Any Any -> (Values (Matrix Flonum) (Listof Index)))
> ((Matrix Flonum) Any Any Pivoting -> (Values (Matrix
> Flonum) (Listof Index)))))
> (define (flonum-matrix-gauss-elim M [jordan? #f] [unitize-pivot? #f]
> [pivoting 'partial])
> (define-values (m n) (matrix-shape M))
> (define rows (matrix->vector* M))
> (let loop ([#{i : Nonnegative-Fixnum} 0]
> [#{j : Nonnegative-Fixnum} 0]
> [#{without-pivot : (Listof Index)} empty])
> (cond
> [(j . fx>= . n)
> (values (vector*->matrix rows)
> (reverse without-pivot))]
> [(i . fx>= . m)
> (values (vector*->matrix rows)
> ;; None of the rest of the columns can have pivots
> (let loop ([#{j : Nonnegative-Fixnum} j] [without-pivot
> without-pivot])
> (cond [(j . fx< . n) (loop (fx+ j 1) (cons j without-pivot))]
> [else (reverse without-pivot)])))]
> [else
> (define-values (p pivot)
> (case pivoting
> [(partial) (find-partial-pivot rows m i j)]
> [(first) (find-first-pivot rows m i j)]))
> (cond
> [(zero? pivot) (loop i (fx+ j 1) (cons j without-pivot))]
> [else
> ;; Swap pivot row with current
> (vector-swap! rows i p)
> ;; Possibly unitize the new current row
> (let ([pivot (if unitize-pivot?
> (begin (vector-scale! (unsafe-vector-ref rows i)
> (unsafe-fl/ 1. pivot))
> (unsafe-fl/ pivot pivot))
> pivot)])
> (flonum-elim-rows! rows m i j pivot (if jordan? 0 (fx+ i 1)))
> (loop (fx+ i 1) (fx+ j 1) without-pivot))])])))
>
> (: flonum-elim-rows!
> ((Vectorof (Vectorof Flonum)) Index Index Index Flonum
> Nonnegative-Fixnum -> Void))
> (define (flonum-elim-rows! rows m i j pivot start)
> (define row_i (unsafe-vector-ref rows i))
> (let loop ([#{l : Nonnegative-Fixnum} start])
> (when (l . fx< . m)
> (unless (l . fx= . i)
> (define row_l (unsafe-vector-ref rows l))
> (define x_lj (unsafe-vector-ref row_l j))
> (unless (= x_lj 0)
> (flonum-vector-scaled-add! row_l row_i (fl* -1. (fl/ x_lj pivot)) j)
> ;; Make sure the element below the pivot is zero
> (unsafe-vector-set! row_l j (- x_lj x_lj))))
> (loop (fx+ l 1)))))
>
>
> (: flonum-matrix-solve
> (All (A) (case->
> ((Matrix Flonum) (Matrix Flonum) -> (Matrix Flonum))
> ((Matrix Flonum) (Matrix Flonum) (-> A) -> (U A (Matrix
> Flonum))))))
> (define flonum-matrix-solve
> (case-lambda
> [(M B) (flonum-matrix-solve
> M B (λ () (raise-argument-error 'flonum-matrix-solve
> "matrix-invertible?" 0 M B)))]
> [(M B fail)
> (define m (square-matrix-size M))
> (define-values (s t) (matrix-shape B))
> (cond [(= m s)
> (define-values (IX wps)
> (parameterize ([array-strictness #f])
> (flonum-matrix-gauss-elim (matrix-augment (list M B)) #t #t)))
> (cond [(and (not (empty? wps)) (= (first wps) m))
> (submatrix IX (::) (:: m #f))]
> [else (fail)])]
> [else
> (error 'flonum-matrix-solve
> "matrices must have the same number of rows; given ~e and ~e"
> M B)])]))
>
> (define-syntax-rule (flonum-vector-generic-scaled-add! vs0-expr
> vs1-expr v-expr start-expr + *)
> (let* ([vs0 vs0-expr]
> [vs1 vs1-expr]
> [v v-expr]
> [n (fxmin (vector-length vs0) (vector-length vs1))])
> (let loop ([#{i : Nonnegative-Fixnum} (fxmin start-expr n)])
> (if (i . fx< . n)
> (begin (unsafe-vector-set! vs0 i (+ (unsafe-vector-ref vs0 i)
> (* (unsafe-vector-ref vs1 i) v)))
> (loop (fx+ i 1)))
> (void)))))
>
> (: flonum-vector-scaled-add!
> (case-> ((Vectorof Flonum) (Vectorof Flonum) Flonum -> Void)
> ((Vectorof Flonum) (Vectorof Flonum) Flonum Index -> Void)))
> (define (flonum-vector-scaled-add! vs0 vs1 s [start 0])
> (flonum-vector-generic-scaled-add! vs0 vs1 s start + *))
>
> (: mx Index)
> (define mx 600)
>
> (: r (Index Index -> Flonum))
> (define (r i j) (random))
>
> (: A : (Matrix Flonum))
> (define A (build-matrix mx mx r))
>
> (: sum : Integer Integer -> Flonum)
> (define (sum i n)
> (let loop ((j 0) (acc 0.0))
> (if (>= j mx) acc
> (loop (+ j 1) (+ acc (matrix-ref A i j))) )))
>
> (: b : (Matrix Flonum))
> (define b (build-matrix mx 1 sum))
>
> (time
> (let [(m (flonum-matrix-solve A b))]
> (matrix-ref m 0 0)))
> (time
> (let [(m (matrix-solve A b))]
> (matrix-ref m 0 0)))
>
> (time
> (let [(m (flonum-matrix-solve A b))]
> (matrix-ref m 0 0)))
> (time
> (let [(m (matrix-solve A b))]
> (matrix-ref m 0 0)))
>
> (time
> (let [(m (flonum-matrix-solve A b))]
> (matrix-ref m 0 0)))
> (time
> (let [(m (matrix-solve A b))]
> (matrix-ref m 0 0)))
>
> /Jens Axel
>
>
> 2014-05-11 23:26 GMT+02:00 Eric Dobson <eric.n.dobson at gmail.com>:
>> Where is the time spent in the algorithm? I assume that most of it is
>> in the matrix manipulation work not the orchestration of finding a
>> pivot and reducing based on that. I.e. `elim-rows!` is the expensive
>> part. Given that you only specialized the outer part of the loop, I
>> wouldn't expect huge performance changes.
>>
>> On Sun, May 11, 2014 at 2:13 PM, Jens Axel Søgaard
>> <jensaxel at soegaard.net> wrote:
>>> I tried restricting the matrix-solve and matrix-gauss-elim to (Matrix Flonum).
>>> I can't observe a change in the timings.
>>>
>>> #lang typed/racket
>>> (require math/matrix
>>> math/array
>>> math/private/matrix/utils
>>> math/private/vector/vector-mutate
>>> math/private/unsafe
>>> (only-in racket/unsafe/ops unsafe-fl/)
>>> racket/fixnum
>>> racket/list)
>>>
>>> (define-type Pivoting (U 'first 'partial))
>>>
>>> (: flonum-matrix-gauss-elim
>>> (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Listof Index)))
>>> ((Matrix Flonum) Any -> (Values (Matrix Flonum) (Listof Index)))
>>> ((Matrix Flonum) Any Any -> (Values (Matrix Flonum) (Listof Index)))
>>> ((Matrix Flonum) Any Any Pivoting -> (Values (Matrix
>>> Flonum) (Listof Index)))))
>>> (define (flonum-matrix-gauss-elim M [jordan? #f] [unitize-pivot? #f]
>>> [pivoting 'partial])
>>> (define-values (m n) (matrix-shape M))
>>> (define rows (matrix->vector* M))
>>> (let loop ([#{i : Nonnegative-Fixnum} 0]
>>> [#{j : Nonnegative-Fixnum} 0]
>>> [#{without-pivot : (Listof Index)} empty])
>>> (cond
>>> [(j . fx>= . n)
>>> (values (vector*->matrix rows)
>>> (reverse without-pivot))]
>>> [(i . fx>= . m)
>>> (values (vector*->matrix rows)
>>> ;; None of the rest of the columns can have pivots
>>> (let loop ([#{j : Nonnegative-Fixnum} j] [without-pivot
>>> without-pivot])
>>> (cond [(j . fx< . n) (loop (fx+ j 1) (cons j without-pivot))]
>>> [else (reverse without-pivot)])))]
>>> [else
>>> (define-values (p pivot)
>>> (case pivoting
>>> [(partial) (find-partial-pivot rows m i j)]
>>> [(first) (find-first-pivot rows m i j)]))
>>> (cond
>>> [(zero? pivot) (loop i (fx+ j 1) (cons j without-pivot))]
>>> [else
>>> ;; Swap pivot row with current
>>> (vector-swap! rows i p)
>>> ;; Possibly unitize the new current row
>>> (let ([pivot (if unitize-pivot?
>>> (begin (vector-scale! (unsafe-vector-ref rows i)
>>> (unsafe-fl/ 1. pivot))
>>> (unsafe-fl/ pivot pivot))
>>> pivot)])
>>> (elim-rows! rows m i j pivot (if jordan? 0 (fx+ i 1)))
>>> (loop (fx+ i 1) (fx+ j 1) without-pivot))])])))
>>>
>>> (: flonum-matrix-solve
>>> (All (A) (case->
>>> ((Matrix Flonum) (Matrix Flonum) -> (Matrix Flonum))
>>> ((Matrix Flonum) (Matrix Flonum) (-> A) -> (U A (Matrix
>>> Flonum))))))
>>> (define flonum-matrix-solve
>>> (case-lambda
>>> [(M B) (flonum-matrix-solve
>>> M B (λ () (raise-argument-error 'flonum-matrix-solve
>>> "matrix-invertible?" 0 M B)))]
>>> [(M B fail)
>>> (define m (square-matrix-size M))
>>> (define-values (s t) (matrix-shape B))
>>> (cond [(= m s)
>>> (define-values (IX wps)
>>> (parameterize ([array-strictness #f])
>>> (flonum-matrix-gauss-elim (matrix-augment (list M B)) #t #t)))
>>> (cond [(and (not (empty? wps)) (= (first wps) m))
>>> (submatrix IX (::) (:: m #f))]
>>> [else (fail)])]
>>> [else
>>> (error 'flonum-matrix-solve
>>> "matrices must have the same number of rows; given ~e and ~e"
>>> M B)])]))
>>>
>>> (: mx Index)
>>> (define mx 600)
>>>
>>> (: r (Index Index -> Flonum))
>>> (define (r i j) (random))
>>>
>>> (: A : (Matrix Flonum))
>>> (define A (build-matrix mx mx r))
>>>
>>> (: sum : Integer Integer -> Flonum)
>>> (define (sum i n)
>>> (let loop ((j 0) (acc 0.0))
>>> (if (>= j mx) acc
>>> (loop (+ j 1) (+ acc (matrix-ref A i j))) )))
>>>
>>> (: b : (Matrix Flonum))
>>> (define b (build-matrix mx 1 sum))
>>>
>>> (time
>>> (let [(m (flonum-matrix-solve A b))]
>>> (matrix-ref m 0 0)))
>>> (time
>>> (let [(m (matrix-solve A b))]
>>> (matrix-ref m 0 0)))
>>>
>>> (time
>>> (let [(m (flonum-matrix-solve A b))]
>>> (matrix-ref m 0 0)))
>>> (time
>>> (let [(m (matrix-solve A b))]
>>> (matrix-ref m 0 0)))
>>>
>>> (time
>>> (let [(m (flonum-matrix-solve A b))]
>>> (matrix-ref m 0 0)))
>>> (time
>>> (let [(m (matrix-solve A b))]
>>> (matrix-ref m 0 0)))
>>>
>>> 2014-05-11 21:48 GMT+02:00 Neil Toronto <neil.toronto at gmail.com>:
>>>> The garbage collection time is probably from cleaning up boxed flonums, and
>>>> possibly intermediate vectors. If so, a separate implementation of Gaussian
>>>> elimination for the FlArray type would cut the GC time to nearly zero.
>>>>
>>>> Neil ⊥
>>>>
>>>>
>>>> On 05/11/2014 01:36 PM, Jens Axel Søgaard wrote:
>>>>>
>>>>> Or ... you could take a look at
>>>>>
>>>>>
>>>>> https://github.com/plt/racket/blob/master/pkgs/math-pkgs/math-lib/math/private/matrix/matrix-gauss-elim.rkt
>>>>>
>>>>> at see if something can be improved.
>>>>>
>>>>> /Jens Axel
>>>>>
>>>>>
>>>>> 2014-05-11 21:30 GMT+02:00 Jens Axel Søgaard <jensaxel at soegaard.net>:
>>>>>>
>>>>>> Hi Eduardo,
>>>>>>
>>>>>> The math/matrix library uses the arrays from math/array to represent
>>>>>> matrices.
>>>>>>
>>>>>> If you want to try the same representation as Bigloo, you could try Will
>>>>>> Farr's
>>>>>> matrix library:
>>>>>>
>>>>>>
>>>>>> http://planet.racket-lang.org/package-source/wmfarr/simple-matrix.plt/1/1/planet-docs/simple-matrix/index.html
>>>>>>
>>>>>> I am interested in hearing the results.
>>>>>>
>>>>>> /Jens Axel
>>>>>>
>>>>>>
>>>>>>
>>>>>> 2014-05-11 21:18 GMT+02:00 Eduardo Costa <edu500ac at gmail.com>:
>>>>>>>
>>>>>>> What is bothering me is the time Racket is spending in garbage
>>>>>>> collection.
>>>>>>>
>>>>>>> ~/wrk/scm/rkt/matrix$ racket matrix.rkt
>>>>>>> 0.9999999999967226
>>>>>>> cpu time: 61416 real time: 61214 gc time: 32164
>>>>>>>
>>>>>>> If I am reading the output correctly, Racket is spending 32 seconds out
>>>>>>> of
>>>>>>> 61 seconds in garbage collection.
>>>>>>>
>>>>>>> I am following Junia Magellan's computer language comparison and I
>>>>>>> cannot
>>>>>>> understand why Racket needs the garbage collector for doing Gaussian
>>>>>>> elimination. In a slow Compaq/HP machine, solving a system of 800 linear
>>>>>>> equations takes 17.3 seconds in Bigloo, but requires 58 seconds in
>>>>>>> Racket,
>>>>>>> even after removing the building of the linear system from
>>>>>>> consideration.
>>>>>>> Common Lisp is also much faster than Racket in processing arrays. I
>>>>>>> would
>>>>>>> like to point out that Racket is very fast in general. The only occasion
>>>>>>> that it lags badly behind Common Lisp and Bigloo is when one needs to
>>>>>>> deal
>>>>>>> with arrays.
>>>>>>>
>>>>>>> Basically, Junia is using Rasch method to measure certain latent traits
>>>>>>> of
>>>>>>> computer languages, like productivity and coaching time. In any case,
>>>>>>> she
>>>>>>> needs to do a lot of matrix calculations to invert the Rasch model.
>>>>>>> Since
>>>>>>> Bigloo works with homogeneous vectors, she wrote a few macros to access
>>>>>>> the
>>>>>>> elements of a matrix:
>>>>>>>
>>>>>>> (define (mkv n) (make-f64vector n))
>>>>>>> (define $ f64vector-ref)
>>>>>>> (define $! f64vector-set!)
>>>>>>> (define len f64vector-length)
>>>>>>>
>>>>>>> (define-syntax $$
>>>>>>> (syntax-rules ()
>>>>>>> (($$ m i j) (f64vector-ref (vector-ref m i) j))))
>>>>>>>
>>>>>>> (define-syntax $$!
>>>>>>> (syntax-rules ()
>>>>>>> (($$! matrix row column value)
>>>>>>> ($! (vector-ref matrix row) column value))))
>>>>>>>
>>>>>>> I wonder whether homogeneous vectors would speed up Racket. In the same
>>>>>>> computer that Racket takes 80 seconds to build and invert a system of
>>>>>>> equations, Bigloo takes 17.3 seconds, as I told before. Common Lisp is
>>>>>>> even
>>>>>>> faster. However, if one subtracts the gc time from Racket's total time,
>>>>>>> the
>>>>>>> result comes quite close to Common Lisp or Bigloo.
>>>>>>>
>>>>>>> ~/wrk/bgl$ bigloo -Obench bigmat.scm -o big
>>>>>>> ~/wrk/bgl$ time ./big
>>>>>>> 0.9999999999965746 1.000000000000774 0.9999999999993039
>>>>>>> 0.9999999999982576
>>>>>>> 1.000000000007648 0.999999999996588
>>>>>>>
>>>>>>> real 0m17.423s
>>>>>>> user 0m17.384s
>>>>>>> sys 0m0.032s
>>>>>>> ~/wrk/bgl$
>>>>>>>
>>>>>>> Well, bigloo may perform global optimizations, but Common Lisp doesn't.
>>>>>>> When
>>>>>>> one is not dealing with matrices, Racket is faster than Common Lisp. I
>>>>>>> hope
>>>>>>> you can tell me how to rewrite the program in order to avoid garbage
>>>>>>> collection.
>>>>>>>
>>>>>>> By the way, you may want to know why not use Bigloo or Common Lisp to
>>>>>>> invert
>>>>>>> the Rasch model. The problem is that Junia and her co-workers are using
>>>>>>> hosting services that do not give access to the server or to the
>>>>>>> jailshell.
>>>>>>> Since Bigloo requires gcc based compilation, Junia discarded it right
>>>>>>> away.
>>>>>>> Not long ago, the hosting service stopped responding to the sbcl Common
>>>>>>> Lisp
>>>>>>> compiler for reasons that I cannot fathom. Although Racket 6.0 stopped
>>>>>>> working too, Racket 6.0.1 is working fine. This left Junia, her
>>>>>>> co-workers
>>>>>>> and students with Racket as their sole option. As for myself, I am just
>>>>>>> curious.
>>>>>>>
>>>>>>>
>>>>>>> 2014-05-11 6:23 GMT-03:00 Jens Axel Søgaard <jensaxel at soegaard.net>:
>>>>>>>
>>>>>>>> 2014-05-11 6:09 GMT+02:00 Eduardo Costa <edu500ac at gmail.com>:
>>>>>>>>>
>>>>>>>>> The documentation says that one should expect typed/racket to be
>>>>>>>>> faster
>>>>>>>>> than
>>>>>>>>> racket. I tested the math/matrix library and it seems to be almost as
>>>>>>>>> slow
>>>>>>>>> in typed/racket as in racket.
>>>>>>>>
>>>>>>>>
>>>>>>>> What was (is?) slow was a call in an untyped module A to a function
>>>>>>>> exported
>>>>>>>> from a typed module B. The functions in B must check at runtime that
>>>>>>>> the values coming from A are of the correct type. If the A was written
>>>>>>>> in Typed Racket, the types would be known at compile time.
>>>>>>>>
>>>>>>>> Here math/matrix is written in Typed Racket, so if you are writing an
>>>>>>>> untyped module, you will in general want to minimize the use of,say,
>>>>>>>> maxtrix-ref. Instead operations that works on entire matrices or
>>>>>>>> row/columns are preferred.
>>>>>>>>
>>>>>>>>> (: sum : Integer Integer -> Flonum)
>>>>>>>>> (define (sum i n)
>>>>>>>>> (let loop ((j 0) (acc 0.0))
>>>>>>>>> (if (>= j mx) acc
>>>>>>>>> (loop (+ j 1) (+ acc (matrix-ref A i j))) )))
>>>>>>>>>
>>>>>>>>> (: b : (Matrix Flonum))
>>>>>>>>> (define b (build-matrix mx 1 sum))
>>>>>>>>
>>>>>>>>
>>>>>>>> The matrix b contains the sums of each row in the matrix.
>>>>>>>> Since matrices are a subset of arrays, you can use array-axis-sum,
>>>>>>>> which computes sum along a given axis (i.e. a row or a column when
>>>>>>>> speaking of matrices).
>>>>>>>>
>>>>>>>> (define A (matrix [[0. 1. 2.]
>>>>>>>> [3. 4. 5.]
>>>>>>>> [6. 7. 8.]]))
>>>>>>>>
>>>>>>>>> (array-axis-sum A 1)
>>>>>>>>
>>>>>>>> - : (Array Flonum)
>>>>>>>> (array #[3.0 12.0 21.0])
>>>>>>>>
>>>>>>>> However as Eric points out, matrix-solve is an O(n^3) algorithm,
>>>>>>>> so the majority of the time is spent in matrix-solve.
>>>>>>>>
>>>>>>>> Apart from finding a way to exploit the relationship between your
>>>>>>>> matrix A and the column vector b, I see no obvious way of
>>>>>>>> speeding up the code.
>>>>>>>>
>>>>>>>> Note that when you benchmark with
>>>>>>>>
>>>>>>>> time racket matrix.rkt
>>>>>>>>
>>>>>>>> you will include startup and compilation time.
>>>>>>>> Therefore if you want to time the matrix code,
>>>>>>>> insert a literal (time ...) call.
>>>>>>>>
>>>>>>>> --
>>>>>>>> Jens Axel Søgaard
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> --
>>>>>> --
>>>>>> Jens Axel Søgaard
>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>> ____________________
>>>> Racket Users list:
>>>> http://lists.racket-lang.org/users
>>>
>>>
>>>
>>> --
>>> --
>>> Jens Axel Søgaard
>>>
>>> ____________________
>>> Racket Users list:
>>> http://lists.racket-lang.org/users
>
>
>
-------------- next part --------------
#lang typed/racket
(require math/matrix
math/array
math/flonum
math/private/vector/vector-mutate
racket/fixnum
racket/flonum
racket/list
racket/unsafe/ops)
(define-type Pivoting (U 'first 'partial))
(: matrix->flvector* ((Matrix Flonum) -> (Vectorof FlVector)))
(define (matrix->flvector* M)
(list->vector (map list->flvector (matrix->list* M))))
(: flvector*->matrix ((Vectorof FlVector) -> (Matrix Flonum)))
(define (flvector*->matrix rows)
(vector*->matrix (vector-map flvector->vector rows)))
(: flonum-matrix-gauss-elim
(case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Listof Index)))
((Matrix Flonum) Any -> (Values (Matrix Flonum) (Listof Index)))
((Matrix Flonum) Any Any -> (Values (Matrix Flonum) (Listof Index)))
((Matrix Flonum) Any Any Pivoting -> (Values (Matrix Flonum) (Listof Index)))))
(define (flonum-matrix-gauss-elim M [jordan? #f] [unitize-pivot? #f] [pivoting 'partial])
(define-values (m n) (matrix-shape M))
(define rows (matrix->flvector* M))
(let loop ([#{i : Nonnegative-Fixnum} 0]
[#{j : Nonnegative-Fixnum} 0]
[#{without-pivot : (Listof Index)} empty])
(cond
[(j . fx>= . n)
(values (flvector*->matrix rows)
(reverse without-pivot))]
[(i . fx>= . m)
(values (flvector*->matrix rows)
;; None of the rest of the columns can have pivots
(let loop ([#{j : Nonnegative-Fixnum} j] [without-pivot
without-pivot])
(cond [(j . fx< . n) (loop (fx+ j 1) (cons j without-pivot))]
[else (reverse without-pivot)])))]
[else
(define-values (p pivot)
(case pivoting
[(partial) (flonum-find-partial-pivot rows m i j)]
[(first) (flonum-find-first-pivot rows m i j)]))
(cond
[(zero? pivot) (loop i (fx+ j 1) (cons j without-pivot))]
[else
;; Swap pivot row with current
(vector-swap! rows i p)
;; Possibly unitize the new current row
(let ([pivot (if unitize-pivot?
(begin (flvector-scale! (unsafe-vector-ref rows i)
(/ 1.0 pivot))
1.0)
pivot)])
(flonum-elim-rows! rows m i j pivot (if jordan? 0 (fx+ i 1)))
(loop (fx+ i 1) (fx+ j 1) without-pivot))])])))
(: flvector-scale! (FlVector Flonum -> Void))
(define (flvector-scale! vs x)
(define n (flvector-length vs))
(let loop ([i : Nonnegative-Fixnum 0])
(when (< i n)
(unsafe-flvector-set! vs i (* x (unsafe-flvector-ref vs i)))
(loop (+ i 1)))))
(: unsafe-flvector2d-ref ((Vectorof FlVector) Index Index -> Flonum))
(define (unsafe-flvector2d-ref vss i j)
(unsafe-flvector-ref (unsafe-vector-ref vss i) j))
(: flonum-find-partial-pivot ((Vectorof FlVector) Index Index Index -> (Values Index Flonum)))
;; Find the element with maximum magnitude in a column
(define (flonum-find-partial-pivot rows m i j)
(define l (fx+ i 1))
(define pivot (unsafe-flvector2d-ref rows i j))
(define mag-pivot (magnitude pivot))
(let loop ([#{l : Nonnegative-Fixnum} l] [#{p : Index} i] [pivot pivot] [mag-pivot mag-pivot])
(cond [(l . fx< . m)
(define new-pivot (unsafe-flvector2d-ref rows l j))
(define mag-new-pivot (magnitude new-pivot))
(cond [(mag-new-pivot . > . mag-pivot) (loop (fx+ l 1) l new-pivot mag-new-pivot)]
[else (loop (fx+ l 1) p pivot mag-pivot)])]
[else (values p pivot)])))
(: flonum-find-first-pivot ((Vectorof FlVector) Index Index Index -> (Values Index Flonum)))
;; Find the first nonzero element in a column
(define (flonum-find-first-pivot rows m i j)
(define pivot (unsafe-flvector2d-ref rows i j))
(if ((magnitude pivot) . > . 0)
(values i pivot)
(let loop ([#{l : Nonnegative-Fixnum} (fx+ i 1)])
(cond [(l . fx< . m)
(define pivot (unsafe-flvector2d-ref rows l j))
(if ((magnitude pivot) . > . 0) (values l pivot) (loop (fx+ l 1)))]
[else
(values i pivot)]))))
(: flonum-elim-rows! ((Vectorof FlVector) Index Index Index Flonum Nonnegative-Fixnum
-> Void))
(define (flonum-elim-rows! rows m i j pivot start)
(define row_i (unsafe-vector-ref rows i))
(let loop ([#{l : Nonnegative-Fixnum} start])
(when (l . fx< . m)
(unless (l . fx= . i)
(define row_l (unsafe-vector-ref rows l))
(define x_lj (unsafe-flvector-ref row_l j))
(unless (= x_lj 0)
(flonum-vector-scaled-add! row_l row_i (fl* -1. (fl/ x_lj pivot)) j)
;; Make sure the element below the pivot is zero
(unsafe-flvector-set! row_l j (- x_lj x_lj))))
(loop (fx+ l 1)))))
(: flonum-matrix-solve
(All (A) (case->
((Matrix Flonum) (Matrix Flonum) -> (Matrix Flonum))
((Matrix Flonum) (Matrix Flonum) (-> A) -> (U A (Matrix
Flonum))))))
(define flonum-matrix-solve
(case-lambda
[(M B) (flonum-matrix-solve
M B (λ () (raise-argument-error 'flonum-matrix-solve
"matrix-invertible?" 0 M B)))]
[(M B fail)
(define m (square-matrix-size M))
(define-values (s t) (matrix-shape B))
(cond [(= m s)
(define-values (IX wps)
(parameterize ([array-strictness #f])
(flonum-matrix-gauss-elim (matrix-augment (list M B)) #t #t)))
(cond [(and (not (empty? wps)) (= (first wps) m))
(submatrix IX (::) (:: m #f))]
[else (fail)])]
[else
(error 'flonum-matrix-solve
"matrices must have the same number of rows; given ~e and ~e"
M B)])]))
(define-syntax-rule (flonum-vector-generic-scaled-add! vs0-expr vs1-expr v-expr start-expr + *)
(let* ([vs0 vs0-expr]
[vs1 vs1-expr]
[v v-expr]
[n (fxmin (flvector-length vs0) (flvector-length vs1))])
(let loop ([#{i : Nonnegative-Fixnum} (fxmin start-expr n)])
(if (i . fx< . n)
(begin (unsafe-flvector-set! vs0 i (+ (unsafe-flvector-ref vs0 i)
(* (unsafe-flvector-ref vs1 i) v)))
(loop (fx+ i 1)))
(void)))))
(: flonum-vector-scaled-add!
(case-> (FlVector FlVector Flonum -> Void)
(FlVector FlVector Flonum Index -> Void)))
(define (flonum-vector-scaled-add! vs0 vs1 s [start 0])
(flonum-vector-generic-scaled-add! vs0 vs1 s start + *))
(: mx Index)
(define mx 600)
(: r (Index Index -> Flonum))
(define (r i j) (random))
(: A : (Matrix Flonum))
(define A (build-matrix mx mx r))
(: sum : Integer Integer -> Flonum)
(define (sum i n)
(let loop ((j 0) (acc 0.0))
(if (>= j mx) acc
(loop (+ j 1) (+ acc (matrix-ref A i j))) )))
(: b : (Matrix Flonum))
(define b (build-matrix mx 1 sum))
(time
(let [(m (flonum-matrix-solve A b))]
(matrix-ref m 0 0)))
(time
(let [(m (matrix-solve A b))]
(matrix-ref m 0 0)))
(time
(let [(m (flonum-matrix-solve A b))]
(matrix-ref m 0 0)))
(time
(let [(m (matrix-solve A b))]
(matrix-ref m 0 0)))
(time
(let [(m (flonum-matrix-solve A b))]
(matrix-ref m 0 0)))
(time
(let [(m (matrix-solve A b))]
(matrix-ref m 0 0)))
![[logo]](/logo.png){kind=logo}