#include <mpl.h>
#include <stdio.h>
#include <stdlib.h>
#include "mpl_blas.h"
#include <sys/time.h> /* For test purposes only*/
#define PRECISION float


/***************************************************************************
*                                                                          *
*   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:                                                             *
*                                                                          *
***************************************************************************/
/*  The following routine performs lower triangular solve with multiple
 *  right hand sides using blocking on the Maspar MP-1.
 *  Each block is square and has sides == nxproc. On rectangular machines
 *  it will therefore use two grids per block.
 *  The routine performs triangular solves on the diagonal blocks in the
 *  A matrix and the B-blocks level to them, and updates the B-blocks below
 *  using matrix multiplication.
 *  The routine expects data to be stored in the Fortran 90 fashion with
 *  cut and stack.
 *
 *  Last bug fix: July 6. -92, TS
 *
 */
#ifdef __STDMPL__
void mpl_strsm_lon(
	int int_diag,
        int m,
        int n,
	plural PRECISION  *a,
	int lda,
	plural PRECISION  *b
	)
#else
void mpl_strsm_lon(int_diag,m,n,a,lda,b)
	int           int_diag,
		      m,n;
	plural PRECISION  *a;
	int           lda;
	plural PRECISION  *b;
#endif
{
 char        diag;


/* 
 * Local variables
 */
 plural PRECISION       *b0_ps,*b1_ps,a0_ps,a1_ps;
 register int        nx=nxproc,ny=nyproc,nbxA,nbyA,nbxB,nbyB,i,j,k,bn,am,t;
 register plural int ix=ixproc,iy=iyproc;
 int l;

 if (int_diag)
   diag = 'u';
 else
   diag = 'n';

 /* 
  * NB ! Below it is meant by "number of blocks", the number of grids 
  * used to store the arrays. The following four variables that appear 
  * on the lhs's are used to keep track of the stacking of the A and B
  * matrices on the dpu.
  */
 nbxA = NBX(m); /* Number of blocks in row direction on A matrix */
 nbyA = NBY(m); /* Number of blocks in column direction on A matrix */
 nbxB = NBX(m); /* Number of blocks in row direction on B matrix */
 nbyB = NBY(n); /* Number of blocks in column direction on B matrix */


 /* Variables b0_ps and b1_ps are temporaries used to hold preskewed  blocks
  * of the B matrix to be used in matrix multiplications. They are arrays
  * and hold multiple blocks. Each block is divided into two parts such that
  * the first half is stored in b0_ps[k] and the second is stored b1_ps[k].
  * The same goes for the variables a0_ps and a1_ps, except that they are not
  * arrays, and they are used for the A matrix blocks.
  */
 b0_ps = (plural PRECISION *) p_malloc(((nbyB+1)/2)*sizeof(PRECISION));

 /* initializing boundary blocks to zero on processors off bounds */

 bn = n&(ny-1); 
 am = m&(nx-1);

 if ((ix >= am)&&(am>0)){ 
   for (i=1;i<=nbyB; i++)
	b[i*nbxB-1] = 0.0;
 }
 if ((iy >= bn)&&(bn>0)) { 
   for (i=0;i<nbxB; i++)
	b[i+(nbyB-1)*nbxB] = 0.0;
 }
  
 /* Start processing. Using seperate code for square and rectangular
    MP-1 */

 if (SQUARE) 
   {              
    for (i=0;i<nbyA;i++)
      {
       /* am is the dimension of the current block of the A matrix */
       if ((i == nbxA-1) && (m&(nx-1)))
	 am = m&(nx-1);
       else
	 am = nx;


       /* Performing solve on diagonal block A(i,i) with all corresponding
	* rhs blocks ( B(i,k) blocks ) and preskew the solved
	* rhs blocks.
        */
       for(k=0;k<nbyB;k++)
         {
	  if ((k == nbyB-1) && (n&(nx-1)))
	    bn = n&(nx-1); /* n%nxproc */
          else
	    bn = nx;
	  mpl_stri_sln(diag,am,bn,a[i*nbxA+i],&(b[i+nbxB*k]));


          /* preskew b blocks */
	  mpl_sq_spsW(b[i+nbxB*k],&(b0_ps[k]),bn,am); 
         }

        /* Updating blocks below current diagonal block using matrix-
         * multiplication.
	 * mpl_sq_smul takes care of active set */                       

	for(j=i+1;j<nbyA;j++)
          {
	   if ((j == nbyA-1) && (m&(ny-1)))
	     am = m&(ny-1);
           else
	     am = ny;
	   mpl_sq_spsN(a[j+nbxA*i],&a0_ps,nx,am);
           a0_ps *= -1;
           for(k=0;k<nbyB;k++)
	     {
	      if ((k == nbyB-1) && (n&(ny-1)))
	        bn = n&(ny-1); /* n%nyproc */
              else
	        bn = ny;
	      mpl_sq_smul('n',a0_ps,b0_ps[k],&b[j+nbxB*k],am,bn,nx);
             }
          }
      }
   }
 else  /**************** if machine is rectangular ****************/
   {
    b1_ps = (plural PRECISION *) p_malloc(((nbyB+1)/2)*sizeof(PRECISION));
    for (i=0;i<(nbyA+1)/2;i++)
      {
       /* 
	* am is the dimension of the current block of the A matrix,
	* and bn is # columns in current block of the B matrix.
	*/
       if ((i == (nbyA+1)/2-1) && (m&(nx-1)))
         am = m&(nx-1);
       else
         am = nx;

       /* Performing solve on diagonal block A(i,i) with all the rhs blocks
	* ( B(i,k) blocks ) and preskewing the solved
	* rhs blocks.
        */
       for(k=0;k<(nbyB+1)/2;k++)
         {
	  if ((k == (nbyB+1)/2-1) && (n&(nx-1)))
	    bn = n&(nx-1); /* n%nxproc */
          else
            bn = nx;
          mpl_stri_rln(diag,am,bn,a[i*nbxA*2+i],a[(2*i+1)*nbxA+i],
		     &(b[i+nbxB*2*k]),&(b[(2*k+1)*nbxB+i]));

          /* preskew b-blocks */
	  if (i< (nbyA+1)/2-1)
	    mpl_rec_spsW(b[i+nbxB*2*k],b[(2*k+1)*nbxB+i],&(b0_ps[k]),
		    &(b1_ps[k]),bn,am);
         }

        /* Updating blocks below current diagonal block using matrix-
         * multiplication.
	 * Padding with zeroes in boundary blocks are done in the
	 * multiplication routine */

       for(j=i+1;j<nbxA;j++)
         {
          if ((j == nbxA-1) && (m&(nx-1)))
            am = m&(nx-1);
          else
            am = nx;
	  mpl_rec_spsN(a[j+nbxA*2*i],a[(2*i+1)*nbxA+j],&a0_ps,&a1_ps,nx,am);
          a0_ps *= -1;
          a1_ps *= -1;
	  for(k=0;k<nbyB/2;k++)
	    {
	     if ((k == (nbyB+1)/2-1) && (n&(nx-1)))
	       bn = n&(nx-1); /* n%nxproc */
             else
               bn = nx;
	     mpl_rec_smul('n',a0_ps,a1_ps,b0_ps[k],b1_ps[k],
			&b[k*nbxB*2+j],&b[(2*k+1)*nbxB+j],nx,bn,am);
	    }
          if (nbyB & 1){
            k = (nbyB-1)/2;
	    if ( n&(nx-1))
	       bn = n&(nx-1); /* n%nxproc */
            else
               bn = nx;
            mpl_rec_smul('n',a0_ps,a1_ps,b0_ps[k],b1_ps[k],
                    &b[j+nbxB*2*k],b1_ps,nx,bn,am);
            }
         }
      }
   p_free(b1_ps); 
   }
 p_free(b0_ps);
}