#include <mpl.h> 
#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:                                                             *
*                                                                          *
***************************************************************************/
/* An upper triangularsolver for nxproc * nxproc Matrix
 on a rectangular MasPar MP-1 (nyproc = nxproc/2).
 
The physical lay-out on the
processor array, transposes the matrix.

The algorithm is the straightforward column update algorithm,
(but remember: Matrix columns = Physical rows of processors,
Thus we update data on rows of processors)

For 64 <m < 128 The non-zero elements are distributed as follows:

	| *  ...           * 0 ... 0 |
    	| 0 *  ..........  * 0 ... 0 |
        | : :              * 0 ... 0 |  64
        | 0 ...  0 *  .... * 0 ... 0 |
        | ===========================| ===
        | 0 ......0 *  ... * 0 ... 0 |
        | :        :       : 0 ... 0 |  m - 64
        | 0  ............0 *  0 ... 0| 
	| 0  .......................0|
        |       m           | 128 -m |



The two parts are stack upon each other. Thus
l(i,j) is a resident of processor (j,i(mod 64)).

The elements of b, b(i), rest on (0,j) and travels downwards
as the update take place. When they arrive at the diagonal
the division takes place and they take the value of the unknown
(x), which is return in the final column.

To improve performence, we actually compute LD' = b; where L is lower 
unit triangular. We then don't have to do the division of the pivot
in n step, but in one single in the beginning. Since division is slow
on MasPar we expect this to give great improvement in the performence

*/

#ifdef __STDMPL__
void mpl_sursolve(
	register int diag,
	register int m,
	register plural PRECISION a0,
	register plural PRECISION a1,
	plural PRECISION *b0,
	plural PRECISION *b1
	)
#else
void mpl_sursolve(diag,m,a0,a1,b0,b1)
	register int diag;
	register int m;
	register plural PRECISION a0,a1;
	plural PRECISION *b0,*b1;
#endif

{
register int i,ii,j,nx= nxproc,ny= nyproc,iter,mmod;
register plural int ix = ixproc, iy = iyproc,iy2;
register plural PRECISION btemp,temp;

/* Putting the righthand side into correct starting position */
mmod = m%ny-1;
if (mmod>=0)
	{if (iy== mmod ) btemp = xnetpN[mmod].*b0; else btemp = 0.0;}
else
 	{if (iy== ny-1 ) btemp = xnetS[1].*b0; else btemp = 0.0;}

iy2 = iy + ny;


/* scaling by the diagonal */
if (ix == iy2) xnetcW[m].temp = a1;
if (ix < iy2 && iy < m-ny) a1 = a1/temp;
   
/* first half */
for (i=m-1; i>=ny; i--) {
   ii = i-ny;
   if (iy == ii) {
	if(ix ==i) xnetcW[i].temp = btemp; 
   	if (ix < i) {
		btemp -= a1*temp;
		xnetN[1].btemp = btemp;
	}
   }
}

/* if scaled, scale the results, spread it out and put into the pointer */

temp = btemp;
if (diag != 'U' && diag != 'u')	if ((ix == iy2)&&(iy<m-ny)) btemp = btemp/a1;
if (ix == iy2) xnetcW[nx-1].btemp = btemp;
*b1 = btemp;
btemp = temp;
/* second half */

/* scaling by the diagonal  */
if (ix == iy) xnetcW[m].temp = a0;
if (ix < iy && iy < m) a0 = a0/temp;
   
iter = ny;
if (m<ny) iter = m;
for (i=iter-1; i>=0; i--) {
 if (iy == i){ 
 	if( ix  ==i) 	xnetcW[i].temp = btemp;
 	if (ix < i){
       		btemp -= a0*temp;
       		xnetN[1].btemp = btemp;
 	}
 }
}

/* if scaled, scale the results */

if (diag != 'U' && diag != 'u')	if ((ix == iy)&&(ix<m))  btemp = btemp/a0;
if (ix == iy)xnetcW[nx-1].btemp = btemp;
*b0 = btemp;

}