/*  pdnmesh - a 2D finite element solver
    Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>  
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.
  
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  $Id: solve.c,v 1.42 2005/02/17 03:43:16 sarod Exp $
*/



#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h> /* for sqrt() */
#include <string.h> /* for memcpy */
#include "types.h"


MY_DOUBLE  g_maxpot,g_minpot;

/* type of equation POISSON, HELMHOLTZ */
eq_type solve_equation;

/* a function to solve a matrix equation by symmetric LU */
/* decomposition or cholevsky decomposition */
static void 
solve_cholevsky( MY_DOUBLE **a, MY_DOUBLE *x, MY_DOUBLE *y, int N) 
{
  /* solves a system in the form [a][x] = [y] order N  where [a] is symmetric*/
  /* first cholevsky decomposition */
  int i,j,k;
  MY_DOUBLE temp;

#ifdef DEBUG
  for ( i = 0; i < N ; i++ ) {
    for ( j = 0 ; j <= i; j++ )
      printf(MDF" ",a[i][j]);
    printf("\n");
  }
#endif

  update_status_bar("Begin solving...");
  for ( j = 0; j < N; j++) {
    /* first diagonal element */
    temp = 0;
    for ( k = 0; k < j; k++) 
      temp += a[j][k]*a[j][k];
    a[j][j] = sqrtf((float)( a[j][j] - temp) );
#ifdef DEBUG
    printf("solve_chol: col %d diag "MDF"\n",j,a[j][j]);
#endif
    /* check for singularities */
    if ( a[j][j] <= 0.0f ) {
      a[j][j] = TOL; /* a small value */
      fprintf(stderr,"warning: solve_cholevsky: singular matrix\n");
    }
    /* now handle the i th column */
    for( i = j+1; i < N; i++) {
	  temp = 0;
	for ( k = 0; k < j; k++)
	  temp += a[i][k]*a[j][k];
	a[i][j]  = ( a[i][j] - temp )/a[j][j];
#ifdef DEBUG
	printf("solve_chol: row %d elem "MDF"\n",i,a[i][j]);
#endif
    }
  }
  
#ifdef DEBUG
  for ( i =0; i < N ; i++ ) {
    for ( j = 0 ; j <= i; j++ )
      printf(MDF" ",a[i][j]);
    printf("\n");
  }
#endif
 
  /* forward elimination */
  for ( i = 0; i < N; i++) {
    temp = 0;
    for ( k = 0; k < i; k++)
      temp += a[i][k]*y[k];
	 /* since we have eliminated singularities, no need to check a[i][i]*/
      y[i] = ( y[i] - temp )/a[i][i];
  }

  /* backward substitution */
  for ( i = N-1; i >= 0; i--) {
    temp = 0;
    for ( j = i+1; j< N; j++)
      temp += x[j]*a[j][i];
	 /* since we have eliminated singularities, no need to check a[i][i]*/
      x[i] = ( y[i] - temp )/a[i][i];
  }

  update_status_bar("Done solving.");
}


/* a function to build the global matrix or the matrix */
/* equation for the system - Poisson Equation */
static void 
build( MY_DOUBLE **P, MY_DOUBLE *q, MY_DOUBLE *x, MY_DOUBLE *y, MY_DOUBLE *phi, int NUN, int NE, int NN, MY_INT *renumber, MY_INT *re_renumber) 
{
  /* infd - input file descriptor                  */
  /* P[][] - global matrix to build   NUNxNUN        */
  /* q[] - global array to build     NUNx1          */
  /* x[] - x coordinates             NNx1          */
  /* y[] - y coordinates             NNx1          */
  /* phi[] - potentials              NNx1          */
  /* NUN - no of unknown potential nodes             */
  /* NE - no of elements                            */
  /* NN - no of nodes                               */

  /* now the running variables for each triangular element */
  int v[3]; /* to keep vertex number of each node */
  MY_DOUBLE pl[3][3];  /* local version of P[][] above */
  MY_DOUBLE ql[3];	/* local version of q[] above */
  MY_DOUBLE rho[3]; /* rho values */ /* conductivity */
  triangle *t;
 
  int    i,j;
  
  MY_DOUBLE 
    epsilon, /* permittivity */
    A, /* area of a triangulr element */
    k, /* value of the determinanat */
    yy0, yy1, yy2;
  
  MY_DOUBLE b[3]; /*local arrays */
  MY_DOUBLE c[3];

  

  update_status_bar("Begin building...");
	/* now read node data with unknown potentials */
  for ( i = 0; i < NUN; i++) {
    x[i]= (Mx(re_renumber[i],M))*g_xscale-g_xoff;
    y[i]= (My(re_renumber[i],M))*g_yscale-g_yoff;
  }
	
  /* now read node data with known potentials , which comes last */
  for ( i = NUN; i < NN; i++ ) { 
    x[i] = (Mx(re_renumber[i],M))*g_xscale-g_xoff;
    y[i] = (My(re_renumber[i],M))*g_yscale-g_yoff;
    phi[i] = (Mz(re_renumber[i],M));
  }
	
  /* now initialize P[][] and q[] */
/*  for ( i = 0; i < NUN; i++) {
    q[i] = 0;
    for ( j = 0; j <= i; j++) 
      P[i][j] = 0;
  } */
  
#ifdef DEBUG
	for  ( i = 0; i< NN; i++)
	    printf("build: %d "MDF", "MDF": "MDF"\n", i, x[i], y[i], phi[i]); 
#endif
	
	/* now read element data, whilst building the global matrices */
	DAG_traverse_list_reset(&dt);
	t=DAG_traverse_prune_list(&dt);
	while(t) {
	  if ( t &&t->status==LIVE ) {
	    
	    v[0] = renumber[t->p0];
	    v[1] = renumber[t->p1];
	    v[2] = renumber[t->p2]; /* ccw direction */
			if (bound.poly_array[t->boundary].rhoexp) {
	      rho[0] = ex(bound.poly_array[t->boundary].rhoexp,t->p0);
	      rho[1] = ex(bound.poly_array[t->boundary].rhoexp,t->p1);
	      rho[2] = ex(bound.poly_array[t->boundary].rhoexp,t->p2);
			} else {
	      rho[0] = bound.poly_array[t->boundary].rho;
	      rho[1] = bound.poly_array[t->boundary].rho;
	      rho[2] = bound.poly_array[t->boundary].rho;
			}

	    epsilon = 1/bound.poly_array[t->boundary].mu;
	    

	    /* now have read one triangle , add it to global */
	    /* first the coefficients */
	    /* some temporary values, neede several times */
	    yy0 = y[v[1]]-y[v[2]];
	    yy1 = y[v[2]]-y[v[0]];
	    yy2 = y[v[0]]-y[v[1]];
	    
	    /* area of the triangle */
	    A = 0.5*( x[v[0]]*yy0 + x[v[1]]*yy1 + x[v[2]]*yy2 ); 
	 /*   A = ABS(A); */ /* points are in ccw order to ensure this */
	    /* another constant */
	    k = 0.25/A;
	    /* now calculate b's and c's */
#ifdef DEBUG
	  printf("adding element with nodes %d %d and %d...\n\n", v[0], v[1], v[2]);
#endif
	  b[0]= yy0;
#ifdef DEBUG
	  printf("b0="MDF"  ", b[0]);
#endif
	  b[1]= yy1;
#ifdef DEBUG
	  printf("b1="MDF"  ", b[1]);
#endif
	  b[2]= yy2;
#ifdef DEBUG
	  printf("b2="MDF" \n", b[2]);
#endif	  
	  c[0] = ( x[v[2]]-x[v[1]] ); 
#ifdef DEBUG
	  printf("c0="MDF"  ", c[0]); 
#endif
	  c[1] = ( x[v[0]]-x[v[2]] );
#ifdef DEBUG
	  printf("c1="MDF"  ", c[1]);
#endif
	  c[2] = ( x[v[1]]-x[v[0]] );
#ifdef DEBUG
	  printf("c2="MDF" \n", c[2]);
#endif
	  
	  /* now build pl[] and ql */
	  for ( i = 0; i < 3; i++ ) 	
	    for ( j = 0; j < 3; j++ )
	      pl[i][j] = k*epsilon*(b[i]*b[j]+c[i]*c[j]); 
	  
	  ql[0] = A*(2*rho[0]+rho[1]+rho[2])*ONE_OVER_TWELVE;
	  ql[1] = A*(rho[0]+2*rho[1]+rho[2])*ONE_OVER_TWELVE;
	  ql[2] = A*(rho[0]+rho[1]+2*rho[2])*ONE_OVER_TWELVE;
#ifdef DEBUG 
	  printf("local p and q ..\n");
	  for ( i = 0; i <3; i++) {
	    printf("| ");
	    for ( j = 0; j < 3; j++)
	      printf(" "MDF" ", pl[i][j]);
	    printf("|| "MDF"  | \n", ql[i]); 
	  }	
#endif 
	  
	  
	  /* better to use multiplication than division */
	  /* now augment the global data set */
	  for ( i = 0; i < 3; i++)
	    if ( v[i] < NUN ) { /* we index from 0, not from 1 */
	      q[v[i]] = q[v[i]] + ql[i];
	      for ( j = 0; j < 3; j++ ) {
		if ( v[j] < NUN )
		  { 
					 if ( v[j] <= v[i] )
					 P[v[i]][v[j]] = P[v[i]][v[j]] + pl[i][j];
					 /* we only build half since we have half only */
		  }
		else 
		  { 
					 q[v[i]] = q[v[i]] - (pl[i][j])*phi[v[j]];
		  } 
	      }
	    } /*else printf("rejecting %d\n", v[i]);*/
#ifdef DEBUG
	  printf("building global matrix ..\n");
#endif
	  }
			t=DAG_traverse_prune_list(&dt);
	}
	
#ifdef DEBUG	  
	  for ( i = 0; i < NUN; i++){
	    printf("|  ");
	    for( j = 0; j <= i; j++)
	      printf(" "MDF"  ", P[i][j]);
	    printf("=  "MDF"  |\n", q[i]);
	    } 
	   for ( i = 0; i < NN; i++) 
		printf("node %d potential "MDF"\n", i, phi[i]); 
#endif
  update_status_bar("Done building.");
}

