chemometrics.f
1 ! @section Synopsis
2 !
3 ! Machine Learning subprograms for VIF. Needs LAPACK.
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-2025 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 ml colmean(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 ml ols(a, m, n, x, y, a inv, 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), a inv(n, m), tol
54 ! Solve A x = b.
55 CALL ml pinv (a, m, n, a inv, u, vt, s, tol)
56 CALL dgemv ('n', n, m, 1.0d0, a inv, n, y, 1, 0.0d0, x, 1)
57 RETURN
58 END
59 !
60 SUBROUTINE ml pinv(a, m, n, a inv, 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), a inv(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 lwork = min(lwmax, int(work(1)))
82 ! compute svd.
83 CALL dgesdd('a', m, n, a, m, s, u, m, vt, n, work, lwork, iwork, info)
84 IF (info > 0) THEN
85 write (*, *) 'failed to converge'
86 END IF
87 ! compute pseudo inverse V * S^+ * U^T
88 ! blas has no routines for multiplying diagonal matrices.
89 s1 = s(1)
90 DO i = 1, n
91 IF (s(i) > tol * s1) THEN
92 s(i) = 1.0d0 / s(i)
93 ELSE
94 s(i) = 0.0d0
95 END IF
96 END DO
97 DO i = 1, n
98 DO j = 1, m
99 sum = 0
100 DO k = 1, m
101 sum = sum + vt(k, i) * s(k) * u (j, k)
102 END DO
103 a inv(i, j) = dble (sum)
104 END DO
105 END DO
106 END
107
108 SUBROUTINE ml testsignal(a, m, n, y)
109 ! make a dummy test signal.
110 IMPLICIT none
111 INTEGER j, k, l, m, n
112 REAL*8 a(m, n), y(n), x, pi
113 CALL pi8(pi)
114 DO j = 1, m
115 y(j) = j
116 DO k = 1, n
117 a(j, k) = j * sin ((k - 0.5) / n * pi)
118 END DO
119 END DO
120 RETURN
121 END
122
123 SUBROUTINE ml svddd(a, m, n, u, vt, s, iwork)
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 IMPLICIT none
144 REAL*8 work
145 INTEGER iwork, lwmax
146 parameter (lwmax = 100000)
147 common /work/ work(lwmax), iwork(lwmax)
148 INTEGER info, lda, ldu, ldvt, lwork, m, n
149 ! work dimension should be at least 8 * min(m, n)
150 REAL*8 a(m, n), u(m, m), vt(n, n), s(n)
151 lda = m
152 ldu = m
153 ldvt = n
154 ! find optimal workspace.
155 lwork = -1
156 CALL dgesdd('a', m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
157 lwork = min(lwmax, int(work(1)))
158 ! compute svd.
159 CALL dgesdd('a', m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
160 IF (info > 0) THEN
161 write (*, *) 'failed to converge'
162 END IF
163 END
164
165 SUBROUTINE ml wmatrow(a, l, n, m)
166 IMPLICIT none
167 INTEGER k, l, m, n
168 REAL*8 a(n, m)
169 ! Print row n of a having width m.
170 IF (m <= 6) THEN
171 write (*, *) (sngl(a(l, k))i, ' ', k = 1, m)
172 ELSE
173 write (*, *) (sngl(a(l, k)), k = 1, 3) , ' ... ', (sngl(a(l, k)), k = m - 2, m)
174 END IF
175 END
176
177 SUBROUTINE ml wmatrix(caption, a, n, m)
178 IMPLICIT none
179 INTEGER k, m, n
180 ! ML Write MATRIX.
181 REAL*8 a(n, m)
182 character*(*) caption
183 write (*, *) caption, ': ', n, 'x', m, ' matrix'
184 IF (n <= 20) THEN
185 DO k = 1, n
186 CALL ml wmatrow(a, k, n, m)
187 END DO
188 ELSE
189 DO k = 1, 10
190 CALL ml wmatrow(a, k, n, m)
191 END DO
192 write (*, *) ' ... '
193 DO k = n - 10, n
194 CALL ml wmatrow(a, k, n, m)
195 END DO
196 END IF
197 END
198
199 SUBROUTINE ml wcolvec(caption, v, n)
200 ! ML Write COLumn VECtor.
201 CALL ml wmatrix(caption, v, 1, n)
202 END
203
204 SUBROUTINE mlwrowvec(caption, v, m)
205 ! ML Write ROW VECtor.
206 CALL ml wmatrix(caption, v, m, 1)
207 END
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
