suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added back partial pivoting to Gaussian elimination

Browse Source
This commit is contained in:
Neil Toronto 2013-01-01 12:11:49 -07:00
parent 0097d59635
commit c29c4056a1
4 changed files with 58 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
collects/math/matrix.rkt
View File

@ -36,16 +36,19 @@
) )
(require/untyped-contract (require/untyped-contract
(begin (require "private/matrix/matrix-types.rkt")) (begin (require "private/matrix/matrix-types.rkt"
"private/matrix/matrix-gauss-elim.rkt"))
"private/matrix/matrix-gauss-elim.rkt" "private/matrix/matrix-gauss-elim.rkt"
[matrix-gauss-elim [matrix-gauss-elim
(case-> ((Matrix Number) -> (Values (Matrix Number) (Listof Index))) (case-> ((Matrix Number) -> (Values (Matrix Number) (Listof Index)))
((Matrix Number) Any -> (Values (Matrix Number) (Listof Index))) ((Matrix Number) Any -> (Values (Matrix Number) (Listof Index)))
((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index))))] ((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index)))
((Matrix Number) Any Any Pivoting -> (Values (Matrix Number) (Listof Index))))]
[matrix-row-echelon [matrix-row-echelon
(case-> ((Matrix Number) -> (Matrix Number)) (case-> ((Matrix Number) -> (Matrix Number))
((Matrix Number) Any -> (Matrix Number)) ((Matrix Number) Any -> (Matrix Number))
((Matrix Number) Any Any -> (Matrix Number)))]) ((Matrix Number) Any Any -> (Matrix Number))
((Matrix Number) Any Any Pivoting -> (Matrix Number)))])
(require/untyped-contract (require/untyped-contract
(begin (require "private/matrix/matrix-types.rkt")) (begin (require "private/matrix/matrix-types.rkt"))

38
collects/math/private/matrix/matrix-gauss-elim.rkt
View File

@ -9,17 +9,22 @@
"../vector/vector-mutate.rkt") "../vector/vector-mutate.rkt")
(provide (provide
Pivoting
matrix-gauss-elim matrix-gauss-elim
matrix-row-echelon) matrix-row-echelon)
(define-type Pivoting (U 'first 'partial))
(: matrix-gauss-elim (: matrix-gauss-elim
(case-> ((Matrix Real) -> (Values (Matrix Real) (Listof Index))) (case-> ((Matrix Real) -> (Values (Matrix Real) (Listof Index)))
((Matrix Real) Any -> (Values (Matrix Real) (Listof Index))) ((Matrix Real) Any -> (Values (Matrix Real) (Listof Index)))
((Matrix Real) Any Any -> (Values (Matrix Real) (Listof Index))) ((Matrix Real) Any Any -> (Values (Matrix Real) (Listof Index)))
((Matrix Number) -> (Values (Matrix Number) (Listof Index))) ((Matrix Real) Any Any Pivoting -> (Values (Matrix Real) (Listof Index)))
((Matrix Number) Any -> (Values (Matrix Number) (Listof Index))) ((Matrix Number) -> (Values (Matrix Number) (Listof Index)))
((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index))))) ((Matrix Number) Any -> (Values (Matrix Number) (Listof Index)))
(define (matrix-gauss-elim M [jordan? #f] [unitize-pivot? #f]) ((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index)))
((Matrix Number) Any Any Pivoting -> (Values (Matrix Number) (Listof Index)))))
(define (matrix-gauss-elim M [jordan? #f] [unitize-pivot? #f] [pivoting 'partial])
(define-values (m n) (matrix-shape M)) (define-values (m n) (matrix-shape M))
(define rows (matrix->vector* M)) (define rows (matrix->vector* M))
(let loop ([#{i : Nonnegative-Fixnum} 0] (let loop ([#{i : Nonnegative-Fixnum} 0]
@ -36,7 +41,10 @@
(cond [(j . fx< . n) (loop (fx+ j 1) (cons j without-pivot))] (cond [(j . fx< . n) (loop (fx+ j 1) (cons j without-pivot))]
[else (reverse without-pivot)])))] [else (reverse without-pivot)])))]
[else [else
(define-values (p pivot) (find-partial-pivot rows m i j)) (define-values (p pivot)
(case pivoting
[(partial) (find-partial-pivot rows m i j)]
[(first) (find-first-pivot rows m i j)]))
(cond (cond
[(zero? pivot) (loop i (fx+ j 1) (cons j without-pivot))] [(zero? pivot) (loop i (fx+ j 1) (cons j without-pivot))]
[else [else
@ -51,12 +59,14 @@
(loop (fx+ i 1) (fx+ j 1) without-pivot))])]))) (loop (fx+ i 1) (fx+ j 1) without-pivot))])])))
(: matrix-row-echelon (: matrix-row-echelon
(case-> ((Matrix Real) -> (Matrix Real)) (case-> ((Matrix Real) -> (Matrix Real))
((Matrix Real) Any -> (Matrix Real)) ((Matrix Real) Any -> (Matrix Real))
((Matrix Real) Any Any -> (Matrix Real)) ((Matrix Real) Any Any -> (Matrix Real))
((Matrix Number) -> (Matrix Number)) ((Matrix Real) Any Any Pivoting -> (Matrix Real))
((Matrix Number) Any -> (Matrix Number)) ((Matrix Number) -> (Matrix Number))
((Matrix Number) Any Any -> (Matrix Number)))) ((Matrix Number) Any -> (Matrix Number))
(define (matrix-row-echelon M [jordan? #f] [unitize-pivot? jordan?]) ((Matrix Number) Any Any -> (Matrix Number))
(let-values ([(M _) (matrix-gauss-elim M jordan? unitize-pivot?)]) ((Matrix Number) Any Any Pivoting -> (Matrix Number))))
(define (matrix-row-echelon M [jordan? #f] [unitize-pivot? #f] [pivoting 'partial])
(let-values ([(M _) (matrix-gauss-elim M jordan? unitize-pivot? pivoting)])
M)) M))

15
collects/math/private/matrix/utils.rkt
View File

@ -83,6 +83,21 @@
[else (loop (fx+ l 1) p pivot mag-pivot)])] [else (loop (fx+ l 1) p pivot mag-pivot)])]
[else (values p pivot)]))) [else (values p pivot)])))
(: find-first-pivot
(case-> ((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real))
((Vectorof (Vectorof Number)) Index Index Index -> (Values Index Number))))
;; Find the first nonzero element in a column
(define (find-first-pivot rows m i j)
(define pivot (unsafe-vector2d-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-vector2d-ref rows l j))
(if ((magnitude pivot) . > . 0) (values l pivot) (loop (fx+ l 1)))]
[else
(values i pivot)]))))
(: elim-rows! (: elim-rows!
(case-> ((Vectorof (Vectorof Real)) Index Index Index Real Nonnegative-Fixnum -> Void) (case-> ((Vectorof (Vectorof Real)) Index Index Index Real Nonnegative-Fixnum -> Void)
((Vectorof (Vectorof Number)) Index Index Index Number Nonnegative-Fixnum -> Void))) ((Vectorof (Vectorof Number)) Index Index Index Number Nonnegative-Fixnum -> Void)))

10
collects/math/tests/matrix-tests.rkt
View File

@ -781,7 +781,15 @@
(check-equal? (matrix-row-echelon (matrix [[ 2 1 -1 8] (check-equal? (matrix-row-echelon (matrix [[ 2 1 -1 8]
[-3 -1 2 -11] [-3 -1 2 -11]
[-2 1 2 -3]]) [-2 1 2 -3]])
#t) #t #t 'partial)
(matrix [[1 0 0 2]
[0 1 0 3]
[0 0 1 -1]]))
(check-equal? (matrix-row-echelon (matrix [[ 2 1 -1 8]
[-3 -1 2 -11]
[-2 1 2 -3]])
#t #t 'first)
(matrix [[1 0 0 2] (matrix [[1 0 0 2]
[0 1 0 3] [0 1 0 3]
[0 0 1 -1]])) [0 0 1 -1]]))