ml.f
1 ! @section Synopsis
2 !
3 ! Machine Learning subprograms for VIF.
4 !
5 ! @author J. Marcel van der Veer
6 !
7 ! @section copyright
8 !
9 ! This file is part of VIF - vintage fortran compiler.
10 ! Copyright 2020-2022 J. Marcel van der Veer <algol68g@xs4all.nl>.
11 !
12 ! @section license
13 !
14 ! This program is free software; you can redistribute it and/or modify it
15 ! under the terms of the gnu general public license as published by the
16 ! free software foundation; either version 3 of the license, or
17 ! (at your option) any later version.
18 !
19 ! This program is distributed in the hope that it will be useful, but
20 ! without any warranty; without even the implied warranty of merchantability
21 ! or fitness for a particular purpose. See the GNU general public license for
22 ! more details. You should have received a copy of the GNU general public
23 ! license along with this program. If not, see <http://www.gnu.org/licenses/>.
24 !
25 subroutine mlcolmean(a, m, n)
26 implicit none
27 integer i, j, m, n
28 real*8 a(m, n), mean
29 real*16 sum
30 do i = 1, n
31 sum = 0.0q0
32 do j = 1, m
33 sum = sum + a(j, i)
34 end do
35 mean = dble (sum / m)
36 do j = 1, m
37 a(j, i) = a(j, i) - mean
38 end do
39 end do
40 end
41 !
42 subroutine mlols(a, m, n, x, y, ainv, u, vt, s, tol)
43 ! Solve 'A' 'x' = 'y' by computing the pseudo-inverse of 'A'.
44 ! 'A' is a 'm'x'n' matrix which is overwritten.
45 ! 'Ainv' will contain the 'n'x'm' pseudo-inverse.
46 ! 'U' will contain U, and 'VT' will contain V^T.
47 ! 'S' will contain reciprocal singular values.
48 ! 'Tol' is the threshold ratio relative to the largest singular value.
49 ! 'Info' is a status indicator. Zero on exit is good.
50 implicit none
51 integer i, info, j, k, lmax, lwork, m, n
52 ! Work dimension should be at least 8 * min(m, n)
53 real*8 a(m, n), x(m), y(m), u(m, m), vt(n, n), s(n), ainv(n, m), tol
54 ! Solve A x = b.
55 call mlpinv (a, m, n, ainv, u, vt, s, tol)
56 call dgemv ('n', n, m, 1.0d0, ainv, n, y, 1, 0.0d0, x, 1)
57 return
58 end
59 !
60 subroutine mlpinv(a, m, n, ainv, u, vt, s, tol)
61 ! Compute pseudo-inverse of 'A'.
62 ! 'A' is a 'm'x'n' matrix which is overwritten.
63 ! 'Ainv' will contain the 'n'x'm' pseudo-inverse.
64 ! 'U' will contain U, and 'VT' will contain V^T.
65 ! 'S' will contain reciprocal singular values.
66 ! 'Tol' is the threshold ratio relative to the largest singular value.
67 ! 'Info' is a status indicator. Zero on exit is good.
68 implicit none
69 real*8 work
70 integer iwork, lwmax
71 parameter (lwmax = 100000)
72 common /work/ work(lwmax), iwork(lwmax)
73 integer i, info, j, k, lmax, lwork, m, n
74 ! Work dimension should be at least 8 * min(m, n)
75 real*8 a(m, n), u(m, m), vt(n, n), s(n), ainv(n, m)
76 real*8 s1, tol
77 real*16 sum
78 ! find optimal workspace.
79 lwork = -1
80 call dgesdd('a', m, n, a, m, s, u, m, vt, n, work, lwork, iwork, info)
81 write (*,*) lwork
82 lwork = min(lwmax, int(work(1)))
83 write (*,*) lwork, work(1)
84 ! compute svd.
85 call dgesdd('a', m, n, a, m, s, u, m, vt, n, work, lwork, iwork, info)
86 if (info > 0) then
87 write (*, *) 'failed to converge'
88 end if
89 ! compute pseudo inverse V * S^+ * U^T
90 ! blas has no routines for multiplying diagonal matrices.
91 s1 = s(1)
92 do i = 1, n
93 if (s(i) > tol * s1) then
94 s(i) = 1.0d0 / s(i)
95 else
96 s(i) = 0.0d0
97 end if
98 end do
99 do i = 1, n
100 do j = 1, m
101 sum = 0
102 do k = 1, m
103 sum = sum + vt(k, i) * s(k) * u (j, k)
104 end do
105 ainv(i, j) = dble (sum)
106 enddo
107 enddo
108 end
109
110 subroutine mltestsignal(a, m, n, y)
111 ! make a dummy test signal.
112 implicit none
113 integer j, k, l, m, n
114 real*8 a(m, n), y(n), x, pi
115 call pi8(pi)
116 do j = 1, m
117 y(j) = j
118 do k = 1, n
119 a(j, k) = j * sin ((k - 0.5) / n * pi)
120 end do
121 end do
122 return
123 end
124
125 ! compute the SVD of rectangular matrix using a divide and conquer algorithm.
126 ! the svd reads
127 !
128 ! a = u * sigma * vt
129 !
130 ! where sigma is an m x n matrix which is zero except for its min(m, n)
131 ! diagonal elements, u is an m x m orthogonal matrix and vt (v transposed)
132 ! is an n x n orthogonal matrix.
133 !
134 ! the diagonal elements of sigma are the real non - negative singular values
135 ! of a- negative, in descending order.
136 ! the first min(m, n) columns of u and v are the left and right singular vectors of a.
137 !
138 ! left singular vectors are stored columnwise, right singular vectors are stored rowwise.
139 ! note that the routine returns vt, not v.
140 !
141 ! iwork dimension should be at least 8 * min(m, n)
142 !
143 subroutine mlsvddd(a, m, n, u, vt, s, iwork)
144 implicit none
145 real*8 work
146 integer iwork, lwmax
147 parameter (lwmax = 100000)
148 common /work/ work(lwmax), iwork(lwmax)
149 integer info, lda, ldu, ldvt, lwork, m, n
150 ! work dimension should be at least 8 * min(m, n)
151 real*8 a(m, n), u(m, m), vt(n, n), s(n)
152 lda = m
153 ldu = m
154 ldvt = n
155 ! find optimal workspace.
156 lwork = -1
157 write (*, *) info
158 call dgesdd('a', m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
159 lwork = min(lwmax, int(work(1)))
160 write (*, *) work(1), lwork
161 ! compute svd.
162 call dgesdd('a', m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
163 if (info > 0) then
164 write (*, *) 'failed to converge'
165 end if
166 end
167
168 subroutine mlwmatrow(a, l, n, m)
169 implicit none
170 integer k, l, m, n
171 real*8 a(n, m)
172 ! Print row n of a having width m.
173 if (m <= 6) then
174 write (*, *) (sngl(a(l, k))i, ' ', k = 1, m)
175 else
176 write (*, *) (sngl(a(l, k)), k = 1, 3) , ' ... ', (sngl(a(l, k)), k = m - 2, m)
177 endif
178 end
179
180 subroutine mlwmatrix(capt, a, n, m)
181 implicit none
182 integer k, m, n
183 ! ML Write MATRIX.
184 real*8 a(n, m)
185 character*(*) capt
186 write (*, *) capt, ': ', n, 'x', m, ' matrix'
187 if (n <= 20) then
188 do k = 1, n
189 call mlwmatrow(a, k, n, m)
190 end do
191 else
192 do k = 1, 10
193 call mlwmatrow(a, k, n, m)
194 end do
195 write (*, *) ' ... '
196 do k = n - 10, n
197 call mlwmatrow(a, k, n, m)
198 end do
199 endif
200 end
201
202 subroutine mlwcolvec(capt, v, n)
203 ! ML Write COLumn VECtor.
204 call mlwmatrix(capt, v, 1, n)
205 end
206
207 subroutine mlwrowvec(capt, v, m)
208 ! ML Write ROW VECtor.
209 call mlwmatrix(capt, v, m, 1)
210 end
© 2002-2024 J.M. van der Veer (jmvdveer@xs4all.nl)
|