Permalink
Browse files

Merge pull request #5 from soegaard/master

New matrix operations
  • Loading branch information...
2 parents b44b883 + 5ebb616 commit d78af9105d9cc0b7a396df71cdaa5c883b1dfc8a @ntoronto ntoronto committed Jul 14, 2012
@@ -3,11 +3,17 @@
(require "private/matrix/matrix-pointwise.rkt"
"private/matrix/matrix-multiply.rkt"
"private/matrix/matrix-constructors.rkt"
+ "private/matrix/matrix-operations.rkt"
+ "private/matrix/matrix-expt.rkt"
+ "private/matrix/matrix-types.rkt"
"private/matrix/utils.rkt")
(provide (all-from-out "private/matrix/matrix-pointwise.rkt"
"private/matrix/matrix-multiply.rkt"
- "private/matrix/matrix-constructors.rkt")
+ "private/matrix/matrix-constructors.rkt"
+ "private/matrix/matrix-operations.rkt"
+ "private/matrix/matrix-expt.rkt"
+ "private/matrix/matrix-types.rkt")
;; From "utils.rkt"
array-matrix?
;; would also like matrix? : (Any -> Boolean : (Array Any)), but we can't have one until we
@@ -1,12 +1,16 @@
-#lang typed/racket/base
+#lang typed/racket
-(require "../../array.rkt"
+(require racket/unsafe/ops
+ math/array
"matrix-types.rkt")
(provide identity-matrix flidentity-matrix
matrix->list list->matrix fllist->matrix
matrix->vector vector->matrix flvector->matrix
- const-matrix)
+ const-matrix
+ matrix-row
+ matrix-col
+ submatrix)
(: identity-matrix (Integer -> (Result-Matrix Real)))
(define (identity-matrix size) (diagonal-array 2 size 1 0))
@@ -41,3 +45,21 @@
(: matrix->vector : (All (A) (Matrix A) -> (Vectorof* A)))
(define (matrix->vector a)
(array->vector a))
+
+(: submatrix : (Matrix Number) (Sequenceof Index) (Sequenceof Index) -> (Result-Matrix Number))
+(define (submatrix a row-range col-range)
+ ((inst array-slice Number) a (list row-range col-range)))
+
+(: matrix-row : (Matrix Number) Index -> (Result-Matrix Number))
+(define (matrix-row a i)
+ (define ds (matrix-dimensions a))
+ (define col-dim (vector-ref ds 1))
+ ((inst array-slice Number) a (list (in-range i (add1 i)) (in-range col-dim))))
+
+(: matrix-col : (Matrix Number) Index -> (Result-Matrix Number))
+(define (matrix-col a j)
+ (define ds (matrix-dimensions a))
+ (define row-dim (vector-ref ds 0))
+ ((inst array-slice Number) a (list (in-range row-dim) (in-range j (add1 j)))))
+
+
@@ -0,0 +1,263 @@
+#lang typed/racket
+
+(require math/array)
+(require "matrix-types.rkt"
+ "matrix-constructors.rkt")
+
+(provide matrix-scale-row
+ matrix-swap-rows
+ matrix-add-scaled-row
+ matrix-gauss-eliminate
+ matrix-gauss-jordan-eliminate
+ matrix-row-echelon-form
+ matrix-reduced-row-echelon-form
+ matrix-ref)
+
+(define-syntax (inline-matrix-scale-row stx)
+ (syntax-case stx ()
+ [(_ i c)
+ (syntax/loc stx
+ (λ (arr)
+ (let ([arr (array-lazy arr)])
+ (define ds (unsafe-array-shape arr))
+ (define g (unsafe-array-proc arr))
+ (cond
+ [(< i 0)
+ (error 'matrix-scale-row "row index must be non-negative, got ~a" i)]
+ [(not (< i (vector-ref ds 0)))
+ (error 'matrix-scale-row "row index must be smaller than the number of rows, got ~a" i)]
+ [else
+ (unsafe-lazy-array ds (λ: ([js : (Vectorof Index)])
+ (if (= i (vector-ref js 0))
+ (* c (g js))
+ (g js))))]))))]))
+
+(: matrix-scale-row : (Matrix Number) Integer Number -> (Result-Matrix Number))
+(define (matrix-scale-row a i c)
+ ((inline-matrix-scale-row i c) a))
+
+(define-syntax (inline-matrix-swap-rows stx)
+ (syntax-case stx ()
+ [(_ i j)
+ (syntax/loc stx
+ (λ (arr)
+ (let ([arr (array-lazy arr)])
+ (define ds (unsafe-array-shape arr))
+ (define g (unsafe-array-proc arr))
+ (cond
+ [(< i 0)
+ (error 'matrix-swap-rows "row index must be non-negative, got ~a" i)]
+ [(< j 0)
+ (error 'matrix-swap-rows "row index must be non-negative, got ~a" j)]
+ [(not (< i (vector-ref ds 0)))
+ (error 'matrix-swap-rows "row index must be smaller than the number of rows, got ~a" i)]
+ [(not (< j (vector-ref ds 0)))
+ (error 'matrix-swap-rows "row index must be smaller than the number of rows, got ~a" j)]
+ [else
+ (unsafe-lazy-array ds (λ: ([js : (Vectorof Index)])
+ (cond
+ [(= i (vector-ref js 0))
+ (g (vector j (vector-ref js 1)))]
+ [(= j (vector-ref js 0))
+ (g (vector i (vector-ref js 1)))]
+ [else
+ (g js)])))]))))]))
+
+(: matrix-swap-rows : (Matrix Number) Integer Integer -> (Result-Matrix Number))
+(define (matrix-swap-rows a i j)
+ ((inline-matrix-swap-rows i j) a))
+
+(define-syntax (inline-matrix-add-scaled-row stx)
+ (syntax-case stx ()
+ [(_ i c j)
+ (syntax/loc stx
+ (λ (arr)
+ (let ([arr (array-lazy arr)])
+ (define ds (unsafe-array-shape arr))
+ (define g (unsafe-array-proc arr))
+ (cond
+ [(< i 0)
+ (error 'matrix-add-scaled-row "row index must be non-negative, got ~a" i)]
+ [(< j 0)
+ (error 'matrix-add-scaled-row "row index must be non-negative, got ~a" j)]
+ [(not (< i (vector-ref ds 0)))
+ (error 'matrix-add-scaled-row
+ "row index must be smaller than the number of rows, got ~a" i)]
+ [(not (< j (vector-ref ds 0)))
+ (error 'matrix-add-scaled-row
+ "row index must be smaller than the number of rows, got ~a" j)]
+ [else
+ (unsafe-lazy-array ds (λ: ([js : (Vectorof Index)])
+ (if (= i (vector-ref js 0))
+ (+ (g js) (* c (g (vector j (vector-ref js 1)))))
+ (g js))))]))))]))
+
+(: matrix-add-scaled-row : (Matrix Number) Integer Number Integer -> (Result-Matrix Number))
+(define (matrix-add-scaled-row a i c j)
+ ((inline-matrix-add-scaled-row i c j) a))
+
+(: flmatrix-add-scaled-row : (Matrix Flonum) Index Flonum Index -> (Result-Matrix Flonum))
+(define (flmatrix-add-scaled-row a i c j)
+ ((inline-matrix-add-scaled-row i c j) a))
+
+
+;;; GAUSS ELIMINATION / ROW ECHELON FORM
+
+(: matrix-ref : (Matrix Number) Integer Integer -> Number)
+(define (matrix-ref M i j)
+ ((inst array-ref Number) M (list i j)))
+
+(: matrix-gauss-eliminate :
+ (case-> ((Matrix Number) Boolean Boolean -> (Values (Result-Matrix Number) (Listof Integer)))
+ ((Matrix Number) Boolean -> (Values (Result-Matrix Number) (Listof Integer)))
+ ((Matrix Number) -> (Values (Result-Matrix Number) (Listof Integer)))))
+(define (matrix-gauss-eliminate M [unitize-pivot-row? #f] [partial-pivoting? #t])
+ (define dims (matrix-dimensions M))
+ (define m (vector-ref dims 0))
+ (define n (vector-ref dims 1))
+ (: loop : (Integer Integer (Matrix Number) Integer (Listof Integer)
+ -> (Values (Matrix Number) (Listof Integer))))
+ (define (loop i j ; i from 0 to m
+ M
+ k ; count rows without pivot
+ without-pivot)
+ (cond
+ [(or (= i m) (= j n)) (values M without-pivot)]
+ [else
+ ; find row to become pivot
+ (define p
+ (if partial-pivoting?
+ ; find element with maximal absolute value
+ (let: max-loop : (U Boolean Integer)
+ ([l : Integer i] ; i<=l<m
+ [max-current : Real -inf.0]
+ [max-index : Integer i])
+ (cond
+ [(= l m) max-index]
+ [else
+ (let ([v (magnitude (matrix-ref M l j))])
+ (if (> (magnitude (matrix-ref M l j)) max-current)
+ (max-loop (+ l 1) v l)
+ (max-loop (+ l 1) max-current max-index)))]))
+ ; find non-zero element in column
+ (let: first-loop : (U Boolean Integer)
+ ([l : Integer i]) ; i<=l<m
+ (cond
+ [(= l m) #f]
+ [(not (zero? (matrix-ref M l j))) l]
+ [else (first-loop (+ l 1))]))))
+ (cond
+ [(boolean? p) ; p=#f
+ ; no pivot found
+ (loop i (+ j 1) M (+ k 1) (cons j without-pivot))]
+ [(integer? p)
+ ; swap if neccessary
+ (let* ([M (if (= i p) M (matrix-swap-rows M i p))]
+ ; now we now (i,j) is a pivot
+ [M ; maybe scale row
+ (if unitize-pivot-row?
+ (let ([pivot (matrix-ref M i j)])
+ (if (zero? pivot)
+ M
+ (matrix-scale-row M i (/ pivot))))
+ M)])
+ (let ([pivot (matrix-ref M i j)])
+ ; remove elements below pivot
+ (let l-loop ([l (+ i 1)] [M M])
+ (if (= l m)
+ (loop (+ i 1) (+ j 1) M k without-pivot)
+ (let ([x_lj (matrix-ref M l j)])
+ (l-loop (+ l 1)
+ (if (zero? x_lj)
+ M
+ (matrix-add-scaled-row M l (- (/ x_lj pivot)) i))))))
+ ))])]))
+ (let-values ([(M without) (loop 0 0 M 0 '())])
+ (values (array-lazy M) without)))
+
+(: matrix-row-echelon-form :
+ (case-> ((Matrix Number) Boolean -> (Result-Matrix Number))
+ ((Matrix Number) Boolean -> (Result-Matrix Number))
+ ((Matrix Number) -> (Result-Matrix Number))))
+(define (matrix-row-echelon-form M [unitize-pivot-row? #f])
+ (let-values ([(M wp) (matrix-gauss-eliminate M unitize-pivot-row?)])
+ M))
+
+(: matrix-gauss-jordan-eliminate :
+ (case-> ((Matrix Number) Boolean Boolean -> (Values (Result-Matrix Number) (Listof Integer)))
+ ((Matrix Number) Boolean -> (Values (Result-Matrix Number) (Listof Integer)))
+ ((Matrix Number) -> (Values (Result-Matrix Number) (Listof Integer)))))
+(define (matrix-gauss-jordan-eliminate M [unitize-pivot-row? #f] [partial-pivoting? #t])
+ (define dims (matrix-dimensions M))
+ (define m (vector-ref dims 0))
+ (define n (vector-ref dims 1))
+ (: loop : (Integer Integer (Matrix Number) Integer (Listof Integer)
+ -> (Values (Matrix Number) (Listof Integer))))
+ (define (loop i j ; i from 0 to m
+ M
+ k ; count rows without pivot
+ without-pivot)
+ (cond
+ [(or (= i m) (= j n)) (values M without-pivot)]
+ [else
+ ; find row to become pivot
+ (define p
+ (if partial-pivoting?
+ ; find element with maximal absolute value
+ (let: max-loop : (U Boolean Integer)
+ ([l : Integer i] ; i<=l<m
+ [max-current : Real -inf.0]
+ [max-index : Integer i])
+ (cond
+ [(= l m) max-index]
+ [else
+ (let ([v (magnitude (matrix-ref M l j))])
+ (if (> (magnitude (matrix-ref M l j)) max-current)
+ (max-loop (+ l 1) v l)
+ (max-loop (+ l 1) max-current max-index)))]))
+ ; find non-zero element in column
+ (let: first-loop : (U Boolean Integer)
+ ([l : Integer i]) ; i<=l<m
+ (cond
+ [(= l m) #f]
+ [(not (zero? (matrix-ref M l j))) l]
+ [else (first-loop (+ l 1))]))))
+ (cond
+ [(boolean? p) ; p=#f
+ ; no pivot found - this implies the matrix is singular (not invertible)
+ (loop i (+ j 1) M (+ k 1) (cons j without-pivot))]
+ [(integer? p)
+ ; swap if neccessary
+ (let* ([M (if (= i p) M (matrix-swap-rows M i p))]
+ ; now we now (i,j) is a pivot
+ [M ; maybe scale row
+ (if unitize-pivot-row?
+ (let ([pivot (matrix-ref M i j)])
+ (if (zero? pivot)
+ M
+ (matrix-scale-row M i (/ pivot))))
+ M)])
+ (let ([pivot (matrix-ref M i j)])
+ ; remove elements above and below pivot
+ (let l-loop ([l 0] [M M])
+ (cond
+ [(= l m) (loop (+ i 1) (+ j 1) M k without-pivot)]
+ [(= l i) (l-loop (+ l 1) M)]
+ [else
+ (let ([x_lj (matrix-ref M l j)])
+ (l-loop (+ l 1)
+ (if (zero? x_lj)
+ M
+ (matrix-add-scaled-row M l (- (/ x_lj pivot)) i))))]))))])]))
+ (let-values ([(M without) (loop 0 0 M 0 '())])
+ (values (array-lazy M) without)))
+
+(: matrix-reduced-row-echelon-form :
+ (case-> ((Matrix Number) Boolean -> (Result-Matrix Number))
+ ((Matrix Number) Boolean -> (Result-Matrix Number))
+ ((Matrix Number) -> (Result-Matrix Number))))
+(define (matrix-reduced-row-echelon-form M [unitize-pivot-row? #f])
+ (let-values ([(M wp) (matrix-gauss-jordan-eliminate M unitize-pivot-row?)])
+ M))
+
+
Oops, something went wrong.

0 comments on commit d78af91

Please sign in to comment.