#!clisp
(load "szcom")
(import '(szcom:PMAT szcom:PSMAT szcom:RMAT))


(defun shooting (A LU)
  (prog
   ((n  (car (array-dimensions A))))
   (loop for k from 0 below n
	 do
	 (setf (aref LU k 0) (aref A k 0))
	 (setf (aref LU k 1)
	       ( - (aref A k 1) 
		   (if (> k 0)
		       (* (aref LU k 0)
			  (aref LU (- k 1) 2)
			  )
		     0)
		   )
	       )
	 (if (< k n)
	     (setf (aref LU k 2)
		   (/ (aref A k 2)
		      (aref LU k 1)
		      )
		   )
	   )
	 )
   )
  )

(defun shoot (LU B X)
  (prog
   ((n  (car (array-dimensions LU))))
   (loop for k from 0 below n
	 do
	 (setf (aref X k 0)
	       (/
		(- (aref B k 0)
		   (if (> k 0)  
		       (* (aref LU k 0) (aref X (- k 1) 0))
		     0 
		     )
		   )
		(aref LU k 1)
		)
	       )
	 )
   (pprint X)
   (loop for k from (- n 1) downto 0
	 do
	 (decf (aref X k 0)
	   (if (< k (- n 1)) 
	       (* (aref LU k 2) (aref X (+ k 1) 0))
	     0)
	   )
	 )
   )
  )


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


(format t "~%Result of shooting: ~%")
(setq LU (make-array (array-dimensions A)))
(shooting A LU)
(PSMat LU)
(format t "~%In Pairs: ~%" )
(pprint LU)

(format t "~%Result of shooting result: ~%")
(setq X (make-array (array-dimensions B)))
(shoot LU B X)
(PMat X)
(pprint X)

;;; (load "lincom")
;;; (import '(lincom:SOLVETD))
;;; (setq X (make-array (array-dimensions B)))
;;; (pprint A)
;;; (SOLVETD A B X)
;;; (pprint X)