/* a function to build the global matrix or the matrix */
/* equation for the system - Helmholtz Equation */
/* FIXME - need to solve only homogeneous case. dont need q */
static void 
build_helmholtz( MY_DOUBLE **P, MY_DOUBLE **T,  MY_DOUBLE *x, MY_DOUBLE *y, int NUN, int NE, int NN, MY_INT *renumber, MY_INT *re_renumber) 
{
  /* infd - input file descriptor                  */
  /* P[][] - global matrix to build   NUNxNUN        */
  /* T[][] - global matrix to build   NUNxNUN        */
  /* x[] - x coordinates             NNx1          */
  /* y[] - y coordinates             NNx1          */
  /* NUN - no of unknown potential nodes             */
  /* NE - no of elements                            */
  /* NN - no of nodes                               */

  /* now the running variables for each triangular element */
  int v[3]; /* to keep vertex number of each node */
  MY_DOUBLE pl[3][3];  /* local version of P[][] above */
  MY_DOUBLE tl[3][3];  /* local version of T[][] above */
  /*MY_DOUBLE ql[3];*/	/* local version of q[] above */
  /*long double rho[3];*/ /* rho values */ /* conductivity */
  triangle *t;
 
  int    i,j;
  
  MY_DOUBLE 
    epsilon, /* permittivity */
    A, /* area of a triangulr element */
    k, /* value of the determinanat */
    yy0, yy1, yy2;
  
  MY_DOUBLE b[3]; /*local arrays */
  MY_DOUBLE c[3];

  

 update_status_bar("Start building Eigenproblem.");
	/* now read node data with unknown potentials */
  for ( i = 0; i < NUN; i++) {
    x[i]= (Mx(re_renumber[i],M))*g_xscale-g_xoff;
    y[i]= (My(re_renumber[i],M))*g_yscale-g_yoff;
  }
	
  /* now read node data with known potentials , which comes last */
  for ( i = NUN; i < NN; i++ ) { 
    x[i] = (Mx(re_renumber[i],M))*g_xscale-g_xoff;
    y[i] = (My(re_renumber[i],M))*g_yscale-g_yoff;
  }
	
#ifdef DEBUG
	for  ( i = 0; i< NN; i++)
	    printf("build: %d->%d "MDF", "MDF"\n", i,re_renumber[i], x[i], y[i]); 
#endif
	
	/* now read element data, whilst building the global matrices */
	DAG_traverse_list_reset(&dt);
	t=DAG_traverse_prune_list(&dt);
	while(t) {
	  if ( t &&t->status==LIVE ) {
	    
	    v[0] = renumber[t->p0];
	    v[1] = renumber[t->p1];
	    v[2] = renumber[t->p2]; /* ccw direction */
	   /* rho[0] = bound.poly_array[t->boundary].rho;
	    rho[1] = bound.poly_array[t->boundary].rho;
	    rho[2] = bound.poly_array[t->boundary].rho; */
	    epsilon = 1/bound.poly_array[t->boundary].mu; 
	    

	    /* now have read one triangle , add it to global */
	    /* first the coefficients */
	    /* some temporary values, neede several times */
	    yy0 = y[v[1]]-y[v[2]];
	    yy1 = y[v[2]]-y[v[0]];
	    yy2 = y[v[0]]-y[v[1]];
	    
	    /* area of the triangle */
	    A = 0.5*( x[v[0]]*yy0 + x[v[1]]*yy1 + x[v[2]]*yy2 ); 
	    /*A = ABS(A); *//* we keep points (0,1,2) in ccw order, A > 0 */
	    /* another constant */
	    k = 0.25/(A*A);
	    /* now calculate b's and c's */
#ifdef DEBUG
	  printf("adding element with nodes %d %d and %d area "MDF"...?("MDF","MDF")("MDF","MDF")("MDF","MDF")\n\n", v[0], v[1], v[2],A,x[v[0]],y[v[0]],x[v[1]],y[v[1]],x[v[2]],y[v[2]]);
#endif
	  b[0]= yy0;
#ifdef DEBUG
	  printf("b0="MDF"  ", b[0]);
#endif
	  b[1]= yy1;
#ifdef DEBUG
	  printf("b1="MDF"  ", b[1]);
#endif
	  b[2]= yy2;
#ifdef DEBUG
	  printf("b2="MDF" \n", b[2]);
#endif	  
	  c[0] = ( x[v[2]]-x[v[1]] ); 
#ifdef DEBUG
	  printf("c0="MDF"  ", c[0]); 
#endif
	  c[1] = ( x[v[0]]-x[v[2]] );
#ifdef DEBUG
	  printf("c1="MDF"  ", c[1]);
#endif
	  c[2] = ( x[v[1]]-x[v[0]] );
#ifdef DEBUG
	  printf("c2="MDF" \n", c[2]);
#endif
	  
	  /* now build pl[] and ql */
	  for ( i = 0; i < 3; i++ ) 	
	    for ( j = 0; j < 3; j++ )
	      pl[i][j] = k*epsilon*(b[i]*b[j]+c[i]*c[j]); 
	  
			 /* build tl */
			 tl[0][0] = tl[1][1] = tl[2][2] = A*ONE_OVER_SIX;
				tl[0][1] = tl[0][2] = tl[1][0] = tl[1][2] = tl[2][0] = tl[2][1] = A*ONE_OVER_TWELVE; 
	 /* ql[0] = A*(2*rho[0]+rho[1]+rho[2])/12.0;
	  ql[1] = A*(rho[0]+2*rho[1]+rho[2])/12.0;
	  ql[2] = A*(rho[0]+rho[1]+2*rho[2])/12.0; */
#ifdef DEBUG 
	  printf("local p and q and t..\n");
	  for ( i = 0; i <3; i++) {
	    printf("| ");
	    for ( j = 0; j < 3; j++)
	      printf(" "MDF" ", pl[i][j]);
	    printf("||"); 
      for ( j = 0; j < 3; j++)
	      printf(" "MDF" ", tl[i][j]);
	    printf("|\n"); 
	  }	
#endif 
	  
	  
	  /* better to use multiplication than division */
	  /* now augment the global data set */
	  for ( i = 0; i < 3; i++)
	    if ( v[i] < NUN ) { /* we index from 0, not from 1 */
	    /*  q[v[i]] = q[v[i]] + ql[i]; */
	      for ( j = 0; j < 3; j++ ) {
		if ( v[j] < NUN )
		  { 
					 if ( v[j] <= v[i] ) { /* matricex P, T are symmetric */
					 P[v[i]][v[j]] = P[v[i]][v[j]] + pl[i][j];
					 T[v[i]][v[j]] = T[v[i]][v[j]] + tl[i][j];
						}
					 /* we only build half since we have half only */
		  }/* else { 
					 q[v[i]] = q[v[i]] - (pl[i][j])*phi[v[j]];
		  } */
	      }
	    } /*else printf("rejecting %d\n", v[i]);*/
#ifdef DEBUG
	  printf("building global matrix ..\n");
    printf("global matrices P and T\n");
	  for ( i = 0; i < NUN; i++){
	    printf("|  ");
	    for( j = 0; j <= i; j++)
	      printf(" "MDF"  ", P[i][j]);
	    printf("|\n");
	    } 
     for ( i = 0; i < NUN; i++){
	    printf("|  ");
	    for( j = 0; j <= i; j++)
	      printf(" "MDF"  ", T[i][j]);
	    printf("|\n");
	    } 

#endif
	  }
			t=DAG_traverse_prune_list(&dt);
	}
	
#ifdef DEBUG	  
	printf("global matrices P and T\n");
	  for ( i = 0; i < NUN; i++){
	    printf("|  ");
	    for( j = 0; j <= i; j++)
	      printf(" "MDF"  ", P[i][j]);
	    printf("|\n");
	    } 
     for ( i = 0; i < NUN; i++){
	    printf("|  ");
	    for( j = 0; j <= i; j++)
	      printf(" "MDF"  ", T[i][j]);
	    printf("|\n");
	    } 
#endif

  update_status_bar("Done building Eigenproblem.");
}

