Topic: Can't link openblas or lapack

I can't link to blas or lapack. In fact, I can't even find lapack or blas.

2 (edited by guitarfdm Yesterday 13:45:01)

Re: Can't link openblas or lapack

I tried a reinstall as admin with virus software turned off, and it is still not here. There is an openblas_config.h in the C:\Program Files (x86)\Simply Fortran 3\mingw-w64\x86_64-w64-mingw32\include but no
libopenblas.a, libopenblas.dll.a, openblas.dll.

Microsoft Windows 11 Professional (x64) Build 26200.8655 (25H2)

It tries to auto-link and fails.

here is the compiler output:
Generating fortran-library.dll
build\burgMemcof.o: In function `linRegMA':
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C:\Users\Admin\Fortran\Burg/./burgMemcof.f90:84: undefined reference to `dgels'
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C:\Users\Admin\Fortran\Burg/./burgMemcof.f90:91: undefined reference to `dgels'
collect2.exe: error: ld returned 1 exit status
Error: Last command making (fortran-library.dll) returned a bad status
Error: Make execution terminatedcompiler output.

Re: Can't link openblas or lapack

Here is the module:

module burgMemcof
    use iso_c_binding, only: C_DOUBLE, C_INT
    implicit none

   
!    interface
!        subroutine dgels(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info) bind(C, name="dgels_")
!            ! Pull the types directly into the interface block itself
!            use, intrinsic :: iso_c_binding, only: C_CHAR, C_INT, C_DOUBLE
!            implicit none
!           
!            character(kind=C_CHAR) :: trans
!            integer(kind=C_INT)    :: m, n, nrhs, lda, ldb, lwork, info
!            real(kind=C_DOUBLE)    :: a(*), b(*), work(*)
!        end subroutine dgels
!    end interface

   
contains

subroutine linRegMA(residual, arCoefNum, maCoefNum, sampleSize, linRegEst) BIND(C, NAME="linRegMA")
    use iso_c_binding, only: C_DOUBLE, C_INT
    implicit none
    !GCC$ ATTRIBUTES DLLEXPORT :: linRegMA

    ! Inputs
    real(C_DOUBLE), dimension(*), intent(in) :: residual
    integer(C_INT), intent(in) :: arCoefNum, maCoefNum, sampleSize

    ! Output
    real(C_DOUBLE), dimension(*), intent(out) :: linRegEst

    ! Local variables
    integer(C_INT) :: m, n
    integer(C_INT) :: i, j
    real(C_DOUBLE), allocatable :: designMatrix(:,:), targetSet(:)
    real(C_DOUBLE), allocatable :: work(:)
    integer(C_INT) :: lwork, info
   
    external dgels
    ! -----------------------------
    ! Compute dimensions
    ! -----------------------------
    m = sampleSize - (arCoefNum + maCoefNum)
    n = maCoefNum + 1      ! +1 for intercept column

    ! -----------------------------
    ! Allocate arrays
    ! -----------------------------
    allocate(designMatrix(m, n))
    allocate(targetSet(m))

    ! -----------------------------
    ! Build targetSet
    ! Mathematica: residual[[arCoefNum + maCoefNum + 1 ;; sampleSize]]
    ! -----------------------------
    do i = 1, m
        targetSet(i) = residual(arCoefNum + maCoefNum + i)
    end do

    ! -----------------------------
    ! Build design matrix
    ! Column 1 = intercept = 1.0
    ! Columns 2..n = lagged residual  values
    ! -----------------------------
    designMatrix(:,1) = 1.0_C_DOUBLE

    do j = 1, maCoefNum
        do i = 1, m
            designMatrix(i, j+1) = residual(arCoefNum + j + (i-1))
        end do
    end do

    ! -----------------------------
    ! Solve Least Squares using DGELS
    ! A = designMatrix (m x n)
    ! B = targetSet (m x 1)
    ! DGELS overwrites B with solution coefficients
    ! -----------------------------
    lwork = -1
    allocate(work(1))

    ! Workspace query
    call dgels('N', m, n, 1, designMatrix, m, targetSet, m, work, lwork, info)

    lwork = int(work(1))
    deallocate(work)
    allocate(work(lwork))

    ! Actual solve
    call dgels('N', m, n, 1, designMatrix, m, targetSet, m, work, lwork, info)

    if (info /= 0) then
        print *, "DGELS failed with INFO=", info
    end if

    ! -----------------------------
    ! Extract coefficients
    ! DGELS stores solution in first n entries of targetSet
    ! -----------------------------
    linRegEst(1:n) = targetSet(1:n)

    ! Cleanup
    deallocate(designMatrix, targetSet, work)
end subroutine linRegMA


subroutine memcof(tseries, n, m, xms, d) BIND(C, NAME="memcof")
    implicit none
    !GCC$ ATTRIBUTES DLLEXPORT :: memcof
   
    ! Argument Types
    integer(C_INT), value :: n, m
    real(C_DOUBLE), intent(out) :: xms
    real(C_DOUBLE), dimension(*), intent(in) :: tseries
    real(C_DOUBLE), dimension(*), intent(out) :: d

    ! Local variables - Must match C_DOUBLE for precision
    integer(C_INT) :: i, j, k
    real(C_DOUBLE) :: p, pneum, denom
    real(C_DOUBLE), allocatable :: wk1(:), wk2(:), wkm(:)

    allocate(wk1(n), wk2(n), wkm(m))

    ! Initial power estimate
    p = 0.0_C_DOUBLE
    do j = 1, n
        p = p + tseries(j)**2
    end do
    xms = p / n

    ! Initialize forward/backward prediction arrays
    wk1(1:n) = tseries(1:n)
    wk2(1:n-1) = tseries(2:n)

    ! Burg recursion
    do k = 1, m
        pneum = 0.0_C_DOUBLE
        denom = 0.0_C_DOUBLE
       
        do j = 1, n - k
            pneum = pneum + wk1(j) * wk2(j)
            denom = denom + wk1(j)**2 + wk2(j)**2
        end do
       
        d(k) = 2.0_C_DOUBLE * pneum / denom
        xms = xms * (1.0_C_DOUBLE - d(k)**2)

        if (k > 1) then
            do i = 1, k-1
                d(i) = wkm(i) - d(k) * wkm(k-i)
            end do
        end if

        if (k == m) exit

        wkm(1:k) = d(1:k)

        do j = 1, n - k - 1
            wk1(j) = wk1(j) - wkm(k) * wk2(j)
            wk2(j) = wk2(j+1) - wkm(k) * wk1(j+1)
        end do
    end do

    deallocate(wk1, wk2, wkm)
end subroutine memcof
   
   
subroutine SymGFilter(tseries, n, lag, width, fout) BIND(C, NAME="SymGFilter")
  implicit none
  !GCC$ ATTRIBUTES DLLEXPORT :: SymGFilter
 
  integer(C_INT), value :: n, lag
  real(C_DOUBLE), intent(in) :: width
  real(C_DOUBLE), dimension(*), intent(in) :: tseries
  real(C_DOUBLE), dimension(*), intent(out) :: fout

  integer(C_INT) :: i, k, fil_len
  real(C_DOUBLE) :: sigma, norm, pi
  real(C_DOUBLE), allocatable :: kernel(:)
  pi = 3.141592653589793_C_DOUBLE

  ! Filter length = 2*lag - 1
  fil_len = 2*lag - 1

  allocate(kernel(fil_len))

  ! Compute sigma
  sigma = real(lag, kind = C_DOUBLE)/width

  ! Build Gaussian kernel centered at index lag
  do k = 1, fil_len
     kernel(k) = exp( - (real(k - lag, kind = C_DOUBLE))**2 / (2.0_C_DOUBLE * sigma**2) ) &
                 / (sqrt(2.0_C_DOUBLE * pi) * sigma)
  end do

  ! Normalize kernel
  norm = 0.0_C_DOUBLE
  do k = 1, fil_len
     norm = norm + kernel(k)
  end do
  kernel = kernel / norm

  ! Initialize output
  do i = 1, n
      fout(i) = 0.0_C_DOUBLE
  end do

  ! Convolution: valid region only
  do i = lag, n - lag + 1
     fout(i + lag - 1) = sum( kernel(:) * tseries(i - lag + 1 : i + lag - 1) )
  end do

  deallocate(kernel)
end subroutine SymGFilter
   
end module burgMemcof

4

Re: Can't link openblas or lapack

In Project Options under the "Linker" category, have you checked the box "Link BLAS and LAPACK Libraries" to enable linking? 

The library is present as libsfopenblas.a in the mingw-w64\x86_64-w64-mingw32\lib directory.  It might be helpful if you select "View Makefile" in the Project menu and share what Simply Fortran generated.  There might be a bug linking to the OpenBLAS library with DLL projects.

Jeff Armstrong
Approximatrix, LLC