Changeset 5643 in project


Ignore:
Timestamp:
08/23/07 09:59:49 (12 years ago)
Author:
Ivan Raikov
Message:

Added procedure make-V

File:
1 edited

Legend:

Unmodified
Added
Removed
  • nordsieck/trunk/nordsieck.scm

    r5626 r5643  
    3232(define-extension nordsieck)
    3333
    34 (declare (export nordsieck-form
    35                  nordsieck-step
    36                  nordsieck-hscale))
     34;;(declare (export nordsieck-form
     35;;               nordsieck-step
     36;;               nordsieck-hscale))
    3737
    3838
     
    5151(define-matrix-map f64)
    5252
    53 (define-record nordsieck-def order Ainv A B C V d dhist dhistlen p1 p2 m n)
     53(define-record nordsieck-def order Ainv A B C V Econst d dhist dhistlen p1 p2 m n)
    5454
    5555(define-record-printer (nordsieck-def x out)
     
    149149
    150150
    151 ;; let ((sum-sqr (lambda (x result) (+ result (* x x)))))
    152 ;;         (fold-ec 0 (:range i 10) i sum-sqr) )
    153 
    154151(define (make-Ainverse order p1 p2)
    155152  (define N (fx+ 2 (fx+ p1 p2)))
    156153  (define N-1 (fx- N 1))
    157154
    158   (let ((A-1  (make-f64vector (fx* N N)  0.0))
     155  (let ((A-1  (matrix-zeros N N))
    159156        (bot  (fx* N (fx+ 1 p1))))
    160157
     
    180177
    181178    A-1))
     179
     180
     181
     182(define (make-V order p1 p2)
     183  (define N (fx+ 2 (fx+ p1 p2)))
     184 
     185  (define (f1 i n)
     186    (if (> i 0) (* n (f1 (- i 1) (- n 1))) 1))
     187
     188  (let ((V (matrix-zeros N N)))
     189    (fill-matrix! order N N V
     190                  (lambda (i j ax)
     191                    (values (cond ((fx= i 0)  (+ 1 j))
     192                                  ((fx= i j)  1.0)
     193                                  ((fx> i j)  0.0)
     194                                  (else (* (/ 1 (f1 i (+ 1 i))) (f1 (+ 1 i) (+ 1 j))))) #f))
     195                  #f)))
     196                                 
    182197
    183198
     
    199214              (let ((D  (let ((K (matrix-ones n n order)))
    200215                          (blas:dgemm! order blas:NoTrans blas:NoTrans n n n 1.0 A BAinv 2.0 K)))
    201                     (d     (let ((v (make-f64vector n 0.0)))
     216                    (d     (let ((v (matrix-zeros n 1)))
    202217                             (blas:dgemv! order blas:NoTrans n n 1.0 A c 0.0 v))))
    203218                  (print "Ainv = " Ainv)
     
    208223                  (let ((C   (let ((K (matrix-ones n n order)))
    209224                               (matrix-map round (blas:daxpy! (fx* n n) -2.0 K D))))
    210                         (V   (matrix-zeros n n)))
     225                        (V   (make-V order p1 p2)))
    211226                    (print "C = " C)
     227                    (print "V = " V)
    212228                    (make-nordsieck-def order Ainv A B C V d (list) dhistlen p1 p2 m n))))))))))
    213229 
     
    218234        (C     (nordsieck-def-C def))
    219235        (d     (nordsieck-def-d def)))
    220   (let ((y0  (let ((v (make-f64vector n 0.0)))
     236  (let ((y0  (let ((v (matrix-zeros n 1)))
    221237               (blas:dgemv! blas:NoTrans n n 1.0 C prev 0.0 v))))
    222238    (let loop ((i 0) (m m) (ym (list y0)))
     
    246262            (let ((DKV  (let ((R (matrix-zeros n n)))
    247263                          (blas:dgemm! order blas:Trans blas:Trans n n n 1.0 D KV 0.0 R))))
    248               (let ((d1  (let ((v (make-f64vector n)))
     264              (let ((d1  (let ((v (matrix-zeros n 1)))
    249265                           (blas:dgemv! order blas:NoTrans n n 1.0 DKV gamma 0.0 v))))
    250                 (let ((Ainv (nordsieck-def-Ainv def))
    251                       (A (nordsieck-def-A def))
    252                       (B (nordsieck-def-B def))
    253                       (C (nordsieck-def-C def))
    254                       (dhist (nordsieck-def-dhist def))
    255                       (dhistlen (nordsieck-def-dhistlen def))
    256                       (p1 (nordsieck-def-p1 def))
    257                       (p2 (nordsieck-def-p2 def))
    258                       (m  (nordsieck-def-m def)))
    259                   (let ((dhist  (if (fx>= (length dhist) dhistlen)
    260                                     (take dhist (fx- dhistlen 1)) dhist)))
    261                     (make-nordsieck-def order Ainv A B C V d1 (cons d dhist) dhistlen p1 p2 m n)))))))))))
     266                (match def
     267                       (($ nordsieck-def order Ainv A B C V d dhist dhistlen p1 p2 m n)
     268                        (let ((dhist  (if (fx>= (length dhist) dhistlen)
     269                                          (take dhist (fx- dhistlen 1)) dhist)))
     270                          (make-nordsieck-def order Ainv A B C V d1 (cons d dhist) dhistlen p1 p2 m n))))))))))))
    262271       
Note: See TracChangeset for help on using the changeset viewer.