Diff for "Hoxide/Lisp/TLU.lisp"

三对角阵的追赶法

;; 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))