/*-----------------------------------------------------------------------------
 * File:  transarray.c
 *
 * Purpose:
 *    General purpose utility routines
 *
 * Invokes:
 *
 * Public functions:
 *    transar - transpose the array and reverse the dimension sizes
 *
 * Lower level functions:
 *
 * Private functions:
 *
 * Remarks:
 *
 *--------------------------------------------------------------------------*/


/*----------------------------------------------------------------------------
 * Name:    transar
 * Purpose: Transpose in memory a data array and reverse its dimension
 *          sizes accordingly. 
 *          It can be used by C programs to transpose
 *          SDS's which were written by Fortran programs using HDF library
 *          previous to HDF3.2r1 (or HDF3.2beta3)
 * Inputs:  rankp: pointer to the rank of the data array
 *          indims: dimension sizes of input (source) data array
 *          outdims: dimension sizes of output (destination) data array
 *          eltsizep: pointer to the size of the data 
 *          indata: source data array
 *          outdata: destination data array --transposation of indata
 *          isfortranp: pointer to the indicator of the caller,
 *			*(isfortranp)=0 if the caller is a C program
 * Returns: 0 on success, -1 on failure 
 * Outputs: the transposed array and reversed dimension sizes
 * Users:   a general utility function. can be used by C programs to transpose
 *	    SDS's  which were written by Fortran programs using HDF library
 *          previous to HDF3.2r1 (or HDF3.2beta3)
 * Invokes:
 * Method:  reverse dimensions. transpose the data array. 
 * Remarks: HDF C programs may call DFKNTsize(numbertype | DFNT_NATIVE) to get 
 *          the memory size of the number type. The Fortran interface of 
 *          DFKNTsize needs to be written. 
 *          
 **********************************************************************/


#ifdef PROTOTYPE
int transar(int *rankp, int *indims, int *outdims,
		  int *eltsizep, char *indata, char *outdata, int *isfortranp)
#else
int transar(rankp, indims, outdims, eltsizep, indata, outdata,isfortranp)

    int     *rankp;	/* number of dimensions of the data array  */
    int     *indims,    /* array containing the coordinates of the input */
                        /* data array  */
            *outdims,	/* array to store the dimensions of the output */
                        /* data array  */
            *eltsizep;   /* size(number of bytes) of the element */
    char    *indata,	/* the input data array */
            *outdata;   /* the transposed data array -- output data array  */
    int     *isfortranp; /* 1--called from Fortran program, 0--from C prog. */

#endif /* PROTOTYPE */

{
    int
        done,           /* true if we are at the end of the slice */
        i, j,           /* temporary loop index */
        /* temp,		 used for reverse of the  dims   */
        rank,           /* rank of the array               */
        NTsize,         /* size of the data                */
        leastsig,        /* rank-1  */
/*        localNTsize,   */
                         /* size of this NT as it occurs in this machine */
        outstride,      /* byte distance,in outdata,of two adjacent elements */
			/* in indata[]           */
        *offset,        /* array for accessing the next element in data[] */
        *outoffset,     /* array for accessing the next element in outdata[] */
        *dims,          /* dimsizes used in transposition   */
        *dimsleft;      /* array for tracking the current position in data[] */
/*        isnative,       */
/*        machinetype;    */
                        /* assigned DF_MT.  used for debugging */
/*    uint8               */
/*        platnumsubclass, */
                        /* class of this NT for this platform */
     unsigned char
        *outdatap,      /* ptr into outdata[] at starting offset */
                           /* of current block  */
        *outdp,		/* ptr into outdata[] at an element of the current row */
        *datap,         /* ptr into data[] at starting offset */
                            /* of current block */
        *dp;            /* ptr into data[] at an element of the current row */
    /* char *FUNC="transarray"; */


    if (!indata || !outdata ) {
        return (-1); 
    }

/*    isnative = DFNT_NATIVE;
    machinetype = DF_MT;
    platnumsubclass = DFKgetPNSC(numtype, DF_MT);
    localNTsize = DFKNTsize(numtype | isnative);
*/
    rank = *rankp;
    NTsize = *eltsizep;
    /* reverse dimensions first  */
    for (i=0; i<rank; i++)    {
        outdims[i] = indims[rank-i-1];
    }
    if (*isfortranp == 1) dims=outdims;  /* use outdims, slowest dim first, to make */
                                        /* the dims consistant with the data in C order */
    else dims=indims;
    /*
     * To transpose the data we must work on a row by row
     * basis and cannot collapse dimensions.
     */
    
    leastsig = rank-1;		/* which is least sig dim */

    /* allocate buffers */
    offset = (int *) malloc(3 * rank * sizeof(int));
    if (!offset) {
        return (-1);
    } 
    outoffset = offset + rank;
    dimsleft = outoffset + rank;

        /* compute initial position in the data */
    for (i=leastsig; i>=0; i--)
        dimsleft[i] = dims[i];

        /* compute offsets for the source array */
    offset[leastsig] = 1*NTsize;
    for (i = leastsig; i>0; i--)
        offset[i-1] = offset[i] * dims[i];
    
        /* compute offsets for the destination array */
    outoffset[0] = 1*NTsize;
    for (i=0; i<leastsig; i++)
        outoffset[i+1] = outoffset[i] * outdims[leastsig - i];
    outstride = outoffset[leastsig];
    datap = (unsigned char *)indata;
    outdatap = (unsigned char *)outdata;
    done = 0;
        /* -- now read in the data */
    do {
            /* move to the next row in source array */
            /* scatter out the elements of one row */
#ifdef UNICOS
#pragma ivdep
#endif
        for (dp=datap, outdp = outdatap, i=0; i<dims[leastsig]; i++) {
             for (j=0; j<NTsize; j++)
                        *(outdp +j) = *(dp +j);
             dp += NTsize;
             outdp += outstride;
        }

            /*
             * Find starting place of the next row/block.
             * Note that all moves are relative:
             *   this preserves the starting offsets of each dimension
             */
        for (i=leastsig-1; i>=0; i--) {
            if (--dimsleft[i] > 0) {
                    /* move to next element in the current dimension */
                datap += offset[i];
                outdatap += outoffset[i];
                break;
            } else {
                dimsleft[i] = dims[i];
                   /*
                     * Note that we are still positioned at the beginning of
                     * the last element in the current dimension
                     */
                    /* move back to the beginning of dimension i */
                datap -= offset[i] * (dims[i]-1);
                    /* move back to beginning outdata position of dimension i */
                outdatap -= outoffset[i] * (outdims[i]-1);
                if (i==0) done = 1;
            }
        }
    } while (!done && leastsig > 0);

    free((char *)offset);
    return(0);
}   /* end transarray */
