// NewtonMethod.cpp: implementation of the CNewtonMethod class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "DiagramEng.h"
#include "NewtonMethod.h"
#include "DiagramEngDoc.h"	// rqd. since GetDocumentGlobalFn() called within funcv()
#include <math.h>			// rqd. for fabs() in fdjac()
#include "assert.h"			// rqd. for assert() used in newt()-based code.

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CNewtonMethod::CNewtonMethod()
{
	// Init the fvec that is "communicated" by the fns within the class.
	fvec = NULL;
}

CNewtonMethod::CNewtonMethod(int t_cnt, double t, double delta_t)
{
	// Init time-based parameters
	m_t_cnt = t_cnt;
	m_t = t;
	m_delta_t = delta_t;

	// Init the fvec that is "communicated" by the fns within the class.
	fvec = NULL;
}

CNewtonMethod::~CNewtonMethod()
{
	// Delete member variable fvec. However, it should already have been deleted in memory_delete_newt().
	if(fvec != NULL)
	{
		delete [] fvec;
		fvec = NULL;
	}
}

// Finite Difference Jacobian Method: Numerical Recipes in C: Section 9.7
//void CNewtonMethod::fdjac(int n, float x[], float fvec[], float **df, void (*vecfunc)(int, float [], float []))
void CNewtonMethod::fdjac(int n, double *x, double *fvec, double **df)
{
	int i;
	int j;
	double h;
	double temp;
	double *f = NULL;

	// MEMORY ALLOC
	f = memory_dvector(n);

	for(j=0; j<n; j++)
	{
		temp = x[j];
		h = EPS*fabs(temp);

		if(h == 0.0)
		{
			h = EPS;			// original NR code
			//h = (double)EPS;	// BF altered
		}

		x[j] = temp + h;	// trick to reduce finite precision error
		h = x[j] - temp;

//		(*vecfunc)(n,x,f);
		funcv(n, x, f);

		x[j] = temp;

		for(i=0; i<n; i++)
		{
			df[i][j] = (f[i] - fvec[i])/h;	// forward difference formula
		}

	}// end for j

	// FREE MEMORY
	if(f != NULL)
	{
		delete [] f;	
		f = NULL;
	}
}

// Function Min: Numerical Recipes in C: Section 9.7
double CNewtonMethod::fmin(double *x, int n)
{
	int i;
	double sum;

//	(CNewtonMethod::*nrfuncv)(n, x, fvec);
	funcv(n, x, fvec);

	for(sum = 0.0, i=0; i<n; i++)
	{
		sum += SQR(fvec[i]);
	}

	return 0.5*sum;
}

