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

visible void 
      mpl_dgemv1(m, n, alpha, a, x, beta, y)
/***************************************************************************
*                                                                          *
*   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:                                                             *
*                                                                          *
***************************************************************************/

      PRECISION alpha, beta;
      int m, n;                    
      register  plural PRECISION *a,*x,*y;


	/* 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, ytmp;
      register plural PRECISION *xtemp;
      register int i,j; 
      register int mask;
      int nbx,nby,ncols;
      register plural int ix=ixproc, iy=iyproc;

      /* put pointer arguments into local variables */

#ifdef DEBUG
   extern show();
#endif

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

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


/*
*        Form  A*x 
*/

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

#ifdef DEBUG
              printf(" Parameters to mpl_dgemv1: m= %d n= %d \n", m,n);

           { plural double debug = *x;
              show(&debug,nyproc,nxproc,"x");
            }
#endif

            /* first "transpose" x  (towards SW) */
            if (iy == 0){
                xtmp = *x;
                xnetcSW[n-1].xtmp = xtmp;
            }

#ifdef DEBUG
            { plural double debug = xtmp;
              show(&debug,nyproc,nxproc,"x'");
            }
#endif

            sum = zero;

            /* distributing x */

            if (ix == 0){
                xnetcE[m-1].xtmp = xtmp;
            }

#ifdef DEBUG
            { plural double debug = xtmp;
              show(&debug,nyproc,nxproc,"x' after distribution");
            }
#endif
            /* multiplying */

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

#ifdef DEBUG
            { plural double debug = sum;
              show(&debug,nyproc,nxproc,"sum = Ax before summation");
            }
#endif

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

            for (mask=1; mask <n; mask<<=1){
                if (!(iy&((mask<<1)-1)))
                   sum += xnetpS[mask].sum;
            }

#ifdef DEBUG
            { plural double debug = sum;
              show(&debug,nyproc,nxproc,"sum = Ax after summation");
            }
#endif
        
            if (iproc < m) {

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

               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;
            }

#ifdef DEBUG
            { plural double debug = *y;
              show(&debug,nyproc,nxproc,"y");
            }
#endif
         }

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

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


         /* "transpose" x  (towards SW) */

         xtmp = *x;
         ytmp = router[ix*nyproc + iy].xtmp; 

 
         /* distribute x */

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

         for (i=0; i<nby; i++){
                if (ix == i){
                    xtmp = ytmp;
                    xnetcE[nxproc-1].xtmp = xtmp;
                }
                xtemp[i] = xtmp;
         }

         for (j=0; j<nbx; j++){
             sum = zero;
             ncols = m-j*nxproc;
             for (i=0; i<nby; i++){

#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< n-i*nyproc)
                    && (ix< ncols))  
                       sum += a[i*nbx+j]*xtemp[i];
                             
             }

#ifdef DEBUG
             { plural double debug = sum;
               show(&debug,nyproc,nxproc,"A*x");
             }
#endif


           /* summing up columns in log(nyproc) steps using xnetp */
  
             for (mask=1; mask <nyproc; mask<<=1){
                 if(!((iy+nyproc-j)&((mask<<1)-1)))
          	        sum += xnetpS[mask].sum;
             }

           /* store the results */

           if (iy == j)
              ytmp = sum;

#ifdef DEBUG
             { plural double debug = ytmp;
               show(&debug,nyproc,nxproc,"ytmp");
             }
#endif

		}

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

        if (iproc < m) {
           if (beta == zero)
              *y = zero;
           else if (beta == -one)
              *y = - *y;
           else if (beta != one)
              *y *= beta;

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

        p_free(xtemp);

       }  /* end array > PE grid  */

/* 
*     End of mpl_dgemv1 .
*/
}