temperature_mobile.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef _MSC_VER
#define strcasecmp    _stricmp
#define strncasecmp   _strnicmp
#else
#include <strings.h>
#endif
#include <math.h>

#include "temperature_mobile.h"
#include "flp.h"
#include "util.h"

/* constructors	*/
blist_t *new_blist(int idx, double occupancy)
{
	blist_t *ptr = (blist_t *) calloc (1, sizeof(blist_t));
	if (!ptr)
		fatal("memory allocation error\n");
	ptr->idx = idx;
	ptr->occupancy = occupancy;
	ptr->next = NULL;
	return ptr;
}

blist_t ***new_b2gmap(int rows, int cols)
{
	int i;
	blist_t ***b2gmap;

	b2gmap = (blist_t ***) calloc (rows, sizeof(blist_t **));
	b2gmap[0] = (blist_t **) calloc (rows * cols, sizeof(blist_t *));
	if (!b2gmap || !b2gmap[0])
		fatal("memory allocation error\n");

	for(i=1; i < rows; i++) 
		b2gmap[i] = b2gmap[0] + cols * i;

	return b2gmap;	
}

/* destructor	*/
void delete_b2gmap(blist_t ***b2gmap, int rows, int cols)
{
	int i, j;
	blist_t *ptr, *temp;

	/* free the linked list	*/
	for(i=0; i < rows; i++)
		for(j=0; j < cols; j++) {
			ptr = b2gmap[i][j];
			while(ptr) {
				temp = ptr->next;
				free(ptr);
				ptr = temp;
			}
		}

	/* free the array space	*/
	free(b2gmap[0]);
	free(b2gmap);
}

/* re-initialize */
void reset_b2gmap(grid_model_t *model, layer_t *layer)
{
	int i, j;
	blist_t *ptr, *temp;

	/* free the linked list	*/
	for(i=0; i < model->rows; i++)
		for(j=0; j < model->cols; j++) {
			ptr = layer->b2gmap[i][j];
			while(ptr) {
				temp = ptr->next;
				free(ptr);
				ptr = temp;
			}
			layer->b2gmap[i][j] = NULL;
		}
}

/* create a linked list node and append it at the end	*/
void blist_append(blist_t *head, int idx, double occupancy)
{
	blist_t *tail = NULL;
	
	if(!head)
		fatal("blist_append called with empty list\n");

	/* traverse till the end	*/
	for(; head; head = head->next)
		tail = head;

	/* append */
	tail->next =  new_blist(idx, occupancy);
}

/* compute the power/temperature average weighted by occupancies	*/
double blist_avg(blist_t *ptr, flp_t *flp, double *v, int type)
{
	double  val = 0.0;
	
	for(; ptr; ptr = ptr->next) {
		if (type == V_POWER)
			val += ptr->occupancy * v[ptr->idx] / (flp->units[ptr->idx].width * 
				   flp->units[ptr->idx].height);
		else if (type == V_TEMP)		   
			val += ptr->occupancy * v[ptr->idx];
		else
			fatal("unknown vector type\n");
	}		

	return val;		   
}

void debug_print_blist(blist_t *head, flp_t *flp);
/* test the block-grid map data structure	*/
void test_b2gmap(grid_model_t *model, layer_t *layer)
{
	int i, j;
	blist_t *ptr;
	double sum;

	/* a correctly formed b2gmap should have the 
	 * sum of occupancies in each linked list
	 * to be equal to 1.0
	 */
	for (i=0; i < model->rows; i++)
		for(j=0; j < model->cols; j++) {
			sum = 0.0;
			for(ptr = layer->b2gmap[i][j]; ptr; ptr = ptr->next)
				sum += ptr->occupancy;
			if (!eq(floor(sum*1e5 + 0.5)/1e5, 1.0)) {
				fprintf(stdout, "i: %d\tj: %d\n", i, j);
				debug_print_blist(layer->b2gmap[i][j], layer->flp);
				fatal("erroneous b2gmap data structure. invalid floorplan?\n");
			}	
		}
}

/* setup the block and grid mapping data structures	*/
void set_bgmap(grid_model_t *model, layer_t *layer)
{
	/* i1, i2, j1 and j2 are indices of the boundary grid cells	*/
	int i, j, u, i1, i2, j1, j2;

	/* shortcuts for cell width(cw) and cell height(ch)	*/
	double cw = model->width / model->cols;
	double ch = model->height / model->rows;

	/* initialize	*/
	reset_b2gmap(model, layer);

	/* for each functional unit	*/
	for(u=0; u < layer->flp->n_units; u++) {
		/* shortcuts for unit boundaries	*/
		double lu = layer->flp->units[u].leftx;
		double ru = lu + layer->flp->units[u].width;
		double bu = layer->flp->units[u].bottomy;
		double tu = bu + layer->flp->units[u].height;

		/* top index (lesser row) = rows - ceil (topy / cell height)	*/
		i1 = model->rows - tolerant_ceil(tu/ch);
		/* bottom index (greater row) = rows - floor (bottomy / cell height)	*/
		i2 = model->rows - tolerant_floor(bu/ch);
		/* left index = floor (leftx / cell width)	*/
		j1 = tolerant_floor(lu/cw);
		/* right index = ceil (rightx / cell width)	*/
		j2 = tolerant_ceil(ru/cw);
		/* sanity check	*/
		if((i1 < 0) || (j1 < 0))
			fatal("negative grid cell start index!\n");
		if((i2 > model->rows) || (j2 > model->cols))
			fatal("grid cell end index out of bounds!\n");
		if((i1 >= i2) || (j1 >= j2))
			fatal("invalid floorplan spec or grid resolution\n");

		/* setup g2bmap	*/
		layer->g2bmap[u].i1 = i1;
		layer->g2bmap[u].i2 = i2;
		layer->g2bmap[u].j1 = j1;
		layer->g2bmap[u].j2 = j2;

		/* setup b2gmap	*/
		/* for each grid cell in this unit	*/
		for(i=i1; i < i2; i++)
			for(j=j1; j < j2; j++)
				/* grid cells fully overlapped by this unit	*/
				if ((i > i1) && (i < i2-1) && (j > j1) && (j < j2-1)) {
					/* first unit in the list	*/
					if (!layer->b2gmap[i][j])
						layer->b2gmap[i][j] = new_blist(u, 1.0);
					else {
					/* this should not occur since the grid cell is 
					 * fully covered and hence, no other unit should 
					 * be sharing it
					 */
						blist_append(layer->b2gmap[i][j], u, 1.0);
						warning("overlap of functional blocks?\n");
					}
				/* boundary grid cells partially overlapped by this unit	*/
				} else {
					/* shortcuts for cell boundaries	*/
					double lc = j * cw, rc = (j+1) * cw;
					double tc = model->height - i * ch;
					double bc = model->height - (i+1) * ch;
					
					/* shortcuts for overlap width and height	*/
					double oh = (MIN(tu, tc) - MAX(bu, bc));
					double ow = (MIN(ru, rc) - MAX(lu, lc));
					double occupancy;
			
					/* overlap tolerance	*/
					if (eq(oh/ch, 0))
						oh = 0;
					else if (eq(oh/ch, 1))
						oh = ch;

					if (eq(ow/cw, 0))
						ow = 0;
					else if (eq(ow/cw, 1))
						ow = cw;

					occupancy = (oh * ow) / (ch * cw);
					if (oh < 0 || ow < 0)
						fatal("negative overlap!\n");

					/* first unit in the list	*/
					if (!layer->b2gmap[i][j])
						layer->b2gmap[i][j] = new_blist(u, occupancy);
					else
					/* append at the end	*/
						blist_append(layer->b2gmap[i][j], u, occupancy);
				}
	}

	/* 
	 * sanity check	
	test_b2gmap(model, layer);
	 */
}