void CNewtonMethod::funcv(int n, double *x, double *f)
{
	int i;
	int valid;
	double *F_vec = NULL;	// vector that receives F(x) of F(x) = x - f(x), when working with CSystemModel fns.
	double *f_vec = NULL;	// vector that receives f(x) of F(x) = x - f(x), when working with CSystemModel fns.
	double *x_vec = NULL;	// vector that receives x of F(x) = x - f(x), when working with CSystemModel fns.
	CString	sMsg;			// string to be displayed
	CString sMsgTemp;		// temp msg.
	UINT nType = MB_OK;		// style of msg. box
	UINT nIDhelp = 0;		// help context ID for the msg.
	CDiagramEngDoc *pDoc = GetDocumentGlobalFn();

	// MEMORY ALLOC
	F_vec = memory_dvector(n);	// F_vec = F(x) = x - f(x) in the CSystemModel centric code, size = nrows
	f_vec = memory_dvector(n);	// f_vec = f(x), of F(x) = x - f(x), in the CSystemModel centric code, size = nrows
	x_vec = memory_dvector(n);	// x_vec = x, of F(x) = x - f(x), in the CSystemModel centric code, size = nrows

	// Assign incoming x vector of the CNewtonMethod centric code, to x_vec of the CSystemModel centric code.
	for(i=0; i<n; i++)	
	{
		x_vec[i] = x[i];
	}

	// DETERMINE SYSTEM EQUATION VECTOR f_vec (f(x) of CSystemModel centric code)
	// Retrieve the functional operation values, given the initial output signal vector used by the block OperateOnData() fns.
	// f_vec = f(x) = [f0(x), f1(x), ... , fs-1(x)], where s = no. of blocks that generate output signals.
	valid = pDoc->GetSystemModel().DetermineSystemEquationsFunctionOperationVector(m_t_cnt, m_t, m_delta_t, x_vec, f_vec); 
	if(valid == 0)
	{
		sMsgTemp.Format("\n CNewtonMethod::funcv()\n");
		sMsg+=sMsgTemp;
		sMsgTemp.Format("\n WARNING! Unsuccessful determination of f(x)!\n");
		sMsg+=sMsgTemp;
		AfxMessageBox(sMsg, nType, nIDhelp);	
	}
		
	// BUILD F(x) = x - f(x) (= 0)
	for(i=0; i<n; i++)
	{
		F_vec[i] = x_vec[i] - f_vec[i];
	}

	// CONVERT F(x) of the CSystemModel centric code, to f(x) of the CNewtonSolver centric code
	for(i=0; i<n; i++)
	{
		f[i] = F_vec[i];
	}
	
	// Print vectors to check values
	/*sMsgTemp.Format("\n Printing x_vec \n");
	sMsg+=sMsgTemp;
	AfxMessageBox(sMsg, nType, nIDhelp);
	PrintVector(x_vec, n);
	
	sMsgTemp.Format("\n Printing f_vec \n");
	sMsg+=sMsgTemp;
	AfxMessageBox(sMsg, nType, nIDhelp);
	PrintVector(f_vec, n);

	sMsgTemp.Format("\n Printing F_vec \n");
	sMsg+=sMsgTemp;
	AfxMessageBox(sMsg, nType, nIDhelp);
	PrintVector(F_vec, n);
	*/
		
	// MEMORY DELETE
	if(F_vec != NULL)
	{
		delete [] F_vec;
		F_vec = NULL;
	}
	if(f_vec != NULL)
	{
		delete [] f_vec;
		f_vec = NULL;
	}
	if(x_vec != NULL)
	{
		delete [] x_vec;
		x_vec = NULL;
	}
}

