1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
;; author: ismdeep
;; blog: https://ismdeep.com
;; blog entry: https://ismdeep.com/238
;; doolittle algorithm to decompose a square matrix

(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))

;; Using doolittle algorithm to decomposes a square matrix A for LU = A and LDV = A
(defun doolittle-lu (A)
(let* ((n (car (array-dimensions A))) (tmp 0)
(L (make-array `(,n ,n) :initial-element 0))
(U (make-array `(,n ,n) :initial-element 0))
(D (make-array `(,n ,n) :initial-element 0))
(V (make-array `(,n ,n) :initial-element 0))
(lu-mmul 0)
(ldv-mmul 0))

(loop for j from 0 to (- n 1) do (setf (aref U 0 j) (aref A 0 j)))
(loop for i from 0 to (- n 1) do (setf (aref L i i) 1))
(loop for i from 1 to (- n 1) do
; generate L
(loop for j from 0 to (- i 1) do
(setf tmp (aref A i j))
(loop for k from 0 to (- j 1) do (setf tmp (- tmp (* (aref L i k) (aref U k j)))))
(setf (aref L i j) (/ tmp (aref U j j)))
)
; generate U
(loop for j from i to (- n 1) do
(setf tmp (aref A i j ))
(loop for k from 0 to (- i 1) do (setf tmp (- tmp (* (aref L i k) (aref U k j)))))
(setf (aref U i j) tmp)
)
)

; output A
(format t "---------- A ------------------~%")
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(format t "~4A " (aref A i j)))
(format t "~%"))

; output L
(format t "~%~%---------- L ------------------~%")
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(format t "~4A " (aref L i j)))
(format t "~%"))

; output U
(format t "~%~%---------- U ------------------~%")
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(format t "~4A " (aref U i j)))
(format t "~%"))

; output L * U
(setf lu-mmul (mmul L U))
(format t "~%~%---------- L * U --------------~%")
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(format t "~4A " (aref lu-mmul i j)))
(format t "~%"))

; generate D
(loop for i from 0 to (- n 1) do
(setf (aref D i i) (aref U i i)))

; generate V
(loop for i from 0 to (- n 1) do
(loop for j from i to (- n 1) do
(setf (aref V i j) (/ (aref U i j) (aref U i i)))))

; output D
(format t "~%~%---------- D ------------------~%")
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(format t "~4A " (aref D i j)))
(format t "~%"))

; output V
(format t "~%~%---------- V ------------------~%")
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(format t "~4A " (aref V i j)))
(format t "~%"))


; output L * D * V
(setf ldv-mmul (mmul (mmul L D) V))
(format t "~%~%---------- L * D * V ----------~%")
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(format t "~4A " (aref ldv-mmul i j)))
(format t "~%"))


;; Return L and U
(values L U)))


(let ((g 0))
(setf g (make-array '(4 4) :initial-contents '((2 1 -5 1) (1 -3 0 -6) (0 2 -1 2) (1 4 -7 6))))
(doolittle-lu g))