|
Size: 1772
Comment:
|
← Revision 4 as of 2009-12-25 07:17:24 ⇥
Size: 1772
Comment: converted to 1.6 markup
|
| No differences found! | |
对称阵的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))