//***********************************************************************
//  pde.c
//  author: Matthew Judd and John Karpovich
//          mrj2p@virginia.edu and jfk3w@Virginia.edu
//          University of Virginia
//          December 7, 1992
//
//	     adapted from code by Ambar Sarkar
//
//  compilation:
//      make serial   - for serial version
//      make all      - for MENTAT parallel version
//
//  description:
//      Implements class definitions for pde and related objects
//
//  user responsibility:
//      Write setRowsAndCols() code - this code initializes the number of rows
//	and columns in the image. It also initializes the window variables to 
//	include the whole image.
//
//      Write getNextStencil() code - this function both returns the next stencil
//      to apply to the matrix, and decides whether or not to stop applying
//      filters.
//
//      Write prepareDest() code - this function prepares the destination file
//      to receive the final image from the pdes.  Any information
//      the destination file needs other than that image must be written here.
//***********************************************************************

// ****************************************************************************
//
//              Copyright (c) 1992
//              University of Virginia
//              All Rights Reserved
//
//              This software is the property of the University of Virginia.
//              This software is provided `as is' and without any express or
//              implied warrenties.
//
//              Send problems & suggestions to:
//              Andrew Grimshaw (grimshaw@cs.virginia.edu)
//
// **************************************************************************


#include <memory.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pde.h"




//************************************************************************
//************************************************************************
//
//  		USER DEFINED pde FUNCTIONS
//
//************************************************************************
//************************************************************************



void pde::workers_alloc(int pieces)
{
#ifndef SERIAL
    pde  prototype;   // a convolve handle used to create the
			   // vector of workers needed. 
    int  N;
    
    workers = (Stenciler *) new pde[pieces];
    
    for (N=0; N<pieces; N++)
	workers[N].create(N);
    
#else
    if (pieces >0) 
	Stenciler::workers = (Stenciler *) new pde[pieces];
#endif
}
#ifdef fred
#endif

void pde::destroy_workers(int pieces)
{
    int i;

#ifndef SERIAL
    for (i = 0; i < pieces; i++)
    {
	workers[i].destroy();
    }
    
#else
    delete workers;
#endif
}


//*************************************************************************
// Name of file with source image is sNm
// 
// This code should determine the number of rows and columns in the data 
// matrix, and set the pde member variables numRows and numCols.  
// Finally, it should set the pde member variables lr_row_window 
// and lr_col_window, which designate the lower right point defining the 
// part of the image which is mutable, to (numRows - 1) and (numCols - 1), 
// respectively.
// Examples of this file for use with .ppx format and user-defined formats
// exist with the convolve and pde examples.
//*************************************************************************
int pde::setRowsCols() 
{
    int szread = 100;

    srcFile = fopen((char *)srcNm, "r");
    if (srcFile != NULL) 
    {
	fseek(srcFile,0,0);
	fread((char *) &numRows, sizeof(int), 1, srcFile);
	szread = fread((char *) &numCols, sizeof(int), 1, srcFile);
    }
    else 
    {
	fprintf(stderr, "Cannot open %s\n", srcNm); 
	fflush(stderr);
	exit (-1);
    }
    lr_row_window = numRows - 1;
    lr_col_window = numCols - 1;

    ulRow = 0;
    ulCol = 0;
    lrRow = numRows;
    lrCol = numCols;

    return (szread);
}





//************************************************************************
//  Determines whether or not to continue applying stencils; if so, 
//  determines what the next stencil is and returns a pointer to it.  If 
//  loop is finished, returns a NULL pointer
//************************************************************************
stencil* pde::getNextStencil()
{
    stencil  *next;
    int i;
    float currConvVal, sum, *temp;

    if ((currIter == 0) && (currentStencil != NULL))
    {
	next = currentStencil->st;
	currIter++;
    }
    else
    {
	temp = new float[myPieces];
	sum = 0.0;

	for (i=0; i<myPieces; i++)
	    temp[i] = workers[i].checkConvergence();

        for (i=0; i<myPieces; i++)
	{
#ifdef STEN_DEBUG
	    fprintf(stderr, "Partial sum from piece %d = %f\n", i, temp[i]);
#endif
	    sum += temp[i];
	}
	
	delete [] temp;
	
	currConvVal = sum;

#ifdef STEN_DEBUG
    fprintf(stderr, "Iteration %d: currConvVal = %f\n\n", currIter, currConvVal);
#endif
	
	
	if (currConvVal <= (convVal*convVal))
	    next = NULL;
	else
	{
	    currIter++;
	    next = currentStencil->st;
	}
    }

    return (next);
}



//************************************************************************
//  This function is used to prepare the destination file to accept the
//  data from the workers.  For example, if the chosen file format requires
//  a header, this function should write that header to the output file.
//  After this function has been called, the output file should need nothing
//  but the actual matrix data to be complete.
//************************************************************************

