29 lines
891 B
Racket
29 lines
891 B
Racket
#lang racket
|
|
(require math/matrix math/array)
|
|
(define-values (T I col size)
|
|
(values ; short names
|
|
matrix-transpose identity-matrix matrix-col matrix-num-rows))
|
|
|
|
(define (scale c A) (matrix-scale A c))
|
|
(define (unit n i) (build-matrix n 1 (λ (j _) (if (= j i) 1 0))))
|
|
|
|
(define (H u)
|
|
(matrix- (I (size u))
|
|
(scale (/ 2 (matrix-dot u u))
|
|
(matrix* u (T u)))))
|
|
|
|
(define (normal a)
|
|
(define a0 (matrix-ref a 0 0))
|
|
(matrix- a (scale (* (sgn a0) (matrix-2norm a))
|
|
(unit (size a) 0))))
|
|
|
|
(define (QR A)
|
|
(define n (size A))
|
|
(for/fold ([Q (I n)] [R A]) ([i (- n 1)])
|
|
(define Hi (H (normal (submatrix R (:: i n) (:: i (+ i 1))))))
|
|
(define Hi* (if (= i 0) Hi (block-diagonal-matrix (list (I i) Hi))))
|
|
(values (matrix* Q Hi*) (matrix* Hi* R))))
|
|
|
|
(QR (matrix [[12 -51 4]
|
|
[ 6 167 -68]
|
|
[-4 24 -41]]))
|