From c29c4056a133c732442ed91805201cebd9706ec5 Mon Sep 17 00:00:00 2001 From: Neil Toronto Date: Tue, 1 Jan 2013 12:11:49 -0700 Subject: [PATCH] Added back partial pivoting to Gaussian elimination --- collects/math/matrix.rkt | 17 +++++---- .../math/private/matrix/matrix-gauss-elim.rkt | 38 ++++++++++++------- collects/math/private/matrix/utils.rkt | 15 ++++++++ collects/math/tests/matrix-tests.rkt | 10 ++++- 4 files changed, 58 insertions(+), 22 deletions(-) diff --git a/collects/math/matrix.rkt b/collects/math/matrix.rkt index b238016fca..601130949f 100644 --- a/collects/math/matrix.rkt +++ b/collects/math/matrix.rkt @@ -36,16 +36,19 @@ ) (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" [matrix-gauss-elim - (case-> ((Matrix Number) -> (Values (Matrix Number) (Listof Index))) - ((Matrix Number) Any -> (Values (Matrix Number) (Listof Index))) - ((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index))))] + (case-> ((Matrix Number) -> (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 Pivoting -> (Values (Matrix Number) (Listof Index))))] [matrix-row-echelon - (case-> ((Matrix Number) -> (Matrix Number)) - ((Matrix Number) Any -> (Matrix Number)) - ((Matrix Number) Any Any -> (Matrix Number)))]) + (case-> ((Matrix Number) -> (Matrix Number)) + ((Matrix Number) Any -> (Matrix Number)) + ((Matrix Number) Any Any -> (Matrix Number)) + ((Matrix Number) Any Any Pivoting -> (Matrix Number)))]) (require/untyped-contract (begin (require "private/matrix/matrix-types.rkt")) diff --git a/collects/math/private/matrix/matrix-gauss-elim.rkt b/collects/math/private/matrix/matrix-gauss-elim.rkt index b13e73d6c2..19719ac7d6 100644 --- a/collects/math/private/matrix/matrix-gauss-elim.rkt +++ b/collects/math/private/matrix/matrix-gauss-elim.rkt @@ -9,17 +9,22 @@ "../vector/vector-mutate.rkt") (provide + Pivoting matrix-gauss-elim matrix-row-echelon) +(define-type Pivoting (U 'first 'partial)) + (: matrix-gauss-elim - (case-> ((Matrix Real) -> (Values (Matrix Real) (Listof Index))) - ((Matrix Real) Any -> (Values (Matrix Real) (Listof Index))) + (case-> ((Matrix Real) -> (Values (Matrix Real) (Listof Index))) + ((Matrix Real) Any -> (Values (Matrix Real) (Listof Index))) ((Matrix Real) Any Any -> (Values (Matrix Real) (Listof Index))) - ((Matrix Number) -> (Values (Matrix Number) (Listof Index))) - ((Matrix Number) Any -> (Values (Matrix Number) (Listof Index))) - ((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index))))) -(define (matrix-gauss-elim M [jordan? #f] [unitize-pivot? #f]) + ((Matrix Real) Any Any Pivoting -> (Values (Matrix Real) (Listof Index))) + ((Matrix Number) -> (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 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 rows (matrix->vector* M)) (let loop ([#{i : Nonnegative-Fixnum} 0] @@ -36,7 +41,10 @@ (cond [(j . fx< . n) (loop (fx+ j 1) (cons j without-pivot))] [else (reverse without-pivot)])))] [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 [(zero? pivot) (loop i (fx+ j 1) (cons j without-pivot))] [else @@ -51,12 +59,14 @@ (loop (fx+ i 1) (fx+ j 1) without-pivot))])]))) (: matrix-row-echelon - (case-> ((Matrix Real) -> (Matrix Real)) - ((Matrix Real) Any -> (Matrix Real)) + (case-> ((Matrix Real) -> (Matrix Real)) + ((Matrix Real) Any -> (Matrix Real)) ((Matrix Real) Any Any -> (Matrix Real)) - ((Matrix Number) -> (Matrix Number)) - ((Matrix Number) Any -> (Matrix Number)) - ((Matrix Number) Any Any -> (Matrix Number)))) -(define (matrix-row-echelon M [jordan? #f] [unitize-pivot? jordan?]) - (let-values ([(M _) (matrix-gauss-elim M jordan? unitize-pivot?)]) + ((Matrix Real) Any Any Pivoting -> (Matrix Real)) + ((Matrix Number) -> (Matrix Number)) + ((Matrix Number) Any -> (Matrix Number)) + ((Matrix Number) Any Any -> (Matrix Number)) + ((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)) diff --git a/collects/math/private/matrix/utils.rkt b/collects/math/private/matrix/utils.rkt index c7c14089bc..1e91557363 100644 --- a/collects/math/private/matrix/utils.rkt +++ b/collects/math/private/matrix/utils.rkt @@ -83,6 +83,21 @@ [else (loop (fx+ l 1) p pivot mag-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! (case-> ((Vectorof (Vectorof Real)) Index Index Index Real Nonnegative-Fixnum -> Void) ((Vectorof (Vectorof Number)) Index Index Index Number Nonnegative-Fixnum -> Void))) diff --git a/collects/math/tests/matrix-tests.rkt b/collects/math/tests/matrix-tests.rkt index 81edb48b33..1a385dfbdf 100644 --- a/collects/math/tests/matrix-tests.rkt +++ b/collects/math/tests/matrix-tests.rkt @@ -781,7 +781,15 @@ (check-equal? (matrix-row-echelon (matrix [[ 2 1 -1 8] [-3 -1 2 -11] [-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] [0 1 0 3] [0 0 1 -1]]))