void pde::prepareDest()
{
    char hdr[512];
    char *tmp;
    int i;
    
    destFile = fopen ((char*) destNm, "w");

    // extra header junk for ppx.

    hdr[0] = 0;
    hdr[1] = 35;
    hdr[2] = 1;

    tmp    = (char *) &numRows;
    hdr[3] = tmp[2];
    hdr[4] = tmp[3];

    tmp    = (char *) &numCols;
    hdr[5] = tmp[2];
    hdr[6] = tmp[3];

    strcpy(&hdr[7],"a00661.hex");
    for (i=17; i<512; i++)
	hdr[i] = 0;

    fwrite (hdr, sizeof(char), 512, destFile);

    fclose (srcFile);
    fclose (destFile);

}




//************************************************************************
//************************************************************************
//
//  		USER DEFINED pde FUNCTIONS
//
//************************************************************************
//************************************************************************




float pde::checkConvergence()
{
    float *temp, sum;
    int i;

    sum = 0.0;

    if (workerType == LEAF)
	sum = checkConvergencePiece();

    else
    {
	temp = new float[myPieces];

	sum = 0.0;

	for (i=0; i<myPieces; i++)
	    temp[i] = workers[i].checkConvergence();

        for (i=0; i<myPieces; i++)
	    sum += temp[i];

	delete [] temp;
    }

#ifndef SERIAL
    mentat_return (sum);
#endif

    return (sum);
}


float pde::checkConvergencePiece()
{
    float temp, sum;
    int i, j, totalCols;
    DATA_TYPE *srcPtr, *oldPtr;

    sum = 0.0;
    totalCols = (maxStencilCols * 2) + numCols;

    srcPtr = srcArray->get_r_ptr(maxStencilRows);
    srcPtr += maxStencilCols;
    oldPtr = destArray->get_r_ptr(maxStencilCols);
    oldPtr += maxStencilCols;

    for (i=maxStencilRows; i < numRows + maxStencilRows; i++)
    {
	for (j=0; j<numCols; j++)
	{
	    temp = srcPtr[j] - oldPtr[j];
	    sum += temp * temp;
	}
	srcPtr += totalCols;
	oldPtr += totalCols;
    }

    return (sum);
}


// **************************************************************************
// Output this section of final image to target file.
//
// This function should open the destination file, copy this workers piece
// into that file in the correct place, and then close the output file.
// Examples for both image convolution and pde's are in the example files;
// both output images into files in .ppx format.
// **************************************************************************
void pde::putMatrixPiece()
{
    DATA_TYPE  *rowPtr;
    int        i, j;
    int        hdrSz, offset, colOffset, totalCols;
    int        ulr, lrr;
    DD_chararray* tmpArray;
    char*	charRow;


    // Open destination file - report any error
    destFile = fopen((char*) destNm, "r+");
    if (destFile == NULL)
    {
	fprintf (stderr, "ERROR!  Cannot open destination file %s\n", 
		 (char*) destNm);
	fflush (stderr);
	exit (-1);
    }

    // Initialize necessary variables
    hdrSz     =  512;   // The size of the header
    offset    =  hdrSz + (((matrixCols * ulRow) + ulCol) * sizeof (char));
    colOffset =  matrixCols - numCols;
    totalCols =  (maxStencilCols * 2) + numCols;
    ulr       =  maxStencilRows;	
    lrr       =  ulr + numRows - 1;
    tmpArray  =  new (1,numCols) DD_chararray(1,numCols);
    charRow   =  tmpArray->get_r_ptr(0);
    rowPtr    =  srcArray->get_r_ptr(ulr);
    rowPtr    += maxStencilCols;

    //  		START OUPUTTING PIECE
    // Find beginning of piece
    fseek (destFile, offset, 0);

    for (i = ulr; i <= lrr; i++)
    {
	// copy/convert floats to chars 
	for (j=0;j<numCols;j++)
	    charRow[j] = (char) rowPtr[j];

	if (0 == (j=fwrite (charRow, sizeof (char), numCols, destFile)))
	{
	    fprintf (stderr, "ERROR!  Cannot fwrite to %s\n", (char*) destNm);
	    fflush (stderr);
	    exit (-1);
	}

	// Find start points for next row
	rowPtr += totalCols;    
	fseek (destFile, colOffset, 1);
    }
    fclose (destFile);
    delete tmpArray;
}	