/* populate default set of layers	*/ 
void populate_default_layers(grid_model_t *model, flp_t *flp_default)
{
	/* silicon	*/
	model->layers[LAYER_SI].no = LAYER_SI;
	model->layers[LAYER_SI].has_lateral = TRUE;
	model->layers[LAYER_SI].has_power = TRUE;
	model->layers[LAYER_SI].k = model->config.k_chip;
	model->layers[LAYER_SI].thickness = model->config.t_chip;
	model->layers[LAYER_SI].sp = model->config.p_chip;
	model->layers[LAYER_SI].flp = flp_default;
	model->layers[LAYER_SI].b2gmap = new_b2gmap(model->rows, model->cols);
	model->layers[LAYER_SI].g2bmap = (glist_t *) calloc(flp_default->n_units, sizeof(glist_t));
	if (!model->layers[LAYER_SI].g2bmap)
		fatal("memory allocation error\n");

	/* interface material	*/
	model->layers[LAYER_INT].no = LAYER_INT;
	model->layers[LAYER_INT].has_lateral = TRUE;
	model->layers[LAYER_INT].has_power = FALSE;
	model->layers[LAYER_INT].k = model->config.k_interface;
	model->layers[LAYER_INT].thickness = model->config.t_interface;
	model->layers[LAYER_INT].sp = model->config.p_interface;
	model->layers[LAYER_INT].flp = flp_default;
	model->layers[LAYER_INT].b2gmap = model->layers[LAYER_SI].b2gmap;
	model->layers[LAYER_INT].g2bmap = model->layers[LAYER_SI].g2bmap;
}

