I solved it by not using `for/sum' and writing this ridiculous function
for the recursive base case and the initial values in `x':
(: zero-of (case-> (Real -> Real)
(Number -> Number)))
(define (zero-of x) 0)
Fortunately, it should get inlined. I also renamed `U' to `V', because
it was overshadowing the type name.
(define matrix-solve-upper
; solve the equation Ux=b
; using back substitution
(case-lambda
[(V b)
(define (default-fail)
(raise-argument-error
'matrix-solve-upper
"The upper triangular matrix is not invertible." 0 V))
(matrix-solve-upper V b default-fail)]
[(V b fail)
(define zero (zero-of (matrix-ref V 0 0)))
(define m (matrix-num-rows V))
(define x (make-vector m zero))
(for ([i (in-range (- m 1) -1 -1)])
(define bi (matrix-ref b i 0))
(define Vii (matrix-ref V i i))
(when (zero? Vii) (fail))
(let loop ([j (+ i 1)] [s zero])
(cond [(j . < . m)
(loop (+ j 1) (+ s (* (matrix-ref V i j)
(vector-ref x j))))]
[else
(vector-set! x i (/ (- bi s) Vii))])))
(vector->matrix m 1 x)]))
Neil ⊥
On 01/03/2013 10:47 AM, Jens Axel Søgaard wrote:
Ignore the previous example. Here is the example again, now
with correct usage of case-lambda. The for/sum problem remains.
/Jens Axel
#lang typed/racket
(require math)
(: matrix-solve-upper
(All (A) (case->
((Matrix Real) (Matrix Real) -> (Matrix Real))
((Matrix Real) (Matrix Real) (-> A) -> (U A (Matrix Real)))
((Matrix Number) (Matrix Number) -> (Matrix Number))
((Matrix Number) (Matrix Number) (-> A) -> (U A (Matrix
Number))))))
(define matrix-solve-upper
; solve the equation Ux=b
; using back substitution
(case-lambda
[(U b)
(define (default-fail)
(raise-argument-error
'matrix-solve-upper
"The upper triangular matrix is not invertible." 0 U))
(matrix-solve-upper U b default-fail)]
[(U b fail)
(define m (matrix-num-rows U))
(define x (make-vector m 0))
(for ([i (in-range (- m 1) -1 -1)])
(define bi (matrix-ref b i 0))
(define Uii (matrix-ref U i i))
(when (zero? Uii) (fail))
(define x.Ui (for/sum ([j (in-range (+ i 1) m)])
(* (matrix-ref U i j) (vector-ref x j))))
(vector-set! x i (/ (- bi x.Ui) Uii)))
(vector->matrix m 1 x)]))
(define U (matrix [[4 -1 2 3]
[0 -2 7 -4]
[0 0 6 5]
[0 0 0 3]]))
(define b (col-matrix [20 -7 4 6]))
b ; expected
(define x (solve-upper-triangular U b))
(matrix* U x)
_________________________
Racket Developers list:
http://lists.racket-lang.org/dev
_________________________
Racket Developers list:
http://lists.racket-lang.org/dev
