#include <mpl.h> 
#include "mpl_blas.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:                                                             *
*                                                                          *
***************************************************************************/


/* A solver for lower-triangular system, for nxproc*nxproc  Matrix 
 on a rectangular MasPar (nyproc = nxproc/2). The matrix is suppose
to be lower triangular, but 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 date 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 avoid the division in the inner loop,
we actually compute LDx = b; where L is lower 
unit triangular.

*/
#ifdef __STDMPL__

void mpl_slrsolve(int diag,
                  register int m,
                  register plural PRECISION a0,
                  register plural PRECISION a1,
                  plural PRECISION *b0,
                  plural PRECISION *b1)
#else

void mpl_slrsolve(diag,m,a0,a1,b0,b1)

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;
register plural int ix = ixproc, iy = iyproc,iy2;
register plural PRECISION btemp,temp;

if (iy==0) btemp = *b0; else btemp = 0.0;
iy2 = iy + ny;
iter = ny;
if (m<ny) iter = m;


/* scaling by the diagonal, using x as a temporary variable */
if (ix == iy) xnetcE[m].temp = a0;
if (ix > iy && ix < m) a0 = a0/temp;
   
/* first half */
for (i=0; i< iter; i++) {
   if (iy == i) {
	if(ix ==i) xnetcE[m-i].temp = btemp; 
   	if (ix > i) {
		btemp -= a0*temp;
		xnetS[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 == iy)&&(ix<m)) btemp = btemp/a0;
if (ix == iy) xnetcE[nx-1].btemp = btemp;
*b0 = btemp;
btemp = temp;

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

/* if scaled, scale the results */

if (diag != 'U' && diag != 'u')	if ((ix == iy2)&&(ix<m))  btemp = btemp/a1;
if (ix == iy2)xnetcE[nx-1].btemp = btemp;
*b1 = btemp;

}