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.
>
> Thanks,
>
> - Michel Salim

Unless the documentation of a library routine such as DGEMM explicitly
declares it to be safe, it is risky to have input and output array
arguments overlap or coincide. Elements of the input array may get
clobbered before they are used.