|
Size: 270
Comment:
|
← Revision 4 as of 2009-12-25 07:10:00 ⇥
Size: 270
Comment: converted to 1.6 markup
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 3: | Line 3: |
| [[Include(Hoxide/Lisp/szcom.lisp)]] | <<Include(Hoxide/Lisp/szcom.lisp)>> |
| Line 6: | Line 6: |
| [[Include(Hoxide/Lisp/LU.lisp)]] | <<Include(Hoxide/Lisp/LU.lisp)>> |
| Line 8: | Line 8: |
| [[Include(Hoxide/Lisp/TLU.lisp)]] | <<Include(Hoxide/Lisp/TLU.lisp)>> |
| Line 10: | Line 10: |
| [[Include(Hoxide/Lisp/LLT.lisp)]] | <<Include(Hoxide/Lisp/LLT.lisp)>> |
| Line 12: | Line 12: |
| [[Include(Hoxide/Lisp/LDLT.lisp)]] | <<Include(Hoxide/Lisp/LDLT.lisp)>> |
| Line 14: | Line 14: |
| [[Include(Hoxide/Lisp/SOR.lisp)]] | <<Include(Hoxide/Lisp/SOR.lisp)>> |
基础工具集 szcom.lisp
;; szcom.lisp
(provide 'szcom)
(defpackage "szcom"
(:use "COMMON-LISP" "EXT")
(:export "RMAT" "PMAT" "PDMAT" "PSMAT" )
(:nicknames "SZCOM")
)
(in-package "szcom")
(defun RMAT (&optional input-stream)
(prog
((m (read input-stream))
(n (read input-stream))
tn
)
(if (not (and (integerp m)
(integerp n)))
(error "RMAT: m or n must be intergers.~%")
)
(setq A (make-array (list m n) :initial-element 0))
(loop for i from 0 below m
do
(loop for j from 0 below n
do
(setq tn (read input-stream))
(cond ((rationalp tn)
(setf (aref A i j) tn))
((null tn) (error "RMAT: unexpected steam END.~%"))
(T (error "RMAT: the member of Array must be numbers~%"))
)
)
)
(return A)
)
)
(defun PMAT (A)
(prog
((m (car (array-dimensions A)))
(n (cadr (array-dimensions A))))
(format t "~%")
(loop for i from 0 to (- m 1)
do
(loop for j from 0 to (- n 1)
do
(format t "~12,6G" (aref A i j))
)
(format t "~%" )
)
)
)
(defun PSMat (A)
(prog
((n (decf (car (array-dimensions A)))))
(loop for i from 0 to n
do
(loop for j from 0 to n
do
(format t "~12,G~^"
(let ((ind (+ (- j i) 1)))
(if (and (<= 0 ind) (< ind 3))
(aref A i ind)
0
)
)
)
)
(format t "~%")
)
)
)
(defun PDMAT (A)
(prog
((n (car (array-dimensions A))))
(format t "~%")
(loop for i from 0 to (- n 1)
do
(loop for j from 0 to (- n 1)
do
(format t "~12,6G~^" (if (= i j) (aref A i) 0)
)
)
(format t "~%" )
)
)
)
LU 分解法解线性方程组
;; 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))
三对角阵的追赶法
;; 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))
对称阵的LLT 分解法
;; LLT.lisp
(load "szcom")
(import '(szcom:RMAT szcom:PMAT))
(defun LLT (A)
(prog
((n (car (array-dimensions A))))
(loop for k from 0 below n
do
(setf (aref A k k)
(sqrt
(- (aref A k k)
(loop for j from 0 to (- k 1)
sum (* (aref A k j)
(aref A j k)
)
)
)
)
)
(loop for i from (+ k 1) below n
do
(setf (aref A i k)
(/ (- (aref A i k)
(loop for j from 0 to (- k 1)
sum (* (aref A i j)
(aref A k j)
)
)
)
(aref A k k)
)
)
(setf (aref A k i ) (aref A i k))
)
)
)
)
(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PMat A)
(pprint A)
(pprint B)
(format t "~%Result of LLT: ~%")
(LLT A)
(PMat A)
(format t "~%In Pairs: ~%" )
(pprint A)
例子
bash-3.00$ cat LLT.dat
4 4
10 7 8 7
7 5 6 5
8 6 10 9
7 5 9 10
4 1
1
2
3
4bash-3.00$ ./LLT.lisp < LLT.dat
A :
10.0000 7.00000 8.00000 7.00000
7.00000 5.00000 6.00000 5.00000
8.00000 6.00000 10.0000 9.00000
7.00000 5.00000 9.00000 10.0000
#2A((10 7 8 7) (7 5 6 5) (8 6 10 9) (7 5 9 10))
#2A((1) (2) (3) (4))
Result of LLT:
3.16228 2.21359 2.52982 2.21359
2.21359 0.316228 1.26491 0.316228
2.52982 1.26491 1.41421 2.12132
2.21359 0.316228 2.12132 0.707108
In Pairs:
#2A((3.1622777 2.2135944 2.529822 2.2135944)
(2.2135944 0.3162276 1.2649105 0.3162276)
(2.529822 1.2649105 1.414214 2.1213198)
(2.2135944 0.3162276 2.1213198 0.70710814))
对称阵的LDLT分解法
;; LDLT.lisp
(load "szcom")
(import '(szcom:RMAT szcom:PMAT szcom:PDMAT))
(defun LDLT (A)
(prog
((n (car (array-dimensions A))))
(setq L (make-array (list n n) :initial-element 0)
D (make-array (list n) :initial-element 0))
(loop for k from 0 below n
do
(setf (aref D k)
(- (aref A k k)
(loop for j from 0 to (- k 1)
sum (* (aref L k j)
(aref L k j)
(aref D j)
)
)
)
)
(setf (aref L k k) 1)
(loop for i from (+ k 1) below n
do
(setf (aref L i k)
(/ (- (aref A i k)
(loop for j from 0 to (- k 1)
sum (* (aref L i j)
(aref L k j)
(aref D j)
)
))
(aref D k)
)
)
)
)
(return (cons L D))
)
)
(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PMat A)
(pprint A)
(pprint B)
(format t "~%Result of LDLT: ~%")
(setq res (LDLT A))
(PMat (car res))
(PDMAT (cdr res))
(format t "~%In Pairs: ~%" )
(pprint res)例子:
bash-3.00$ cat LLT.dat 4 4 10 7 8 7 7 5 6 5 8 6 10 9 7 5 9 10 4 1 1 2 3 4 bash-3.00$ ./LDLT.lisp < LLT.dat A : 10.0000 7.00000 8.00000 7.00000 7.00000 5.00000 6.00000 5.00000 8.00000 6.00000 10.0000 9.00000 7.00000 5.00000 9.00000 10.0000 #2A((10 7 8 7) (7 5 6 5) (8 6 10 9) (7 5 9 10)) #2A((1) (2) (3) (4)) Result of LDLT: 1.00000 0.000000 0.000000 0.000000 0.700000 1.00000 0.000000 0.000000 0.800000 4.00000 1.00000 0.000000 0.700000 1.00000 1.50000 1.00000 10.0000 0.000000 0.000000 0.000000 0.000000 0.100000 0.000000 0.000000 0.000000 0.000000 2.00000 0.000000 0.000000 0.000000 0.000000 0.500000 In Pairs: (#2A((1 0 0 0) (7/10 1 0 0) (4/5 4 1 0) (7/10 1 3/2 1)) . #(10 1/10 2 1/
SOR 迭代法解线性方程组
;; SOR.lisp
(load "szcom")
(import '(szcom:PMAT szcom:RMAT))
(defun SOR (A B X &key ((:weight w) 1.0) ((:ep ep) 1.e-8 ))
(prog ((n (car (array-dimensions B)))
(k 0) (e ep) tx)
(loop until (< e ep)
do
(setf e 0)
(loop for k from 0 below n
do
(setf tx
(/
(- (aref B k 0)
(loop for j from 0 below n
sum
(if (= j k) 0
(* (aref A k j)
(aref X j 0) )
)))
(aref A k k)
)
)
(setf tx (+ (* w tx) (* (- 1 w) (aref X k 0)) ))
(setf e (max e (abs (- (aref X k 0) tx))))
(setf (aref X k 0) tx)
)
(incf k)
;(pprint X)
)
(return 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 SOR result: ~%")
(setq X (make-array (array-dimensions B) :initial-element 0))
(setq k (SOR A B X :weight 1.3 :ep 1e-5))
(format t "~% Loop ~D time(s)" k)
(PMat X)
(pprint X)例子:
bash-3.00$ cat SOR.dat 4 4 -4 1 1 1 1 -4 1 1 1 1 -4 1 1 1 1 -4 4 1 1 1 1 1 bash-3.00$ ./SOR.lisp < SOR.dat A : -4.00000 1.00000 1.00000 1.00000 1.00000 -4.00000 1.00000 1.00000 1.00000 1.00000 -4.00000 1.00000 1.00000 1.00000 1.00000 -4.00000 #2A((-4 1 1 1) (1 -4 1 1) (1 1 -4 1) (1 1 1 -4)) B: 1.00000 1.00000 1.00000 1.00000 #2A((1) (1) (1) (1)) Result of SOR result: Loop 12 time(s) -1.00000 -.999999 -1.00000 -1.00000 #2A((-1.0000014) (-0.99999917) (-1.0) (-1.0000004))