#include <mpl.h>
#include <stdio.h>
#include <stdlib.h>
#include "mpl_blas.h"
#define PRECISION float
/* #define DEBUG */

/***************************************************************************
*                                                                          *
*   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:                                                             *
*                                                                          *
***************************************************************************/
#ifdef __STDMPL__ 
visible void 
      mpl_sgemv2(int m, int n, 
                 PRECISION alpha, 
                 plural PRECISION *a, 
                 plural PRECISION *x, 
                 PRECISION beta,  
                 plural PRECISION *y)
#else
visible void 
      mpl_sgemv2(m, n, alpha, a, x, beta, y)
      PRECISION  alpha, beta;
      int m, n;                    
      register  plural PRECISION *a,*x,*y;
#endif

	/* MasPar Blas 2 routine

		y := beta*y + alpha*A'*x
	
        written July 20, 1990 
        revised November 1, 1990
        NEW VERSION Dec 12, 1990  (array A has columns in the x-direction)
    
        by Erik Boman, University of Bergen 


     * Comments:

     * works only for m and n < nproc 
     * the use of register variables should be optimized

    */

{
      register plural PRECISION sum,xtmp;
      register plural PRECISION xtrans,ytrans;
      register plural PRECISION *xtemp;
      register int i,j; 
      register int mask;
      int nbx,nby,nrows;
      register plural int ix=ixproc, iy=iyproc;

      /* put pointer arguments into local variables */

	if (m>nproc) {
		printf("Error in mpl_sgemv2 : m > nproc\n");
		exit(1);
	}

	if (n>nproc) {
		printf("Error in mpl_sgemv2: n > nproc\n");
		exit(1);
	}


/*
*        Form  A'*x
*/

         if ((m <= nxproc) && (n <= nyproc)){

            sum = zero;

            /* distributing x */

            if (iy == 0){
                xtmp = *x;
                xnetcS[n-1].xtmp = xtmp;
            }

            /* multiplying */

            if (iy<n && ix<m)
                sum = *a * xtmp; 

            /* summing up in log(m) steps */

            for (mask=1; mask<m; mask<<=1){
                if (!(ix&((mask<<1)-1)))
                   sum += xnetpE[mask].sum;
            }
        
			/* new faster "transpose" */

			if (ix ==0)
				xnetcNE[n-1].sum = sum;

/*
*              y = beta*y + alpha*(A'*x)
*/

            if (iproc < n) {

               if (beta == zero)
                  *y = zero;
               else if (beta == -one)
                  *y = - *y;
               else if (beta != one)
                  *y *= beta;

               if (alpha == one)
                  *y += sum;
               else if (alpha == -one)
                  *y -= sum;
               else
                  *y += alpha*sum;
            }
         }

         else { /* more than one block */
            
/*
*        Compute number of blocks
*/

         nbx = NBX(m);
         nby = NBY(n);

         /* first distribute x */

         xtemp = (plural PRECISION *)p_malloc(nbx*sizeof(PRECISION));

         for (j=0; j<nbx; j++){
                if (iy == j){
                    xtmp = *x;
                    xnetcS[nyproc-1].xtmp = xtmp;
                }
                xtemp[j] = xtmp;
         }

         for (i=0; i<nby; i++){
             sum = zero;
             nrows = n-i*nyproc;
             for (j=0; j<nbx; j++){
/*
#ifdef DEBUG
                printf("i= %d j=%d \n",i,j);
                show(&a[i*nbx+j],nyproc,nxproc,"matrix A'");
#endif
*/
                /* do a block matrix-vector multiplication */

            
                /* then multiply and add */
            
                if ((iy<nrows)
                    && (ix< m-j*nxproc))  
                       sum += a[i*nbx+j]*xtemp[j];
/*
#ifdef DEBUG
             { plural double debug = sum;
               printf("j = %2d\n",j);
               show(&xtemp[j],nyproc,nxproc,"x in loop");
               show(&debug,nyproc,nxproc,"cumulative A'*x in loop");
             }
#endif
*/
             }


           /* summing up rows in log(nxproc) steps using xnetp */
  
             for (mask=1; mask <nxproc; mask<<=1){
                 if(!((ix+nxproc-i)&((mask<<1)-1)))
          	        sum += xnetpE[mask].sum;
             }

           /* store the result in column i */

           if (ix == i)
              ytrans = sum;
/*
#ifdef DEBUG
             { plural double debug = ytrans;
               show(&debug,nyproc,nxproc,"ytrans");
             }
#endif
*/
		}

     	/* transpose (towards NE) to put A'*x in the right place ! */

        router[ix*nyproc + iy].sum = ytrans; 
/*
#ifdef DEBUG
             { plural double debug = sum;
               show(&debug,nyproc,nxproc,"sum");
             }
#endif
*/ 
        /*
        *     y = beta*y + alpha*(A'*x)
        */


            if (iproc < n) {

               if (beta == zero)
                  *y = zero;
               else if (beta == -one)
                  *y = - *y;
               else if (beta != one)
                  *y *= beta;

               if (alpha == one)
                  *y += sum;
               else if (alpha == -one)
                  *y -= sum;
               else
                  *y += alpha*sum;
            }

        p_free(xtemp);

       }  /* end array > PE grid  */

/* 
*     End of mpl_sgemv2 .
*/
}