/* a function to build the global matrix or the matrix */
/* equation for the system - Helmholtz Equation */
static void 
build_helmholtz_sparse( SPMat *P, SPMat *T,  MY_DOUBLE *x, MY_DOUBLE *y, int NUN, int NE, int NN, MY_INT *renumber, MY_INT *re_renumber) 
{
  /* infd - input file descriptor                  */
  /* P[][] - global matrix to build   NUNxNUN        */
  /* T[][] - global matrix to build   NUNxNUN        */
  /* x[] - x coordinates             NNx1          */
  /* y[] - y coordinates             NNx1          */
  /* NUN - no of unknown potential nodes             */
  /* NE - no of elements                            */
  /* NN - no of nodes                               */

  /* now the running variables for each triangular element */
  int v[3]; /* to keep vertex number of each node */
  MY_DOUBLE pl[3][3];  /* local version of P[][] above */
  MY_DOUBLE tl[3][3];  /* local version of T[][] above */
  /*MY_DOUBLE ql[3];*/	/* local version of q[] above */
  /*long double rho[3];*/ /* rho values */ /* conductivity */
  triangle *t;
 
  int    i,j;
  
  MY_DOUBLE 
    epsilon, /* permittivity */
    A, /* area of a triangulr element */
    k, /* value of the determinanat */
    yy0, yy1, yy2;
  
  MY_DOUBLE b[3]; /*local arrays */
  MY_DOUBLE c[3];

  

 update_status_bar("Start building Eigenproblem.");
	/* now read node data with unknown potentials */
  for ( i = 0; i < NUN; i++) {
    x[i]= (Mx(re_renumber[i],M))*g_xscale-g_xoff;
    y[i]= (My(re_renumber[i],M))*g_yscale-g_yoff;
  }
	
  /* now read node data with known potentials , which comes last */
  for ( i = NUN; i < NN; i++ ) { 
    x[i] = (Mx(re_renumber[i],M))*g_xscale-g_xoff;
    y[i] = (My(re_renumber[i],M))*g_yscale-g_yoff;
  }
	
#ifdef DEBUG
	for  ( i = 0; i< NN; i++)
	    printf("build: %d->%d "MDF", "MDF"\n", i,re_renumber[i], x[i], y[i]); 
#endif
	
	/* now read element data, whilst building the global matrices */
	DAG_traverse_list_reset(&dt);
	t=DAG_traverse_prune_list(&dt);
	while(t) {
	  if ( t &&t->status==LIVE ) {
	    
	    v[0] = renumber[t->p0];
	    v[1] = renumber[t->p1];
	    v[2] = renumber[t->p2]; /* ccw direction */
	   /* rho[0] = bound.poly_array[t->boundary].rho;
	    rho[1] = bound.poly_array[t->boundary].rho;
	    rho[2] = bound.poly_array[t->boundary].rho; */
	    epsilon = 1/bound.poly_array[t->boundary].mu; 
	    

	    /* now have read one triangle , add it to global */
	    /* first the coefficients */
	    /* some temporary values, neede several times */
	    yy0 = y[v[1]]-y[v[2]];
	    yy1 = y[v[2]]-y[v[0]];
	    yy2 = y[v[0]]-y[v[1]];
	    
	    /* area of the triangle */
	    A = 0.5*( x[v[0]]*yy0 + x[v[1]]*yy1 + x[v[2]]*yy2 ); 
	    /*A = ABS(A); *//* we keep points (0,1,2) in ccw order, A > 0 */
	    /* another constant */
	    k = 0.25/(A*A);
	    /* now calculate b's and c's */
#ifdef DEBUG
	  printf("adding element with nodes %d %d and %d area "MDF"...?("MDF","MDF")("MDF","MDF")("MDF","MDF")\n\n", v[0], v[1], v[2],A,x[v[0]],y[v[0]],x[v[1]],y[v[1]],x[v[2]],y[v[2]]);
#endif
	  b[0]= yy0;
#ifdef DEBUG
	  printf("b0="MDF"  ", b[0]);
#endif
	  b[1]= yy1;
#ifdef DEBUG
	  printf("b1="MDF"  ", b[1]);
#endif
	  b[2]= yy2;
#ifdef DEBUG
	  printf("b2="MDF" \n", b[2]);
#endif	  
	  c[0] = ( x[v[2]]-x[v[1]] ); 
#ifdef DEBUG
	  printf("c0="MDF"  ", c[0]); 
#endif
	  c[1] = ( x[v[0]]-x[v[2]] );
#ifdef DEBUG
	  printf("c1="MDF"  ", c[1]);
#endif
	  c[2] = ( x[v[1]]-x[v[0]] );
#ifdef DEBUG
	  printf("c2="MDF" \n", c[2]);
#endif
	  
	  /* now build pl[] and ql */
	  for ( i = 0; i < 3; i++ ) 	
	    for ( j = 0; j < 3; j++ )
	      pl[i][j] = k*epsilon*(b[i]*b[j]+c[i]*c[j]); 
	  
			 /* build tl */
			 tl[0][0] = tl[1][1] = tl[2][2] = A*ONE_OVER_SIX;
				tl[0][1] = tl[0][2] = tl[1][0] = tl[1][2] = tl[2][0] = tl[2][1] = A*ONE_OVER_TWELVE; 
	 /* ql[0] = A*(2*rho[0]+rho[1]+rho[2])/12.0;
	  ql[1] = A*(rho[0]+2*rho[1]+rho[2])/12.0;
	  ql[2] = A*(rho[0]+rho[1]+2*rho[2])/12.0; */
#ifdef DEBUG 
	  printf("local p and q and t..\n");
	  for ( i = 0; i <3; i++) {
	    printf("| ");
	    for ( j = 0; j < 3; j++)
	      printf(" "MDF" ", pl[i][j]);
	    printf("||"); 
      for ( j = 0; j < 3; j++)
	      printf(" "MDF" ", tl[i][j]);
	    printf("|\n"); 
	  }	
#endif 
	  
	  
	  /* better to use multiplication than division */
	  /* now augment the global data set */
	  for ( i = 0; i < 3; i++)
	    if ( v[i] < NUN ) { /* we index from 0, not from 1 */
	    /*  q[v[i]] = q[v[i]] + ql[i]; */
	      for ( j = 0; j < 3; j++ ) {
		if ( v[j] < NUN )
		  { 
					 if ( v[j] <= v[i] ) { /* matricex P, T are symmetric */
					 /*P[v[i]][v[j]] = P[v[i]][v[j]] + pl[i][j]; */
           SPM_add(P,v[i],v[j],pl[i][j]);
					 /*T[v[i]][v[j]] = T[v[i]][v[j]] + tl[i][j]; */
           SPM_add(T,v[i],v[j],tl[i][j]);
						}
					 /* we only build half since we have half only */
		  }
	   }
	  } /*else printf("rejecting %d\n", v[i]);*/
#ifdef DEBUG
	  printf("building global matrix ..\n");
    printf("global matrices P and T\n");
	  for ( i = 0; i < NUN; i++){
	    printf("|  ");
	    for( j = 0; j <= i; j++)
	      printf(" "MDF"  ", SPM_e(P,i,j));
	    printf("|\n");
	    } 
     for ( i = 0; i < NUN; i++){
	    printf("|  ");
	    for( j = 0; j <= i; j++)
	      printf(" "MDF"  ", SPM_e(T,i,j));
	    printf("|\n");
	    } 

#endif
	  }
			t=DAG_traverse_prune_list(&dt);
	}
	
#ifdef DEBUG	  
	printf("global matrices P and T\n");
	  for ( i = 0; i < NUN; i++){
	    printf("|  ");
	    for( j = 0; j <= i; j++)
	      printf(" "MDF"  ", SPM_e(P,i,j));
	    printf("|\n");
	    } 
     for ( i = 0; i < NUN; i++){
	    printf("|  ");
	    for( j = 0; j <= i; j++)
	      printf(" "MDF"  ", SPM_e(T,i,j));
	    printf("|\n");
	    } 
#endif

  update_status_bar("Done building Eigenproblem.");
}