// Line Search Method: Numerical Recipes in C: Section 9.7
//void CNewtonMethod::lnsrch(int n, float xold[], float fold, float g[], float p[], float x[], float *f, float stpmax, int *check, float (*func)(float [], int n))
int CNewtonMethod::lnsrch(int n, double *xold, double fold, double *g, double *p, double *x, double *f, double stpmax, int *check)
{
	int i;
	int nr_status = 1;	// value indicating the nr method status: -2 (numerical error), -1 (memory allocation error), 1 (no error)
	double a;
	double alam;
	double alam2;
	double alamin;
	double b;
	double disc;
	double f2;
	double rhs1;
	double rhs2;
	double slope;
	double sum;
	double temp;
	double test;
	double tmplam;

	*check = 0;
	for(sum=0.0, i=0; i<n; i++)
	{
		sum += p[i]*p[i];
	}

	// Scale if attempted step is too big.
	sum = sqrt(sum);
	if(sum > stpmax)
	{
		for(i=0; i<n; i++)
		{
			p[i] *= stpmax/sum;
		}
	}

	for(slope = 0.0, i=0; i<n; i++)
	{
		slope += g[i]*p[i];
	}

	if(slope >= 0.0)
	{
		// DiagramEng-altered
		nr_status = -2;
		nrerror("Roundoff problem in lnsrch().", nr_status);
		return nr_status;
	}

	// Compute lambda min.
	test = 0.0;
	for(i=0; i<n; i++)
	{
		temp = fabs(p[i])/FMAX(fabs(xold[i]), 1.0);

		if(temp > test)
		{
			test = temp;
		}
	}

	// Always try full Newton step first
	// Assertion (BF altered)
	assert(test != 0.0);		
	alamin = TOLX/test;
	alam = 1.0;

	// Start of iteration loop
	for(;;)
	{
		for(i=0; i<n; i++)
		{
			x[i] = xold[i] + alam*p[i];
		}
//		*f = (*func)(x, n);
		*f = fmin(x, n);

		// Convergence on delta_x. For zero finding, the calling program should verify cgce.
		if(alam < alamin)
		{
			for(i=0; i<n; i++)
			{
				x[i] = xold[i];
			}
			*check = 1;

			// DiagramEng-altered
			return nr_status;
		}	
		else if(*f <= fold + ALF*alam*slope)	// Sufficient fn decrease.
		{
			// DiagramEng-altered
			return nr_status;
		}
		else									// Backtrack. 
		{	
			if(alam == 1.0)
			{
				// Assertion (BF altered)
				assert((*f - fold - slope) != 0.0);
				tmplam = -slope/(2.0*(*f - fold - slope));
			}
			else
			{
				// Assertion (BF altered)
				assert(alam != 0.0);
				assert(alam2 != 0.0);
				assert((alam - alam2) != 0.0);

				rhs1 = *f - fold - alam*slope;
				rhs2 = f2 - fold - alam2*slope;
				a = (rhs1/(alam*alam) - rhs2/(alam2*alam2))/(alam - alam2);
				b = (-alam2*rhs1/(alam*alam) + alam*rhs2/(alam2*alam2))/(alam - alam2);

				if(a == 0.0)
				{
					// Assertion (BF altered)
					assert(b != 0.0);		
					tmplam = -slope/(2.0*b);
				}
				else
				{
					disc = b*b - 3.0*a*slope;
					if(disc < 0.0)
					{
						tmplam = 0.5*alam;
					}
					else if(b <= 0.0)
					{
						// Assertion (BF altered)
						assert(a != 0.0);
						tmplam = (-b + sqrt(disc))/(3.0*a);
					}
					else
					{	
						// Assertion (BF altered)
						assert((b + sqrt(disc)) != 0.0);
						tmplam = -slope/(b + sqrt(disc));
					}
				}// end else

				if(tmplam > 0.5*alam)
				{
					tmplam = 0.5*alam;		// lambda <= 0.5*lambda_1
				}

			}// end else

		}// end else

		alam2 = alam;
		f2 = *f;
		alam = FMAX(tmplam, 0.1*alam);		// lambda >= 0.1*lambda_1

	}// end for

	// DiagramEng-altered
	return nr_status;
}

// Lower Upper Decomposition with Back Substitution: Numerical Recipes in C: Section 2.3
void CNewtonMethod::lubksb(double **a, int n, int *indx, double *b)
{
	int i;
	int ii = 0;
	int ip;
	int j;
	double sum;

	for(i=0; i<n; i++)
	{
		ip = indx[i];
		sum = b[ip];
		b[ip] = b[i];

		if(ii != 0)
		{
			for(j=ii-1; j<i; j++)
			{
				sum -= a[i][j]*b[j];
			}
		}
		else if(sum != 0.0)
		{
			ii = i + 1;
		}
		b[i] = sum;
	}

	for(i=n-1; i>=0; i--)
	{
		sum = b[i];
		for(j=i+1; j<n; j++)
		{
			sum -= a[i][j]*b[j];
		}
		
		// Assertion (BF altered)
		assert(a[i][i] != 0.0);	
		AssertManual(a[i][i] != 0.0);
		b[i] = sum/a[i][i];
	}
}

