Solution at the end
I’m having a little bit of trouble with Blas’ band matrix operation subroutine DGBMV.
From the man page:
DGBMV performs one of the matrix-vector operations y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y, where alpha and beta are scalars, x and y are vectors and A is an m by n band matrix, with kl sub-diagonals and ku super-diagonals.
And so I’m using this to compute a matrix-vector multiplication, therefore I’m setting beta=0.0 and alpha=1.0. I wrote a little function to help with the usage:
function lmm(A, x, ninj, nsup, nsub) result(Cout) implicit none real(8), intent(in) :: A(:,:), x(: ) real(8) :: Cout(size(x)) real(8) :: AC(size(A,1),size(A,2)), xC(size(x)) integer, intent(in) :: nsub, nsup, ninj AC(:,:) = A(:,:) xC(:) = x(:) Cout = 0.0 call dgbmv('n', ninj, ninj, nsup, nsub, 1.0,& AC, size(AC,1), xC, 1, 0.0, Cout, 1)
end function lmm
Everything works but the results are completely wrong, however when I use matmul the results are as expected. I’ve also noticed that if I set Cout to another value before calling the subroutine I get different results.
I’ve fount in the IBM docs that the vector y should not have any element in common with both the matrix A and the vector x, maybe the problem is related to that but to be honest I couldn’t figure it out.
To convert the matrix from the dense to sparse format I’m using the following function:
function denseToBand(A, nsup, nsub) result(AB) implicit none integer, intent(in) :: nsub, nsup real(8), intent(in) :: A(:,: ) integer :: j, i real(8) :: AB(2*nsub+nsup+1, size(A,1)) do j=1, size(A,2) do i = max(1, j-nsup), min(size(A,1), j+nsub) AB(nsup+1-j+i,j) = A(i,j) enddo enddo end function denseToBand
The Case I’m working on:(following the advice of Arjen)
I’m currently solving a 2D differential equation and and the coefficient matrix A (in Ax=b of the linear system) has the following visual representation:
In which the green are zeros. The diagonal is composed by columns of the domain grid, as it can be seen at the start and end of each of these sections the value is different corresponding to a boundary condition. This implies that the x vector (in Ax=b) is of the form:
|x_1,1| |x_2,1| ... |x_n-1,n| |x_n,n|
which essentially are the columns of the domain grid stacked on top of each other.
The explicit solution involves multiplying A by the previous solution so: matmul(A, x). When I use matmul I get the following result:
However when I use the Blas subroutine I get the following:
Which in this case means that the subroutine is outputting a 0 valued vector, i.e. the multiplication is outputting 0. And the only lines in the code I changed were from this:
RHS = dertg*dt*matmul(A, phi) + dt*sig + phi
phiCopy = phi RHS = dertg*dt*lmm(Adense, phiCopy, ninj, nsup, nsub) + dt*sig + phi
Thank you so much.
Following Ivan advice I set the parameters alpha and beta to double precision and it worked perfectly. It should also be noted that my original formulation for the banded matrix was wrong, which is the appropriate for Blas’ DGBSV (linear system solver), or DGBTRF as pointed out by Ivan. For the DGBMV the sparse formulation is the one provided in the docs and pointed out by Ivan.
- Use double precision in both alpha and beta parameters
- Use appropriate sparse matrix formulation provided by the docs.
I cannot thank enough Ivan for his thorough replies !
Edit 1-3: Formatting
Edit 4: changed the function used to convert the dense to band matrix since I’m now using the one presented in the DGBMV documentation and I’m still getting the same results. This is just to clarify following the suggestion of Ivan.
Edit 5: Supplying the example.
Edit 6: Solution provided at the end the problem was the precision of alpha and beta.
Edit 7: Typos