三对角阵的追赶法
;; TLU.lisp
(load "szcom")
(import '(szcom:PMAT szcom:PSMAT szcom:RMAT))
(defun shooting (A LU)
(prog
((n (car (array-dimensions A))))
(loop for k from 0 below n
do
(setf (aref LU k 0) (aref A k 0))
(setf (aref LU k 1)
( - (aref A k 1)
(if (> k 0)
(* (aref LU k 0)
(aref LU (- k 1) 2)
)
0)
)
)
(if (< k n)
(setf (aref LU k 2)
(/ (aref A k 2)
(aref LU k 1)
)
)
)
)
)
)
(defun shoot (LU B X)
(prog
((n (car (array-dimensions LU))))
(loop for k from 0 below n
do
(setf (aref X k 0)
(/
(- (aref B k 0)
(if (> k 0)
(* (aref LU k 0) (aref X (- k 1) 0))
0
)
)
(aref LU k 1)
)
)
)
(pprint X)
(loop for k from (- n 1) downto 0
do
(decf (aref X k 0)
(if (< k (- n 1))
(* (aref LU k 2) (aref X (+ k 1) 0))
0)
)
)
)
)
(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PSMat A)
(pprint A)
(format t "~% B: ~%")
(PMat B)
(pprint B)
(format t "~%Result of shooting: ~%")
(setq LU (make-array (array-dimensions A)))
(shooting A LU)
(PSMat LU)
(format t "~%In Pairs: ~%" )
(pprint LU)
(format t "~%Result of shooting result: ~%")
(setq X (make-array (array-dimensions B)))
(shoot LU B X)
(PMat X)
(pprint X)例子:
bash-3.00$ cat ex3.dat
5 3
0 13 12
11 10 9
8 7 6
5 4 3
2 1 0
5 1
3
0
-2
6
8
bash-3.00$ ./ex3.lisp < ex3.dat
A :
13. 12. 0.0 0.0 0.0
11. 10. 9. 0.0 0.0
0.0 8. 7. 6. 0.0
0.0 0.0 5. 4. 3.
0.0 0.0 0.0 2. 1.
#2A((0 13 12) (11 10 9) (8 7 6) (5 4 3) (2 1 0))
B:
3.00000
0.000000
-2.00000
6.00000
8.00000
#2A((3) (0) (-2) (6) (8))
Result of shooting:
13. .9230769 0.0 0.0 0.0
11. -.15384616 -58.5 0.0 0.0
0.0 8. 475. 1.263157900E-2 0.0
0.0 0.0 5. 3.9368422 .7620321
0.0 0.0 0.0 2. -.5240642
In Pairs:
#2A((0 13 12/13)
(11 -2/13 -117/2)
(8 475 6/475)
(5 374/95 285/374)
(2 -98/187 0))
Result of shooting result:
#2A((3/13) (33/2) (-134/475) (32/17) (-396/49))
5.71837
-5.94490
-.383673
8.04082
-8.08163
#2A((1401/245) (-2913/490) (-94/245) (394/49) (-396/49))