// Lower Uppder Decomposition: Numerical Recipes in C: Section 2.3
int CNewtonMethod::ludcmp(double **a, int n, int *indx, double *d)
{
	int i;
	int imax;
	int j;
	int k;
	int nr_status = 1; // value indicating the nr method status: -2 (numerical error), -1 (memory allocation error), 1 (no error)
	double big;
	double dum;
	double sum;
	double temp;
	double *vv = NULL;

	// Memory alloc
	vv = memory_dvector(n);

	*d = 1.0;
	for(i=0; i<n; i++)
	{
		big = 0.0;
		for(j=0; j<n; j++)
		{
			if((temp = fabs(a[i][j])) > big)
			{
				big = temp;
			}
		}
		if(big == 0.0)
		{
			// DiagramEng-altered
			nr_status = -2;
			nrerror("Singular matrix in routine ludcmp().", nr_status);
			if(vv != NULL)
			{
				delete [] vv;
				vv = NULL;
			}
			return nr_status;
		}
		vv[i] = 1.0/big;
	}// end for i

	for(j=0; j<n; j++)
	{
		for(i=0; i<j; i++)
		{
			sum = a[i][j];
			for(k=0; k<i; k++)
			{
				sum -= a[i][k]*a[k][j];
			}
			a[i][j] = sum;
		}

		big = 0.0;
		for(i=j; i<n; i++)
		{
			sum = a[i][j];
			for(k=0; k<j; k++)
			{
				sum -= a[i][k]*a[k][j];
			}
			a[i][j] = sum;
			if((dum = vv[i]*fabs(sum)) >= big)
			{
				big = dum;
				imax = i;
			}
		}

		if(j != imax)
		{
			for(k=0; k<n; k++)
			{
				dum = a[imax][k];
				a[imax][k] = a[j][k];
				a[j][k] = dum;
			}
			*d = -(*d);
			vv[imax] = vv[j];
		}
		indx[j] = imax;

		if(a[j][j] == 0.0)
		{
			a[j][j] = TINY;				// original NR code
			//a[j][j] = (double)TINY;	// BF altered.
		}

		if(j != (n-1))
		{
			dum = 1.0/(a[j][j]);
			for(i=j+1; i<n; i++)
			{
				a[i][j] *= dum;
			}
		}
	}// end for j

	// Memory delete
	if(vv != NULL)
	{
		delete [] vv;
		vv = NULL;
	}

	// DiagramEng-altered
	return nr_status;
}

