### 1. 对称阵的LDLT分解法

```;; LDLT.lisp

(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/```

Hoxide/Lisp/LDLT.lisp (last edited 2009-12-25 07:16:25 by localhost)

• Page.execute = 0.006s
• getACL = 0.026s
• init = 0.001s