Changeset 5643 in project
 Timestamp:
 08/23/07 09:59:49 (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

nordsieck/trunk/nordsieck.scm
r5626 r5643 32 32 (defineextension nordsieck) 33 33 34 (declare (export nordsieckform35 nordsieckstep36 nordsieckhscale))34 ;;(declare (export nordsieckform 35 ;; nordsieckstep 36 ;; nordsieckhscale)) 37 37 38 38 … … 51 51 (definematrixmap f64) 52 52 53 (definerecord nordsieckdef order Ainv A B C V d dhist dhistlen p1 p2 m n)53 (definerecord nordsieckdef order Ainv A B C V Econst d dhist dhistlen p1 p2 m n) 54 54 55 55 (definerecordprinter (nordsieckdef x out) … … 149 149 150 150 151 ;; let ((sumsqr (lambda (x result) (+ result (* x x)))))152 ;; (foldec 0 (:range i 10) i sumsqr) )153 154 151 (define (makeAinverse order p1 p2) 155 152 (define N (fx+ 2 (fx+ p1 p2))) 156 153 (define N1 (fx N 1)) 157 154 158 (let ((A1 (ma kef64vector (fx* N N) 0.0))155 (let ((A1 (matrixzeros N N)) 159 156 (bot (fx* N (fx+ 1 p1)))) 160 157 … … 180 177 181 178 A1)) 179 180 181 182 (define (makeV 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 (matrixzeros N N))) 189 (fillmatrix! 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 182 197 183 198 … … 199 214 (let ((D (let ((K (matrixones n n order))) 200 215 (blas:dgemm! order blas:NoTrans blas:NoTrans n n n 1.0 A BAinv 2.0 K))) 201 (d (let ((v (ma kef64vector n 0.0)))216 (d (let ((v (matrixzeros n 1))) 202 217 (blas:dgemv! order blas:NoTrans n n 1.0 A c 0.0 v)))) 203 218 (print "Ainv = " Ainv) … … 208 223 (let ((C (let ((K (matrixones n n order))) 209 224 (matrixmap round (blas:daxpy! (fx* n n) 2.0 K D)))) 210 (V (ma trixzeros n n)))225 (V (makeV order p1 p2))) 211 226 (print "C = " C) 227 (print "V = " V) 212 228 (makenordsieckdef order Ainv A B C V d (list) dhistlen p1 p2 m n)))))))))) 213 229 … … 218 234 (C (nordsieckdefC def)) 219 235 (d (nordsieckdefd def))) 220 (let ((y0 (let ((v (ma kef64vector n 0.0)))236 (let ((y0 (let ((v (matrixzeros n 1))) 221 237 (blas:dgemv! blas:NoTrans n n 1.0 C prev 0.0 v)))) 222 238 (let loop ((i 0) (m m) (ym (list y0))) … … 246 262 (let ((DKV (let ((R (matrixzeros n n))) 247 263 (blas:dgemm! order blas:Trans blas:Trans n n n 1.0 D KV 0.0 R)))) 248 (let ((d1 (let ((v (ma kef64vector n)))264 (let ((d1 (let ((v (matrixzeros n 1))) 249 265 (blas:dgemv! order blas:NoTrans n n 1.0 DKV gamma 0.0 v)))) 250 (let ((Ainv (nordsieckdefAinv def)) 251 (A (nordsieckdefA def)) 252 (B (nordsieckdefB def)) 253 (C (nordsieckdefC def)) 254 (dhist (nordsieckdefdhist def)) 255 (dhistlen (nordsieckdefdhistlen def)) 256 (p1 (nordsieckdefp1 def)) 257 (p2 (nordsieckdefp2 def)) 258 (m (nordsieckdefm def))) 259 (let ((dhist (if (fx>= (length dhist) dhistlen) 260 (take dhist (fx dhistlen 1)) dhist))) 261 (makenordsieckdef order Ainv A B C V d1 (cons d dhist) dhistlen p1 p2 m n))))))))))) 266 (match def 267 (($ nordsieckdef 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 (makenordsieckdef order Ainv A B C V d1 (cons d dhist) dhistlen p1 p2 m n)))))))))))) 262 271
Note: See TracChangeset
for help on using the changeset viewer.