Algol 68 Genie

ml.f

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

This website is archived by the National Library of the Netherlands.