1. LU 分解法解线性方程组

#!lisp
;; LU.lisp
(load "szcom")
(import '(szcom:PMAT szcom:RMAT))

(defun LU-S (A LU)
  (prog
   ((n  (car (array-dimensions A)) ))
   (loop for k from 0 to n
         do
         (loop for j from k  below n
               do
               (setf (aref LU k j) 
                     (- (aref A k j) 
                        (loop for q from 0 below k
                              sum (* (aref LU k q) (aref LU q j))
                              )
                        )
                     )
               )
         (loop for i from (+ k 1) below n
               do
               (setf (aref LU i k)
                     (/ 
                      (- (aref A i k) 
                         (loop for q from 0 below k
                               sum (* (aref LU i q) (aref LU q k))
                               )
                         )
                      (aref LU k k)
                      )
                     )
               )
         )
   )
  )

(defun LU-C (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)
                    (loop for i from 0 below k
                          sum (* (aref X i 0) (aref LU k i)))
                    )
                )
          )
    (loop for k from (- n 1) downto 0
          do
          (setf (aref X k 0)
                (/
                 (- (aref X k 0)
                    (loop for j from (+ k 1) below n
                          sum (* (aref X j 0) (aref LU k j))
                          )
                    )
                 (aref LU k k)
                 )
                )
          )
    )
  )

(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PMat A)
(pprint A)
(format t "~% B: ~%")
(PMat B)
(pprint B)

(format t "~%Result of LU: ~%")
;(setq LU (make-array (array-dimensions A)))
(setq LU A)
(LU-S A LU)
(PMat LU)
(format t "~%In Pairs: ~%" )
(pprint LU)

(format t "~%Result of LUS: ~%")
;(setq X (make-array (array-dimensions B)))
(setq X B)
(LU-C LU B X)
(PMat X)
(pprint X)

执行结果:

bash-3.00$ cat ex2.dat
4 4
1  4 -2 3
2  2  0 4
3  0 -1 2
1  2  2 -3

4 1
6
2
1
8

bash-3.00$ ./ex2.lisp < ex2.dat

 A :

 1.00000     4.00000    -2.00000     3.00000
 2.00000     2.00000    0.000000     4.00000
 3.00000    0.000000    -1.00000     2.00000
 1.00000     2.00000     2.00000    -3.00000

#2A((1 4 -2 3) (2 2 0 4) (3 0 -1 2) (1 2 2 -3))
 B:

 6.00000
 2.00000
 1.00000
 8.00000

#2A((6) (2) (1) (8))
Result of LU:

 1.00000     4.00000    -2.00000     3.00000
 2.00000    -6.00000     4.00000    -2.00000
 3.00000     2.00000    -3.00000    -3.00000
 1.00000    0.333333    -.888889    -8.00000

In Pairs:

#2A((1 4 -2 3) (2 -6 4 -2) (3 2 -3 -3) (1 1/3 -8/9 -8))
Result of LUS:

 1.00000
 2.00000
0.000000
-1.00000

#2A((1) (2) (0) (-1))