/* parse the layer file open for reading	*/
void parse_layer_file(grid_model_t *model, FILE *fp)
{
	char line[LINE_SIZE], *ptr, cval;
	int count, i = 0, field = LCF_SNO, ival;
	double dval;

	fseek(fp, 0, SEEK_SET);
	count = 0;
	while (!feof(fp) && count < (model->n_layers * LCF_NPARAMS)) {
		fgets(line, LINE_SIZE, fp);
		if (feof(fp))
			break;

		/* ignore comments and empty lines	*/
		ptr = strtok(line, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;
			
		switch (field) 
		{
			case LCF_SNO:
						if (sscanf(ptr, "%d", &ival) != 1)
							fatal("invalid layer number\n");
						if(ival >= model->n_layers || ival < 0)
							fatal("layer number must be >= 0 and < no. of layers\n");
						if (model->layers[ival].no != 0)
							fatal("layer numbers must be unique\n");
						i = ival;
						model->layers[i].no = ival;
						break;
			case LCF_LATERAL:
						if (sscanf(ptr, "%c", &cval) != 1)
							fatal("invalid layer heat flow indicator\n");
						if (cval == 'Y' || cval == 'y')
							model->layers[i].has_lateral = TRUE;
						else if (cval == 'N' || cval == 'n')			
							model->layers[i].has_lateral = FALSE;
						else			
							fatal("invalid layer heat flow indicator\n");
						break;
			case LCF_POWER:
						if (sscanf(ptr, "%c", &cval) != 1)
							fatal("invalid layer power dissipation indicator\n");
						if (cval == 'Y' || cval == 'y')
							model->layers[i].has_power = TRUE;
						else if (cval == 'N' || cval == 'n')			
							model->layers[i].has_power = FALSE;
						else			
							fatal("invalid layer power dissipation indicator\n");
						break;
			case LCF_SP:
						if (sscanf(ptr, "%lf", &dval) != 1)
							fatal("invalid specific heat\n");
						model->layers[i].sp = dval;
						break;
			case LCF_RHO:
						if (sscanf(ptr, "%lf", &dval) != 1)
							fatal("invalid resistivity\n");
						model->layers[i].k = 1.0 / dval;
						break;
			case LCF_THICK:
						if (sscanf(ptr, "%lf", &dval) != 1)
							fatal("invalid thickness\n");
						model->layers[i].thickness = dval;
						break;
			case LCF_FLP:
						model->layers[i].flp = read_flp(ptr, FALSE);
						/* first layer	*/
						if (count < LCF_NPARAMS) {
							model->width = get_total_width(model->layers[i].flp);
							model->height = get_total_height(model->layers[i].flp);
						} else if(!eq(model->width, get_total_width(model->layers[i].flp)) || 
								  !eq(model->height, get_total_height(model->layers[i].flp)))
							fatal("width and height differ across layers\n");
						break;
			default:
						fatal("invalid field id\n");
						break;
		}
		field = (field + 1) % LCF_NPARAMS;
		count++;
	}

	/* allocate the block-grid maps */
	for(i=0; i < model->n_layers; i++) {
		model->layers[i].b2gmap = new_b2gmap(model->rows, model->cols);
		model->layers[i].g2bmap = (glist_t *) calloc(model->layers[i].flp->n_units, 
								   sizeof(glist_t));
		if (!model->layers[i].g2bmap)
			fatal("memory allocation error\n");
	}
}

/* populate layer info either from the default floorplan or from
 * the layer configuration file (lcf)
 */
void populate_layers_grid(grid_model_t *model, flp_t *flp_default)
{
	char str[STR_SIZE];
	FILE *fp = NULL;

	/* lcf file specified	*/
	if (model->has_lcf) {
		if (!strcasecmp(model->config.grid_layer_file, "stdin"))
			fp = stdin;
		else
			fp = fopen (model->config.grid_layer_file, "r");
		if (!fp) {
			sprintf(str, "error opening file %s\n", model->config.grid_layer_file);
			fatal(str);
		}
	}

	/* compute the no. of layers	*/
	if (model->has_lcf) {
		model->n_layers = count_significant_lines(fp);
		if (model->n_layers % LCF_NPARAMS)
			fatal("wrong no. of lines in layer file\n");
		model->n_layers /= LCF_NPARAMS;
	/* default no. of layers when lcf file is not specified	*/	
	} else
		model->n_layers = DEFAULT_CHIP_LAYERS;

	/* allocate initial memory */
	model->layers = (layer_t *) calloc (model->n_layers, sizeof(layer_t));
	if (!model->layers)
		fatal("memory allocation error\n");

	/* read in values from the lcf when specified	*/
	if (model->has_lcf) {
		parse_layer_file(model, fp);
		warning("layer configuration file specified. overriding default floorplan with those in lcf file\n");
	/* default set of layers	*/
	} else
		populate_default_layers(model, flp_default);
	
	if (model->has_lcf && fp != stdin)
		fclose(fp);
}

/* constructor	*/ 
grid_model_t *alloc_grid_model(thermal_config_t *config, flp_t *flp_default)
{
	int i;
	grid_model_t *model;

	if (config->grid_rows & (config->grid_rows-1) ||
		config->grid_cols & (config->grid_cols-1))
		fatal("grid rows and columns should both be powers of two\n");

	model = (grid_model_t *) calloc (1, sizeof(grid_model_t));
	if (!model)
		fatal("memory allocation error\n");
	model->config = *config;
	model->rows = config->grid_rows;
	model->cols = config->grid_cols;

	if(!strcasecmp(model->config.grid_map_mode, GRID_AVG_STR))
		model->map_mode = GRID_AVG;
	else if(!strcasecmp(model->config.grid_map_mode, GRID_MIN_STR))
		model->map_mode = GRID_MIN;
	else if(!strcasecmp(model->config.grid_map_mode, GRID_MAX_STR))
		model->map_mode = GRID_MAX;
	else if(!strcasecmp(model->config.grid_map_mode, GRID_CENTER_STR))
		model->map_mode = GRID_CENTER;
	else
		fatal("unknown mapping mode\n");

	/* layer configuration file specified?	*/
	if(strcmp(model->config.grid_layer_file, NULLFILE))
		model->has_lcf = TRUE;
	else {
		model->has_lcf = FALSE;
		model->base_n_units = flp_default->n_units;
	}

	/* get layer information	*/
	populate_layers_grid(model, flp_default);

	/* count the total no. of blocks */
	model->total_n_blocks = 0;
	for(i=0; i < model->n_layers; i++)
		model->total_n_blocks += model->layers[i].flp->n_units;

	/* allocate internal state	*/
	model->last_steady = new_grid_model_vector(model);
	model->last_trans = new_grid_model_vector(model);

	return model;
}

void populate_R_model_grid(grid_model_t *model, flp_t *flp)
{
	int i;
	double cw, ch;

	/* setup the block-grid maps; flp parameter is ignored */
	if (model->has_lcf)
		for(i=0; i < model->n_layers; i++)
			set_bgmap(model, &model->layers[i]);
	/* only the silicon layer has allocated space for the maps. 
	 * all the rest just point to it. so it is sufficient to
	 * setup the block-grid map for the silicon layer alone.
	 * further, for default layer configuration, the `flp' 
	 * parameter should be the same as that of the silicon 
	 * layer. finally, the chip width and height information 
	 * need to be populated for default layer configuration
	 */
	else {
		if (flp != model->layers[LAYER_SI].flp)
			fatal("mismatch between the floorplan and the thermal model\n");
		model->width = get_total_width(flp);
		model->height = get_total_height(flp);
		set_bgmap(model, &model->layers[LAYER_SI]);
	}

	/* shortcuts for cell width(cw) and cell height(ch)	*/
	cw = model->width / model->cols;
	ch = model->height / model->rows;

	/* layer specific resistances	*/
	for(i=0; i < model->n_layers; i++) {
		if (model->layers[i].has_lateral) {
			model->layers[i].rx =  getr(model->layers[i].k, cw, ch * model->layers[i].thickness);
			model->layers[i].ry =  getr(model->layers[i].k, ch, cw * model->layers[i].thickness);
		} else {
			/* positive infinity	*/
			model->layers[i].rx = LARGENUM;
			model->layers[i].ry = LARGENUM;
		}
		model->layers[i].rz =  getr(model->layers[i].k, model->layers[i].thickness, cw * ch);
		/* last layer is connected to ambient. divide r_convec proportional to cell area */
		if (i == model->n_layers - 1)
			model->layers[i].rz += model->config.r_convec * 
								   (model->width * model->height) / (cw * ch);
	}

	/* done	*/
	model->r_ready = TRUE;
}

void populate_C_model_grid(grid_model_t *model, flp_t *flp)
{
	int i;

	/* shortcuts for cell width(cw) and cell height(ch)	*/
	double cw = model->width / model->cols;
	double ch = model->height / model->rows;

	if (!model->r_ready)
		fatal("R model not ready\n");
	if (!model->has_lcf && flp != model->layers[LAYER_SI].flp)
		fatal("different floorplans for R and C models!\n");

	/* layer specific capacitances	*/
	for(i=0; i < model->n_layers; i++)
		model->layers[i].c =  getcap(model->layers[i].sp, model->layers[i].thickness, cw * ch);

	/* last layer is connected to the ambient. 
	 * divide c_convec proportional to cell area 
	 */
	model->layers[model->n_layers-1].c += C_FACTOR * model->config.c_convec * (cw * ch) /
										  (model->width * model->height);

	/* done	*/	
	model->c_ready = TRUE;
}

/* destructor	*/
void delete_grid_model(grid_model_t *model)
{
	int i;

	if (model->has_lcf) 
		for(i=0; i < model->n_layers; i++) {
			delete_b2gmap(model->layers[i].b2gmap, model->rows, model->cols);
			free(model->layers[i].g2bmap);
			free_flp(model->layers[i].flp, FALSE);
		}
	/* only the silicon layer has allocated space for the maps. 
	 * all the rest just point to it. also, its floorplan was
	 * allocated elsewhere. so, we don't need to deallocate those.
	 */
	else {
		delete_b2gmap(model->layers[LAYER_SI].b2gmap, model->rows, model->cols);
		free(model->layers[LAYER_SI].g2bmap);
	}
	
	free_grid_model_vector(model->last_steady);
	free_grid_model_vector(model->last_trans);
	free(model->layers);
	free(model);
}

/* differs from 'dvector()' in that memory for internal nodes is also allocated	*/
double *hotspot_vector_grid(grid_model_t *model)
{
	if (model->total_n_blocks <= 0)
		fatal("total_n_blocks is not greater than zero\n");
	return dvector(model->total_n_blocks);
}

/* copy 'src' to 'dst' except for a window of 'size'
 * elements starting at 'at'. useful in floorplan
 * compaction. can be used only with default layer
 * configuration as all layers should have the same
 * floorplan. incompatible with layer configuration
 * file.
 */
void trim_hotspot_vector_grid(grid_model_t *model, double *dst, double *src, 
						 	  int at, int size)
{
	int i;
	if (model->has_lcf)
		fatal("trim_hotspot_vector_grid called with lcf file\n");
	for (i=0; i < at && i < model->total_n_blocks; i++)
		dst[i] = src[i];
	for(i=at+size; i < model->total_n_blocks; i++)
		dst[i-size] = src[i];
}

/* update the model corresponding to floorplan compaction	*/						 
void resize_thermal_model_grid(grid_model_t *model, int n_units)
{
	int i;

	if (model->has_lcf)
		fatal("resize_thermal_model_grid called with lcf file\n");
	if (n_units > model->base_n_units)
		fatal("resizing grid model to more than the allocated space\n");

	/* count the total no. of blocks again */
	model->total_n_blocks = 0;
	for(i=0; i < model->n_layers; i++)
		model->total_n_blocks += model->layers[i].flp->n_units;

	/* nothing more needs to be done because the only data structure
	 * that is dependent on flp->n_units is g2bmap (others are
	 * dependent on 'grid size' which does not change because
	 * of resizing). g2bmap is a 1-d array and needs no reallocation
	 */
}

/* sets the temperature of a vector 'temp' allocated using 'hotspot_vector'	*/
void set_temp_grid(grid_model_t *model, double *temp, double val)
{
	int i;
	if (model->total_n_blocks <= 0)
		fatal("total_n_blocks is not greater than zero\n");
	for(i=0; i < model->total_n_blocks; i++)
		temp[i] = val;
}

/* dump the top layer grid temperatures of 'temp' onto 'file'	*/
void dump_top_layer_temp_grid(grid_model_t *model, char *file, 
							  grid_model_vector_t *temp)
{
	int i, j;
	char str[STR_SIZE];
	FILE *fp;

	if (!model->r_ready)
		fatal("R model not ready\n");

	if (!strcasecmp(file, "stdout"))
		fp = stdout;
	else if (!strcasecmp(file, "stderr"))
		fp = stderr;
	else 	
		fp = fopen (file, "w");

	if (!fp) {
		sprintf (str,"error: %s could not be opened for writing\n", file);
		fatal(str);
	}

	for(i=0;  i < model->rows; i++)
		for(j=0;  j < model->cols; j++)
			fprintf(fp, "%d\t%.2f\n", i*model->cols+j, temp->cuboid[0][i][j]); 
		
	if(fp != stdout && fp != stderr)
		fclose(fp);	
}

/* dump the steady state grid temperatures of the top layer onto 'file'	*/
void dump_steady_temp_grid (grid_model_t *model, char *file)
{
	/* top layer of the most-recently computed steady state temperature	*/
	dump_top_layer_temp_grid(model, file, model->last_steady);
}

/* dump temperature vector alloced using 'hotspot_vector' to 'file' */ 
void dump_temp_grid(grid_model_t *model, double *temp, char *file)
{
	int i, n, base = 0;
	char str[STR_SIZE];
	FILE *fp;

	if (!strcasecmp(file, "stdout"))
		fp = stdout;
	else if (!strcasecmp(file, "stderr"))
		fp = stderr;
	else 	
		fp = fopen (file, "w");

	if (!fp) {
		sprintf (str,"error: %s could not be opened for writing\n", file);
		fatal(str);
	}

	/* layer temperatures	*/
	for(n=0; n < model->n_layers; n++) {
		/* default set of layers	*/
		if (!model->has_lcf) {
			switch(n)
			{
				case LAYER_SI:
						strcpy(str,"");
						break;
				case LAYER_INT:
						strcpy(str,"iface_");
						break;
				default:
						fatal("unknown layer\n");
						break;
			}
		/* layer configuration file	*/
		} else {
			sprintf(str,"layer_%d_", n);
		}

		for(i=0; i < model->layers[n].flp->n_units; i++)
			fprintf(fp, "%s%s\t%.2f\n", str, 
					model->layers[n].flp->units[i].name, temp[base+i]);
		base += model->layers[n].flp->n_units;	
	}

	if (base != model->total_n_blocks)
		fatal("total_n_blocks failed to tally\n");

	if(fp != stdout && fp != stderr)
		fclose(fp);	
}

void copy_temp_grid(grid_model_t *model, double *dst, double *src)
{
	copy_dvector(dst, src, model->total_n_blocks);
}

/* 
 * read temperature vector alloced using 'hotspot_vector' from 'file'
 * which was dumped using 'dump_temp'. values are clipped to thermal
 * threshold based on 'clip'
 */ 
void read_temp_grid(grid_model_t *model, double *temp, char *file, int clip)
{
	int i, n, idx, base = 0;
	double max=0, val;
	char *ptr, str1[LINE_SIZE], str2[LINE_SIZE];
	char name[STR_SIZE], format[STR_SIZE];
	FILE *fp;

	if (!strcasecmp(file, "stdin"))
		fp = stdin;
	else
		fp = fopen (file, "r");

	if (!fp) {
		sprintf (str1,"error: %s could not be opened for reading\n", file);
		fatal(str1);
	}	

	/* temperatures of the different layers	*/
	for (n=0; n < model->n_layers; n++) {
		/* default set of layers	*/
		if (!model->has_lcf) {
			switch(n)
			{
				case LAYER_SI:
						strcpy(format,"%s%lf");
						break;
				case LAYER_INT:
						strcpy(format,"iface_%s%lf");
						break;
				default:
						fatal("unknown layer\n");
						break;
			}
		/* layer configuration file	*/
		} else {
			sprintf(format,"layer_%d_%%s%%lf", n);
		}

		for (i=0; i < model->layers[n].flp->n_units; i++) {
			fgets(str1, LINE_SIZE, fp);
			if (feof(fp))
				fatal("not enough lines in temperature file\n");
			strcpy(str2, str1);
			/* ignore comments and empty lines	*/
			ptr = strtok(str1, " \r\t\n");
			if (!ptr || ptr[0] == '#') {
				i--;
				continue;
			}
			if (sscanf(str2, format, name, &val) != 2)
				fatal("invalid temperature file format\n");
			idx = get_blk_index(model->layers[n].flp, name);
			if (idx >= 0)
				temp[base+idx] = val;
			else	/* since get_blk_index calls fatal, the line below cannot be reached	*/
				fatal ("unit in temperature file not found in floorplan\n");

			/* find max temp on the top layer 
			 * (silicon for the default set of layers)
			 */
			if (n == 0 && temp[idx] > max)
				max = temp[idx];
		}
		base += model->layers[n].flp->n_units;
	}

	if (base != model->total_n_blocks)
		fatal("total_n_blocks failed to tally\n");

	fgets(str1, LINE_SIZE, fp);
	if (!feof(fp))
		fatal("too many lines in temperature file\n");

	if(fp != stdin)
		fclose(fp);	

	/* clipping	*/
	if (clip && (max > model->config.thermal_threshold)) {
		/* if max has to be brought down to thermal_threshold, 
		 * (w.r.t the ambient) what is the scale down factor?
		 */
		double factor = (model->config.thermal_threshold - model->config.ambient) / 
						(max - model->config.ambient);
	
		/* scale down all temperature differences (from ambient) by the same factor	*/
		for (i=0; i < model->total_n_blocks; i++)
			temp[i] = (temp[i]-model->config.ambient)*factor + model->config.ambient;
	}
}

/* dump power numbers to file	*/
void dump_power_grid(grid_model_t *model, double *power, char *file)
{
	int i, n, base = 0;
	char str[STR_SIZE];
	FILE *fp;

	if (!strcasecmp(file, "stdout"))
		fp = stdout;
	else if (!strcasecmp(file, "stderr"))
		fp = stderr;
	else 	
		fp = fopen (file, "w");
	if (!fp) {
		sprintf (str,"error: %s could not be opened for writing\n", file);
		fatal(str);
	}

	/* dump values only for the layers dissipating power	*/
	for(n=0; n < model->n_layers; n++) {
		if (model->layers[n].has_power) {
			for(i=0; i < model->layers[n].flp->n_units; i++)
				if (model->has_lcf)
					fprintf(fp, "layer_%d_%s\t%.3f\n", n, 
							model->layers[n].flp->units[i].name, power[base+i]);
				else 
					fprintf(fp, "%s\t%.3f\n", 
							model->layers[n].flp->units[i].name, power[base+i]);
		}				
		base += model->layers[n].flp->n_units;		
	}

	if(fp != stdout && fp != stderr)
		fclose(fp);	
}

/* 
 * read power vector alloced using 'hotspot_vector' from 'file'
 * which was dumped using 'dump_power'. 
 */ 
void read_power_grid (grid_model_t *model, double *power, char *file)
{
	int i, idx, n, base = 0;
	double val;
	char *ptr, str1[LINE_SIZE], str2[LINE_SIZE]; 
	char name[STR_SIZE], format[STR_SIZE];
	FILE *fp;

	if (!strcasecmp(file, "stdin"))
		fp = stdin;
	else
		fp = fopen (file, "r");
	if (!fp) {
		sprintf (str1,"error: %s could not be opened for reading\n", file);
		fatal(str1);
	}

	/* lcf file could potentially specify more than one power dissipating 
	 * layer. hence, units with zero power within a layer cannot be left
	 * out in the power file.
	 */
	if (model->has_lcf) {
		for(n=0; n < model->n_layers; n++) {
			if (model->layers[n].has_power)
				for(i=0; i < model->layers[n].flp->n_units; i++) {
					fgets(str1, LINE_SIZE, fp);
					if (feof(fp))
						fatal("not enough lines in power file\n");
					strcpy(str2, str1);

					/* ignore comments and empty lines	*/
					ptr = strtok(str1, " \r\t\n");
					if (!ptr || ptr[0] == '#') {
						i--;
						continue;
					}

					sprintf(format, "layer_%d_%%s%%lf", n);
		  			if (sscanf(str2, format, name, &val) != 2)
						fatal("invalid power file format\n");
		  			idx = get_blk_index(model->layers[n].flp, name);
		  			if (idx >= 0)
						power[base+idx] = val;
					/* since get_blk_index calls fatal, the line below cannot be reached	*/
		  			else
						fatal ("unit in power file not found in floorplan\n");
				}
			base += model->layers[n].flp->n_units;	
		}
		fgets(str1, LINE_SIZE, fp);
		if (!feof(fp))
			fatal("too many lines in power file\n");
	/* default layer configuration. so only one layer
	 * has power dissipation. units with zero power 
	 * can be omitted in the power file
	 */
	} else {
		while(!feof(fp)) {
			fgets(str1, LINE_SIZE, fp);
			if (feof(fp))
				break;
			strcpy(str2, str1);
	
			/* ignore comments and empty lines	*/
			ptr = strtok(str1, " \r\t\n");
			if (!ptr || ptr[0] == '#')
				continue;
	
			if (sscanf(str2, "%s%lf", name, &val) != 2)
				fatal("invalid power file format\n");
			idx = get_blk_index(model->layers[LAYER_SI].flp, name);
			if (idx >= 0)
				power[idx] = val;
			else	/* since get_blk_index calls fatal, the line below cannot be reached	*/
				fatal ("unit in power file not found in floorplan\n");
		}
	}

	if(fp != stdin)
		fclose(fp);
}

double find_max_temp_grid(grid_model_t *model, double *temp)
{
	int i;
	double max = 0.0;
	/* max temperature occurs on the top-most layer	*/
	for(i=0; i < model->layers[0].flp->n_units; i++) {
		if (temp[i] < 0)
			fatal("negative temperature!\n");
		else if (max < temp[i])
			max = temp[i];
	}

	return max;
}

double find_avg_temp_grid(grid_model_t *model, double *temp)
{
	int i, n, base = 0, count = 0; 
	double sum = 0.0;
	/* average temperature of all the power dissipating blocks	*/
	for(n=0; n < model->n_layers; n++) {
		if (model->layers[n].has_power) {
			for(i=0; i < model->layers[n].flp->n_units; i++) {
				if (temp[base+i] < 0)
					fatal("negative temperature!\n");
				else 
					sum += temp[base+i];
			}
			count += model->layers[n].flp->n_units;
		}	
		base += model->layers[n].flp->n_units;
	}

	if (!count)
		fatal("no power dissipating units?!\n");
	return (sum / count);
}

/* grid_model_vector routines	*/

/* constructor	*/
grid_model_vector_t *new_grid_model_vector(grid_model_t *model)
{
	grid_model_vector_t *v;

	v = (grid_model_vector_t *) calloc (1, sizeof(grid_model_vector_t));
	if (!v)
		fatal("memory allocation error\n");

	v->cuboid = dcuboid_tail(model->rows, model->cols, model->n_layers, 0);
	return v;
}

/* destructor	*/
void free_grid_model_vector(grid_model_vector_t *v)
{
	free_dcuboid(v->cuboid);
	free(v);
}

/* translate power/temperature between block and grid vectors	*/
void xlate_vector_b2g(grid_model_t *model, double *b, grid_model_vector_t *g, int type)
{
	int i, j, n, base = 0;
	double area;

	/* area of a single grid cell	*/
	area = (model->width * model->height) / (model->cols * model->rows);

	for(n=0; n < model->n_layers; n++) {
		for(i=0; i < model->rows; i++)
			for(j=0; j < model->cols; j++) {
				/* for each grid cell, the power density / temperature are 
				 * the average of the power densities / temperatures of the 
				 * blocks in it weighted by their occupancies
				 */
				/* convert power density to power	*/ 
				if (type == V_POWER)
					g->cuboid[n][i][j] = blist_avg(model->layers[n].b2gmap[i][j], 
										 model->layers[n].flp, &b[base], type) * area;
				/* no conversion necessary for temperature	*/ 
				else if (type == V_TEMP)
					g->cuboid[n][i][j] = blist_avg(model->layers[n].b2gmap[i][j], 
										 model->layers[n].flp, &b[base], type);
				else
					fatal("unknown vector type\n");
			}
		/* keep track of the beginning address of this layer in the 
		 * block power vector
		 */
		base += model->layers[n].flp->n_units;							 
	}

}

/* translate temperature between grid and block vectors	*/
void xlate_temp_g2b(grid_model_t *model, double *b, grid_model_vector_t *g)
{
	int i, j, n, u, base = 0, count;
	int i1, j1, i2, j2, ci1, cj1, ci2, cj2;
	double min, max, avg;

	for(n=0; n < model->n_layers; n++) {
		for(u=0; u < model->layers[n].flp->n_units; u++) {
			/* extent of this unit in grid cell units	*/
			i1 = model->layers[n].g2bmap[u].i1;
			j1 = model->layers[n].g2bmap[u].j1;
			i2 = model->layers[n].g2bmap[u].i2;
			j2 = model->layers[n].g2bmap[u].j2;

			/* map the center grid cell's temperature to the block	*/
			if (model->map_mode == GRID_CENTER) {
				/* center co-ordinates	*/	
				ci1 = (i1 + i2) / 2;
				cj1 = (j1 + j2) / 2;
				/* in case of even no. of cells, center 
				 * is the average of two central cells
				 */
				/* ci2 = ci1-1 when even, ci1 otherwise	*/  
				ci2 = ci1 - !((i2-i1) % 2);
				/* cj2 = cj1-1 when even, cj1 otherwise	*/  
				cj2 = cj1 - !((j2-j1) % 2);

				b[base+u] = (g->cuboid[n][ci1][cj1] + g->cuboid[n][ci2][cj1] + 
						    g->cuboid[n][ci1][cj2] + g->cuboid[n][ci2][cj2]) / 4;
				continue;
			}
		
			/* find the min/max/avg temperatures of the 
			 * grid cells in this block
			 */
			avg = 0.0;
			count = 0;
			min = max = g->cuboid[n][i1][j1];
			for(i=i1; i < i2; i++)
				for(j=j1; j < j2; j++) {
					avg += g->cuboid[n][i][j];
					if (g->cuboid[n][i][j] < min)
						min = g->cuboid[n][i][j];
					if (g->cuboid[n][i][j] > max)
						max = g->cuboid[n][i][j];
					count++;
				}

			/* map to output accordingly	*/
			switch (model->map_mode)
			{
				case GRID_AVG:
					b[base+u] = avg / count;
					break;
				case GRID_MIN:
					b[base+u] = min;
					break;
				case GRID_MAX:
					b[base+u] = max;
					break;
				/* taken care of already	*/	
				case GRID_CENTER:
					break;
				default:
					fatal("unknown mapping mode\n");
					break;
			}
		}
		/* keep track of the beginning address of this layer in the 
		 * block power vector
		 */
		base += model->layers[n].flp->n_units;							 
	}

}

/* set up initial temperatures for the steady state solution
 * heuristically (ignoring the lateral resistances)
 */
void set_heuristic_temp(grid_model_t *model, grid_model_vector_t *power, 
						grid_model_vector_t *temp)
{
	int n, i, j, nl, nr, nc;
	double **sum;

	/* shortcuts	*/
	nl = model->n_layers;
	nr = model->rows;
	nc = model->cols;

	/* layer temperatures	*/
	/* add up power for each grid cell across all layers */
	sum = dmatrix(nr, nc);
	for(n=0; n < nl; n++)
		scaleadd_dvector(sum[0], sum[0], power->cuboid[n][0], nr*nc, 1.0);

	/* last layer	*/
	for(i=0; i < nr; i++)
		for(j=0; j < nc; j++)
			temp->cuboid[nl-1][i][j] = model->config.ambient + sum[i][j] * 
									 model->layers[nl-1].rz;
	/* subtract away the layer's power	*/			
	scaleadd_dvector(sum[0], sum[0], power->cuboid[nl-1][0], nr*nc, -1.0);
			
	/* go from last-1 to first	*/
	for(n=nl-2; n >= 0; n--) {
		/* nth layer temp is n+1th temp + cumul_power * rz of the nth layer	*/
		scaleadd_dvector(temp->cuboid[n][0], temp->cuboid[n+1][0], sum[0], 
						 nr*nc, model->layers[n].rz);
		/* subtract away the layer's power	*/			
		scaleadd_dvector(sum[0], sum[0], power->cuboid[n][0], nr*nc, -1.0);
	}
	
	free_dmatrix(sum);
}

/* macros for calculating conductances	*/
/* conductance to the next cell north. zero if on northern boundary	*/
# define NC(l,n,i,j,nl,nr,nc)		((i > 0) ? (1.0/l[n].ry) : 0.0)
/* conductance to the next cell south. zero if on southern boundary	*/
# define SC(l,n,i,j,nl,nr,nc)		((i < nr-1) ? (1.0/l[n].ry) : 0.0)
/* conductance to the next cell east. zero if on eastern boundary	*/
# define EC(l,n,i,j,nl,nr,nc)		((j < nc-1) ? (1.0/l[n].rx) : 0.0)
/* conductance to the next cell west. zero if on western boundary	*/
# define WC(l,n,i,j,nl,nr,nc)		((j > 0) ? (1.0/l[n].rx) : 0.0)
/* conductance to the next cell below. zero if on bottom face		*/
# define BC(l,n,i,j,nl,nr,nc)		((n < nl-1) ? (1.0/l[n].rz) : 0.0)
/* conductance to the next cell above. zero if on top face			*/
# define AC(l,n,i,j,nl,nr,nc)		((n > 0) ? (1.0/l[n-1].rz) : 0.0)

/* macros for calculating weighted temperatures	*/
/* weighted T of the next cell north. zero if on northern boundary	*/
# define NT(l,v,n,i,j,nl,nr,nc)		((i > 0) ? (v[n][i-1][j]/l[n].ry) : 0.0)
/* weighted T of the next cell south. zero if on southern boundary	*/
# define ST(l,v,n,i,j,nl,nr,nc)		((i < nr-1) ? (v[n][i+1][j]/l[n].ry) : 0.0)
/* weighted T of the next cell east. zero if on eastern boundary	*/
# define ET(l,v,n,i,j,nl,nr,nc)		((j < nc-1) ? (v[n][i][j+1]/l[n].rx) : 0.0)
/* weighted T of the next cell west. zero if on western boundary	*/
# define WT(l,v,n,i,j,nl,nr,nc)		((j > 0) ? (v[n][i][j-1]/l[n].rx) : 0.0)
/* weighted T of the next cell below. zero if on bottom face		*/
# define BT(l,v,n,i,j,nl,nr,nc)		((n < nl-1) ? (v[n+1][i][j]/l[n].rz) : 0.0)
/* weighted T of the next cell above. zero if on top face			*/
# define AT(l,v,n,i,j,nl,nr,nc)		((n > 0) ? (v[n-1][i][j]/l[n-1].rz) : 0.0)

/* single steady state iteration of grid solver - silicon part */
double single_iteration_steady_grid(grid_model_t *model, grid_model_vector_t *power,
									grid_model_vector_t *temp)
{
	int n, i, j;
	double prev, delta, max = 0;
	/* sum of the conductances	*/
	double csum;
	/* weighted sum of temperatures	*/
	double wsum;

	/* shortcuts	*/
	double ***v = temp->cuboid;
	thermal_config_t *c = &model->config;
	layer_t *l = model->layers;
	int nl = model->n_layers;
	int nr = model->rows;
	int nc = model->cols;

	/* for each grid cell	*/
	for(n=0; n < nl; n++)
		for(i=0; i < nr; i++)
			for(j=0; j < nc; j++) {
				/* sum the conductances to cells north, south, 
				 * east, west, above and below
				 */
				csum = NC(l,n,i,j,nl,nr,nc) + SC(l,n,i,j,nl,nr,nc) + 
					   EC(l,n,i,j,nl,nr,nc) + WC(l,n,i,j,nl,nr,nc) + 
					   AC(l,n,i,j,nl,nr,nc) + BC(l,n,i,j,nl,nr,nc);

				/* sum of the weighted temperatures of all the neighbours	*/
				wsum = NT(l,v,n,i,j,nl,nr,nc) + ST(l,v,n,i,j,nl,nr,nc) + 
					   ET(l,v,n,i,j,nl,nr,nc) + WT(l,v,n,i,j,nl,nr,nc) + 
					   AT(l,v,n,i,j,nl,nr,nc) + BT(l,v,n,i,j,nl,nr,nc);

				/* last layer */
				if (n == nl-1) {
					/* all nodes are connected to the ambient	*/
					csum += 1.0/l[n].rz;
					wsum += c->ambient/l[n].rz;
				}

				/* update the current cell's temperature	*/	   
				prev = v[n][i][j];
				v[n][i][j] = (power->cuboid[n][i][j] + wsum) / csum;

				/* compute maximum delta	*/
				delta =  fabs(prev - v[n][i][j]);
				if (delta > max)
					max = delta;
			}
	return max;
}

/* restriction operator for multigrid solver. given a power vector
 * corresponding to a fine grid, outputs a vector corresponding
 * to one level coarser grid (half the no. of rows and cols)
 * NOTE: model->rows and model->cols denote the size of the
 * coarser grid
 */
void multigrid_restrict_power(grid_model_t *model, grid_model_vector_t *dst, 
							  grid_model_vector_t *src)
{
	/* coarse grid indices	*/
	int n, i, j;

	/* grid cells - add the four nearest neighbours	*/
	for(n=0; n < model->n_layers; n++)
		for(i=0; i < model->rows; i++)
			for(j=0; j < model->cols; j++)
			{
				dst->cuboid[n][i][j] = (src->cuboid[n][2*i][2*j] +
										src->cuboid[n][2*i+1][2*j] +
										src->cuboid[n][2*i][2*j+1] +
										src->cuboid[n][2*i+1][2*j+1]);
			}
}

/* prolongation(interpolation) operator for multigrid solver. 
 * given a temperature vector corresponding to a coarse grid,
 * outputs a (bi)linearly interpolated vector corresponding 
 * to one level finer grid (twice the no. of rows and cols)
 * NOTE: model->rows and model->cols denote the size of the
 * coarser grid
 */
void multigrid_prolong_temp(grid_model_t *model, grid_model_vector_t *dst, 
							grid_model_vector_t *src)
{
	/* coarse grid indices	*/
	int n, i, j;

	/* shortcuts	*/
	int nr = model->rows;
	int nc = model->cols;
	double ***d = dst->cuboid;
	double ***s = src->cuboid;

	/* For the fine grid cells not on the boundary,
	 * we want to linearly interpolate in the region 
	 * surrounded by the coarse grid cells (i,j),
	 * (i+1,j) (i,j+1) and (i+1, j+1). The fine
	 * grid cells in that region are (2i+1, 2j+1),
	 * (2i+2, 2j+1), (2i+1, 2j+2) and (2i+2, 2j+2).
	 * To interpolate bilinearly, we first interpolate
	 * along one axis and then along the other. In
	 * each axis, due to the proximity of the fine
	 * grid cell co-ordinates to one of the coarse
	 * grid cells, the weights are 3/4 and 1/4 (not
	 * 1/2 and 1/2) So, repeating them along the  
	 * other axis also results in weights of 9/16,
	 * 3/16, 3/16 and 1/16
	 */ 
	for(n=0; n < model->n_layers; n++)
		for(i=0; i < nr-1; i++)
			for(j=0; j < nc-1; j++) {
				d[n][2*i+1][2*j+1] = 9.0/16.0 * s[n][i][j] +
									 3.0/16.0 * s[n][i+1][j] +
									 1.0/16.0 * s[n][i+1][j+1] +
									 3.0/16.0 * s[n][i][j+1];
				d[n][2*i+2][2*j+1] = 3.0/16.0 * s[n][i][j] +
									 9.0/16.0 * s[n][i+1][j] +
									 3.0/16.0 * s[n][i+1][j+1] +
									 1.0/16.0 * s[n][i][j+1];
				d[n][2*i+2][2*j+2] = 1.0/16.0 * s[n][i][j] +
									 3.0/16.0 * s[n][i+1][j] +
									 9.0/16.0 * s[n][i+1][j+1] +
									 3.0/16.0 * s[n][i][j+1];
				d[n][2*i+1][2*j+2] = 3.0/16.0 * s[n][i][j] +
									 1.0/16.0 * s[n][i+1][j] +
									 3.0/16.0 * s[n][i+1][j+1] +
									 9.0/16.0 * s[n][i][j+1];
			}

	/* for the cells on the boundary, we perform
	 * a zeroth order interpolation. i.e., copy
	 * the nearest coarse grid cell's value
	 * as it is
	 */
	for(n=0; n < model->n_layers; n++) {
		for(i=0; i < nr; i++) {
			d[n][2*i][0] = d[n][2*i+1][0] = s[n][i][0];
			d[n][2*i][2*nc-1] = d[n][2*i+1][2*nc-1] = s[n][i][nc-1];
		}	
		for(j=0; j < nc; j++) {
			d[n][0][2*j] = d[n][0][2*j+1] = s[n][0][j];
			d[n][2*nr-1][2*j] = d[n][2*nr-1][2*j+1] = s[n][nr-1][j];
		}
	}
}

/* recursive multigrid solver. it uses the Gauss-Seidel (GS) iterative 
 * solver to solve at a particular grid granularity. Although GS removes
 * high frequency errors in a solution estimate pretty quickly, it 
 * takes a lot of time to eliminate the low frequency ones. These 
 * low frequency errors can be eliminated easily by coarsifying the
 * grid. This is the core principle of mulrigrid solvers. The version here
 * is called "nested iteration". It solves the problem (iteratively using GS)
 * first in the coarsest granularity and then utilizes that solution to 
 * estimate the solution in the next finer grid. This is repeated until
 * the solution is found in the finest desired grid. For details, take 
 * a look at Numerical Recipes in C (2nd edition), sections 19.5 and 19.6
 * (http://www.nrbook.com/a/bookcpdf/c19-5.pdf and
 * http://www.nrbook.com/a/bookcpdf/c19-6.pdf). Also refer to Prof.
 * Jayathi Murthy's ME 608 notes from Purdue - Chapter 8, sections 8.7-8.9.
 * (http://meweb.ecn.purdue.edu/~jmurthy/me608/main.pdf - pg. 175-192)
 */ 
void recursive_multigrid(grid_model_t *model, grid_model_vector_t *power,
						 grid_model_vector_t *temp)
{
	double delta;
	#if VERBOSE > 1
	unsigned int i = 0;
	#endif
	grid_model_vector_t *coarse_power, *coarse_temp;
	int n;

	/* setup heuristic initial temperatures	at the coarsest level*/
	if (model->rows <= 1 || model->cols <= 1) {
		set_heuristic_temp(model, power, temp);
	
	/* for finer grids. use coarser solutions as estimates	*/
	} else {
		/* make the grid coarser	*/
		model->rows /= 2;
		model->cols /= 2;
		for(n=0; n < model->n_layers; n++) {
			/* only rz's and c's change. rx's and 
			 * ry's remain the same	
			 */
			model->layers[n].rz /= 4;
			if (model->c_ready)
				model->layers[n].c *= 4;
		}

		/* vectors for the coarse grid	*/
		coarse_power = new_grid_model_vector(model);
		coarse_temp = new_grid_model_vector(model);

		/* coarsen the power vector	*/
		multigrid_restrict_power(model, coarse_power, power);

		/* solve recursively	*/
		recursive_multigrid(model, coarse_power, coarse_temp);

		/* interpolate the solution to the current fine grid	*/
		multigrid_prolong_temp(model, temp, coarse_temp);

		/* cleanup	*/
		free_grid_model_vector(coarse_power);
		free_grid_model_vector(coarse_temp);
		
		/* restore the grid */
		model->rows *= 2;
		model->cols *= 2;
		for(n=0; n < model->n_layers; n++) {
			model->layers[n].rz *= 4;
			if (model->c_ready)
				model->layers[n].c /= 4;
		}
	}
	/* refine solution iteratively till convergence	*/
	do {
		delta = single_iteration_steady_grid(model, power, temp);
		#if VERBOSE > 1
			i++;
		#endif
	} while (!eq(delta, 0));
	#if VERBOSE > 1
	fprintf(stdout, "no. of iterations for steady state convergence (%d x %d grid): %d\n", 
					model->rows, model->cols, i);
	#endif
}

void steady_state_temp_grid(grid_model_t *model, double *power, double *temp)
{
	grid_model_vector_t *p;

	if (!model->r_ready)
		fatal("R model not ready\n");

	p = new_grid_model_vector(model);

	/* map the block power numbers to the grid	*/
	xlate_vector_b2g(model, power, p, V_POWER);

	/* solve recursively. use grid model's internal 
	 * state vector to store the grid temperatures
	 */ 
	recursive_multigrid(model, p, model->last_steady);

	/* map the temperature numbers back	*/
	xlate_temp_g2b(model, temp, model->last_steady);

	free_grid_model_vector(p);
}

/* function to access a 1-d array as a 3-d matrix	*/
#define A3D(array,n,i,j,nl,nr,nc)		(array[(n)*(nr)*(nc) + (i)*(nc) + (j)])

/* macros for calculating currents(power values)	*/
/* current(power) from the next cell north. zero if on northern boundary	*/
# define NP(l,v,n,i,j,nl,nr,nc)		((i > 0) ? ((A3D(v,n,i-1,j,nl,nr,nc)-A3D(v,n,i,j,nl,nr,nc))/l[n].ry) : 0.0)
/* current(power) from the next cell south. zero if on southern boundary	*/
# define SP(l,v,n,i,j,nl,nr,nc)		((i < nr-1) ? ((A3D(v,n,i+1,j,nl,nr,nc)-A3D(v,n,i,j,nl,nr,nc))/l[n].ry) : 0.0)
/* current(power) from the next cell east. zero if on eastern boundary	*/
# define EP(l,v,n,i,j,nl,nr,nc)		((j < nc-1) ? ((A3D(v,n,i,j+1,nl,nr,nc)-A3D(v,n,i,j,nl,nr,nc))/l[n].rx) : 0.0)
/* current(power) from the next cell west. zero if on western boundary	*/
# define WP(l,v,n,i,j,nl,nr,nc)		((j > 0) ? ((A3D(v,n,i,j-1,nl,nr,nc)-A3D(v,n,i,j,nl,nr,nc))/l[n].rx) : 0.0)
/* current(power) from the next cell below. zero if on bottom face		*/
# define BP(l,v,n,i,j,nl,nr,nc)		((n < nl-1) ? ((A3D(v,n+1,i,j,nl,nr,nc)-A3D(v,n,i,j,nl,nr,nc))/l[n].rz) : 0.0)
/* current(power) from the next cell above. zero if on top face			*/
# define AP(l,v,n,i,j,nl,nr,nc)		((n > 0) ? ((A3D(v,n-1,i,j,nl,nr,nc)-A3D(v,n,i,j,nl,nr,nc))/l[n-1].rz) : 0.0)

/* compute the slope vector for the grid cells. the transient
 * equation is CdV + sum{(T - Ti)/Ri} = P 
 * so, slope = dV = [P + sum{(Ti-T)/Ri}]/C
 */
void slope_fn_grid(grid_model_t *model, double *v, grid_model_vector_t *p, double *dv)
{
	int n, i, j;
	/* sum of the currents(power values)	*/
	double psum;

	/* shortcuts	*/
	thermal_config_t *c = &model->config;
	layer_t *l = model->layers;
	int nl = model->n_layers;
	int nr = model->rows;
	int nc = model->cols;

	/* for each grid cell	*/
	for(n=0; n < nl; n++)
		for(i=0; i < nr; i++)
			for(j=0; j < nc; j++) {
				/* sum the currents(power values) to cells north, south, 
				 * east, west, above and below
				 */
				psum = NP(l,v,n,i,j,nl,nr,nc) + SP(l,v,n,i,j,nl,nr,nc) + 
					   EP(l,v,n,i,j,nl,nr,nc) + WP(l,v,n,i,j,nl,nr,nc) + 
					   AP(l,v,n,i,j,nl,nr,nc) + BP(l,v,n,i,j,nl,nr,nc);

				/* last layer	*/
				if (n == nl-1) {
					/* all nodes are connected to the ambient	*/
					psum += (c->ambient - A3D(v,n,i,j,nl,nr,nc))/l[n].rz;
				}

				/* update the current cell's temperature	*/	   
				A3D(dv,n,i,j,nl,nr,nc) = (p->cuboid[n][i][j] + psum) / l[n].c;
			}
}

void compute_temp_grid(grid_model_t *model, double *power, double *temp, double time_elapsed)
{
	double t, h, new_h;
	grid_model_vector_t *p;
	#if VERBOSE > 1
	unsigned int i = 0;
	#endif
	
	if (!model->r_ready || !model->c_ready)
		fatal("grid model not ready\n");

	p = new_grid_model_vector(model);

	/* map the block power/temp numbers to the grid	*/
	xlate_vector_b2g(model, power, p, V_POWER);

	/* if temp is NULL, re-use the temperature from the
	 * last call. otherwise, translate afresh and remember 
	 * the grid and block temperature arrays for future use
	 */
	if (temp != NULL) {
		xlate_vector_b2g(model, temp, model->last_trans, V_TEMP);
		model->last_temp = temp;
	}	

	/* Obtain temp at time (t+time_elapsed). 
	 * Instead of getting the temperature at t+time_elapsed directly, we
	 * do it in multiple steps with the correct step size at each time 
	 * provided by rk4. 
	 */
	for (t = 0, new_h = MIN_STEP; t + new_h <= time_elapsed; t+=h) {
		h = new_h;
		/* pass the entire grid as a 1-d array */
		new_h = rk4(model, model->last_trans->cuboid[0][0],  p, 
				 /* array size	*/
				 model->rows * model->cols * model->n_layers, h,
				 model->last_trans->cuboid[0][0], 
				 /* the slope function callback is typecast accordingly */
				 (slope_fn_ptr) slope_fn_grid);
		#if VERBOSE > 1
			i++;
		#endif	
	}

	/* remainder	*/
	if (time_elapsed > t)
		rk4(model, model->last_trans->cuboid[0][0],  p, 
			model->rows * model->cols * model->n_layers,
			/* the slope function callback is typecast accordingly */
			time_elapsed - t, model->last_trans->cuboid[0][0], 
			(slope_fn_ptr) slope_fn_grid);

	#if VERBOSE > 1
	fprintf(stdout, "no. of rk4 calls during compute_temp: %d\n", i+1);
	#endif

	/* map the temperature numbers back	*/
	xlate_temp_g2b(model, model->last_temp, model->last_trans);

	free_grid_model_vector(p);
}

/* debug print	*/
void debug_print_blist(blist_t *head, flp_t *flp)
{
	blist_t *ptr;
	fprintf(stdout, "printing blist information...\n");
	for(ptr = head; ptr; ptr = ptr->next) {
		fprintf(stdout, "unit: %s\n", flp->units[ptr->idx].name);
		fprintf(stdout, "occupancy: %f\n", ptr->occupancy);
	}
}

void debug_print_glist(glist_t *array, flp_t *flp)
{
	int i;
	fprintf(stdout, "printing glist information...\n");
	for(i=0; i < flp->n_units; i++)
		fprintf(stdout, "unit: %s\tstartx: %d\tendx: %d\tstarty: %d\tendy: %d\n",
				flp->units[i].name, array[i].j1, array[i].j2, array[i].i1, array[i].i2);
}

void debug_print_layer(grid_model_t *model, layer_t *layer)
{
	int i, j;
	fprintf(stdout, "printing layer information...\n");
	fprintf(stdout, "no: %d\n", layer->no);
	fprintf(stdout, "has_lateral: %d\n", layer->has_lateral);
	fprintf(stdout, "has_power: %d\n", layer->has_power);
	fprintf(stdout, "k: %f\n", layer->k);
	fprintf(stdout, "thickness: %f\n", layer->thickness);
	fprintf(stdout, "sp: %f\n", layer->sp);
	fprintf(stdout, "rx: %f\try: %f\trz: %f\tc: %f\n", 
			layer->rx, layer->ry, layer->rz, layer->c);

	fprintf(stdout, "printing b2gmap information...\n");
	for(i=0; i < model->rows; i++)
		for(j=0; j < model->cols; j++) {
			fprintf(stdout, "row: %d, col: %d\n", i, j);
			debug_print_blist(layer->b2gmap[i][j], layer->flp);
		}

	fprintf(stdout, "printing g2bmap information...\n");
	debug_print_glist(layer->g2bmap, layer->flp);
}

void debug_print_grid_model_vector(grid_model_vector_t *v, int nl, int nr, int nc)
{
	int n;

	fprintf(stdout, "printing cuboid information...\n");
	for(n=0; n < nl; n++)
		dump_dmatrix(v->cuboid[n], nr, nc);
}

void debug_print_grid(grid_model_t *model)
{
	int i;
	fprintf(stdout, "printing grid model information...\n");
	fprintf(stdout, "rows: %d\n", model->rows);
	fprintf(stdout, "cols: %d\n", model->cols);
	fprintf(stdout, "width: %f\n", model->width);
	fprintf(stdout, "height: %f\n", model->height);

	fprintf(stdout, "total_n_blocks: %d\n", model->total_n_blocks);
	fprintf(stdout, "map_mode: %d\n", model->map_mode);
	fprintf(stdout, "r_ready: %d\n", model->r_ready);
	fprintf(stdout, "c_ready: %d\n", model->c_ready);
	fprintf(stdout, "has_lcf: %d\n", model->has_lcf);

	fprintf(stdout, "printing last_steady information...\n");
	debug_print_grid_model_vector(model->last_steady, model->n_layers, 
								  model->rows, model->cols);

	fprintf(stdout, "printing last_trans information...\n");
	debug_print_grid_model_vector(model->last_trans, model->n_layers,
								  model->rows, model->cols);

	fprintf(stdout, "printing last_temp information...\n");
	if (model->last_temp)
		dump_dvector(model->last_temp, model->total_n_blocks);
	else
		fprintf(stdout, "(null)\n");

	for(i=0; i < model->n_layers; i++)
		debug_print_layer(model, &model->layers[i]);

	fprintf(stdout, "base_n_units: %d\n", model->base_n_units);
}