Re: [Help-gsl] gsl blas idamax

[Rhys Ulerich]

Hi Andrew,

> If I call gsl blas idamax as follows:
>
>    double x[3] = { 11, 33, 22 };
>    . . .
>    i = cblas_idamax(3, x, 1);
>
> then which returned value is correct, 1 or 2 ? Gsl blas idamax returns
> 1 that corresponds to 0-base indexing. However, AFAIK, some C blas
> implementations (e.g. IBM ESSL) follow Fortran conventions; otherwise,
> the standard blas tests would fail.

Looking into cblas_idamax.c from
http://www.netlib.org/blas/blast-forum/cblas.tgz shows the reference
implementation as follows...

/*
 * cblas_idamax.c
 *
 * The program is a C interface to idamax.
 * It calls the fortran wrapper before calling idamax.
 *
 * Written by Keita Teranishi.  2/11/1998
 *
 */
#include "cblas.h"
#include "cblas_f77.h"
CBLAS_INDEX cblas_idamax( const int N, const double *X, const int incX)
{
   int iamax;
#ifdef F77_INT
   F77_INT F77_N=N, F77_incX=incX;
#else
   #define F77_N N
   #define F77_incX incX
#endif
   F77_idamax_sub( &F77_N, X, &F77_incX, &iamax);
   return iamax ? iamax-1 : 0;
}

where the reference Fortran routine is
http://www.netlib.org/blas/idamax.f.  So, I expect 0-based indexing
from GSL because the reference implementation uses ' iamax ? iamax-1 :
0' to convert to 0-based indexing.

> I also would like to note that in case of n = 0 gsl blas idamax
> (see blas/source_iamax_r.h) returns 0 rather than -1. Probably it is
> a bug, since the reference Fortran version of idamax returns 0 in this
> case.

This is odd but not a bug, it seems.  When N = 0
http://www.netlib.org/blas/idamax.f returns 0 and 'iamax ? iamax-1 :
0' evaluates to 0.  Or, rather, it's an API bug but GSL is consistent
with the reference implementation given the behavior you report.

- Rhys