#include <stdio.h>
#include <math.h>
#include "pvm3.h"

/*
    PVM 3.0 version
	Node program for column Cholesky factorization of
	a symmetric matrix partitioned by vertical strips
	using synchronous fan-out broadcast communication
*/

main()
/* node process */
{
	int mytid, tids[100];
	long ti[8] ;
	int n, i, ncols, p, me, alg, *map, *mycols, info[3] ;
	float *b, **col, tm[4], err, maxerr ;

    /* enroll in pvm */
    mytid = pvm_mytid();

    /* receive tids array */
    pvm_recv(-1, 16000);
    pvm_upkint(&p,1,1);
    pvm_upkint(tids, p, 1);
    for( i=0; i<p ; i++ )
       if( mytid == tids[i] ){ me = i; break; }

	pvm_recv(-1,16001);
	ti[6] = time(0);
	pvm_upkint(info,3,1);

	n = info[0] ;
	alg = info[1] ;
	p = info[2];
	if (n == 0 || n < p) {
		fprintf(stderr,"Node no %d : n=%d; p=%d\n",me,n,p);
		exit(1);
	}
	fprintf(stderr,"Node no %d : n=%d; p=%d\n",me,n,p);
	map = (int *)malloc(n*sizeof(int)) ;

	pvm_recv(-1,16002);
	pvm_upkint(map,n,1);
	fprintf(stderr, "Received 16002\n");

	ncols = 0 ;
	for ( i = 0 ; i < n; i++ ) if (map[i] == me) ncols++ ;

	/* allocate storage */
	mycols = (int *)malloc(ncols*sizeof(int)) ;
	b = (float *)malloc(ncols*sizeof(float)) ;
	col = (float * *)malloc(ncols*sizeof(float *)) ;

	fprintf( stderr, "Storage Allocated\n" );

	/* set up data structures for problem to be solved */
	prob1(me, n, ncols, map, mycols, col, b) ;

	fprintf( stderr, "after prob1\n" );

	/* Cholesky factorization - fan-out synchronous algorithm */
	ti[0] = time(0);
	cholfo(n, ncols, map, mycols, col, me, p, tids) ;
	ti[1] = time(0);

	fprintf( stderr, "after cholfo\n" );

	/* forward substitution */
	ti[2] = time(0);
	forsub(n, ncols, map, mycols, col, b, me, p, tids) ;
	ti[3] = time(0);

	fprintf( stderr, "after forsub\n" );

	/* back substitution */
	ti[4] = time(0);
	backsub(n, ncols, map, mycols, col, b, me, p, tids) ;
	ti[5] = time(0);

	fprintf( stderr, "after backsub\n" );

	/* compute maximum error in solution and send to cholhost */
	maxerr = 0 ;
	for ( i = 0 ; i < ncols ; i++ ) {
		err = b[i] ;
		err -= mycols[i]+1 ;
		maxerr = maxerr < fabs(err) ? fabs(err) : maxerr ;
	}
	pvm_initsend(PvmDataRaw);
	pvm_pkfloat(&maxerr,1,1);
	pvm_send(pvm_parent(),16003);

	/* send local times to host */
	ti[7] = time(0);
	pvm_initsend(PvmDataRaw);
	for ( i = 0 ; i < 4 ; i++ ) {
		tm[i] = ti[2*i+1]-ti[2*i] ;
		pvm_pkfloat(&tm[i],1,1);
	}
	pvm_send(pvm_parent(),16004);

	/* free storage */
	for ( i = 0 ; i < ncols ; i++ ) free(col[i]) ;
	free(col) ;
	free(b) ;
	free(mycols) ;
	free(map) ;
	
	pvm_exit();
	exit(0);
}


prob1 ( me, n, ncols, map, mycols, col, b )
	int me, n, ncols, *map, *mycols ;
	float **col, *b ;

