reference common lisp lup-decomposition implementation, Common Lisp

reference common lisp lup-decomposition implementation

(defun mmul (A B)
  (let* ((m (car (array-dimensions A)))
         (n (cadr (array-dimensions A)))
         (l (cadr (array-dimensions B)))
         (C (make-array `(,m ,l) :initial-element 0)))
    (loop for i from 0 to (- m 1) do
              (loop for k from 0 to (- l 1) do
                    (setf (aref C i k)
                          (loop for j from 0 to (- n 1)
                                sum (* (aref A i j)
                                       (aref B j k))))))
    C))

;; Creates a nxn identity matrix.
(defun eye (n)
  (let ((I (make-array `(,n ,n) :initial-element 0)))
    (loop for j from 0 to (- n 1) do
          (setf (aref I j j) 1))
    I))
 
;; Swap two rows l and k of a mxn matrix A, which is a 2D array.
(defun swap-rows (A l k)
  (let* ((n (cadr (array-dimensions A)))
         (row (make-array n :initial-element 0)))
    (loop for j from 0 to (- n 1) do
          (setf (aref row j) (aref A l j))
          (setf (aref A l j) (aref A k j))
          (setf (aref A k j) (aref row j)))))
 
;; Creates the pivoting matrix for A.
(defun pivotize (A)
  (let* ((n (car (array-dimensions A)))
         (P (eye n)))
    (loop for j from 0 to (- n 1) do
          (let ((max (aref A j j))
                (row j))
            (loop for i from j to (- n 1) do
                  (if (> (aref A i j) max)
                      (setq max (aref A i j)
                            row i)))
            (if (not (= j row))
                (swap-rows P j row))))
 
  ;; Return P.
  P))
 
;; Decomposes a square matrix A by PA=LU and returns L, U and P.
(defun lu (A)
  (let* ((n (car (array-dimensions A)))
         (L (make-array `(,n ,n) :initial-element 0))
         (U (make-array `(,n ,n) :initial-element 0))
         (P (pivotize A))
         (A (mmul P A)))
 
    (loop for j from 0 to (- n 1) do
          (setf (aref L j j) 1)
          (loop for i from 0 to j do
                (setf (aref U i j)
                      (- (aref A i j)
                         (loop for k from 0 to (- i 1)
                               sum (* (aref U k j)
                                      (aref L i k))))))
          (loop for i from j to (- n 1) do
                (setf (aref L i j)
                      (/ (- (aref A i j)
                            (loop for k from 0 to (- j 1)
                                  sum (* (aref U k j)
                                         (aref L i k))))
                         (aref U j j)))))
 
  ;; Return L, U and P.
  (values L U P)))

(setf h (make-array '(4 4) :initial-contents '((11 9 24 2)(1 5 2 6)(3 17 18 1)(2 5 7 1))))

(multiple-value-bind (L U P) (lu h)
    (format t "~A~%" L)
    (format t "~A~%" U)
    (format t "~A~%" P))


run \| edit \| history \| help	0

Please log in to post a comment.