#include <mpl.h>
#include "mpl_blas.h"
#define PRECISION float

visible void

mpl_isamax(n, x, max_index)
/***************************************************************************
*                                                                          *
*   DATA PARALLEL BLAS based on MPL                                        *
*                                                                          *
*   Internal routine, this routine is not supposed to be                   *
*                     called by user programs.                             *
*                                                                          *
*   Version 1.0   1/4-92 ,                                                 *
*   For MasPar MP-1 computers                                              *
*                                                                          *
*   para//ab, University of Bergen, NORWAY                                 *
*                                                                          *
*   The calling sequence may be changed in a future version.               *
*   Please report any BUGs, ideas for improvement or other                 *
*   comments to                                                            *
*                    adm@parallab.uib.no                                   *
*                                                                          *
*   Future versions may then reflect your suggestions.                     *
*   The most current version of this software is available                 *
*   from netlib@nac.no , send the message `send index from maspar'         *
*                                                                          *
*   REVISIONS:                                                             *
*                                                                          *
***************************************************************************/
/*************************************************************************
**                                                                      **
**   Find the index of the largest element in a vector                 **
**   If there are several max elements then isamax returns              **
**   with the smallest index that "points" to a max element.            **
**                                                                      **
**************************************************************************/

int n,*max_index;
plural PRECISION *x;
{
    plural int    tmp_index;
    plural PRECISION tmp_max; 
    PRECISION max;
    int    nfb, nb;
    int    i, return_index;

    nfb = NFB(n);
    nb  = NB(n);

    /*
     *    There are two cases.: 
     *    - Array with less than or equal to nproc elements.
     *    - Array with more than nproc elements.
     */
    if( nb == 1 ) {
        if ( iproc < n ) {
			tmp_index = iproc;
                        max = reduceMaxd(ABS(*x));
			if(max == ABS(*x))
			return_index = reduceMin64(tmp_index);
        }
		return_index++;
	}
	else
	{
        /*
         *    Find local max og max index for each processor     
         *    for arrays bigger than nproc.                     
	 *    The biggest element for each processor are stored
	 *    in tmp_max
         */
        tmp_max = x[0];
        tmp_index = iproc;

        for(i=1 ; i < nfb ; i++ ) {
            if(ABS(x[i]) > ABS(tmp_max)) {
                tmp_max  = x[i];
                tmp_index = i*nproc+iproc;
            }
        }
        /*
         *    The rest block
         */
        if(ENDVEC(n))  
            if(ABS(x[nb-1]) > ABS(tmp_max)) {
                tmp_max = x[nb-1];
                tmp_index = (nb-1)*nproc+iproc;
            }       
        max = reduceMaxd(ABS(tmp_max)); 
        if(max == ABS(tmp_max)) 
            return_index = reduceMin64(tmp_index); 
            return_index++;
    }

   *max_index = return_index;
}