// Newton Method: Numerical Recipes in C: Section 9.7
//void CNewtonMethod::newt(float x[], int n, int *check, void(*vecfunc)(int, float[], float[]))
int CNewtonMethod::newt(double *x, int n, int *check)
{
	int i;
	int its;
	int j;
	int nr_status = 1; // value indicating the nr method status: -2 (numerical error), -1 (memory allocation error), 1 (no error)
	int *indx = NULL;
	double d;
	double den;
	double f;
	double fold;
	double stpmax;
	double sum;
	double temp;
	double test;
	double **fjac = NULL;
	double *g = NULL;
	double *p = NULL;
	double *xold = NULL;

	// MEMORY ALLOC
	indx = memory_ivector(n);
	g = memory_dvector(n);
	p = memory_dvector(n);
	xold = memory_dvector(n);
	fvec = memory_dvector(n);	// global var used (globally) in newt() and fmin(), and passed as an arg to fdjac() and funcv().
	fjac = memory_dmatrix(n,n);
	
	// define global variable
//	nrfuncv = vecfunc;

	f = fmin(x, n);				// fvec is also computed by this call

	test = 0.0;					// test fo initial guess being a root. Use a more stringent test than simply TOLF.

	for(i=0; i<n; i++)
	{
		if(fabs(fvec[i]) > test)
		{
			test = fabs(fvec[i]);
		}
	}

	if(test < 0.01*TOLF)
	{
		*check = 0;
		memory_delete_newt(indx, g, p, xold, fjac, n);
		return nr_status;
	}

	for(sum=0.0, i=0; i<n; i++)
	{
		sum += SQR(x[i]);		// calculate stpmax for line searches
	}
	stpmax = STPMX*FMAX(sqrt(sum), (double)n);

	// Start of iteration loop
	for(its = 0; its < MAXITS; its++)
	{
//		fdjac(n, x, fvec, fjac, vecfunc);
		fdjac(n, x, fvec, fjac);

		// -- DiagramEng-altered test of fjac
		nr_status = CheckJacobian(fjac, n, n);
		if(nr_status == -2)
		{
			nrerror("Numerical problem concerning fjac() detected in newt().", nr_status);
			memory_delete_newt(indx, g, p, xold, fjac, n);
			return nr_status;
		}

		for(i=0; i<n; i++)
		{
			for(sum=0.0, j=0; j<n; j++)
			{
				sum += fjac[j][i]*fvec[j];
			}
			g[i] = sum;
		}

		for(i=0; i<n; i++)
		{
			xold[i] = x[i];		// store x
		}
		fold = f;				// store f

		for(i=0; i<n; i++)
		{
			p[i] = -fvec[i];	// rhs for linear equs
		}

		// -- DiagramEng-altered
		nr_status = ludcmp(fjac, n, indx, &d);	// solve linear equs by LU decomp
		if(nr_status == -2)
		{
			nrerror("Numerical problem concerning ludcmp() detected in newt().", nr_status);
			memory_delete_newt(indx, g, p, xold, fjac, n);
			return nr_status;
		}

		lubksb(fjac, n, indx, p);

		// lnsrch returns new x and f. It also calculates fvec at the new x when it calls fmin.
//		lnsrch(n, xold, fold, g, p, x, &f, stpmax, check, fmin);

		// -- DiagramEng-altered
		nr_status = lnsrch(n, xold, fold, g, p, x, &f, stpmax, check);
		if(nr_status == -2)
		{
			nrerror("Numerical problem concerning lnsrch() detected in newt().", nr_status);
			memory_delete_newt(indx, g, p, xold, fjac, n);
			return nr_status;
		}

		test = 0.0;		// test for cgce on fn vals

		for(i=0; i<n; i++)
		{
			if(fabs(fvec[i]) > test)
			{
				test = fabs(fvec[i]);
			}
		}

		if(test < TOLF)
		{
			*check = 0;
			memory_delete_newt(indx, g, p, xold, fjac, n);
			return nr_status;
		}

		// Check for gradient of f zero, i.e. spurious cgce.
		if(*check)
		{
			test = 0.0;
			den = FMAX(f, 0.5*n);
			for(i=0; i<n; i++)
			{
				temp = fabs(g[i])*FMAX(fabs(x[i]), 1.0)/den;
				if(temp > test)
				{
					test = temp;
				}
			}
			*check = (test < TOLMIN ? 1 : 0);
			memory_delete_newt(indx, g, p, xold, fjac, n);
			return nr_status;
		}

		// Test for cgce on delta x
		test = 0.0;
		for(i=0; i<n; i++)
		{
			temp = (fabs(x[i] - xold[i]))/FMAX(fabs(x[i]), 1.0);
			if(temp > test)
			{
				test = temp;
			}
		}

		if(test < TOLX)
		{
			memory_delete_newt(indx, g, p, xold, fjac, n);
			return nr_status;
		}

	}// end "for its"

	// DiagramEng-altered
	nr_status = -2;
	nrerror("MAXITS exceeded in newt().", nr_status);
	memory_delete_newt(indx, g, p, xold, fjac, n);
	return nr_status;
}

