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)