#if defined (USE_LAPACK) || defined (USE_MKL)
 /* nothing */
#else
static void
solve_helmholtz(MY_DOUBLE **P,MY_DOUBLE **T, MY_DOUBLE **x,  MY_DOUBLE *ev,  MY_INT  NUN, MY_INT n_to_find)
{
 /* solves the symmetrix, generalized eigenvalue problem 
		* (P-lambda. T)x= 0
		* P, T = NUN by NUN matrices, only lower triangle
		* x = eigenvectors to be found size (m x n_to_find)
    * ev = eigenvalues to be found, n_to_find by 1 array
		*/
  MY_INT i,j,k;
  MY_DOUBLE **L,**Pt;
	MY_DOUBLE temp;
		/* first find the cholesky factorization of T */
		/* T=L L' */
  if((L= (MY_DOUBLE**) calloc(NUN, sizeof(MY_DOUBLE*)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
	for (i = 0; i <NUN; ++i) {
	  if((L[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 memcpy((void*)L[i],(void*)T[i],(size_t)(i+1)*sizeof(MY_DOUBLE));
 }

 for ( j = 0; j < NUN; j++) {
    /* first diagonal element */
    temp = 0;
    for ( k = 0; k < j; k++) 
      temp += L[j][k]*L[j][k];
    L[j][j] = sqrtf((float)( L[j][j] - temp) );
#ifdef DEBUG
    printf("solve_helmholtz: col %d diag "MDF"\n",j,L[j][j]);
#endif
    /* check for singularities */
    if ( L[j][j] <= 0.0f ) {
      L[j][j] = TOL; /* a small value */
      fprintf(stderr,"warning: solve_cholevsky: singular matrix\n");
    }
    /* now handle the i th column */
    for( i = j+1; i < NUN; i++) {
	  temp = 0;
	for ( k = 0; k < j; k++)
	  temp += L[i][k]*L[j][k];
	L[i][j]  = ( L[i][j] - temp )/L[j][j];
#ifdef DEBUG
	printf("solve_helm: row %d elem "MDF"\n",i,L[i][j]);
#endif
    }
  }

#ifdef DEBUG
  for ( i =0; i < NUN ; i++ ) {
    for ( j = 0 ; j <= i; j++ )
      printf(""MDF" ",L[i][j]);
    printf("\n");
  }
#endif
 
		/* now we solve the system
			* (L^{-1} P L^{-T}-lambda I)y =0 
			* and Ly =x
			* to get x
			*/
		/* first: create a full matrix for P */
   if((Pt= (MY_DOUBLE**) calloc(NUN, sizeof(MY_DOUBLE*)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
	for (i = 0; i <NUN; ++i) {
	  if((Pt[i] = ( MY_DOUBLE * ) calloc( (size_t)NUN, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 memcpy((void*)Pt[i],(void*)P[i],(size_t)(i+1)*sizeof(MY_DOUBLE));
	/* fill up the rest of Pt */
			for (j=0;j<i;j++) {
         Pt[j][i]=P[i][j];
		  }
 }


		/* calculation of L^{-1}P */
	/* column by column */
  /* forward elimination */
		for (k=0;k<NUN;k++) {
  for (i=0; i<NUN; i++) {
						temp=0;
						for (j=0; j<i;j++) {
													temp+=L[i][j]*Pt[j][k];
						}
					 Pt[i][k]=(Pt[i][k]-temp)/L[i][i];
		}
		}

		/* calculation of P L^{-T} */
		/* row by row */
		for (k=0;k<NUN;k++) {
		for (i=0;i<NUN;i++) {
            temp=0;
						for (j=0; j<i;j++) {
													temp+=L[i][j]*Pt[k][j];
						}
            Pt[k][i]=(Pt[k][i]-temp)/L[i][i];
		 }
		}
#ifdef DEBUG
  for ( i =0; i < NUN ; i++ ) {
    for ( j = 0 ; j <NUN; j++ )
      printf(""MDF" ",Pt[i][j]);
    printf("\n");
  }
#endif

  /* find the eigenvector */
  /*find_eigenvectors(Pt, x, ev, NUN-2, NUN); */
	/*x -  m x n_to_find */
  find_eigenvectors(Pt, x, ev, NUN, n_to_find);

	/* now we have L^{T} x = eigenvector
	 * hence solve to get true eigenvector */
	for (k=0;k<n_to_find;k++) {
	for (i=NUN-1;i>=0;i--) {
            temp=0;
						for (j=NUN-1; j>i;j--) {
													temp+=L[j][i]*x[i][k];
						}
            x[i][k]=(x[i][k]-temp)/L[i][i];

#ifdef DEBUG
					printf("x[%d][%d]="MDF" ",i,k,x[i][k]);
#endif
  }
	}
#ifdef DEBUG
    printf("\n");
#endif
	for (i=0; i<NUN; i++)  {
								free(L[i]);
								free(Pt[i]);
	}
	free(L);
	free(Pt);
}

static void
solve_helmholtz_sparse(SPMat *P,SPMat *T, MY_DOUBLE **x,  MY_DOUBLE *ev,  MY_INT  NUN, MY_INT n_to_find)
{
 /* solves the symmetrix, generalized eigenvalue problem 
		* (P-lambda. T)x= 0
		* P, T = NUN by NUN matrices, only lower triangle
		* x = eigenvectors to be found size (m x n_to_find)
    * ev = eigenvalues to be found, n_to_find by 1 array
		*/
  MY_INT i,j,k;
  MY_DOUBLE **L,**Pt;
	MY_DOUBLE temp;
		/* first find the cholesky factorization of T */
		/* T=L L' */
  if((L= (MY_DOUBLE**) calloc(NUN, sizeof(MY_DOUBLE*)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
	for (i = 0; i <NUN; ++i) {
	  if((L[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 //memcpy((void*)L[i],(void*)T[i],(size_t)(i+1)*sizeof(MY_DOUBLE));
  for (j=0; j<=i; j++)
   L[i][j]=SPM_e(T,i,j);
 }

 for ( j = 0; j < NUN; j++) {
    /* first diagonal element */
    temp = 0;
    for ( k = 0; k < j; k++) 
      temp += L[j][k]*L[j][k];
    L[j][j] = sqrtf((float)( L[j][j] - temp) );
#ifdef DEBUG
    printf("solve_helmholtz: col %d diag "MDF"\n",j,L[j][j]);
#endif
    /* check for singularities */
    if ( L[j][j] <= 0.0f ) {
      L[j][j] = TOL; /* a small value */
      fprintf(stderr,"warning: solve_cholevsky: singular matrix\n");
    }
    /* now handle the i th column */
    for( i = j+1; i < NUN; i++) {
	  temp = 0;
	for ( k = 0; k < j; k++)
	  temp += L[i][k]*L[j][k];
	L[i][j]  = ( L[i][j] - temp )/L[j][j];
#ifdef DEBUG
	printf("solve_helm: row %d elem "MDF"\n",i,L[i][j]);
#endif
    }
  }

#ifdef DEBUG
  for ( i =0; i < NUN ; i++ ) {
    for ( j = 0 ; j <= i; j++ )
      printf(""MDF" ",L[i][j]);
    printf("\n");
  }
#endif
 
		/* now we solve the system
			* (L^{-1} P L^{-T}-lambda I)y =0 
			* and Ly =x
			* to get x
			*/
		/* first: create a full matrix for P */
   if((Pt= (MY_DOUBLE**) calloc(NUN, sizeof(MY_DOUBLE*)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
	for (i = 0; i <NUN; ++i) {
	  if((Pt[i] = ( MY_DOUBLE * ) calloc( (size_t)NUN, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 //memcpy((void*)Pt[i],(void*)P[i],(size_t)(i+1)*sizeof(MY_DOUBLE));
 for(j=0;j<=i;j++)
  Pt[j][i]=Pt[i][j]=SPM_e(P,i,j);
 }


	/* calculation of L^{-1}P */
	/* column by column */
  /* forward elimination */
		for (k=0;k<NUN;k++) {
  for (i=0; i<NUN; i++) {
						temp=0;
						for (j=0; j<i;j++) {
							temp+=L[i][j]*Pt[j][k];
						}
					 Pt[i][k]=(Pt[i][k]-temp)/L[i][i];
		}
		}

		/* calculation of P L^{-T} */
		/* row by row */
		for (k=0;k<NUN;k++) {
		for (i=0;i<NUN;i++) {
            temp=0;
						for (j=0; j<i;j++) {
								temp+=L[i][j]*Pt[k][j];
						}
            Pt[k][i]=(Pt[k][i]-temp)/L[i][i];
		 }
		}
#ifdef DEBUG
  for ( i =0; i < NUN ; i++ ) {
    for ( j = 0 ; j <NUN; j++ )
      printf(""MDF" ",Pt[i][j]);
    printf("\n");
  }
#endif

  /* find the eigenvector */
  /*find_eigenvectors(Pt, x, ev, NUN-2, NUN); */
	/*x -  m x n_to_find */
  find_eigenvectors(Pt, x, ev, NUN, n_to_find);

	/* now we have L^{T} x = eigenvector
	 * hence solve to get true eigenvector */
	for (k=0;k<n_to_find;k++) {
	for (i=NUN-1;i>=0;i--) {
            temp=0;
						for (j=NUN-1; j>i;j--) {
													temp+=L[j][i]*x[i][k];
						}
            x[i][k]=(x[i][k]-temp)/L[i][i];

#ifdef DEBUG
					printf("x[%d][%d]="MDF" ",i,k,x[i][k]);
#endif
  }
	}
#ifdef DEBUG
    printf("\n");
#endif
	for (i=0; i<NUN; i++)  {
								free(L[i]);
								free(Pt[i]);
	}
	free(L);
	free(Pt);
}

#endif /* !USE_LAPACK || USE_MKL */

/* after genarating the mesh, need to re-number points */
/* so that unknown points are first */
/* then we solve the problem -poisson equation*/
int 
re_number_and_solve_poisson(void)
{
  unsigned MY_INT i,j;
 MY_INT * renumber; /* re-numbering array */
 MY_INT * re_renumber; /* re-renumbering array */
 /* used to map potentials back to the original points */


  unsigned MY_INT unknowns=0; /* unknown points */
  unsigned MY_INT total = 0; /* total points */
  MY_INT triangles = 0; /* known points */


  triangle *t;
  
  /* first the global node data, all total size */
  MY_DOUBLE *x ; /* x coordinates */
  MY_DOUBLE *y ; /* y coordinates */
  MY_DOUBLE *phi ; /* potentials at each node */
  
  /* the global variables */
  MY_DOUBLE **P ; /* size unknownsxunknowns note we do not need a NN x NN matrix, only the unknowns
		     */
  MY_DOUBLE *q ; /* unknowns sizesame here */


  /* allocate memory for arrays */
  if ((renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT )))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((re_renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  /* need to solve for only points referred by triangles */
  /* SO -- */
		DAG_traverse_list_reset(&dt);
		t=DAG_traverse_prune_list(&dt);
  while ( t ) {
    if ( t && t->status==LIVE ) {
      /* mark the points in this triangle */
      /* use re_renumber[] array for this */
      re_renumber[t->p0] = 1;
      re_renumber[t->p1] = 1;
      re_renumber[t->p2] = 1;

      triangles++;
    }
		t=DAG_traverse_prune_list(&dt);
  }

  /* now find unknowns, or INSIDE points */
  j = 0;
  for ( i = 0; i < M.count; i++ )
    if (  Mval(i,M) == VAR  && re_renumber[i] == 1) {
      renumber[i] = j++; 
#ifdef DEBUG
						printf("renumber_and_solve: renumbering variable point %d as %d\n",i,renumber[i]);
#endif
      unknowns++;
    }
  /* now renumber the known points */
  total = unknowns;
  for ( i = 0; i < M.count; i++ )
    if ( Mval(i,M) == FX  && re_renumber[i]==1 ) {
      renumber[i] = j++; 
#ifdef DEBUG
						printf("renumber_and_solve: renumbering fixed point %d as %d\n",i,renumber[i]);
#endif
      total++;
    }
 /* finally give numbers to points 0,1,2 */
	/* we will not use these points */
		renumber[0]=j++;
		renumber[1]=j++;
		renumber[2]=j++;

/* now construct re-renmbering array */
		  i = 0;
				/* do an insertion sort */
				while ( i < M.count ) {
						for ( j = 0; j < M.count; j++ ) {
						   if ( i==(unsigned)renumber[j] )
									  break;
						}
										re_renumber[i] = j;
										i++;
				}

#ifdef DEBUG
printf("build: renumbering %d points with %d unknowns out of %d total\n",total,unknowns,M.count);
#endif
  /*renumber_points(renumber, re_renumber, total, unknowns); */
  
#ifdef DEBUG
  printf("i,  renum, rerenum\n");
  for ( i = 0; i < total; i++ )
    printf("%d %d %d\n",i,renumber[i],re_renumber[i]);
  for ( i = 0; i < total; i++ )
    printf("%d  "MDF"  "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
#endif

  
  /* we free these two arrays only after solving the problem */
  /* now solve the problem */
  /* first allocate memory for the arrays */
  if((x=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((y=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((phi=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((q=(MY_DOUBLE*)calloc(unknowns, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

	
	/* now the matrix */
	/* STEP 1: SET UP THE ROWS. */
	if((P= (MY_DOUBLE**) calloc(unknowns , sizeof(MY_DOUBLE*)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
	/* STEP 2: SET UP THE COLUMNS. */
	for (i = 0; i < unknowns; ++i)
			  /* variable length columns */
	  if((P[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

	/* Poisson equation */
	/* now the memory allocation is over */
  build( P, q, x, y, phi, unknowns, triangles, total, renumber, re_renumber);
  /* if user has chosen so, solve by conj. grad method*/
  /* by default solve by cholevsky decomposition */
  solve_cholevsky(P,phi, q, unknowns);
  
	/* after solution, array phi[] will contain the new potentials */
  /* now put them back to the point data */
  /* we also note the maximum and minimum values for contour plotting */
  g_maxpot = -1000;
  g_minpot = 1000;
  for ( i = 0; i < total; i++ ) {
    Mz(re_renumber[i],M) = phi[i];
    if ( phi[i] > g_maxpot )
      g_maxpot = phi[i];
    if ( phi[i] < g_minpot )
      g_minpot = phi[i];
  }



  /* now we can free the renumber[] array */
  free(renumber);
  /* now free this array too */
  free(re_renumber);
  
  /* free everything else */
  
  /* first the arrays */
  free(x);
  free(y);
  free(phi);
  free(q);
  
  /* next the matrix */
  for ( i = 0; i < unknowns ; i++)
    free(P[i]);
  free(P);
 
  /* now printout a summery */
#ifdef DEBUG
  for ( i = 0; i < total; i++ )
    printf("renumber_and_solve: %d  "MDF"  "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
#endif


			 return(0);
}

/* solve helmholtz equation */
int 
re_number_and_solve_helmholtz(void)
{
  unsigned MY_INT i,j;
  MY_INT k;
 MY_INT * renumber; /* re-numbering array */
 MY_INT * re_renumber; /* re-renumbering array */
/* used to map potentials back to the original points */


  unsigned MY_INT unknowns=0; /* unknown points */
  unsigned MY_INT total = 0; /* total points */
  MY_INT triangles = 0; /* known points */


  triangle *t;
  
  /* first the global node data, all total size */
  MY_DOUBLE *x ; /* x coordinates */
  MY_DOUBLE *y ; /* y coordinates */
  MY_DOUBLE **v; /* eigenvector matrix */
  
  MY_DOUBLE *max_value_array;

  /* the global variables */
  MY_DOUBLE **P ; /* size unknowns x unknowns note we do not need a NN x NN matrix, only the unknowns
		     */
		MY_DOUBLE **T=0; /* used in Helmholtz equation */


  /* allocate memory for arrays */
  if ((renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT )))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((re_renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  /* need to solve for only points referred by triangles */
  /* SO -- */
		DAG_traverse_list_reset(&dt);
		t=DAG_traverse_prune_list(&dt);
  while ( t ) {
    if ( t && t->status==LIVE ) {
      /* mark the points in this triangle */
      /* use re_renumber[] array for this */
      re_renumber[t->p0] = 1;
      re_renumber[t->p1] = 1;
      re_renumber[t->p2] = 1;

      triangles++;
    }
		t=DAG_traverse_prune_list(&dt);
  }

  /* now find unknowns, or INSIDE points */
  j = 0;
  for ( i = 0; i < M.count; i++ )
    if (  Mval(i,M) == VAR  && re_renumber[i] == 1) {
      renumber[i] = j++; 
#ifdef DEBUG
						printf("renumber_and_solve: renumbering variable point %d as %d\n",i,renumber[i]);
#endif
      unknowns++;
    }
  /* now renumber the known points */
  total = unknowns;
  for ( i = 0; i < M.count; i++ )
    if ( Mval(i,M) == FX  && re_renumber[i]==1 ) {
      renumber[i] = j++; 
#ifdef DEBUG
						printf("renumber_and_solve: renumbering fixed point %d as %d\n",i,renumber[i]);
#endif
      total++;
    }
 /* finally give numbers to points 0,1,2 */
	/* we will not use these points */
		renumber[0]=j++;
		renumber[1]=j++;
		renumber[2]=j++;

/* now construct re-renmbering array */
		  i = 0;
				/* do an insertion sort */
				while ( i < M.count ) {
						for ( j = 0; j < M.count; j++ ) {
						   if ( renumber[j] == (int)i )
									  break;
						}
										re_renumber[i] = j;
										i++;
				}

#ifdef DEBUG
printf("build: renumbering %d points with %d unknowns out of %d total\n",total,unknowns,M.count);
#endif
  /*renumber_points(renumber, re_renumber, total, unknowns); */
  
#ifdef DEBUG
  printf("i,  renum, rerenum\n");
  for ( i = 0; i < total; i++ )
    printf("%d %d %d\n",i,renumber[i],re_renumber[i]);
  for ( i = 0; i < total; i++ )
    printf("%d  "MDF"  "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
#endif

  
  /* we free these two arrays only after solving the problem */
  /* now solve the problem */
  /* first allocate memory for the arrays */
  if((x=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((y=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

	/* structure to store all eigenvectors - matrix */
  if((v=(MY_DOUBLE **)calloc(unknowns, sizeof(MY_DOUBLE*)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 for (i = 0; i <unknowns ; i++)
	  if((v[i] = ( MY_DOUBLE * ) calloc( (size_t) degree_of_freedom, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

	
	/* now the matrix */
	/* STEP 1: SET UP THE ROWS. */
	if((P= (MY_DOUBLE**) calloc(unknowns , sizeof(MY_DOUBLE*)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
	/* STEP 2: SET UP THE COLUMNS. */
	for (i = 0; i < unknowns; ++i)
			  /* variable length columns */
	  if((P[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((T= (MY_DOUBLE**) calloc(unknowns , sizeof(MY_DOUBLE*)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
	for (i = 0; i < unknowns; ++i)
	  if((T[i] = ( MY_DOUBLE * ) calloc( (size_t)i+1, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
#ifdef DEBUG	
	printf("solving problem with %d unknowns, %d triangles\n",unknowns,triangles);
#endif

/*build_helmholtz( MY_DOUBLE **P, MY_DOUBLE **T,  MY_DOUBLE *x, MY_DOUBLE *y, int NUN, int NE, int NN)  */
  build_helmholtz( P, T, x, y, unknowns, triangles, total, renumber, re_renumber);

/*solve_helmholtz(MY_DOUBLE **P,MY_DOUBLE **T, MY_DOUBLE **x,  MY_INT  NUN, MY_INT n_to_find) */
#if defined (USE_LAPACK) || defined (USE_MKL)
	solve_helmholtz_lapack(P,T,v,eig_array,unknowns,degree_of_freedom);
#else
	solve_helmholtz(P,T,v,eig_array,unknowns,degree_of_freedom);
#endif /* ! USE_LAPACK || USE_MKL */

	/* all eigenvectors are normalized, hence normlization has to be done individually */
  	/* we copy back elements of eigenvectors to z arrays of each point */
	/* in order to do that, we need to change size of those z arrays */
  /* find ranges for potentials */
	g_maxpot = -1.0;
  g_minpot = 1.0;
  /* normalize each eigenvector such that the absolute max
    value is 1 */
	 if((max_value_array=( MY_DOUBLE * ) calloc( (size_t) degree_of_freedom, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
   }

  for(i=0; i<unknowns;i++) {
  	for (k=0;k< degree_of_freedom;k++) {
			  if ( ABS(v[i][k])> max_value_array[k]) max_value_array[k]=ABS(v[i][k]);
		}
  } 
 /* find reciprocals of max values for division */
  for(k=0;k<degree_of_freedom;k++) {
    if(!IS_ZERO(max_value_array[k]))
    max_value_array[k]=1/max_value_array[k];
  }

  /* after finding the max values, divide each column 
    of v by its max abs value */
  for(i=0; i<unknowns;i++) {
  	for (k=0;k< degree_of_freedom;k++) {
			  v[i][k] *=max_value_array[k];
		}
  } 

  /* reallocate memory for z array and copy all eigenvectors */
  for ( i = 0; i < unknowns; i++ ) {
		if (((M.Narray[re_renumber[i]])->z=(MY_DOUBLE *)realloc((void*)(M.Narray[re_renumber[i]])->z,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
     fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
		 exit(2);
		}
    memcpy((void*)(M.Narray[re_renumber[i]])->z,(void*)v[i],(size_t)degree_of_freedom*sizeof(MY_DOUBLE));

  }
 /* now consider fixed points: fill their z arrays with zeros */
  for ( i = unknowns; i < total; i++ ) {
		if (((M.Narray[re_renumber[i]])->z=(MY_DOUBLE *)realloc((void*)(M.Narray[re_renumber[i]])->z,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
     fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
		 exit(2);
		}
   /* fill zeros */
    memset((void*)(M.Narray[re_renumber[i]])->z,0,(size_t)degree_of_freedom*sizeof(MY_DOUBLE));
  }



  /* now we can free the renumber[] array */
  free(renumber);
  /* now free this array too */
  free(re_renumber);
  
  /* free everything else */
  
  /* first the arrays */
  free(x);
  free(y);
  
  /* next the matrix */
  for ( i = 0; i < unknowns ; i++) {
    free(P[i]);
		free(v[i]);
    free(T[i]);
	}
  free(P);
	free(v);
  free(max_value_array);
  free(T);
 
  /* now printout a summery */
#ifdef DEBUG
  for ( i = 0; i < total; i++ )
    printf("renumber_and_solve: %d  "MDF"  "MDF" "MDF"\n",i,Mx(i,M),My(i,M),Mz(i,M));
#endif

			 return(0);
}

/* solve helmholtz equation using sparce matrices */
int 
re_number_and_solve_helmholtz_sparse(Mesh Ms, DAG_tree Dt)
{
 unsigned MY_INT i,j;
 MY_INT k;
 MY_INT * renumber; /* re-numbering array */
 MY_INT * re_renumber; /* re-renumbering array */
/* used to map potentials back to the original points */


  unsigned MY_INT unknowns=0; /* unknown points */
  unsigned MY_INT total = 0; /* total points */
  MY_INT triangles = 0; /* known points */


  triangle *t;
  
  /* first the global node data, all total size */
  MY_DOUBLE *x ; /* x coordinates */
  MY_DOUBLE *y ; /* y coordinates */
  MY_DOUBLE **v; /* eigenvector matrix */
  
  MY_DOUBLE *max_value_array;
  /* the global variables */
  SPMat P ; /* size unknowns x unknowns note we do not need a NN x NN matrix, only the unknowns
		     */
  SPMat T; /* used in Helmholtz equation */


  /* allocate memory for arrays */
  if ((renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT )))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((re_renumber = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  /* need to solve for only points referred by triangles */
  /* SO -- */
		DAG_traverse_list_reset(&Dt);
		t=DAG_traverse_prune_list(&Dt);
  while ( t ) {
    if ( t && t->status==LIVE ) {
      /* mark the points in this triangle */
      /* use re_renumber[] array for this */
      re_renumber[t->p0] = 1;
      re_renumber[t->p1] = 1;
      re_renumber[t->p2] = 1;

      triangles++;
    }
		t=DAG_traverse_prune_list(&Dt);
  }

  /* now find unknowns, or INSIDE points */
  j = 0;
  for ( i = 0; i < Ms.count; i++ )
    if (  Mval(i,Ms) == VAR  && re_renumber[i] == 1) {
      renumber[i] = j++; 
#ifdef DEBUG
		printf("renumber_and_solve: renumbering variable point %d as %d\n",i,renumber[i]);
#endif
      unknowns++;
    }
  /* now renumber the known points */
  total = unknowns;
  for ( i = 0; i < Ms.count; i++ )
    if ( Mval(i,Ms) == FX  && re_renumber[i]==1 ) {
      renumber[i] = j++; 
#ifdef DEBUG
		printf("renumber_and_solve: renumbering fixed point %d as %d\n",i,renumber[i]);
#endif
      total++;
    }
 /* finally give numbers to points 0,1,2 */
	/* we will not use these points */
		renumber[0]=j++;
		renumber[1]=j++;
		renumber[2]=j++;

/* now construct re-renmbering array */
		  i = 0;
				/* do an insertion sort */
				while ( i < Ms.count ) {
						for ( j = 0; j < Ms.count; j++ ) {
						   if ( renumber[j] == (int)i )
									  break;
						}
										re_renumber[i] = j;
										i++;
				}

#ifdef DEBUG
printf("build: renumbering %d points with %d unknowns out of %d total\n",total,unknowns,Ms.count);
#endif
  /*renumber_points(renumber, re_renumber, total, unknowns); */
  
#ifdef DEBUG
  printf("i,  renum, rerenum\n");
  for ( i = 0; i < total; i++ )
    printf("%d %d %d\n",i,renumber[i],re_renumber[i]);
  for ( i = 0; i < total; i++ )
    printf("%d  "MDF"  "MDF" "MDF"\n",i,Mx(i,Ms),My(i,Ms),Mz(i,Ms));
#endif

  
  /* we free these two arrays only after solving the problem */
  /* now solve the problem */
  /* first allocate memory for the arrays */
  if((x=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

  if((y=(MY_DOUBLE *)calloc(total, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

	/* structure to store all eigenvectors - matrix */
  if((v=(MY_DOUBLE **)calloc(unknowns, sizeof(MY_DOUBLE*)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 for (i = 0; i <unknowns ; i++)
	  if((v[i] = ( MY_DOUBLE * ) calloc( (size_t) degree_of_freedom, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }


	
	/* now the matrices */
  SPM_init(&P,unknowns);
  SPM_init(&T,unknowns);
#ifdef DEBUG	
	printf("solving problem with %d unknowns, %d triangles\n",unknowns,triangles);
#endif

/*build_helmholtz( MY_DOUBLE **P, MY_DOUBLE **T,  MY_DOUBLE *x, MY_DOUBLE *y, int NUN, int NE, int NN)  */
  build_helmholtz_sparse(&P,&T,x,y, unknowns, triangles, total, renumber, re_renumber);

/*solve_helmholtz(MY_DOUBLE **P,MY_DOUBLE **T, MY_DOUBLE **x,  MY_INT  NUN, MY_INT n_to_find) */
#if defined (USE_LAPACK) || defined (USE_MKL)
	solve_helmholtz_lapack_sparse(&P,&T,v,eig_array,unknowns,degree_of_freedom);
#else 
	solve_helmholtz_sparse(&P,&T,v,eig_array,unknowns,degree_of_freedom);
#endif /* !USE_LAPACK || USE_MKL */
	/* all eigenvectors are normalized, hence normlization has to be done individually */
  	/* we copy back elements of eigenvectors to z arrays of each point */
	/* in order to do that, we need to change size of those z arrays */
  /* find ranges for potentials */
	g_maxpot = -1.0;
  g_minpot = 1.0;
  /* normalize each eigenvector such that the absolute max
    value is 1 */
	 if((max_value_array=( MY_DOUBLE * ) calloc( (size_t) degree_of_freedom, sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
   }

  for(i=0; i<unknowns;i++) {
  	for (k=0;k< degree_of_freedom;k++) {
			  if ( ABS(v[i][k])> max_value_array[k]) max_value_array[k]=ABS(v[i][k]);
		}
  } 
 /* find reciprocals of max values for division */
  for(k=0;k<degree_of_freedom;k++) {
    if(!IS_ZERO(max_value_array[k]))
    max_value_array[k]=1/max_value_array[k];
  }

  /* reallocate memory for z array and copy all eigenvectors */
  for ( i = 0; i < unknowns; i++ ) {
		if (((Ms.Narray[re_renumber[i]])->z=(MY_DOUBLE *)realloc((void*)(Ms.Narray[re_renumber[i]])->z,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
     fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
		 exit(2);
		}
    memcpy((void*)(Ms.Narray[re_renumber[i]])->z,(void*)v[i],(size_t)degree_of_freedom*sizeof(MY_DOUBLE));

				/* find current limits */
				for (k=0;k< degree_of_freedom;k++) {
										if ( v[i][k]> g_maxpot ) g_maxpot=v[i][k];
										if ( v[i][k]< g_minpot ) g_minpot=v[i][k];
				}
  }
 /* now consider fixed points: fill their z arrays with zeros */
  for ( i = unknowns; i < total; i++ ) {
		if (((Ms.Narray[re_renumber[i]])->z=(MY_DOUBLE *)realloc((void*)(Ms.Narray[re_renumber[i]])->z,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0) {
     fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
		 exit(2);
		}
   /* fill zeros */
    memset((void*)(Ms.Narray[re_renumber[i]])->z,0,(size_t)degree_of_freedom*sizeof(MY_DOUBLE));
  }



  /* now we can free the renumber[] array */
  free(renumber);
  /* now free this array too */
  free(re_renumber);
  
  /* free everything else */
  
  /* first the arrays */
  free(x);
  free(y);
  
  /* next the matrix */
  for ( i = 0; i < unknowns ; i++) {
		free(v[i]);
	}
	free(v);
  free(max_value_array);
  SPM_destroy(&P);
  SPM_destroy(&T);
 
  /* now printout a summery */
#ifdef DEBUG
  for ( i = 0; i < total; i++ )
    printf("renumber_and_solve: %d  "MDF"  "MDF" "MDF"\n",i,Mx(i,Ms),My(i,Ms),Mz(i,Ms));
#endif

  return(0);
}