Diff for "Hoxide/Lisp/LDLT.lisp"

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