Re: Problem using BLAS' xGEMM
Michel Salim wrote:
> I'm definitely doing something wrong here, but could not figure out
> what the error is - could anyone help?
> Given a (sub)matrix A, with nrows number of rows and ncols number of
> columns and leading dimension LDA, I'm trying to convert this
> expression to one that uses DGEMM
> A(k, k+1:ncols) = A(k, k+1:ncols) - A(k, 1:k-1) * A(1:k-1, k+1:ncols)
> DGEMM calculates C = alpha * op(A) * op(B) + beta*C
> DGEMM parameters:
> M = # rows of C = 1
> N = # cols of C = ncols-k
> K = # rows of B = k-1
> alpha = -1.0d0
> beta = 1d0
> A = A(k, 1:k-1)
> B = A(1:k-1, k+1:ncols)
> C = A(k, k+1:ncols)
> resulting in the call
> CALL dgemm('N', 'N', 1, (ncols-k), (k-1), -1.0d0, A(k, 1:k-1), lda,
> A(1:k-1, k+1:ncols), lda, 1.0d0, A(k, k+1:ncols), lda)
> When tested on A =
> 3 -1 -2 -1
> 2 -7 -5 -1
> 4 -21 -9 6
> 4 0 -2 15
> it should result in updating the last two elements of the second row to
> -1 1, but instead the result is -1 -1.
> Any pointer really appreciated.
> - Michel Salim
I guess a temporary copy is done for the noncontiguous arrays, which
means that their leading dimension changes to 1.
btw your case is actually a matrix-vector multiply, it would be better
handled using DGEMV.
Alternatively, you can use my smart BLAS95 interface and call
(the interface calulates correct M,N,K,LDA,LDB and LDC and calls xGEMM)
in which case DGEMV is called instead.