[Help-gsl] blas rotm index computation

[Marco Maggi]

Ciao,

  I see this sequence of nested calls (with non-relevant bits
stripped of):

int
gsl_blas_drotm (gsl_vector * X, gsl_vector * Y, const double P[])
{
  cblas_drotm (INT (X->size),
               X->data, INT (X->stride),
               Y->data, INT (Y->stride), P);
}
void
cblas_drotm (const int N, double *X, const int incX,
             double *Y, const int incY, const double *P)
{
#define BASE double
#include "source_rotm.h"
#undef BASE
}
/* cblas/cblas.h */
#define OFFSET(N, incX) ((incX) > 0 ?  0 : ((N) - 1) * (-(incX)))
/* blas/source_rotmg.h */
{
  ...
  INDEX i = OFFSET(N, incX);
  INDEX j = OFFSET(N, incY);
  ... 
  for (n = 0; n < N; n++) {
    const BASE w = X[i];
    const BASE z = Y[j];
    ...
    i += incX;
    j += incY;
  }
}

so when the vector stride is negative the elements
of the vector are iterated starting from index:

i = (X->size - 1) * -(X->stride);

which equals:

i = (X->size - 1) * abs(X->stride);

no? But is this correct? "X->data" references the
first element in the vector, not the memory slot
in the block with the miminum address value. So
that, as documented in the node "Accessing vector
elements", it is possible to do:

X->data[i * X->stride]

to get element at index "i". What happens when
one vector has positive stride and the other negative
stride?

Notice that the test suite in "cblas/test_rotm.c"
tests vectors of size==1, so this problem does not
show itself.
-- 
Marco Maggi