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

/*
	PVM 3.0 version
	Host program driver for column Cholesky factorization
	of a symmetric matrix partitioned by vertical strips
	USAGE: cholhost <numprocs> <order of matrix> <stripsize>
*/

main(argc,argv)
int argc;
char **argv;
/* host process */
{
	int mytid, tids[100] ;
	float f0, f1[4], err, secs[4], mflops[3], eff[3] ;
	int info[3], *map ;
	int me, p, n, blksize, i, j, nprocs;

	if (argc < 4)  { 
		printf("Usage: cholhost #procs N blksize \n"); 
		exit(0);
	}
	nprocs = atoi(*++argv);
	p = nprocs;
	n = atoi(*++argv);
	if (n < p) {
		printf("Problem size (%d) < nprocs (%d)\n",n,p);
		exit(1);
	}
	blksize = atoi(*++argv);
	if (blksize > n/p) blksize = n/p ;
    if ((mytid = pvm_mytid()) < 0) {
        exit(1);
    }
    printf("i'm t%x\n", mytid);

    /* start up tasks */
    pvm_spawn("cholnode", (char**)0, 0, "", nprocs, tids);

    /* multicast tids array to children */
    pvm_initsend( PvmDataDefault );
    pvm_pkint(&nprocs, 1, 1);
    pvm_pkint(tids, nprocs, 1);
    pvm_mcast(tids, nprocs, 16000);

	pvm_initsend(PvmDataRaw);
	info[0] = n ; 
	info[1] = 0;
	info[2] = p; 
	pvm_pkint(info,3,1);
	if (pvm_mcast(tids,nprocs,16001) < 0) {
		printf("HOST: Broadcast of 16001 failed!\n");
		for (i=0;i<nprocs;i++) pvm_kill(tids[i]);
		pvm_exit();
		exit(1);
	}

	map = (int *)malloc(n*sizeof(int)) ;
	for ( i = 0 ; i < n ; i++ ) map[i] = ((i/blksize)%p);

	pvm_initsend(PvmDataRaw);
	pvm_pkint(map,n,1);
	if (pvm_mcast(tids,nprocs,16002) < 0) {
		printf("HOST: send map (partitioning info) failed!!!\n");
		for (i=0;i<nprocs;i++) pvm_kill(tids[i]);
		pvm_exit();
		exit(1);
	}

	printf("Number of processors = %d, order of matrix = %d\n", p, n) ;
	printf("Cholesky factorization - wrap mapping, natural ordering, ");
	printf("blocksize = %d\n",blksize);
	printf("Full topology, fan-out factorization, synchronous algorithm\n");
	err=0.0;
	
	for (i=0;i<p;i++) {
		pvm_recv(-1,16003);
		pvm_upkfloat(&f0,1,1);
		if (f0 > err) err=f0;
	}
	printf("maxerr\t%e\n", err) ;

	for (i=0;i<4;i++) secs[i]=0.0;
	
	for (i=0;i<p;i++) {
		pvm_recv(-1,16004);
		pvm_upkfloat(f1,4,1); 
		for (j=0;j<4;j++) if (f1[j] > secs[j]) secs[j] = f1[j];
	}

	printf("      \t factor\t fsolve\t bsolve\t total\n") ;
	printf("secs  \t%7.2f\t%7.2f\t%7.2f\t%7.2f\n", secs[0], secs[1], 
		secs[2], secs[3]) ;
	free(map) ;

	pvm_exit();
}