// DiagramEng altered (nr_status)
void CNewtonMethod::nrerror(char *error_text, int nr_status)
{
	// Numerical Recipes standard error handler
	CString sMsg;		// string to be displayed
	CString sMsgTemp;	// temp msg.
	UINT nType = MB_OK;	// style of msg. box
	UINT nIDhelp = 0;	// help context ID for the msg.

	// NR Impl.
	//fprintf(stderr, "Numerical Recipes run-time error ...\n");
	//fprintf(stderr, "%s\n", error_text);
	//fprintf(stderr, "...now exiting to system...\n");

	// DiagramEng-altered
	// Start Message
	sMsgTemp.Format(" Numerical Recipes run-time error ...\n\n");
	sMsg+=sMsgTemp;
	sMsgTemp.Format(" %s\n", error_text);
	sMsg+=sMsgTemp;

	// Filter nr_status flag
	// nr_status is the value indicating the nr method status: -2 (numerical error), -1 (memory allocation error), 1 (no error)
	if(nr_status == -2)				// Print message indicating a numerical error.
	{
		AfxMessageBox(sMsg, nType, nIDhelp);
	}
	else if(nr_status == -1)		// Exit system only if there is a memory deallocation failure
	{
		sMsgTemp.Format(" Now exiting to system. \n");
		sMsg+=sMsgTemp;
		AfxMessageBox(sMsg, nType, nIDhelp);
		exit(1);
	}
	else
	{
		// no action for now
	}
}

int CNewtonMethod::CheckJacobian(double **fjac, int nrows, int ncols)
{
	int i;
	int j;
	int nr_status = 1; // nr method status: -2 (numerical error), -1 (memory allocation error), 1 (no error)
	double big;
	double temp;
	
	// Check for zero row i.e. a singular Jacobian matrix
	for(i=0; i<nrows; i++)
	{
		big = 0.0;
		for(j=0; j<ncols; j++)
		{
			if((temp = fabs(fjac[i][j])) > big)
			{
				big = temp;
			}
		}
		if(big == 0.0)
		{
			PrintMatrix(fjac,nrows,ncols);
			nr_status = -2;
			nrerror("CheckJacobian(): Jacobian fjac is singular.", nr_status);
			break;	// break from for loop since zero row is found
		}
	}

	return nr_status;
}

double* CNewtonMethod::memory_dvector(int n)
{
	double *v = NULL;
	
	// Memory Alloc
	v = new double[n];

	// Check memory
	if(!v)
	{
		//nrerror("memory_dvector(), memory allocation failure.");
		nrerror("memory_dvector(), memory allocation failure.", -1);
	}

	return v;
}

int* CNewtonMethod::memory_ivector(int n)
{
	int *v = NULL;
	
	// Memory Alloc
	v = new int[n];

	// Check memory
	if(!v)
	{
		//nrerror("memory_ivector(), memory allocation failure.");
		nrerror("memory_ivector(), memory allocation failure.", -1);
	}

	return v;
}

double** CNewtonMethod::memory_dmatrix(int nrows, int ncols)
{
	int i;
	double **m = NULL;

	// Memory alloc of primary array
	// Allocate an array of ptrs-to-double and return the & to m.
	// Note: an array is of type ptr, hence an array of ptrs is of type ptr-ptr.
	m = new double *[nrows];		

	// Check new
	if(!m)
	{
		//nrerror("memory_dmatrix(), memory allocation failure (array 1).");
		nrerror("memory_dmatrix(), memory allocation failure (array 1).", -1);
	}

	// Memory alloc of secondary arrays
	// Allocate an array of double for each of the entries (rows of the column vector) of matrix, and return the & to m[i].
	// Note: an array is of type ptr, and its & is stored at location i in the array m[]
	for(i=0;i<nrows;i++)
	{
		m[i] = new double[ncols];	// allocate an array of doubles, return & to m[i]

		// Check new
		if(!m[i])
		{
			//nrerror("memory_dmatrix(), memory allocation failure (array 2).");
			nrerror("memory_dmatrix(), memory allocation failure (array 2).", -1);
		}
	}
	
	return m;
}

