1. 对称阵的LLT 分解法
#!lisp ;; 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))