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