void CNewtonMethod::memory_delete_newt(int *indx, double *g, double *p, double *xold, double **fjac, int nrows)
{
	int i;

	// Memory Delete Vectors
	if(indx != NULL)
	{
		delete [] indx;
		indx = NULL;
	}
	if(g != NULL)
	{
		delete [] g;
		g = NULL;
	}
	if(p != NULL)
	{
		delete [] p;
		p = NULL;
	}
	if(xold != NULL)
	{
		delete [] xold;
		xold = NULL;
	}

	// WARNING! fvec is a class member variable: if it is declared locally in the fn header and fvec passed into the
	// memory_delete_newt() fn from the newt() calling environment, then the assignment fvec = NULL would not be 
	// persistent outside of this fn as it would be a local statement. 
	// Hence, fvec should not be passed into memory_delete_newt(), but rather simply deleted and reassigned to NULL
	// as it is a class member variable and memory_delete_newt() has access to it as a member fn.
	if(fvec != NULL)
	{
		delete [] fvec;
		fvec = NULL;
	}

	// Memory Delete Matrix
	if(fjac != NULL)
	{
		for(i=0; i<nrows; i++)
		{
			delete [] fjac[i];
		}
		delete [] fjac;
		fjac = NULL;
	}

	return;
}

int CNewtonMethod::PrintVector(double *v, int ncols)
{
	int	i;
	CString sMsg;		// main msg string
	CString sMsgTemp;	// temp msg string
	UINT nType = MB_OK;	// style of msg. box
	UINT nIDhelp = 0;	// help context ID for the msg.

	// If the vector is NULL then return.
	if(v == NULL)
	{
		return 1;
	}

	sMsgTemp.Format("\n PrintVector()\n");
	sMsg += sMsgTemp;
	
	// Print Vector
	sMsgTemp.Format("\n");
	sMsg += sMsgTemp;

	for(i=0; i<ncols; i++)
	{
		//sMsgTemp.Format(" %-10.3e", v[i]);	// left adjusted, total width 10, 3 decimal places, exponential notation
		sMsgTemp.Format(" %lf", v[i]);			// standard double formatting
		sMsg += sMsgTemp;
	}
	sMsgTemp.Format("\n");
	sMsg += sMsgTemp;

	// Display msg. box
	AfxMessageBox(sMsg, nType, nIDhelp);

	return 0;
}

int CNewtonMethod::PrintMatrix(double **M, int nrows, int ncols)
{
	int	i;
	int	j;
	CString sMsg;		// main msg string
	CString sMsgTemp;	// temp msg string
	UINT nType = MB_OK;	// style of msg. box
	UINT nIDhelp = 0;	// help context ID for the msg.

	// If the matrix is NULL then return.
	if(M == NULL)
	{
		return 1;
	}

	sMsgTemp.Format("\n PrintMatrix()\n");
	sMsg += sMsgTemp;
	
	// Print Matrix
	sMsgTemp.Format("\n");
	sMsg += sMsgTemp;
	
	for(i=0; i<nrows; i++)
	{
		for(j=0; j<ncols; j++)
		{
			//sMsgTemp.Format(" %-10.3e", M[i][j]); // left adjusted, total width 10, 3 decimal places, exponential notation
			sMsgTemp.Format(" %lf", M[i][j]);		// standard double formatting
			sMsg += sMsgTemp;
		}
		sMsgTemp.Format("\n");
		sMsg += sMsgTemp;
	}
	sMsgTemp.Format("\n");
	sMsg += sMsgTemp;

	// Display msg. box
	AfxMessageBox(sMsg, nType, nIDhelp);

	return 0;
}

void CNewtonMethod::CheckDenomZero(double denom, char *error_text)
{
	int nr_status = 1; // value indicating nr method status: -2 (numerical error), -1 (memory allocation error), 1 (no error)

	// If denominator = 0, return flt_min.
	if(denom == 0.0)
	{
		nr_status = -2;			// numerical error due to an attempted divide by zero
		nrerror(error_text, nr_status);
	}

	return;
}