/*
 *  Set up data structures for problem to be solved
 */
{
	int i, j, k ;
	float *p, t ;

	j = 0 ;
	for ( k = 0 ; k < n ; k++ )
		if ( map[k] == me ) {  
			p = (float *)malloc((n-k)*sizeof(float)) ;
			col[j] = p ;
			mycols[j] = k ;
			j++ ;
			if ( k == 0 ) {
				t = n ;
				*p++ = t*(t+1)/2 ;
				for ( i = 1 ; i < n ; i++ ) *p++ = -1 ;
			}
			else {
				*p++ = 2 ;
				for ( i = 1 ; i < n-k ; i++ ) *p++ = 0 ;
			}
		}
		for ( i = 0 ; i < ncols ; i++ ) b[i] = 2*mycols[i]+1 ;
}

cholfo ( n, ncols, map, mycols, col ,me, nprocs, tids)
	int n, ncols, *map, *mycols ,me;
	float **col ;
	int nprocs, *tids;

/* synchronous Cholesky factorization */
{
	int  i, j, k ;
	float t, *p, *q, *colk ;

	colk = (float *)malloc(n*sizeof(float)) ;
	j = 0 ;
	for ( k = 0 ; k < n ; k++ ) {
		printf("me=%d k=%d\n", me, k);
		if (map[k] == me) {
			t = 1.0/sqrt(*col[j]) ;
			for ( p = col[j] , q = colk ; p < col[j]+n-mycols[j] ;
				p++ , q++ ) {
				*p *= t ;
				*q = *p ;
			}
			j++ ;
		}
		if (map[k] == me) {
			if (nprocs > 1) {
				pvm_initsend(PvmDataRaw);
				pvm_pkfloat(colk,(n-k),1);
				pvm_mcast(tids,nprocs,k); 
			    /* pvm_recv(-1,k); /* mcast sends to itself */
			}
		}
		else {
			pvm_recv(-1,k);
			pvm_upkfloat(colk,(n-k),1);
		}

		for ( i = j ; i < ncols ; i++ ) {
			for ( p = col[i] , q = colk+mycols[i]-k , t = *q ;
				p < col[i]+n-mycols[i] ; p++ , q++ )
					*p -= *q * t ;
		}
	}
	free(colk) ;
}


forsub ( n, ncols, map, mycols, col, b ,me, nprocs, tids)
	int me, n, ncols, *map, *mycols ;
	float **col, *b ;
	int nprocs, *tids;

/* lower triangular forward solve */
/* column-oriented fan-in algorithm */

{
	int  i, j, k ;
	float t ,t2;

	j = 0 ;
	if (map[0] == me) {
		b[j] /= *col[j] ;
		j++ ;
	}
	for ( k = 1 ; k < n ; k++ ) {
		t = 0 ;
		for ( i = 0 ; i < j ; i++ ) t += b[i] * *(col[i]+k-mycols[i]) ;
		if (map[k] == me) {
			for (i=0;i<(nprocs-1);i++) {
				pvm_recv(-1,k);
				pvm_upkfloat(&t2,1,1);
				t+=t2;
			}
		}
		else {
			pvm_initsend(PvmDataRaw);
			pvm_pkfloat(&t,1,1);
			pvm_send(tids[map[k]],k);
		}

		if (map[k] == me) {
			b[j] -= t ;
			b[j] /= *col[j] ;
			j++ ;
		}
	}
}


backsub ( n, nrows, map, myrows, row, b ,me, nprocs, tids)
	int me, n, nrows, *map, *myrows ;
	float **row, *b ;
	int nprocs, *tids;

/* upper triangular backward solve */
/* row-oriented fan-out algorithm */

{
	int host, dir, i, j, k ;
	float x ;

	j = nrows-1 ;
	for ( k = n-1 ; k > 0 ; k-- ) {
		if (map[k] == me) {
			b[j] /= *row[j] ;
			x = b[j] ;
			j-- ;
		}
		if (map[k] == me) {
			pvm_initsend(PvmDataRaw);
			pvm_pkfloat(&x,1,1);
			pvm_mcast(tids,nprocs,k);
			/* pvm_recv(-1,k);  /* mcast sends to itself */
		}
		else {
			pvm_recv(-1,k);
			pvm_upkfloat(&x,1,1);
		}
		for ( i = j ; i >= 0 ; i-- ) b[i] -= x * *(row[i]+k-myrows[i]) ;
	}
	if (map[0] == me) b[j] /= *row[j] ;
}