1. SOR 迭代法解线性方程组
#!lisp ;; 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))