Diff for "Hoxide/Lisp/SOR.lisp"

SOR 迭代法解线性方程组

;; SOR.lisp

(load "szcom")
(import '(szcom:PMAT szcom:RMAT))

(defun SOR (A B X &key ((:weight w) 1.0) ((:ep ep) 1.e-8 ))
  (prog ((n (car (array-dimensions B)))
         (k 0) (e ep) tx)
        (loop until (< e ep)
              do
              (setf e 0)
              (loop for k from 0 below n
                    do
                    (setf tx 
                          (/ 
                           (- (aref B k 0)
                              (loop for j from 0 below n
                                    sum
                                    (if (= j k) 0
                                      (* (aref A k j)
                                         (aref X j 0) )
                                      )))
                           (aref A k k)
                           )
                          )
                    (setf tx (+ (* w tx) (* (- 1 w) (aref X k 0)) ))
                    (setf e (max e (abs (- (aref X k 0) tx)))) 
                    (setf (aref X k 0) tx)
                    )
              (incf k)
              ;(pprint X)
              )
        (return k)
        )
  )

(setq A (RMAT))
(setq B (RMAT))
(format t "~% A : ~%")
(PMat A)
(pprint A)
(format t "~% B: ~%")
(PMat B)
(pprint B)


(format t "~%Result of SOR result: ~%")
(setq X (make-array (array-dimensions B) :initial-element 0))
(setq k (SOR A B X :weight 1.3 :ep 1e-5))
(format t "~% Loop ~D time(s)" k)
(PMat X)
(pprint X)

例子:

bash-3.00$ cat SOR.dat
4 4
-4  1  1  1
 1 -4  1  1
 1  1 -4  1
 1  1  1 -4
4 1
1
1
1
1

bash-3.00$ ./SOR.lisp < SOR.dat

 A :

-4.00000     1.00000     1.00000     1.00000
 1.00000    -4.00000     1.00000     1.00000
 1.00000     1.00000    -4.00000     1.00000
 1.00000     1.00000     1.00000    -4.00000

#2A((-4 1 1 1) (1 -4 1 1) (1 1 -4 1) (1 1 1 -4))
 B:

 1.00000
 1.00000
 1.00000
 1.00000

#2A((1) (1) (1) (1))
Result of SOR result:

 Loop 12 time(s)
-1.00000
-.999999
-1.00000
-1.00000

#2A((-1.0000014) (-0.99999917) (-1.0) (-1.0000004))