对称阵的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/