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

void mpl_daxpy(n, a, x, y)
/***************************************************************************
*                                                                          *
*   DATA PARALLEL BLAS based on MPL                                        *
*                                                                          *
*   Internal routine, this routine is not supposed to be                   *
*                     called by user programs.                             *
*                                                                          *
*   Version 1.0   20/9-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:                                                             *
*                                                                          *
***************************************************************************/
 
        int	n;
	register PRECISION a;
	register plural	PRECISION *x;
	register plural	PRECISION *y;

{
    /* local variables */
    register	plural	PRECISION x1, y1, x2, y2;
    register		int	i, no_block, rest_vec; 
    register plural	PRECISION *d;

/* Preload daxpy. Loads and stores from/to PMEM (the "M-machine") are
 * overlapped with subsequent execution instructions (the "E-machine").
 * By priming the loop with an extra initial load, we can then overlap
 * subsequent loads with calculations required in the main body of the
 * loop. Any operations not involving the registers that are being loaded
 * into will operate concurrently with the memory fetch. If enough operations
 * are carried before the values being loaded are needed, the load is
 * essentially free of cost.
 */


    no_block = NFB(n);   /* number of full blocks */
    rest_vec = n % nproc; 

    if (no_block == 0) 
       if (ENDVEC(n)) *y += a* *x;
    else if ((no_block == 1) && (rest_vec == 0))
       *y += a* *x;
            
    if (no_block > 0){  
    /* Sett pointer for where to store = pointer 
     * for where to fetch
     */

    if (rest_vec > 0) no_block++;  /* Total number of blocks */
    d = y;

    /* Preload first element from x & y arrays.
     */

    x1 = *x;
    y1 = *y;
    x++;
    y++;

    /* Main loop: unroll to level two to avoid having to move x2 & y2
     * into x1 & y1. Loop is terminated one iteration early to avoid
     * loading one element beyond the end of the x & y arrays when n
     * is even.
     */
    for (i=0; i<(no_block-1)>>1; i++) {

	/* Preload even elements.
	 */
	x2 = *x;
	y2 = *y;
	x++;
	y++;

	/*	*d++ = a * x1 + y1;
	 */
	x1 *= a;
	x1 += y1;
	*d = x1;
	d++;

	/* Preload odd elements.
	 */
	x1 = *x;
	y1 = *y;
	x++;
	y++;

	/*	*d++ = a * x2 + y2;
	 */
	x2 *= a;
	x2 += y2;
	*d = x2;
	d++;
    }

    /* Load last element from x & y arrays when n is even; all elements
     * have already been loaded for n odd. Make the two last operations
     * for n even.
     */

    if (!(no_block & 1)) {
	x2 = *x;
	y2 = *y;
        x1 *= a;
        x1 += y1;
        *d = x1;
        d++;   
	if (ENDVEC(n)){
	   x2 *= a;
	   x2 += y2;
	   *d = x2;
	}
    }

    /* Last operation for n odd. 
     */
    else{
	if (ENDVEC(n)){
           x1 *= a;
           x1 += y1;
           *d = x1;
	}
    }
    } 
}

