|
| 1 | +program metapackage_blas |
| 2 | + implicit none |
| 3 | + |
| 4 | + interface |
| 5 | + subroutine dgesv(n, nrhs, a, lda, ipiv, b, ldb, info) |
| 6 | + integer, intent(in) :: n, nrhs, lda, ldb |
| 7 | + double precision, intent(in out) :: a(lda,*), b(ldb,*) |
| 8 | + integer, intent(out) :: ipiv(*), info |
| 9 | + end subroutine dgesv |
| 10 | + end interface |
| 11 | + |
| 12 | + integer, parameter :: dp = kind(1.0d0) |
| 13 | + real(dp), dimension(:,:), allocatable :: a |
| 14 | + real(dp), dimension(:), allocatable :: b |
| 15 | + integer :: info |
| 16 | + |
| 17 | + allocate(a(3,3), b(3)) |
| 18 | + a = reshape([1.0_dp, 2.0_dp, 3.0_dp, & |
| 19 | + 4.0_dp, 5.0_dp, 6.0_dp, & |
| 20 | + 7.0_dp, 8.0_dp, 9.0_dp], [3,3]) |
| 21 | + b = [1.0_dp, 2.0_dp, 3.0_dp] |
| 22 | + |
| 23 | + call solve_eqsys(a, b, info) |
| 24 | + if (info /= 0) error stop |
| 25 | + |
| 26 | + stop 0 |
| 27 | + |
| 28 | + contains |
| 29 | + |
| 30 | + !> simple wrapper for solvers for real system of linear |
| 31 | + !> equations A * X = B |
| 32 | + subroutine solve_eqsys(a, b, info) |
| 33 | + |
| 34 | + real(dp), dimension(:,:), intent(inout) :: a |
| 35 | + real(dp), dimension(:), intent(inout) :: b |
| 36 | + integer, intent(out) :: info |
| 37 | + integer :: i_alloc |
| 38 | + integer :: n, nrhs, lda, ldb |
| 39 | + integer, dimension(:), allocatable :: ipiv |
| 40 | + ! ------------------------------------------------------------------ |
| 41 | + |
| 42 | + lda = size(a,1) |
| 43 | + n = size(a,2) |
| 44 | + ldb = size(b,1) |
| 45 | + nrhs = 1 |
| 46 | + |
| 47 | + allocate(ipiv(n), stat = i_alloc) |
| 48 | + if (i_alloc /= 0) stop 'solve_eqsys: Allocation for array failed!' |
| 49 | + |
| 50 | + call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info) |
| 51 | + |
| 52 | + info = 0 |
| 53 | + |
| 54 | + deallocate(ipiv, stat = i_alloc) |
| 55 | + if (i_alloc /= 0) stop 'solve_eqsys: Deallocation for array failed!' |
| 56 | + |
| 57 | + end subroutine solve_eqsys |
| 58 | +end program metapackage_blas |
0 commit comments