# DGBMV blas issue with 2d array **URL:** \ **Category:** Help\ **Created:** [May 22, 2023, 4:15pm UTC](https://fortran-lang.discourse.group/t/dgbmv-blas-issue-with-2d-array/5842 "2023-05-22T16:15:49Z")\ **Posts on this page:** 8\ **Page:** 1
**Author:** ![sabbraxcaddabra](https://avatars.discourse-cdn.com/v4/letter/s/9fc348/32.png) [@sabbraxcaddabra](https://fortran-lang.discourse.group/u/sabbraxcaddabra)\ **Post date:** [May 22, 2023, 4:15pm UTC](https://fortran-lang.discourse.group/t/dgbmv-blas-issue-with-2d-array/5842/1 "2023-05-22T16:15:49Z")
I learn how to use blas and lapack libs in fortran for operationg with banded matricies. The code below i use to solve 2d heat equation by Peaceman‐Rachford schema, but i cut all not needed code. What i need? I need to multiply banded matrix coefficients on vector of function to get second deriv of function. I create two arrays: test\_array\_col(2, N+1) and test\_array\_row(N+1, 2) and matrix band\_matrix\_y(3, N+1). I init array by 1.0. After calling `dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, test_array_col(1, :), 1, 1.0_dp, test_array_col(1, :), 1)` or `dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, test_array_row(:, 1), 1, 1.0_dp, test_array_row(:, 1), 1)` i expect get array where all elements are equal to 1.0, but i get this results only in first call. The issue is that i need to store 2d array of my function and call dgbmv on each time step first for each row (work) and for each column (dont work). What can i do without transposing my 2d array of function on each time step. Of course i can implement this functionality myself, but i want to explain what happen below The whole code here ```auto program dgbmv_issue implicit none integer, parameter:: N=10, M=10, dp=(kind(0.d0)) double precision, dimension(:, :):: band_matrix_y(3, N+1), test_array_col(2, N+1), test_array_row(N+1, 2) integer:: i double precision:: DX, DY, L_X, L_Y, DT external:: dgbmv test_array_col = 1.0_dp test_array_row = 1.0_dp L_X = 1.0_dp L_Y = 1.0_dp DX = L_X / M DY = L_Y / N DT = 0.01_dp band_matrix_y = 0.0_dp band_matrix_y(1, 3:) = 0.5_dp * DT / DY / DY band_matrix_y(2, 2:N) = -DT / DY / DY band_matrix_y(3, :N-1) = 0.5_dp * DT / DY / DY open(1, file='out.txt') do i = 1, 3 write(1, *) band_matrix_y(i, :) end do call dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, test_array_col(1, :), 1, 1.0_dp, test_array_col(1, :), 1) write(1, *) 'X' write(1, *) test_array_col(1, :) call dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, test_array_row(:, 1), 1, 1.0_dp, test_array_row(:, 1), 1) write(1, *) 'X' write(1, *) test_array_row(:, 1) close(1) end program dgbmv_issue ``` And here’s meson build file ```auto project('dgbmv_issue', 'fortran') lapack = dependency('lapack') blas = dependency('blas') executable('dgbmv_issue', 'dgbmv_issue.f90', dependencies: [lapack, blas]) ``` Can anybody help me? Please🙏 ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [May 22, 2023, 7:44pm UTC](https://fortran-lang.discourse.group/t/dgbmv-blas-issue-with-2d-array/5842/2 "2023-05-22T19:44:45Z")
One possible problem is that the input argument `x(*)` is aliased to the inout argument `y(*)`. This applies to both calls, but apparently one call appears to work alright, while the other doesn’t. The reason it works alright is probably because of copy-in/copy-out argument association, so it is just by luck that it happens to work. The compiler is not required to notify you of this error, but there might be some combination of compiler options where it would, at least if `dgbmv()` had an explicit interface. You are using implicit interfaces, so these kinds of errors are always easy to make and difficult to find. You are using a combination of array slice notation actual argument with an implicit interface. I think you are doing that correctly, but it is confusing nonetheless. It might be simpler for you, and for anyone reading your code, if you were to avoid the array slice actual argument. I don’t know how to suggest that this should be changed because it also depends on how you choose to avoid the aliasing of the actual arguments. ---
**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\ **Post date:** [May 22, 2023, 8:27pm UTC](https://fortran-lang.discourse.group/t/dgbmv-blas-issue-with-2d-array/5842/3 "2023-05-22T20:27:05Z")
The first thing to fix is the aliasing, for both calls. ```auto real, allocatable :: tmp(:) tmp = test_array_col(1, :) ) call dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, tmp, 1, 1.0_dp, test_array_col(1, :), 1) tmp = test_array_row(:,1) call dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, tmp, 1, 1.0_dp, test_array_row(:, 1), 1) ``` Then replace the array sections by using the BLAS classical INC\* arguments. This can work with array sections, but it should be more efficient without. ```auto real, allocatable :: tmp(:) tmp = test_array_col(1, :) ) call dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, tmp, 1, 1.0_dp, test_array_col, 2) tmp = test_array_row(:,1) call dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, tmp, 1, 1.0_dp, test_array_row, 1) ``` ---
**Author:** ![sabbraxcaddabra](https://avatars.discourse-cdn.com/v4/letter/s/9fc348/32.png) [@sabbraxcaddabra](https://fortran-lang.discourse.group/u/sabbraxcaddabra)\ **Post date:** [May 22, 2023, 9:22pm UTC](https://fortran-lang.discourse.group/t/dgbmv-blas-issue-with-2d-array/5842/5 "2023-05-22T21:22:32Z")
Thank you for your reply! The fact is that first dgbmv example call (for row) is stable, in full code for heat equation I call it in cycle and it’s work fine on all iterations, although one array passing as input and output argument. I don’t really know what to do ---
**Author:** ![sabbraxcaddabra](https://avatars.discourse-cdn.com/v4/letter/s/9fc348/32.png) [@sabbraxcaddabra](https://fortran-lang.discourse.group/u/sabbraxcaddabra)\ **Post date:** [May 22, 2023, 9:23pm UTC](https://fortran-lang.discourse.group/t/dgbmv-blas-issue-with-2d-array/5842/6 "2023-05-22T21:23:43Z")
Thank you for your reply, but I need to compute expression y = alpha_A_x + beta\*y where x=y, so I would rather not call vector+vector subroutine after dgbmv. In other words I can either copy column/row to tmp vector or call vector+vector subroutine after dgbmv? Is there no way to do what I try to do without extra memory allocation or subroutine calls? ---
**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\ **Post date:** [May 22, 2023, 9:51pm UTC](https://fortran-lang.discourse.group/t/dgbmv-blas-issue-with-2d-array/5842/7 "2023-05-22T21:51:50Z")
> [@sabbraxcaddabra](#): > > Thank you for your reply! The fact is that first dgbmv example call (for row) is stable As stated by @RonShepard this is by chance. Argument aliasing is illegal wrt the Fortran standard, as soon as one of the arguments is an output one. The behavior is undefined : it can work, or give wrong results, or crash, etc, depending on the compiler, on the compilation options, etc… > [@sabbraxcaddabra](#): > > Is there no way to do what I try to do without extra memory allocation or subroutine calls? No, you can’t. You need a temporary array at some point, either for the input or for the output. I’ve corrected my previous answer according to what you want. ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [May 22, 2023, 10:49pm UTC](https://fortran-lang.discourse.group/t/dgbmv-blas-issue-with-2d-array/5842/8 "2023-05-22T22:49:58Z")
> [@PierU](#): > > > [@sabbraxcaddabra](#): > > > > Is there no way to do what I try to do without extra memory allocation or subroutine calls? > > No, you can’t. You need a temporary array at some point, either for the input or for the output. I’ve corrected my previous answer according to what you want. I also think this is correct. The arguments cannot be aliased, so you need to make a copy somehow. Another way do to that is to turn one of them into an expression in the call statement. The compiler then most likely creates the copy on the stack rather than the heap (where the allocatable array `tmp(:)` probably lives), so there are both advantages and disadvantages to this approach. Modifying the above code would give something like: > [@PierU](#): > > ```auto > call dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, (test_array_col(1, :)), 1, 1.0_dp, test_array_col, 2) > > call dgbmv('N', N+1, N+1, 1, 1, 1.0_dp, band_matrix_y, 3, (test_array_row(:,1)), 1, 1.0_dp, test_array_row, 1) > > ``` It is the extra parentheses around the argument that does the trick. Note that you are still making the copy, you are just letting the compiler do the heavy lifting. If that argument is needed several times, then it is probably best to allocate an assign it once in `tmp(:)` and then reuse it. If it is a one-time call, then maybe the expression is best. If the array slice is long, many megabytes, then the compiler might exhaust stack space, and then the `tmp(:)` array would have some advantages. I think the increment of 1 is correct in the first call above, but that is one of the things that looks confusing when mixing modern fortran conventions with the f77-style blas/lapack conventions. I’m more confident that the increment of 1 is correct for the second call. ---
**Author:** ![PierU](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pieru/32/1848_2.png) [@PierU](https://fortran-lang.discourse.group/u/PierU)\ **Post date:** [May 23, 2023, 6:00am UTC](https://fortran-lang.discourse.group/t/dgbmv-blas-issue-with-2d-array/5842/9 "2023-05-23T06:00:55Z")
> [@RonShepard](#): > > Another way do to that is to turn one of them into an expression in the call statement. … > It is the extra parentheses around the argument that does the trick. Indeed. It’s important to note that this technique must not be used with an output argument ! > [@RonShepard](#): > > I think the increment of 1 is correct in the first call above Yes it is, as the actual argument is a non-contiguous array slice, which is copied to a contiguous array because of the use of an implicit interface that assumes contiguous dummy arguments)