1. 对称阵的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))

Hoxide/Lisp/LLT.lisp (last edited 2009-12-25 07:17:24 by localhost)