// **************************************************************************
// Apply given stencil to workers section of image.  
//
// This function describes the actual work to be done to each piece of the
// input matrix with each stencil.  Examples for image convolution and 
// partial differential equations are in the example files.
// **************************************************************************
void pde::doStencilPiece(stencil *sten)
{
    DATA_TYPE    *srcr, *destr, *currSrcPtr, *currDstPtr;
    int          ulr, ulc, lrr, lrc;
    int          i, j;
    int 	 totalCols;
    DATA_TYPE	 destval;
    MATRIX_TYPE  *tmpArray, *oldArray;


    currStRows = 1;
    currStCols = 1;
 
    // set region to work on in terms of the source array's coords.
    // window coords are relative to "real" data, i.e 0,0 = top left corner of
    // real data, not source array.  That's why we need to add in maxStencilxxx
    // to calc offset.

    ulr = (w_ulRow - ulRow) + maxStencilRows;
    ulc = (w_ulCol - ulCol) + maxStencilCols;
    lrr = (maxStencilRows + numRows - 1) - (lrRow - w_lrRow);
    lrc = (maxStencilCols + numCols - 1) - (lrCol - w_lrCol);

    totalCols = (maxStencilCols * 2) + numCols;

/*
//		Method I
    for (i=ulr; i <= lrr; i++)
    {
	// point srcr to the beginning of where the data for this stencil will 
	// come from (i.e. to the upper left of the current destination point)

	srcr = srcArray->get_r_ptr(i-currStRows);
	srcr += ulc;
	destr = destArray->get_r_ptr(i);
	destr += ulc;
	
        for (j=ulc; j <= lrc; j++)
	{
	    currSrcPtr = srcr;
	    destval = currSrcPtr[0];
	    destval += currSrcPtr[totalCols*2];
	    destval += currSrcPtr[totalCols-1];
	    destval += currSrcPtr[totalCols+1];
	    *destr = (DATA_TYPE) (destval/4.0);
	    destr++;
	    srcr++;
	}
    }
*/

//		Method II
    srcr = srcArray->get_r_ptr(ulr-currStRows);
    srcr += ulc;
    destr = destArray->get_r_ptr(ulr);
    destr += ulc;
 
    for (i=ulr; i <= lrr; i++)
    {
	// point srcr to the beginning of where the data for this stencil will 
	// come from (i.e. to the upper left of the current destination point)

	currSrcPtr = srcr;
	currDstPtr = destr;
        for (j=ulc; j <= lrc; j++)
	{
	    destval = currSrcPtr[0];
	    destval += currSrcPtr[totalCols*2];
	    destval += currSrcPtr[totalCols-1];
	    destval += currSrcPtr[totalCols+1];
	    *currDstPtr = (DATA_TYPE) (destval/4.0);

	    currSrcPtr++;
	    currDstPtr++;
	}
	srcr += totalCols;
	destr += totalCols;
    }

    oldArray = srcArray->extract_region(maxStencilRows,maxStencilCols,
				 maxStencilRows+numRows,maxStencilCols+numCols);
    tmpArray = destArray->extract_region(ulr,ulc,lrr,lrc);
    srcArray->overlay_region(ulr,ulc,tmpArray);
    destArray->overlay_region(maxStencilRows,maxStencilCols,oldArray);
    delete tmpArray;
    delete oldArray;
}



// **************************************************************************
// Read workers section of Matrix from given input file, using values sent
// from pde through init().  Called from the init() function.
//
// This function must get each workers pice from the input file.  It should
// use pde member variables matrixCols, ulRow, ulCol, and numCols
// to determine the offset from the beginning of the file to the start of the
// first piece row, and the distance in the input file from the end of one
// piece row to the beginning of the next.
// **************************************************************************
void pde::getMatrixPiece()
{
    int       hdrSz, offset, colOffset, totalCols;
    int       ulr, lrr;
    int       i, j;
    DATA_TYPE *rowPtr;
    DATA_TYPE *tmp;

    // Open Source File - report any errors
    srcFile = fopen ((char*) srcNm, "r");
    if (srcFile == NULL)
    {
	fprintf (stderr, "ERROR! Cannot open file %s\n", (char*) srcNm);
	fflush (stderr);
	exit (-1);
    }

    // Initialize necessary variables
    hdrSz     =  2*sizeof(int);   // The size of the header
    offset    =  hdrSz + (((matrixCols * ulRow) + ulCol) * sizeof (DATA_TYPE));
    colOffset =  (matrixCols - numCols)*sizeof(DATA_TYPE);
    totalCols =  (maxStencilCols * 2) + numCols;
    ulr       =  maxStencilRows;	
    lrr       =  ulr + numRows - 1;
    rowPtr    =  srcArray->get_r_ptr(ulr);
    rowPtr    += maxStencilCols;

    // Find beginning of matrix piece in source file
    fseek (srcFile, offset, 0);

    // Start reading in matrix piece
    for (i=ulr; i <= lrr; i++)
    {
	if (0 == (j=fread ((char*)rowPtr, sizeof(DATA_TYPE), numCols, srcFile)))
	{
	    fprintf (stderr, "ERROR! Cannot fread %d\n", (char*)srcNm);
	    fflush (stderr);
	    exit(-1);
	}

	// Find starting points for next row
	rowPtr += totalCols;    
	fseek (srcFile, colOffset, 1);
    }
    fclose (srcFile);

    tmp = srcArray->get_r_ptr(maxStencilRows);
#ifndef SERIAL
    tm.start();
#endif    
}