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

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
call gemm(A(k:k,1:k-1),A(1:k-1,k+1:ncols),A(k:k,k+1:ncols))
(the interface calulates correct M,N,K,LDA,LDB and LDC and calls xGEMM)
or
call
gemm(A(1:k-1,k+1:ncols),A(k,1:k-1),A(k,k+1:ncols),transa=blas_trans)
in which case DGEMV is called instead.