/* <OCN_g03.h>
 *
 * Copyright 1996 Numerical Algorithms Group
 *
 * Include file for NAG C Library g03 Chapter
 *
 * Mark 5, 1997.
 * Mark 6 revised. IER-3042 (July 2000).
 */

#ifndef NAGG03
#define NAGG03

//#importdll "ONAG" // NAG DLL prepared by OriginLab
#pragma dll(ONAG)	
#include <NAG\nag_types.h>


/* begin proto */

/** g03aac
		performs a principal component analysis on a data matrix; both the
		principal component loadings and the principal component scores are returned.

Example 1:
	Assume we have four Matrix, "Matrix1"'s size is 10X3 and stores 30 number; "Matrix2"'s
	size is 3X6; "Matrix3"'s size is 3X3, and "Matrix4"'s size is 12X3.  We will store the 
	value to "Matrix2", "Matrix3" and "Matrix4".
	
	matrix p;
	p.SetSize(3,3);
	double s[3];
	matrix e;
	e.SetSize(3,6);
	matrix v;
	v.SetSize(12,3);
	double wt[12];
	double *wtptr;
	wtptr = NULL;
	int isx[3] = {1, 1, 1};
	int nvar = 3;
	int tdx = 3; 
	int tde = 6;
	int tdp = 3;
	int tdv = 3;
	int i, j, m = 3, n = 10;
	Nag_PrinCompMat pcmatrix;
	Nag_PrinCompScores scores;
	Matrix xx("Matrix1");
	matrix x = xx;
	char weight, mat, std;
	mat = 'V';
	std = 'E';
	weight = 'U';

	if (mat == 'C')
		pcmatrix = Nag_MatCorrelation;
	else if (mat == 'S')
		pcmatrix = Nag_MatStandardised;
	else if (mat == 'U')
		pcmatrix = Nag_MatSumSq;
	else
		pcmatrix = Nag_MatVarCovar;
	if (std == 'S')
		scores = Nag_ScoresStand;
	else if (std == 'U')
		scores = Nag_ScoresNotStand;
	else if (std == 'Z')
		scores = Nag_ScoresUnitVar;
	else
		scores = Nag_ScoresEigenval;
	
	nag_mv_prin_comp(pcmatrix, scores, n, m, x, tdx, isx, s, wtptr, nvar,
				e, tde, p, tdp, v, tdv);
	//put the result to the matrix.
	//first column is Eigenvalues, second column is Percentage variation, 
	//third column is Cumulative variation, fourth column is Chisp, 
	//fifth column is DF and sixth column is Sig.	
	Matrix ee("Matrix2");
	ee = e;
	Matrix pp("Matrix3");		//Eigenvalues
	pp = p;
	Matrix vv("Matrix4");		//Principal component scores
	vv = v;		
Example 2:
	A data set is taken from Cooley and Lohnes (1971), it consists of ten observations on 
	three variables.  The unweighted principal components based on the variance-covariance 
	matrix are computed and unstandardised principal component scores requested.

void test_nag_mv_prin_comp()
{
	matrix p;
	p.SetSize(3,3);
	double s[3];
	matrix e;
	e.SetSize(3,6);
	matrix v;
	v.SetSize(12,3);
	matrix x = {{7.0, 4.0, 3.0}, {4.0, 1.0, 8.0}, {6.0, 3.0, 5.0},
	{8.0, 6.0, 1.0}, {8.0, 5.0, 7.0}, {7.0, 2.0, 9.0}, {5.0, 3.0, 3.0}, 
	{9.0, 5.0, 8.0}, {7.0, 4.0, 5.0}, {8.0, 2.0, 2.0}};
	double wt[12];
	double *wtptr;
	int isx[3];
	int nvar = 3;
	int tdx = 3; 
	int tde = 6;
	int tdp = 3;
	int tdv = 3;
	int i, j, m = 3, n = 10;
	Nag_PrinCompMat pcmatrix;
	Nag_PrinCompScores scores;
	char weight, mat, std;
	
 	for(i = 0; i < 3; i++)
 		isx[i] = 1;
 
	mat = 'V';
	std = 'E';
	weight = 'U';

	if (mat == 'C')
		pcmatrix = Nag_MatCorrelation;
	else if (mat == 'S')
		pcmatrix = Nag_MatStandardised;
	else if (mat == 'U')
		pcmatrix = Nag_MatSumSq;
	else
		pcmatrix = Nag_MatVarCovar;
	if (std == 'S')
		scores = Nag_ScoresStand;
	else if (std == 'U')
		scores = Nag_ScoresNotStand;
	else if (std == 'Z')
		scores = Nag_ScoresUnitVar;
	else
		scores = Nag_ScoresEigenval;
	
	wtptr = NULL;
	
	nag_mv_prin_comp(pcmatrix, scores, n, m, x, tdx, isx, s, wtptr, nvar,
					e, tde, p, tdp, v, tdv);
					
	printf("Eigenvalues  Percentage   Cumulative   Chisq   DF   Sig\n");
	printf("              variation   variation\n\n");

	for (i = 0; i < nvar; ++i)
	{
		for (j = 0; j < 6; ++j)
			printf("%11.4f",e[i][j]);
		printf("\n");
	}
	printf("\nEigenvalues\n\n");
	for (i = 0; i < nvar; ++i)
	{
		for (j = 0; j < nvar; ++j)
			printf("%9.4f",p[i][j]);
		printf("\n");
	}
	printf("\nPrincipal component scores \n\n");
	for (i = 0; i < n; ++i)
	{
		printf("%2ld", i+1);
		for (j = 0; j < nvar; ++j)
			printf("%9.3f", v[i][j]);
		printf("\n");	
	}
}


	The output is following:
	
	Eigenvalues		Percentage 		Cumulative 		Chisq 		DF 			Sig
					variation 		variation
	8.2739 			0.6515 			0.6515 			8.6127 		5.0000 		0.1255
	3.6761 			0.2895 			0.9410 			4.1183 		2.0000 		0.1276
	0.7499 			0.0590 			1.0000 			0.0000 		0.0000 		0.0000
	
	Eigenvalues
	
	 0.1376 	0.6990   0.7017
	 0.2505 	0.6609 	-0.7075
	-0.9583 	0.2731 	-0.0842
	
	Principal component scores
	
	1 	 2.151 	 -0.173 	-0.107
	2 	-3.804   -2.887 	-0.510
	3 	-0.153   -0.987 	-0.269
	4 	 4.707 	  1.302 	-0.652
	5 	-1.294 	  2.279 	-0.449
	6 	-4.099 	  0.144 	 0.803
	7 	 1.626 	 -2.232 	-0.803
	8 	-2.114 	  3.251 	 0.168
	9 	 0.235 	  0.373 	-0.275
	10 	 2.746 	 -1.069 	 2.094
Return:
	The function returns NAG error code or 0 if no error.
	
	70:  On entry, parameter pcmatrix had an illegal value.  On entry, parameter scores had an illegal value.
	11:  On entry, m must not be less than 1: m = _value_.  On entry, n must not be less than 2: n = _value_.  On entry, nvar must not be less than 1: nvar = _value_.  On entry, tde must not be less than 6: tde = _value_.
	17:  On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m.  On entry, tdv = _value_ while nvar = _value_. These parameters must satisfy tdv = nvar.  On entry, tdp = _value_ while nvar = _value_. These parameters must satisfy tdp = nvar.  
	19:  On entry, nvar = _value_ while m = _value_. These parameters must satisfy nvar = m.
	20:  On entry, nvar = _value_ while n = _value_. These parameters must satisfy nvar < n.
	500: On entry, wt[_value_] = _value_.  Constraint: when referenced, all elements of wt must be non-negative.
	501: The number of variables, nvar in the analysis = _value_, while the number of variables included in the analysis via array isx = _value_.  Constraint: these two numbers must be the same.
	504: On entry, the standardisation element s[_value_] = _value_, while the variable to be included isx[_value_] = _value_.  Constraint: when a variable is to included, the standardisation element must be positive.  
	503: With weighted data, the e.ective number of observations given by the sum of weights = _value_, while the number of variables included in the analysis, nvar = _value_.  Constraint: effective number of observations > nvar + 1.
	427: The singular value decomposition has failed to converge.  This is an unlikely error exit.
	505: All eigenvalues/singular values are zero.  This will be caused by all the variables being constant.
	73:  Memory allocation failed.
	74:  An internal error has occurred in this function.  Check the function call and array sizes. If the function call is correct, please consult NAG for assistance.
	
	successfully call of the nag_mv_prin_comp function.
*/
int nag_mv_prin_comp(
    Nag_PrinCompMat pcmatrix,
    Nag_PrinCompScores scores,
    int n, // the number of observations
    int m, // the number of variables in the data matrix
    double x[], // contain the ith observation for the jth variable, for i = 1, 2, . . ., n; j = 1, 2, . . .,m.
    int tdx, 	// the last dimension of the array x.
    int isx[], // indicates whether or not the jth variable is to be included in the analysis.
    double s[], // the standardizations to be used, if any.
    double w[], // contain the weights to be used in the principal component analysis. 
    int nvar, 	// the number of variables in the principal component analysis.
    double e[],	// the statistics of the principal component analysis.
    int tde,	// the last dimension of the array e.
    double p[], // the first nvar columns of p contain the principal component loadings. The jth column of p contains the nvar coefficients for the jth principal component.
    int tdp,	// the last dimension of the array p 
    double v[],	// the first nvar columns of v contain the principal component scores. The jth column of v contains the n scores for the jth principal component.
    int tdv		// the last dimension of the array v.
); // Eigenvalues, Percentage Variation, Cumulative Variation, Principal Component Scores......

/** g03acc
		performs a canonical variate (canonical descrimination) analysis.

Example:
	A sample of nine observations, each consisting of three variables plus group indicator, is read in.
	There are three groups. An unweighted canonical variate analysis is performed.
	
void test_nag_mv_canon_var()
{
	
	matrix e;
	e.SetSize(3,6);
	matrix x = {{13.3, 10.6, 21.2}, {13.6, 10.2, 21.0}, {14.2, 10.7, 21.1},
	{13.4, 9.4, 21.0}, {13.2, 9.6, 20.1}, {13.9, 10.4, 19.8}, {12.9, 10.0, 20.5}, 
	{12.2, 9.9, 20.7}, {13.9, 11.0, 19.1}};	
	double wt[9];
	matrix cvm;
	cvm.SetSize(3,3);
	double tol;
	matrix cvx;
	cvx.SetSize(3,3);
	int i, j, m = 3, n = 9;
	int ng = 3;
	int nx = 3;
	int nig[3], ncv;
	int ing[9] = {1, 2, 3, 1, 2, 3, 1, 2, 3};
	int irx, isx[6];
	int tdx = 3;
	int tdc = 3;
	int tde = 6;
	char wtchar = 'U';
	Nag_Weightstype weight;

	weight = Nag_NoWeights; 
	
	for (j = 0; j < m; ++j)
		isx[j] = 1;
	tol = 1e-6;
	
	nag_mv_canon_var(weight, n, m, x, tdx, isx, nx, ing, ng, wt, nig,
					cvm, tdc, e, tde, &ncv, cvx, tdc, tol, &irx);
			
	printf("%s%2ld\n\n","Rank of x = ",irx);
	printf("      Canonical   Eigenvalues Percentage  ChiSQ       DF          SIG\n");
    printf("      Correlations            Variation\n");
	for (i = 0; i < ncv; ++i)
	{
		for (j = 0; j < 6; ++j)
			printf("%12.4f",e[i][j]);
			printf("\n");
	}
	printf("\nCanonical Coefficients for X\n");
	for (i = 0; i < nx; ++i)
	{
		for (j = 0; j < ncv; ++j)
			printf("%9.4f",cvx[i][j]);
		printf("\n");
	}
	printf("\nCanonical variate means\n");
	for (i = 0; i < ng; ++i)
	{
		for (j = 0; j < ncv; ++j)
			printf("%9.4f",cvm[i][j]);
		printf("\n");
	}
}	
	
	
	The output is following:
		
	Rank of x = 3
	Canonical 		Eigenvalues 	Percentage 		CHISQ 		DF 			SIG
	Correlations 					Variation	
	0.8826 			3.5238 			0.9795 			7.9032 		6.0000 		0.2453
	0.2623 			0.0739 			0.0205 			0.3564 		2.0000 		0.8368
	
	Canonical Coefficients for X
	-1.7070 	0.7277
	-1.3481 	0.3138
	 0.9327 	1.2199
	 
	Canonical variate means
	 0.9841 	0.2797
	 1.1805 	-0.2632
	-2.1646 	-0.0164

Return:
	The function returns NAG error code or 0 if no error. 

	70: On entry, parameter weight had an illegal value.
	11: On entry, nx must not be less than 1: nx = _value_.  On entry, ng must not be less than 2: ng = _value_.  On entry, tde must not be less than 6: tde = _value_.
	5:  On entry, tol must not be less than 0.0: tol = _value_.
	17: On entry, m = _value_ while nx = _value_.  These parameters must satisfy m = nx.  On entry, tdx = _value_ while m = _value_.  These parameters must satisfy tdx = m.  On entry, tdcvx = _value_ while ng = _value_.  These parameters must satisfy tdcvx = ng-1.  On entry, tdcvm = _value_ while nx = _value_.  These parameters must satisfy tdcvm = nx.  
	29: On entry, n = _value_, nx = _value_ and ng = _value_.  These parameters must satisfy n = nx + ng.
	112: On entry, ing[_value_] = _value_, ng = _value_.  Constraint: 1 = ing[i - 1] = ng, i = 1, 2, . . ., n.
	510: The wt array argument must not be NULL when the weight argument indicates weights.
	500: On entry, wt[_value_] = _value_.  Constraint: When referenced, all elements of wt must be non-negative.
	501: The number of variables, nx in the analysis = _value_, while number of variables included in the analysis via array isx = _value_.  Constraint: these two numbers must be the same.  
	427: The singular value decomposition has failed to converge.  This is an unlikely error exit.
	506: A canonical correlation is equal to one.  This will happen if the variables provide an exact indication as to which group every observation is allocated.
	507: Either the e.ective number of groups is less than two or the e.ective number of groups plus the number of variables, nx is greater than the the e.ective number of observations.  
	508: The rank of the variables is zero.  This will happen if all the variables are constants.
	73: Memory allocation failed.
	74: An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
		
	successfully call of the nag_mv_canon_var function.	
*/
int nag_mv_canon_var(
    Nag_Weightstype weights, //indicates the type of weights to be used in the analysis.
    int n, // the number of observations.
    int m, // the total number of variables.
    double x[], // contain the ith observation for the jth variable, for i = 1, 2, . . ., n; j = 1, 2, . . .,m.
    int tdx,	// the last dimension of the array x.
    int isx[],	// indicates whether or not the jth variable is to be included in the analysis.
    int nx,		// the number of variables in the analysis.
    int ing[],	// indicates which group the ith observation is in, for i = 1, 2, . . ., n. 
    int ng,		// The number of groups.
    double *wtptr, // if weight = Nag Weightsfreq or Nag Weightsvar then the elements of wt must contain the weights to be used in the analysis.
    int nig[],	// the number of observations in group j, for j = 1, 2, . . . , ng.
    double cvm[], // the mean of the jth canonical variate for the ith group, for i = 1, 2, . . . , ng; j = 1, 2, . . . , l; 
    int tdcvm,	// the last dimension of the array cvm.
    double e[], // the statistics of the canonical variate analysis.
    int tde,	// the last dimension of the array e.
    int *ncv,	// the number of canonical variates.
    double cvx[], // the canonical variate loadings.
    int tdcvx,	// the last dimension of the array cvx.
    double tol, // the value of tol is used to decide if the variables are of full rank.
    int *irankx // the rank of the dependent variables.
); // Canonical Coeffiecients, Canonical Correlation, Eigenvalues, Canonical Variate Means....
  
/** g03adc
		performs canonical correlation analysis upon input data matrices.
	
Example:
	A sample of nine observations with two variables in each set is read in. The second and third
	variables are x variables while the first and last are y variables. Canonical variate analysis is performed.

void test_nag_mv_canon_corr()
{
	
	matrix e, cvx, cvy;
	e.SetSize(2,6);
	cvx.SetSize(2,2);
	cvy.SetSize(2,2);
	matrix z = {{80.0, 58.4, 14.0, 21.0}, {75.0, 59.2, 15.0, 27.0}, {78.0, 60.3, 15.0, 27.0},
	{75.0, 57.4, 13.0, 22.0}, {79.0, 59.5, 14.0, 26.0}, {78.0, 58.1, 14.5, 26.0},
	{75.0, 58.0, 12.5, 23.0}, {64.0, 55.5, 11.0, 22.0}, {80.0, 59.2, 12.5, 22.0}};
	double wt[9];
	double tol;	
	double *wtptr;
	int i, j, m = 4, n = 9;
	int ix = 2, iy = 2;
	int nx, ny;
	int ncv;
	int isz[4] = {-1, 1, 1, -1};
	int tdz= 4;
	int tde=6;
	int tdc=2;
	char weight = 'U';

	wtptr = NULL;
	tol = 1e-6;
	nx = ix;
	ny= iy;
	
	nag_mv_canon_corr(n, m, z, tdz, isz, nx, ny, wtptr, e, tde, &ncv, cvx, tdc, cvy, tdc, tol);
	
	printf("\n%s%2ld%s%2ld\n\n","Rank of x = ",nx," Rank of y= ",ny);
	printf(" Canonical Eigenvalues Percentage Chisq DF Sig\n");
	printf(" correlations          variation\n");
	for (i = 0; i < ncv; ++i)
	{
		for (j = 0; j < 6; ++j)
			printf("%12.4f",e[i][j]);
		printf("\n");
	}
	printf("\nCanonical coefficients for x\n");
	for (i = 0; i < ix; ++i)
	{
		for (j = 0; j < ncv; ++j)
			printf("%9.4f",cvx[i][j]);
		printf("\n");
	}
	printf("\nCanonical coefficients for y\n");
	for (i = 0; i < iy ;++i)
	{
		for (j = 0; j < ncv; ++j)
			printf("%9.4f",cvy[i][j]);
		printf("\n");
	}
}

	The output is following:
	
	Rank of x = 2 Rank of y= 2
	Canonical 		Eigenvalues 	Percentage 	Chisq 		DF 			Sig
	correlations 					variation
	
	0.9570 			10.8916 		0.9863 		14.3914 	4.0000 		0.0061
	0.3624 		 	 0.1512 		0.0137 		 0.7744 	1.0000 		0.3789
	
	Canonical coefficients for x
	-0.4261 		 1.0337
	-0.3444 		-1.1136
	
	Canonical coefficients for y
	-0.1415 		 0.1504
	-0.2384 		-0.3424	

Return:
	The function returns NAG error code or 0 if no error. 

	11: On entry, nx must not be less than 1: nx = _value_.  On entry, ny must not be less than 1: ny = _value_.  On entry, tde must not be less than 6: tde = _value_.
	5: On entry, tol must not be less than 0.0: tol = _value_.
	17: On entry, tdz = _value_ while m = _value_.  These parameters must satisfy tdz = m.
	29: On entry, m = _value_, nx = _value_ and ny = _value_.  These parameters must satisfy m = nx + ny.  On entry, n = _value_, nx = _value_ and ny = _value_.  These parameters must satisfy n > nx + ny.  On entry, tdcvx = _value_, nx = _value_ and ny = _value_.  These parameters must satisfy tdcvx = min(nx,ny).  On entry, tdcvy = _value_, nx = _value_ and ny = _value_.  These parameters must satisfy tdcvy = min(nx,ny).  
	500: On entry, wt[_value_] = _value_.  Constraint: When referenced, all elements of wt must be non-negative.
	501: The number of variables, nx in the analysis = _value_, while the number of x variables included in the analysis via array isz = _value_.  Constraint: these two numbers must be the same.  The number of variables, ny in the analysis = _value_, while the number of y variables included in the analysis via array isz = _value_.  Constraint: these two numbers must be the same.  
	503: With weighted data, the e.ective number of observations given by the sum of weights = _value_, while number of variables included in the analysis, nx + ny = _value_.  Constraint: E.ective number of observations = nx + ny + 1.
	427: The singular value decomposition has failed to converge.  This is an unlikely error exit.  
	506: A canonical correlation is equal to one.  This will happen if the x and y variables are perfectly correlated.
	509: The rank of the x matrix or the rank of the y matrix is zero.  This will happen if all the x and y variables are constants.
	73: Memory allocation failed.
	74: An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
		
	successfully call of the nag_mv_canon_corr function.		
*/
int nag_mv_canon_corr(
	int n, 		// the number of observations.
	int m, 		// the total number of variables.
	double z[], // contain the ith observation for the jth variable, for i = 1, 2, . . ., n; j = 1, 2, . . .,m.
	int tdz,	// the last dimension of the array z.
	int isz[], 	// indicates whether or not the jth variable is to be included in the analysis and to which set of variables it belongs.
	int nx, 	// the number of x variables in the analysis.
	int ny, 	// the number of y variables in the analysis.
	double wt[], // the weights to be used in the analysis.
	double e[],	// the statistics of the canonical variate analysis.
	int tde, 	// the last dimension of the array e.
	int *ncv, 	// the number of canonical correlations.
	double cvx[], // the canonical variate loadings for the x variables.
	int tdcvx,	// the last dimension of the array cvx.
	double cvy[], // the canonical variate loadings for the y variables.
	int tdcvy, 	// the last dimension of the array cvy.
	double tol 
); // Canonical coefficients, Eigenvalues, Percentage variation, Canonical coefficients...


/** g03bac
		computes orthogonal rotations for a matrix of loadings using a generalized orthomax criterion.
		
Example:
	The example is taken from page 75 of Lawley and Maxwell (1971). The results from a 
	factor analysis of ten variables using three factors are input and rotated using 
	varimax rotations without standardising rows.
	
	
void test_nag_mv_orthomax()
{

	double g = 1.0;
	matrix r, flr;
	r.SetSize(3,3);
	flr.SetSize(10,3);
	matrix fl = {{0.788, -0.152, -0.352}, {0.874, 0.381, 0.041}, {0.814, -0.043, -0.213},
	{0.798, -0.170, -0.204}, {0.641, 0.070, -0.042}, {0.755, -0.298, 0.067}, {0.782, -0.221, 0.028},
	{0.767, -0.091, 0.358}, {0.733, -0.384, 0.229}, {0.771, -0.101, 0.071}};
	double acc = 0.00001;	
	int iter, nvar;
	int i, j, k;
	int maxit = 20;
	int tdf = 3;
	int tdr = 3;
	char char_stand;
	Nag_RotationLoading stand;
	
	nvar = 10;
	k = 3;
	char_stand = 'U';
	
	if (char_stand == 'S')
		stand = Nag_RoLoadStand;
	else
		stand = Nag_RoLoadNotStand;

	nag_mv_orthomax(stand, g, nvar, k, fl, tdf, flr, r, tdr, acc, maxit, &iter);
	
	printf("\n Rotated factor loadings\n\n");
	for (i = 0; i < nvar; ++i)
	{
		for (j = 0; j < k ;++j)
			printf(" %8.3f",flr[i][j]);
		printf("\n");
	}
	printf("\n Rotation matrix\n\n");
	for (i = 0; i < k ;++i)
	{
		for (j = 0; j < k ;++j)
			printf(" %8.3f",r[i][j]);
		printf("\n");
	}	
}	
		
	The output is following:
		
	Rotated factor loadings
	
	0.329 	-0.289 	-0.759
	0.849 	-0.273 	-0.340
	0.450 	-0.327 	-0.633
	0.345 	-0.397 	-0.657
	0.453 	-0.276 	-0.370
	0.263 	-0.615 	-0.464
	0.332 	-0.561 	-0.485
	0.472 	-0.684 	-0.183
	0.209 	-0.754 	-0.354
	0.423 	-0.514 	-0.409
	
	Rotation matrix
	
	0.633 	-0.534 	-0.560
	0.758 	 0.573 	 0.311
	0.155 	-0.622 	 0.768

Return:
	The function returns NAG error code or 0 if no error. 
	
	70: On entry, parameter stand had an illegal value.
	11: On entry, k must not be less than 2: k = _value_.
	12: On entry, maxit must not be less than or equal to 0 : maxit = _value_.
	5:  On entry, g must not be less than 0.0: g = _value_.  On entry, acc must not be less than 0.0: acc = _value_.
	17: On entry, nvar = _value_ while k = _value_.  These parameters must satisfy nvar = k.  On entry, tdf = _value_ while k = _value_.  These parameters must satisfy tdf = k.  On entry, tdr = _value_ while k = _value_.  These parameters must satisfy tdr = k.
	427: The singular value decomposition has failed to converge.  This is an unlikely error exit.
	511: The algorithm to .nd R has failed to reach the required accuracy in the given number of iterations, _value_. Try increasing acc or increasing maxit. The returned solution should be a reasonable approximation.
	73: Memory allocation failed.  
	74: An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.  
		
	successfully call of the nag_mv_orthomax function.	
*/
int nag_mv_orthomax(
    Nag_RotationLoading stand,
    double g, // the criterion constant.
    int nvar, // The number of original variables.
    int k,	// The number of derived variates or factors.	  
    double fl[], // the matrix of loadings.
    int tdf, 	// the last dimension of the arrays fl and flr.
    double flr[], // the rotated matrix of loadings.
    double r[], // the matrix of rotations.
    int tdr,	// the last dimension of the array r.
    double acc, // indicates the accuracy required.
    int maxit,	// the maximum number of iterations.
    int *iter	// the number of iterations performed.
); // Rotation factor loading, Rotation matrix.

/** g03bcc
		computes Procrustes rotations in which an orthogonal rotation is found 
		so that a transformed matrix best matches a target matrix.
	
Example:
	Three points representing the vertices of a triangle in two dimensions are input. 
	The points are translated and rotated to match the triangle given by (0,0),(1,0),(0,2) 
	and scaling is applied after rotation. 

void test_nag_mv_procustes()
{
	matrix r, yhat;
	r.SetSize(2,2);
	yhat.SetSize(3,2);
	matrix x = {{0.63, 0.58}, {1.36, 0.39}, {1.01, 1.76}};
	matrix y = {{0.0, 0.0}, {1.0, 0.0}, {0.0, 2.0}};
	double res[3];
	double alpha;
	double rss;
	int i, j, m = 2, n = 3;
	int tdx = 2;
	int tdr = 2;
	int tdy = 2;
	char char_scale;
	char char_stand;
	Nag_TransNorm stand;
	Nag_RotationScale scale;

	char_stand = 'C';
	char_scale = 'S';

	if (char_stand == 'N')
		stand = Nag_NoTransNorm;
	
	else if (char_stand == 'Z')
		stand = Nag_Orig;
	
	else if (char_stand == 'C')
		stand = Nag_OrigCentroid;
	
	else if (char_stand == 'U')
		stand = Nag_Norm;
	
	else if (char_stand == 'S')
		stand = Nag_OrigNorm;
	
	else if (char_stand == 'M')
		stand = Nag_OrigNormCentroid;
	
	if (char_scale == 'S')
		scale = Nag_LsqScale;
	
	else if (char_scale == 'U')
		scale = Nag_NotLsqScale;
	
	nag_mv_procustes(stand, scale, n, m, x, tdx, y, tdy, yhat, r, tdr, &alpha, &rss, res);
	
	printf("\n Rotation matrix\n\n");
	for (i = 0; i < m; ++i)
	{
		for (j = 0; j < m; ++j)
			printf(" %7.3f ",r[i][j]);
		printf("\n");
	}
	if (char_scale == 'S' || char_scale == 's')
		printf("\n%s%10.3f\n"," Scale factor = ",alpha);
	
	printf("\n Target matrix \n\n");
	for (i = 0; i < n; ++i)
	{
		for (j = 0; j < m; ++j)
			printf(" %7.3f ",y[i][j]);
		printf("\n");
	}
	printf("\n Fitted matrix\n\n");
	for (i = 0; i < n; ++i)
	{
		for (j = 0; j < m; ++j)
			printf(" %7.3f ",yhat[i][j]);
		printf("\n");
	}
	printf("\n%s%10.3f\n","RSS = ",rss);	
}	
	
	
	The output is following:
	
	Rotation matrix
	 0.967 		0.254
	-0.254 		0.967
	
	Scale factor = 1.556
	
	Target matrix
	0.000 	0.000
	1.000 	0.000
	0.000 	2.000
	
	Fitted matrix
	-0.093 		0.024
	 1.080 		0.026
	 0.013 		1.950
	 
	RSS = 0.019	

Return:
	The function returns NAG error code or 0 if no error. 
	
	70: On entry, parameter stand had an illegal value.  On entry, parameter pscale had an illegal value.
	11: On entry, n must not be less than 1: n = _value_.  On entry, m must not be less than 1: m = _value_.
	17: On entry, tdx = _value_ while m = _value_.  These parameters must satisfy tdx = m.  On entry, tdy = _value_ while m = _value_.  These parameters must satisfy tdy = m.  On entry, tdr = _value_ while m = _value_.  These parameters must satisfy tdr = m.
	19: On entry, m = _value_ while n = _value_.  These parameters must satisfy m = n.
	427: The singular value decomposition has failed to converge.  This is an unlikely error exit.
	512: On entry, either x or y contains only zero-points (possibly after translation) when normalization is to be applied.
	513: The fitted matrix Y-hat, contains only zero-points when least-squares scaling is applied.
	73: Memory allocation failed.
	74: An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
		
	successfully call of the nag_mv_orthomax function.
*/
int nag_mv_procustes(
    Nag_TransNorm stand,
    Nag_RotationScale pscale,
    int n,		// the number of points. 
    int m,		// the number of dimensions.
    double x[], // the matrix to be rotated.
    int tdx,	// the last dimension of the array x.
    double y[], // the target matrix.
    int tdy,	// the last dimension of the arrays y.
    double yhat[], // the fitted matrix
    double r[],	// the matrix of rotations.
    int tdr,	// the last dimension of the array r.
    double *alpha, // if pscale = Nag LsqScale the scaling factor, a; otherwise alpha is not set
    double *rss, // the residual sum of squares.
    double res[] // the residuals.
); // Rotation matrix, Target matrix, Fitted matrix, RSS, Scale factor.
 
/** g03cac
		computes the maximum likelihood estimates of the parameters of a factor analysis 
		model. Either the data matrix or a correlation/covariance matrix may be input. 
		Factor loadings, communalities and residual correlations are returned.
	
Example:
	The example is taken from Lawley and Maxwell (1971). The correlation matrix for nine variables
	is input and the parameters of a factor analysis model with three factors are estimated and printed.

	// Note: This function needs the support of struct. It has not been done yet!

void test_nag_mv_factor()
{
	
	matrix fl, x;
	fl.SetSize(9,9);
	x.SetSize(9,9);
	double com[9], e[9];	
	double psi[9], res[36], stat[4], wt[9];	
	double *wtptr=NULL;
	double eps;
	int nfac, nvar;
	int i, j, l, m, n;
	int isx[9];
	int tdfl = 9, tdx = 9;
	char weight, char_mat;
	Nag_FacMat mat;
	Nag_E04_Opt options;
	
	Vprintf("g03cac Example Program Results\n\n");
	Vscanf("%*[^\n]");
	Vscanf("%s",char_mat);
	Vscanf("%s",weight);
	Vscanf("%ld",&n);
	Vscanf("%ld",&m);
	Vscanf("%ld",&nvar);
	Vscanf("%ld",&nfac);
	
	if (m <= MMAX && (*char_mat == 'C' || n <= NMAX ))
	{
		if (*char_mat == 'C')
		{
			for (i = 0; i < m; ++i)
			{
				for (j = 0; j < m; ++j)
					Vscanf("%lf",&x[i][j]);
			}
		}
		else
		{
			if (*weight == 'W' || *weight == 'w')
			{
				for (i = 0; i < n; ++i)
				{
					for (j = 0; j < m; ++j)
						Vscanf("%lf",&x[i][j]);
						
					Vscanf("%lf",&wt[i]);
				}
				wtptr = wt;
			}
			else
			{
				for (i = 0; i < n; ++i)
				{
					for (j = 0; j < m; ++j)
						Vscanf("%lf",&x[i][j]);
				}
			}
		}
		for (j = 0; j < m; ++j)
			Vscanf("%ld",&isx[j]);

		3.g03cac.6 [NP3275/5/pdf]
		g03 - Multivariate Methods g03cac
		if (*char_mat == 'D')
		{
			mat = Nag_DataCorr;
		}
		else if (*char_mat == 'S')
		{
			mat = Nag_DataCovar;
		}
		else if (*char_mat == 'C')
		{
			mat = Nag_MatCorr_Covar;
		}
		e04xxc(&options);
		options.max_iter = 500;
		options.optim_tol = 1e-2;
		eps = 1e-5;
		g03cac(mat, n, m, (double *)x, tdx, nvar, isx, nfac, wtptr, e,
		stat, com, psi, res, (double *)fl, tdfl, &options, eps, NAGERR_DEFAULT);
		Vprintf("\nEigenvalues\n\n");
		for (j = 0; j < m; ++j)
		{
			Vprintf(" %12.4e%s",e[j], (j+1)%6==0 ? "\n" : "");
		}
		Vprintf("\n\n%s%6.3f\n"," Test Statistic = ",stat[1]);
		Vprintf("%s%6.3f\n"," df = ",stat[2]);
		Vprintf("%s%6.3f\n\n","Significance level = ",stat[3]);
		Vprintf("Residuals\n\n");
		l = 1;
		for (i = 1; i <= nvar-1; ++i)
		{
			for (j = l; j <= l+i-1; ++j)
				Vprintf(" %8.3f",res[j-1]);
			Vprintf("\n");
			l += i;
		}
		Vprintf("\nLoadings, Communalities and PSI\n\n");
		for (i = 0; i < nvar; ++i)
		{
			for (j = 0; j < nfac; ++j)
				Vprintf(" %8.3f",fl[i][j]);
				Vprintf("%8.3f%8.3f\n", com[i], psi[i]);
		}
		exit(EXIT_SUCCESS);
	}
	else
	{
		Vprintf("Incorrect input value of n or m.\n");
		exit(EXIT_FAILURE);
	}	
}	
	
	The output is following:
		
	Parameters to e04lbc
	--------------------
	Number of variables........... 9
	optim_tol............... 1.00e-02 linesearch_tol.......... 9.00e-01
	step_max................ 2.70e+01 max_iter................ 500
	print_level.........Nag_Soln_Iter machine precision....... 1.11e-16
	deriv_check............. FALSE
	outfile................. stdout
	Memory allocation:
	state................... User
	hesl.................... User hesd................... User
	
	Iterations performed = 0, function evaluations = 1
	Criterion = 8.635756e-02
	
	Variable 		Standardized
					Communalities
	1 				0.5755
	2 				0.5863
	3 				0.4344
	4 				0.7496
	5 				0.6203
	6 				0.7329
	7 				0.6061
	8 				0.4053
	9 				0.7104
	
	Iterations performed = 1, function evaluations = 3
	Criterion = 3.603203e-02
	
	Variable 		Standardized
					Communalities
	1 				0.5517
	2 				0.5800
	3 				0.3936
	4 				0.7926
	5 				0.6140
	6 				0.8254
	7 				0.6052
	8 				0.5076
	9 				0.7569
	
	Iterations performed = 2, function evaluations = 4
	Criterion = 3.502097e-02
	
	Variable 		Standardized
					Communalities
	1 				0.5496
	2 				0.5731
	3 				0.3838
	4 				0.7875
	5 				0.6200
	6 				0.8238
	7 				0.6006
	8 				0.5349
	9 				0.7697
	
	Iterations performed = 3, function evaluations = 5
	Criterion = 3.501729e-02
	
	Variable 		Standardized
					Communalities
	1 				0.5495
	2 				0.5729
	3 				0.3835
	4 				0.7877
	5 				0.6195
	6 				0.8231
	7 				0.6005
	8 				0.5384
	9 				0.7691
	
	Eigenvalues
	1.5968e+01 4.3577e+00 1.8474e+00 1.1560e+00 1.1190e+00 1.0271e+00
	9.2574e-01 8.9508e-01 8.7710e-01
	
	Test Statistic = 7.149
	df = 12.000
	Significance level = 0.848
	
	Residuals
	 0.000
	-0.013  0.022
	 0.011 -0.005  0.023
	-0.010 -0.019 -0.016  0.003
	-0.005  0.011 -0.012 -0.001 -0.001
	 0.015 -0.022 -0.011  0.002  0.029 -0.012
	-0.001 -0.011  0.013  0.005 -0.006 -0.001  0.003
	-0.006  0.010 -0.005 -0.011  0.002  0.007  0.003 -0.001
	
	Loadings, Communalities and PSI
	0.664 -0.321  0.074 0.550 0.450
	0.689 -0.247 -0.193 0.573 0.427
	0.493 -0.302 -0.222 0.383 0.617
	0.837  0.292 -0.035 0.788 0.212
	0.705  0.315 -0.153 0.619 0.381
	0.819  0.377  0.105 0.823 0.177
	0.661 -0.396 -0.078 0.600 0.400
	0.458 -0.296  0.491 0.538 0.462
	0.766 -0.427 -0.012 0.769 0.231	
	
Return:
	The function returns NAG error code or 0 if no error. 

	70: On entry, parameter matrix had an illegal value.
	11: On entry, nfac must not be less than 1: nfac = _value_.  On entry, nvar must not be less than 2: nvar = _value_.
	17: On entry, m = _value_ while nvar = _value_.  These parameters must satisfy m = nvar.  On entry, tdx = _value_ while m = _value_.  These parameters must satisfy tdx = m.  On entry, td. = _value_ while nfac = _value_.  These parameters must satisfy td. = nfac.
	18: On entry, n = _value_ while nvar = _value_.  These parameters must satisfy n > nvar.
	19: On entry, nfac = _value_ while nvar = _value_.  These parameters must satisfy nfac = nvar.
	44: Value _value_ given to eps is not valid.  Correct range is machine precision = eps < 1.0.
	500: On entry, wt[_value_] = _value_.  Constraint: When referenced, all elements of wt must be non-negative.
	501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_.  Constraint: these two numbers must be the same.
	503: With weighted data, the e.ective number of observations given by the sum of weights = _value_, while the number of variables included in the analysis, nvar = _value_.  Constraint: e.ective number of observations > nvar + 1.
	427: A singular value decomposition has failed to converge.  This is a very unlikely error exit.
	284: The conditions for a minimum have not all been satis.ed but a lower point could not be found.  Note that in this case all the results are computed.  See nag opt bounds 2nd deriv (e04lbc) for further details.
	285:The maximum number of iterations, _value_, have been performed.
	514: On entry, matrix = Nag DataCorr or matrix = Nag DataCovar and the data matrix is not of full column rank, or matrix = Nag MatCorr Covar and the input correlation/variancecovariance matrix is not positive-de.nite.  This exit may also be caused by two of the eigenvalues of S* being equal; this is rare (see Lawley and Maxwell (1971)) and may be due to the data/correlation matrix being almost singular.  
	73: Memory allocation failed.
	74: An internal error has occurred in this function.  Check the function call and array sizes. If the function call is correct, please consult NAG for assistance.  Additional error messages will output if the optimisation fails to converge or if the options are set incorrectly. Details of these can be found in the nag opt bounds 2nd deriv (e04lbc) document.  	
	
	successfully call of the nag_mv_factor function.
*/  

/**
  int nag_mv_factor(
    Nag_FacMat mat,
    int n,		// the number of observations .
    int m,		// the number of variables in the data/correlation/variance-covariance matrix.
    double x[], // the input matrix.
    int tdx,	// the last dimension of the array x.
    int nvar,	// the number of variables in the factor analysis.
    int isx[],	// indicates whether or not the jth variable is to be included in the factor analysis.
    int nfac,	// the number of factors.
    double *wtptr, // the weights to be used in the factor analysis.
    double e[],	// the eigenvalues.
    double stat[], // the test statistics.
    double com[], // the communalities.
    double psi[], // the estimates of psi(i), for i = 1, 2, . . . , p.
    double res[], // the residual correlations.
    double fl[], // the factor loadings.
    int tdfl,	// the last dimension of the array fl.
    Nag_E04_Opt *options,
    double eps, // A lower bound for the value of Øi.
  ); // Loadings, Communalities and PSI, Eigenvalues, Standardized Communalities, ....
*/


/** g03ccc
		computes factor score coefficients from the result of fitting a factor analysis 
		model by maximum likelihood as performed by nag mv factor (g03cac).
	
Example:

	The example is taken from Lawley and Maxwell (1971). The correlation matrix 
	for 220 observations on six school subjects is input and a factor analysis 
	model with two factors .tted using nag mv factor (g03cac). The factor score 
	coefficients are computed using the regression method.
	
// Note: This one needs the support of the the results of the g03cac

Return:
	The function returns NAG error code or 0 if no error. 
	
	70:  On entry, parameter method had an illegal value.  On entry, parameter rotate had an illegal value.
	11:  On entry, nfac must not be less than 1: nfac = _value_.
	17:  On entry, nvar = _value_ while nfac = _value_.  These parameters must satisfy nvar = nfac.  On entry, td. = _value_ while nfac = _value_.  These parameters must satisfy td. = nfac.  On entry, tdfs = _value_ while nfac = _value_.  These parameters must satisfy tdfs = nfac.
	109: On entry, tdr = _value_ while nfac = _value_ and rotate = Nag FacRotate.  These parameters must satisfy tdr = nfac when rotate = Nag FacRotate.
	55:  On entry, e[_value_] = _value_.  Constraint: e[_value_] > 1.0.  On entry, psi[_value_] = _value_.  Constraint: psi[_value_] > 0.0.
	73:  Memory allocation failed.
	74:  An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.		
*/

/**  
  int nag_mv_fac_score(
    Nag_FacScoreMethod method,
    Nag_FacRotation rotate,
    int nvar, // the number of observed variables in the factor analysis.
    int nfac, // the number of factors in the factor analysis.
    double fl[], // the matrix of unrotated factor loadings, Ë, as returned by nag_mv_factor (g03cac).
    int tdfl,	// the last dimension of the array fl.
    double psi[], // the diagonal elements of Ø, as returned by nag mv factor (g03cac).
    double e[],	// the eigenvalues of the matrix S*, as returned by nag mv factor (g03cac).
    double r[], // the orthogonal rotation matrix, R, as returned by nag mv orthomax (g03bac) if rotate = Nag FacRotate.
    int tdr, // the last dimension of the array r.
    double fs[], // the matrix of factor score coefficients
    int tdfs // the last dimension of the array fs.
  ); // Loadings, Communalities and PSI, Factor score coefficients, Standardized Communalities ...
*/

/** g03dac 
		computes a test statistic for the equality of within-group covariance matrices 
		and also computes matrices for use in discriminant analysis.
		
Example:
	The data, taken from Aitchison and Dunsmore (1975), is concerned with the diagnosis 
	of three 'types' of Cushing's syndrome. The variables are the logarithms of the 
	urinary excretion rates (mg/24hr) of two steroid metabolites. Observations for a 
	total of 21 patients are input and the statistics computed by nag mv discrim (g03dac). 
	The printed results show that there is evidence that the within-group variance-covariance 
	matrices are not equal.

void test_nag_mv_discrim()
{
	double stat;
	double det[3], gc[12];
	matrix gmean;
	gmean.SetSize(3,2);
	matrix x = {{1.1314, 2.4596}, {1.0986, 0.2624}, {0.6419, -2.3026}, {1.3350, -3.2189},
	{1.4110, 0.0953}, {0.6419, -0.9163}, {2.1163, 0.0000}, {1.3350, -1.6094}, {1.3610, -0.5108}, 
	{2.0541, 0.1823}, {2.2083, -0.5108}, {2.7344, 1.2809}, {2.0412, 0.4700}, {1.8718, -0.9163},
	{1.7405, -0.9163}, {2.6101, 0.4700}, {2.3224, 1.8563}, {2.2192, 2.0669}, {2.2618, 1.1314}, 
	{3.9853, 0.9163}, {2.7600, 2.0281}};	
	double wt[21];
	double df;
	double *wtptr = NULL;
	double sig;
	int nvar = 2;
	int i, j, m = 2, n = 21;
	int isx[2], nig[3];
	int ing[21] = {1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3};
	int ng = 3;
	int tdgmean = 2;
	int tdx = 2;
	char weight = 'U';	

	for (j = 0; j < m; ++j)
		isx[j] = 1;
		
	nag_mv_discrim(n, m, x, tdx, isx, nvar, ing, ng, wtptr, nig, gmean, tdgmean, det, gc, 
				&stat, &df, &sig);
				
	printf("\nGroup means\n\n");
	for (i = 0; i < ng; ++i)
	{
		for (j = 0; j < nvar; ++j)
			printf("%10.4f",gmean[i][j]);
		printf("\n");
	}
	printf("\nLOG of determinants\n\n");
	for (j = 0; j < ng; ++j)
		printf("%10.4f",det[j]);
	printf("\n\n%s%7.4f\n","stat = ",stat);
	printf("%s%7.4f\n"," df = ",df);
	printf("%s%7.4f\n"," sig = ",sig);
}

	The output is following:
	
	Group means

    1.0433   -0.6034
    2.0073   -0.2060
    2.7097    1.5998

	LOG of determinants

   		-0.8273   -3.0460   -2.2877

	stat = 19.2410
 	df =  6.0000
 	sig =  0.0038

Return:
	The function returns NAG error code or 0 if no error. 
	
	11:  On entry, n must not be less than 1: n = _value_.  On entry, nvar must not be less than 1: nvar = _value_.  On entry, ng must not be less than 2: ng = _value_.
	17:  On entry, m = _value_ while nvar = _value_.  These parameters must satisfy m = nvar.  On entry, tdx = _value_ while m = _value_.  These parameters must satisfy tdx = m.  On entry, tdg = _value_ while nvar = _value_.  These parameters must satisfy tdg = nvar.  
	112: On entry, ing[_value_] = _value_, ng = _value_.  Constraint: 1 = ing[i - 1] = ng, i = 1, 2, . . ., n.
	500: On entry, wt[_value_] = _value_. Constraint: when referenced, all elements of wt must be non-negative.
	501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_.  Constraint: these two numbers must be the same.
	515: On entry, group _value_ has _value_ e.ective observations.  Constraint: in each group the e.ective number of observations must be = 1.
	516: On entry, group _value_ has _value_ members, while nvar = _value_.  Constraint: number of members in each group = nvar
	520: The variables in group _value_ are not of full rank.
	519: The variables are not of full rank.
	73:  Memory allocation failed.
	74:  An internal error has occurred in this function.  Check the function call and array sizes. If the function call is correct, please consult NAG for assistance.
	
	successfully call of the nag_mv_discrim function.	
*/
int nag_mv_discrim(
    int n, // the number of observations.
    int m, // the number of variables in the data array x.
    double x[], // contain the kth observation for the lth variable, for k = 1, 2, . . ., n; l = 1, 2, . . ., m .
    int tdx,	// the last dimension of the array x.
    int isx[], 	// indicates whether or not the lth variable in x is to be included in the variance-covariance matrices.
    int nvar, 	// the number of variables in the variance-covariance matrices.
    int ing[],	// indicates to which group the kth observation belongs, for k = 1, 2, . . ., n.
    int ng,		// the number of groups.
    double *wtptr, // the weights to be used in the analysis and the effective number of observations for a group is the sum of the weights of the observations in that group.
    int nig[], 	// the number of observations in the jth group, for j = 1, 2, . . . , ng
    double gmean[], // the jth row of gmean contains the means of the p selected variables for the jth group, for j = 1, 2, . . . , ng.
    int tdg,	// the last dimension of the array gmean.
    double det[], // the logarithm of the determinants of the within-group variance-covariance matrices.
    double gc[], // the first p(p+1)/2 elements of gc contain R and the remaining ng blocks of p(p+1)/2 elements contain the Rj matrices.
    double *stat, // the likelihood-ratio test statistic, G.
    double *df,	// the degrees of freedom for the distribution of G.
    double *sig // the significance level for G.
);
  

/** g03dbc
		computes Mahalanobis squared distances for groups or pooled variance-covariance 
		matrices. It is intended for use after nag mv discrim (g03dac).
		
Example:
	The data, taken from Aitchison and Dunsmore (1975), are concerned with 
	the diagnosis of three 'types' of Cushing's syndrome. The variables are 
	the logarithms of the urinary excretion rates (mg/24hr) of two steroid 
	metabolites. Observations for a total of 21 patients are input and the
	group means and R matrices are computed by nag mv discrim (g03dac). A 
	further six observations of unknown type are input, and the distances 
	from the group means of the 21 patients of known type are computed under 
	the assumption that the within-group variance-covariance matrices are 
	not equal. These results are printed and indicate that the .rst four 
	are close to one of the groups while observations 5 and 6 are some 
	distance from any group.
	
void test_nag_mv_discrim_mahaldist()
{	
	matrix d, gmean;
	d.SetSize(21,3);
	gmean.SetSize(3,2);
	matrix x = {{1.1314, 2.4596}, {1.0986, 0.2624}, {0.6419, -2.3026}, {1.3350, -3.2189},
	{1.4110, 0.0953}, {0.6419, -0.9163}, {2.1163, 0.0000}, {1.3350, -1.6094}, {1.3610, -0.5108}, 
	{2.0541, 0.1823}, {2.2083, -0.5108}, {2.7344, 1.2809}, {2.0412, 0.4700}, {1.8718, -0.9163},
	{1.7405, -0.9163}, {2.6101, 0.4700}, {2.3224, 1.8563}, {2.2192, 2.0669}, {2.2618, 1.1314}, 
	{3.9853, 0.9163}, {2.7600, 2.0281}}; 
	double det[3], gc[12];	
	double	wt[21];
	double stat;
	double df;
	double sig;
	double *wtptr;
	int nobs, nvar;
	int isx[2], nig[3];
	int ing[21] = {1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3};
	int i, j, m = 2, n = 21;
	int ng = 3;
	int tdd = 3;
	int tdgmean = 2;
	int tdx = 2;
	char char_equal;
	char weight = 'U';
	Nag_GroupCovars equal;

	wtptr = NULL;
	nvar = 2;
	
	for (j = 0; j < m; ++j)
		isx[j] = 1;
	
	nag_mv_discrim(n, m, x, tdx, isx, nvar, ing, ng, wtptr, nig, gmean, tdgmean, det, gc, 
				&stat, &df, &sig);
	
	nobs = 6;
	char_equal = 'U';
	
	x[0][0] = 1.6292;
	x[0][1] = -0.9163;
	
	x[1][0] = 2.5572;
	x[1][1] = 1.6094;
	
	x[2][0] = 2.5649;
	x[2][1] = -0.2231;
	
	x[3][0] = 0.9555;
	x[3][1] = -2.3026;
	
	x[4][0] = 3.4012;
	x[4][1] = -2.3026;
	
	x[5][0] = 3.0204;
	x[5][1] = -0.2231;
	
	if(char_equal == 'E')
		equal = Nag_EqualCovar;
	else if (char_equal == 'U')
		equal = Nag_NotEqualCovar;
	
	nag_mv_discrim_mahaldist(equal, Nag_SamplePoints, nvar, ng, gmean, tdgmean, gc,
							nobs, m, isx, x, tdx, d, tdd);

	printf("\n Obs Distances\n\n");
	for (i = 0; i < nobs; ++i)
	{
		printf(" %3ld",i+1);
		for (j = 0; j < ng; ++j)
			printf("%10.3f",d[i][j]);
		printf("\n");
	}
}

	The Output is following:
	
	Obs Distances

   1     3.339     0.752    50.928
   2    20.777     5.656     0.060
   3    21.363     4.841    19.498
   4     0.718     6.280   124.732
   5    55.000    88.860    71.785
   6    36.170    15.785    15.749

Return:
  	The function returns NAG error code or 0 if no error. 
  	
  	70: On entry, parameter equal had an illegal value.  On entry, parameter mode had an illegal value.
  	11: On entry, nvar must not be less than 1: nvar = _value_.  On entry, ng must not be less than 2: ng = _value_.
  	17: On entry, tdg = _value_ while nvar = _value_.  These parameters must satisfy tdg = nvar.  On entry, tdd = _value_ while ng = _value_.  These parameters must satisfy tdd = ng.
  	108: On entry, nobs = _value_ while mode = Nag SamplePoints.  These parameters must satisfy nobs = 1 w hen mode = Nag SamplePoints.
  	109: On entry, m = _value_ while nvar = _value_ and mode = Nag SamplePoints.  These parameters must satisfy m = nvar when mode = Nag SamplePoints.  On entry, tdx = _value_ while m = _value_ and mode = Nag SamplePoints.  These parameters must satisfy tdx = max(1,m) w hen mode = Nag SamplePoints.
  	502: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_.  Constraint: These two numbers must be the same when mode = Nag SamplePoints.
  	517: A diagonal element of R is zero when equal = Nag EqualCovar.
  	518: A diagonal element of R is zero for some j, w hen equal = Nag NotEqualCovar.
	73: Memory allocation failed.
	74: An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.

  	successfully call of the nag_mv_discrim_mahaldist function.	 	
*/
int nag_mv_discrim_mahaldist(
    Nag_GroupCovars equal,
    Nag_MahalDist mode,
    int nvar, // the number of variables, p, in the variance-covariance matrices as specified to nag mv discrim (g03dac).
    int ng,   // the number of groups.
    double gmean[], // the jth row of gmean contains the means of the p selected variables for the jth group, for j = 1, 2, . . . , ng. 
    int tdg, // the last dimension of the array gmean.
    double gc[], // the first p(p + 1)/2 elements of gc should contain the upper triangular matrix R and the next ng blocks of p(p+1)/2 elements should contain the upper triangular matrices Rj.
    int nobs, 	// if mode = Nag SamplePoints the number of sample points in x for which distances are to be calculated.
    int m,		// if mode = Nag SamplePoints the number of variables in the data array x. 
    int isx[], 	// if mode = Nag SamplePoints, isx[l-1] indicates if the lth variable in x is to be included in the distance calculations.
    double x[], // contain the kth sample value for the lth variable for k = 1, 2, . . . ,nobs; l = 1, 2, . . . ,m. 
    int tdx,	// the last dimension of the array x.
    double d[], // the squared distances.
    int tdd		// the last dimension of the array dd.
);

/** g03dcc
		allocates observations to groups according to selected rules. It is intended 
		for use after nag mv discrim (g03dac).
		
Example
	The data, taken from Aitchison and Dunsmore (1975), are concerned with the 
	diagnosis of three 'types' of Cushing's syndrome. The variables are the 
	logarithms of the urinary excretion rates (mg/24hr) of two steroid 
	metabolites. Observations for a total of 21 patients are input and the 
	group means and R matrices are computed by nag mv discrim (g03dac). A 
	further six observations of unknown type are input and allocations made 
	using the predictive approach and under the assumption that the within-group 
	covariance matrices are not equal. The posterior probabilities of group 
	membership, q[j], and the atypicality index are printed along with the 
	allocated group. The atypicality index shows that observations 5 and 6 do not
	seem to be typical of the three types present in the initial 21 observations.

void test_nag_mv_discrim_group()
{
	double stat;
	matrix ati, gmean, p;
	ati.SetSize(21,3);
	gmean.SetSize(3,2);
	p.SetSize(21,3);
	matrix x = {{1.1314, 2.4596}, {1.0986, 0.2624}, {0.6419, -2.3026}, {1.3350, -3.2189},
	{1.4110, 0.0953}, {0.6419, -0.9163}, {2.1163, 0.0000}, {1.3350, -1.6094}, {1.3610, -0.5108}, 
	{2.0541, 0.1823}, {2.2083, -0.5108}, {2.7344, 1.2809}, {2.0412, 0.4700}, {1.8718, -0.9163},
	{1.7405, -0.9163}, {2.6101, 0.4700}, {2.3224, 1.8563}, {2.2192, 2.0669}, {2.2618, 1.1314}, 
	{3.9853, 0.9163}, {2.7600, 2.0281}}; 
	double det[3], gc[12];	
	double prior[3], wt[21];
	double df;
	double sig;
	double *wtptr = NULL;
	int nobs, nvar;
	int i, j, m, n;
	int iag[21], isx[2], nig[3];
	int ing[21] = {1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3};
	int ng = 3;
	int tdgmean = 2;
	int tdp = 3;
	int tdx = 2;
	Boolean atiq = true;
	char char_type;
	char char_equal;
	char weight = 'U';
	Nag_DiscrimMethod type;
	Nag_GroupCovars equal;
	
	n = 21;
	m = 2;
	nvar = 2;

	for (j = 0; j < m; ++j)
		isx[j] = 1;
	
	nag_mv_discrim(n, m, x, tdx, isx, nvar, ing, ng, wtptr, nig, gmean, tdgmean, det, gc, 
				&stat, &df, &sig);
				
	nobs = 6;
	char_equal ='U';
	char_type = 'P';
	x[0][0] = 1.6292;
	x[0][1] = -0.9163;
	
	x[1][0] = 2.5572;
	x[1][1] = 1.6094;
	
	x[2][0] = 2.5649;
	x[2][1] = -0.2231;
	
	x[3][0] = 0.9555;
	x[3][1] = -2.3026;
	
	x[4][0] = 3.4012;
	x[4][1] = -2.3026;
	
	x[5][0] = 3.0204;
	x[5][1] = -0.2231;

	if (char_type == 'E')
		type = Nag_DiscrimEstimate;
	if(char_type == 'P')
		type = Nag_DiscrimPredict;
	
	if(char_equal == 'E')
		equal = Nag_EqualCovar;
	if (char_equal == 'U')
		equal = Nag_NotEqualCovar;
	
	nag_mv_discrim_group(type, equal, Nag_EqualPrior, nvar, ng, nig, gmean,
					tdgmean, gc, det, nobs, m, isx, x, tdx, prior, p, tdp, iag,
					 atiq, ati);
	printf("\n");
	printf(" Obs Posterior Allocated ");
	printf(" Atypicality ");
	printf("\n");
	printf(" probabilities to group index ");
	printf("\n");
	printf("\n");
	for (i = 0; i < nobs; ++i)
	{
		printf(" %6ld ",i+1);
		for (j = 0; j < ng; ++j)
			printf("%6.3f",p[i][j]);
		
		printf(" %6ld ",iag[i]);
		for (j = 0; j < ng; ++j)
			printf("%6.3f",ati[i][j]);
		
		printf("\n");
	}
}

	The output is following:
	
	Obs Posterior Allocated  Atypicality 
 	probabilities to group index 

      1  0.094 0.905 0.002      2  0.596 0.254 0.975
      2  0.005 0.168 0.827      3  0.952 0.836 0.018
      3  0.019 0.920 0.062      2  0.954 0.797 0.912
      4  0.697 0.303 0.000      1  0.207 0.860 0.993
      5  0.317 0.013 0.670      3  0.991 1.000 0.984
      6  0.032 0.366 0.601      3  0.981 0.978 0.887

Return:
	The function returns NAG error code or 0 if no error. 
	
	70:  On entry, parameter type had an illegal value.  On entry, parameter equal had an illegal value.  On entry, parameter priors had an illegal value.  
	11:  On entry, nvar must not be less than 1: nvar = _value_.  On entry, ng must not be less than 2: ng = _value_.  On entry, nobs must not be less than 1: nobs = _value_.
	17:  On entry, m = _value_ while nvar = _value_.  These parameters must satisfy m = nvar.  On entry, tdx = _value_ while m = _value_.  These parameters must satisfy tdx = m.  On entry, tdp = _value_ while ng = _value_.  These parameters must satisfy tdp = ng.  On entry, tdg = _value_ while nvar = _value_.  These parameters must satisfy tdg = nvar.
	501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_.  Constraint: these two numbers must be the same.
	111: On entry, nig[_value_] = _value_.  Constraint: nig[i - 1] > 0, i = 1, 2, . . . ,ng when equal = Nag EqualCovar.
	112: On entry, nig[_value_] = _value_, nvar = _value_.  Constraint: nig[i - 1] > nvar, i = 1, 2, . . . ,ng when equal = Nag NotEqualCovar.
	115: On entry, prior[_value_] = _value_.  Constraint: prior[j - 1] > 0, j=1,2,. . . ,ng when priors = Nag UserPrior.
	522: On entry,summation of prior[j - 1] = _value_.  Constraint: summation of prior[j - 1] must be within 10 × machine precision of 1 when priors = Nag UserPrior.
	521: On entry, summation of nig[j - 1] = _value_, ng = _value_, nvar = _value_.  Constraint: summation of nig[j - 1] > ng + nvar when equal = Nag EqualCovar.
	517: A diagonal element of R is zero when equal = Nag EqualCovar.
	518: A diagonal element of R is zero for some j, when equal = Nag NotEqualCovar
	73:  Memory allocation failed.
	74:  An internal error has occurred in this function.  Check the function call and array sizes. If the function call is correct, please consult NAG for assistance.

	successfully call of the nag_mv_discrim_group function.	
*/
int nag_mv_discrim_group(
    Nag_DiscrimMethod type,
    Nag_GroupCovars equal,
    Nag_PriorProbability priors,
    int nvar, // the number of variables, p, in the variance-covariance matrices as specified to nag mv discrim (g03dac).
    int ng,	// the number of groups.
    int nig[], // the number of observations in each group training set.
    double gmean[], // the jth row of gmean contains the means of the p variables for the jth group, for j = 1, 2,.., nj.
    int tdg, // the last dimension of the array gmean.
    double gc[], // the first p(p + 1)/2 elements of gc should contain the upper triangular matrix R and the next ng blocks of p(p + 1)/2 elements should contain the upper triangular matrices Rj .
    double det[], // if equal = Nag NotEqualCovar the logarithms of the determinants of the within-group variance-covariance matrices as returned by nag mv discrim (g03dac).
    int nobs, // the number of observations in x which are to be allocated.
    int m,	// the number of variables in the data array x.
    int isx[], // indicates if the lth variable in x is to be included in the distance calculations.
    double x[], // contain the kth observation for the lth variable, for k = 1, 2, . . . ,nobs; l = 1, 2, . . . ,m.
    int tdx, // the last dimension of the array x.
    double prior[], // if priors = Nag UserPrior, the prior probabilities for the ng groups.
    double p[], // contains the posterior probability p[k][j] for allocating the kth observation to the jth group, for k = 1, 2, . . . ,nobs; j = 1, 2, . . ., ng.
    int tdp, // the last dimension of the array p and ati.
    int iag[], // the groups to which the observations have been allocated.
    Boolean atiq, // atiq must be TRUE if atypicality indices are required. .
    double ati[] //  the atypicality index for the kth observation with respect to the jth group, 
    			 // for k = 1, 2, . . . ,nobs; j = 1, 2, . . ., ng if atiq is TRUE, 
);

/** g03eac
		computes a distance (dissimilarity) matrix.
		
Example:
	A data matrix of five observations and three variables is read in
	and a distance matrix is calculated from variables 2 and 3 using 
	squared Euclidean distance with no scaling. This matrix is then printed.
	
void test_nag_mv_distance_mat()
{
	double d[45], s[10];
	matrix x;
	x.SetSize(10,10);
	int isx[10];
	int i, j, m, n;
	int tdx = 10;
	char char_scale;
	char char_update;
	char char_dist;
	Nag_MatUpdate update;
	Nag_DistanceType dist;
	Nag_VarScaleType scale;

	n = 5;
	m = 3;
	char_update = 'I';
	char_dist = 'S';
	char_scale = 'U';	
	
	x[0][0] = 1.0;
	x[0][1] = 1.0;
	x[0][2] = 1.0;
	
	x[1][0] = 2.0;
	x[1][1] = 1.0;
	x[1][2] = 2.0;
	
	x[2][0] = 3.0;
	x[2][1] = 6.0;
	x[2][2] = 3.0;
	
	x[3][0] = 4.0; 
	x[3][1] = 8.0;
	x[3][2] = 2.0;
  
	x[4][0] = 5.0;
	x[4][1] = 8.0;
	x[4][2] = 0.0; 
	
	for (i = 0; i < m; ++i)
		isx[i] = 1;
	isx[0] = 0;
	for (i = 0; i < m; ++i)
		s[i] = 1.0;
	
	if (char_update == 'U')
		update = Nag_MatUp;
	if (char_update == 'I')
		update = Nag_NoMatUp;
	if (char_dist == 'A')
		dist = Nag_DistAbs;
	if (char_dist == 'E')
		dist = Nag_DistEuclid;
	if (char_dist == 'S')
		dist = Nag_DistSquared;
	if (char_scale == 'S')
		scale = Nag_VarScaleStd;
	if (char_scale == 'R')
		scale = Nag_VarScaleRange;
	if (char_scale == 'G')
		scale = Nag_VarScaleUser;
	if (char_scale == 'U')
		scale = Nag_NoVarScale;

	nag_mv_distance_mat(update, dist, scale, n, m, x, tdx, isx, s, d);
	
	printf("\n");
	printf(" Distance matrix ");
	printf("\n");
	printf("\n");
	printf(" %s\n","      1	   2    3	 4");
	printf("\n");
	for (i = 2; i <= n; ++i)
	{
		printf("%2ld ",i);
		for (j=(i-1)*(i-2)/2+1; j<=i*(i - 1)/2; ++j)
			printf("%5.2f ",d[j-1]);
		printf("\n");
	}
}

	The output is following:
	
	Distance matrix 

  	1 	   2 	  3 	 4

 2  1.00 
 3 29.00 26.00 
 4 50.00 49.00  5.00 
 5 50.00 53.00 13.00  4.00 

Return:
	The function returns NAG error code or 0 if no error. 
	
	70: On entry, parameter dist had an illegal value.  On entry, parameter update had an illegal value.  On entry, parameter scale had an illegal value.
	11: On entry, n must not be less than 2: n = _value_.
	12: On entry, m must not be less than or equal to 0: m = _value_.
	17: On entry, tdx = _value_ while m = _value_.  These parameters must satisfy tdx = m.
	111: On entry, isx[_value_] = _value_.  Constraint: isx[i - 1] > 0, for at least one i, i = 1, 2, . . . ,m.
	115: On entry, d[_value_] = _value_.  Constraint: d[i - 1] = 0, i = 1, 2, . . . ,n*(n-1)/2 when update = Nag MatUp.  On entry, s[_value_] = _value_.  Constraint: s[j - 1] > 0, j = 1, 2, . . . ,m when scale = Nag VarScaleUser and isx[j - 1] > 0.  
	523: On entry, scale = Nag VarScaleRange or scale = Nag VarScaleStd, and x[i - 1][j - 1] = x[i][j - 1], for i = 1, 2, . . . , n - 1, for some j with isx[i - 1] > 0.
	74: An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
	
	successfully call of the nag_mv_distance_mat function.		
*/  
int nag_mv_distance_mat(
    Nag_MatUpdate update,
    Nag_DistanceType dist,
    Nag_VarScaleType scale,
    int n,		// the number of observations.
    int m,		// the total number of variables in array x.
    double x[], // the value of the jth variable for the ith object, for i = 1, 2, . . . , n; j = 1, 2, . . . ,m.
    int tdx,	// the last dimension of the array x.
    int isx[],	// indicates whether or not the jth variable in x is to be included in the distance computations.
    double s[],	// Input:if scale = Nag VarScaleUser and isx[j - 1] > 0 then s[j - 1] must contain the scaling for variable j, for j = 1, 2, . . . ,m.  Output:if scale = Nag VarScaleStd and isx[j - 1] > 0 then s[j - 1] contains the standard deviation of the variable in the jth column of x. 
    double d[] // Input:if update = Nag MatUp then d must contain the strictly lower triangle of the distance matrix D to be updated.  Output:the strictly lower triangle of the distance matrix D.
);
  
/** g03ecc
		performs hierarchical cluster analysis.
		
Example:
	Data consisting of three variables on five objects are read in. Euclidean squared 
	distances based on two variables are computed using nag mv distance mat (g03eac), 
	the objects are clustered using nag mv hierar cluster analysis and the dendrogram 
	computed using nag mv dendrogram (g03ehc).  The dendrogram is then printed.
	
void test_nag_mv_hierar_cluster_analysis()
{
	double cd[9], d[45], dord[10], s[10];
	matrix x;
	x.SetSize(10,10);
	double dmin_, dstep, ydist;
	int ilc[9], iord[10], isx[10], iuc[9];
	int nsym;
	int i, j;
	int m = 3, n = 5;
	int int_method;
	int tdx = 10;
	char name[10] = "ABCDE";
	char char_dist;
	char **c= NULL;
	char char_scale;
	char char_update;

	Nag_ClusterMethod method;
	Nag_MatUpdate update;
	Nag_DistanceType dist;
	Nag_VarScaleType scale;

	int_method = 5;
	
	if (int_method == 1)
		method = Nag_SingleLink;
	else if (int_method == 2)
		method = Nag_CompleteLink;
	else if (int_method == 3)
		method = Nag_GroupAverage;
	else if (int_method == 4)
		method = Nag_Centroid;
	else if (int_method == 5)
		method = Nag_Median;
	else
		method = Nag_MinVariance;
		
	char_update = 'I';
	
	if (char_update == 'U')
		update = Nag_MatUp;
	else
		update = Nag_NoMatUp;
	char_dist = 'S';
	if (char_dist == 'A')
		dist = Nag_DistAbs;
	else if (char_dist == 'E')
		dist = Nag_DistEuclid;
	else
		dist = Nag_DistSquared;
		
	char_scale = 'U';

	if (char_scale == 'S')
		scale = Nag_VarScaleStd;
	else if (char_scale == 'R')
		scale = Nag_VarScaleRange;
	else if(char_scale == 'G')
		scale = Nag_VarScaleUser;
	else
		scale = Nag_NoVarScale;

	x[0][0] = 1;
	x[0][1] = 5.0;
	x[0][2] = 2.0;
	
	x[1][0] = 2;
	x[1][1] = 1.0;
	x[1][2] = 1.0;
	
	x[2][0] = 3;
	x[2][1] = 4.0;
	x[2][2] = 3.0;
	
	x[3][0] = 4;
	x[3][1] = 1.0;
	x[3][2] = 2.0;
	
	x[4][0] = 5;
	x[4][1] = 5.0;
	x[4][2] = 0.0;

	for (i = 0; i < m; ++i)
		isx[i] = 1;
	isx[0] = 0;
	for (i = 0; i < m; ++i)
		s[i] = 1.0;
	
	nag_mv_distance_mat(update, dist, scale, n, m, x, tdx, isx, s, d);
	
	nag_mv_hierar_cluster_analysis(method, n, d, ilc, iuc, cd, iord, dord);
	
	printf("\n Distance Clusters Joined\n\n");
	for (i = 1; i <= n-1; ++i)
		printf("%10.3f %3c%3c\n",cd[i-1], name[ilc[i-1]-1], name[iuc[i-1]-1]);
		
	nsym = 20;
	dmin_ = 0.0;
	dstep = cd[n - 2] / (1.0 * nsym);

	nag_mv_dendrogram(Nag_DendSouth, n, dord, dmin_, dstep, nsym, &c);
	
	printf("\n");
	printf("Dendrogram ");
	printf("\n");
	printf("\n");
	ydist = cd[n - 2];
	for (i = 0; i < nsym; ++i)
	{
		if ((i+1) % 3 == 1)
		{
			printf("%10.3f%6s",ydist,"");
			printf("%s",c[i]);
			printf("\n");
		}
		else
		{
			printf("%16s%s","", c[i]);
			printf("\n");
		}
		ydist -= dstep;
	}
	printf("\n");
	printf("%14s","");
	for (i = 0; i < n; ++i)
		printf("%3c",name[iord[i]-1]);
	printf("\n");
	nag_mv_dend_free(&c);
}

	The output is as following:
 Distance Clusters Joined

     1.000   B  D
     2.000   A  C
     6.500   A  E
    14.125   A  B

Dendrogram 

    14.125            -------  
                      I     I  
                      I     I  
    12.006            I     I  
                      I     I  
                      I     I  
     9.887            I     I  
                      I     I  
                      I     I  
     7.769            I     I  
                   ---*     I  
                   I  I     I  
     5.650         I  I     I  
                   I  I     I  
                   I  I     I  
     3.531         I  I     I  
                   I  I     I  
                ---*  I     I  
     1.412      I  I  I  ---*  
                I  I  I  I  I  

                A  C  E  B  D

	
Return:
	The function returns NAG error code or 0 if no error. 
	
	70:  On entry, parameter method had an illegal value.
	11:  On entry, n must not be less than 2: n = _value_.
	115: On entry, d[_value_] = _value_.  Constraint: d[i - 1] = 0.0, i = 1, 2, . . . , n * (n - 1)/2.
	524: A true dendrogram cannot be formed because the distances at which clusters have merged are not increasing for all steps, i.e., cd[i - 1] < cd[i - 2] for some i = 2, 3, . . . , n - 1.  This can occur for the Nag Centroid and Nag Median methods.
	73:  Memory allocation failed.
	74:  An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
	
	successfully call of the nag_mv_hierar_cluster_analysis function.	
*/
int nag_mv_hierar_cluster_analysis(
    Nag_ClusterMethod method,
    int n,	// the number of objects.
    double d[], // the strictly lower triangle of the distance matrix. 
    int ilc[], // contains the number, j, of the cluster merged with cluster k (see iuc), j < k, at step l for l = 1, 2, . n - 1.
    int iuc[], // contains the number, k, of the cluster merged with cluster j, j < k, at step l for l = 1, 2, . . . , n - 1.
    double cd[], // the distance djk, between clusters j and k, j < k, merged at step l for l = 1, 2, . . . , n - 1.
    int iord[], // the objects in dendrogram order.
    double dord[] // the clustering distances corresponding to the order in iord. 
);

  
/** g03efc
		performs K-means cluster analysis.
		
Example:
	The data consists of observations of .ve variables on twenty soils 
	(Kendall and Stuart (1976)). The data is read in, the K-means 
	clustering performed and the results printed.
	
void test_nag_mv_kmeans_cluster_analysis()
{
	double css[3], csw[3], wt[20];
	double *wtptr;
	int nvar, i, j, k;
	int m, n;
	int inc[20], isx[5], nic[5];
	int maxit;
	int tdc = 5;
	int tdx = 5;
	char weight;	
	matrix x = {{77.3, 13.0, 9.7, 1.5, 6.4}, {82.5, 10.0, 7.5, 1.5, 6.5}, {66.9, 20.6, 12.5, 2.3, 7.0},
	{47.2, 33.8, 19.0, 2.8, 5.8}, {65.3, 20.5, 14.2, 1.9, 6.9}, {83.3, 10.0, 6.7, 2.2, 7.0},
	{81.6, 12.7, 5.7, 2.9, 6.7}, {47.8, 36.5, 15.7, 2.3, 7.2}, {48.6, 37.1, 14.3, 2.1, 7.2}, 
	{61.6, 25.5, 12.9, 1.9, 7.3}, {58.6, 26.5, 14.9, 2.4, 6.7}, {69.3, 22.3, 8.4, 4.0, 7.0},
	{61.8, 30.8, 7.4, 2.7, 6.4}, {67.7, 25.3, 7.0, 4.8, 7.3}, {57.2, 31.2, 11.6, 2.4, 6.5},
	{67.2, 22.7, 10.1, 3.3, 6.2}, {59.2, 31.2, 9.6, 2.4, 6.0}, {80.2, 13.2, 6.6, 2.0, 5.8}, 
	{82.2, 11.1, 6.7, 2.2, 7.2}, {69.7, 20.7, 9.6, 3.1, 5.9}};	
	matrix cmeans = {{82.5, 10.0, 7.5, 1.5, 6.5}, {47.8, 36.5, 15.7, 2.3, 7.2}, {67.2, 22.7, 10.1, 3.3, 6.2}};
	
	weight = 'U';
	n = 20;
	m = 5;
	nvar = 5;
	k = 3;
	maxit = 10;
	
	for (j = 0; j < m; ++j)
		isx[j] = 1;
	
	nag_mv_kmeans_cluster_analysis(n, m, x, tdx, isx, nvar, k, cmeans, tdc, wtptr, 
								inc, nic, css, csw, maxit);
			
	printf("\nThe cluster each point belongs to\n");
	for (i = 0; i < n; ++i)
	{	
		printf(" %6ld",inc[i]);
		if((i + 1) % 10 == 0)
			printf("\n");
	}
		
	printf("\n\nThe number of points in each cluster\n");
	for (i = 0; i < k; ++i)
		printf(" %6ld",nic[i]);
	printf("\n\nThe within-cluster sum of weights of each cluster\n");
	for (i = 0; i < k; ++i)
		printf(" %9.2f",csw[i]);
	printf("\n\nThe within-cluster sum of squares of each cluster\n\n");
	for (i = 0; i < k; ++i)
		printf(" %13.4f",css[i]);
	printf("\n\nThe final cluster centres\n");
	printf(" 1 2 3 4 5\n");
	for (i = 0; i < k; ++i)
	{
		printf(" %5ld ",i+1);
		for (j = 0; j < nvar; ++j)
			printf("%8.4f",cmeans[i][j]);
		printf("\n");
	}	
}

	The output is following:
	
	The cluster each point belongs to
      1      1      3      2      3      1      1      2      2      3
      3      3      3      3      3      3      3      1      1      3


	The number of points in each cluster
      6      3     11

	The within-cluster sum of weights of each cluster
      6.00      3.00     11.00

	The within-cluster sum of squares of each cluster

      46.5717       20.3800      468.8964

	The final cluster centres
 			1 		2 		3 		4 		5
     1  81.1833 11.6667  7.1500  2.0500  6.6000
     2  47.8667 35.8000 16.3333  2.4000  6.7333
     3  64.0455 25.2091 10.7455  2.8364  6.6545
	
Return:		
	The function returns NAG error code or 0 if no error. 
	
	11: On entry, n must not be less than 2: n = _value_.  On entry, k must not be less than 2: k = _value_.  On entry, nvar must not be less than 1: nvar = _value_.
	12: On entry, maxit must not be less than or equal to 0: maxit = _value_.
	17: On entry, m = _value_ while nvar = _value_.  These parameters must satisfy m = nvar.  On entry, tdx = _value_ while m = _value_.  These parameters must satisfy tdx = m.  On entry, tdc = _value_ while nvar = _value_.  These parameters must satisfy tdc = nvar.
	501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_.  Constraint: these two numbers must be the same.
	500: On entry, wt[_value_] = _value_.  Constraint: When referenced, all elements of wt must be non-negative.
	525: At least two elements of wt must be greater than zero.
	526: At least one cluster is empty after the initial assignment.  Try a di.erent set of initial cluster centres in cmeans and also consider decreasing the value of k. The empty clusters may be found by examining the values in nic.
	3: Too many iterations (_value_).  Convergence has not been achieved within the maximum number of iterations given by maxit.  Try increasing maxit and, if possible, use the returned values in cmeans as the initial cluster centres.  
	73: Memory allocation failed.  
	74: An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
	
	successfully call of the nag_mv_kmeans_cluster_analysis function.
*/
int nag_mv_kmeans_cluster_analysis(
    int n, // the number of observations.
    int m, // the number of variables in the array x.
    double x[], // the value of jth variable for the ith object for i = 1, 2, . . ., n; j = 1, 2, . . . ,m.
    int tdx, // the last dimension of the array x.
    int isx[], // indicates whether or not the jth variable is to be included in the analysis.
    int nvar, // the number of variables included in the sum of squares calculations.
    int k,	// the number of clusters.
    double cmeans[], // Input: contain the value of the jth variable for the ith initial cluster centre, for i = 1, 2, . . .,K; j = 1, 2, . . . , p.  Output: contains the value of the jth variable for the ith computed cluster centre, for i = 1, 2, . . .,K; j = 1, 2, . . . , p.
    int tdc, // the last dimension of the array cmeans.
    double *wtptr, // contain the weights to be used in the analysis.
    int inc[],	// the cluster to which the ith object has been allocated, for i = 1, 2, . . . , n.
    int nic[], 	// the number of objects in the ith cluster, for i = 1, 2, . . .,K. css[k]
    double css[], // the within-cluster (weighted) sum of squares of the ith cluster, for i = 1, 2, . . .,K.
    double csw[], // s the within-cluster sum of weights of the ith cluster, for i = 1, 2, . . .,K. 
    int maxit 	  // the maximum number of iterations allowed in the analysis.	
);
 
/** g03ehc
		produces a dendrogramfrom the results of nag mv hierar cluster analysis.
		
Example:
	Data consisting of three variables on five objects are read in. Euclidean squared 
	distances are computed using nag mv distance mat (g03eac) and median clustering 
	performed by nag mv hierar cluster analysis (g03ecc). nag mv dendrogram (g03ehc) 
	is used to produce a dendrogram with orientation east and a dendrogram with 
	orientation south. The two dendrograms are printed.  Note the use of 
	nag mv dend free (g03xzc)to free the memory allocated internally to the 
	character array pointed to by c.
	
void test_nag_mv_dendrogram()
{
	double cd[9], d[45], dord[10], s[10];
	matrix x;
	x.SetSize(10,10);
	double dmin_;
	double dstep;
	int ilc[9], iord[10], isx[10], iuc[9];
	int  nsym;
	int i, j;
	int m = 3, n = 5;
	int int_method;
	int tdx;
	char char_dist;
	char char_scale;
	char char_update;
	char **c = NULL;
	Nag_ClusterMethod method;
	Nag_MatUpdate update;
	Nag_DistanceType dist;
	Nag_VarScaleType scale;
	tdx = 10;

	x[0][0] = 1.0;
	x[0][1] = 1.0;
	x[0][2] = 1.0;
	
	x[1][0] = 2.0;
	x[1][1] = 1.0;
	x[1][2] = 2.0;
	
	x[2][0] = 3.0;
	x[2][1] = 6.0;
	x[2][2] = 3.0;
	
	x[3][0] = 4.0; 
	x[3][1] = 8.0;
	x[3][2] = 2.0;
  
	x[4][0] = 5.0;
	x[4][1] = 8.0;
	x[4][2] = 0.0; 

	for (i = 0; i < m; ++i)
		isx[i] = 1;
	isx[0] = 0;
	for (i = 0; i < m; ++i)
		s[i] = 1;
	
	int_method = 5;
	if (int_method == 1)
		method = Nag_SingleLink;
	else if (int_method == 2)
		method = Nag_CompleteLink;
	else if (int_method == 3)
		method = Nag_GroupAverage;
	else if (int_method == 4)
		method = Nag_Centroid;
	else if (int_method == 5)
		method = Nag_Median;
	else
		method = Nag_MinVariance;
		
	char_update = 'I';
	char_dist = 'S';
	char_scale = 'U';

	if (char_update == 'U')
		update = Nag_MatUp;
	else
		update = Nag_NoMatUp;
	if (char_dist == 'A')
		dist = Nag_DistAbs;
	else if (char_dist == 'E')
		dist = Nag_DistEuclid;
	else
		dist = Nag_DistSquared;
	
	if (char_scale == 'S')
		scale = Nag_VarScaleStd;
	else if (char_scale == 'R')
		scale = Nag_VarScaleRange;
	else if (char_scale == 'G')
		scale = Nag_VarScaleUser;
	else
		scale = Nag_NoVarScale;

	dmin_ = 0.0;
	dstep = 1.1;
	nsym = 40;
	
	nag_mv_distance_mat(update, dist, scale, n, m, x, tdx, isx, s, d);
	
	nag_mv_hierar_cluster_analysis(method, n, d, ilc, iuc, cd, iord, dord);

	nag_mv_dendrogram(Nag_DendEast, n, dord, dmin_, dstep, nsym, &c);
	
	printf("\nDendrogram, Orientation East\n\n");
	for (i=0; i< n; i++)
	 	printf("%s\n", c[i]);
	dmin_ = 0.0;
	dstep = 1.0;
	nsym = 40;
	nag_mv_dend_free(&c);
	nag_mv_dendrogram(Nag_DendSouth, n, dord, dmin_, dstep, nsym, &c);
	printf("\n\n Dendrogram, Orientation South\n\n");
	for (i=0; i< nsym; i++)
	{
		printf("%s\n", c[i]);
	}	
	nag_mv_dend_free(&c);
}	

	The output is as following:
Dendrogram, Orientation East

                                        
        ...............................(
       (                         .......
       (                        (    ...
       (........................(...(...


 Dendrogram, Orientation South

               
               
               
   ----------  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I        I  
   I  ------*  
   I  I     I  
   I  I     I  
   I  I     I  
   I  I  ---*  
   I  I  I  I  
   I  I  I  I  
---*  I  I  I  
I  I  I  I  I  

Return:	
	The function returns NAG error code or 0 if no error. 
	
	70: On entry, parameter orient had an illegal value.  
	11: On entry, n must not be less than 2: n = _value_.  On entry, nsym must not be less than 1: nsym = _value_.
	5:  On entry, dmin must not be less than 0.0: dmin = _value_.
	6:  On entry, dstep must not be less than or equal to 0.0: dstep = _value_.
	527: On entry, n = _value_, dord[_value_] = _value_.  Constraint: dord[n-1] = dord[i - 1], i = 1, 2, . . . ,n-1.
	74: An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
	
	successfully call of the nag_mv_dendrogram function.	
*/
int nag_mv_dendrogram(
    Nag_DendOrient orient,
    int n,	// the number of objects in the cluster analysis.
    double dord[], // contains the distances, in dendrogram order, at which clustering takes place.
    double dmin_, // the clustering distance at which the dendrogram begins.
    double dstep, // the distance represented by one symbol of the dendrogram.
    int nsym, // the number of character positions used in the dendrogram. Hence the clustering distance at which the dendrogram terminates is given by dmin + nsym × dstep.
    char ***c // a pointer to an array of character pointers, containing consecutive lines of the dendrogram. 
);

  
/** g03ejc
  		computes a cluster indicator variable from the results of nag mv hierar cluster analysis (g03ecc).
  		
Example:
	Data consisting of three variables on five objects are input. Euclidean 
	squared distances are computed using nag mv distance mat (g03eac) and 
	median clustering performed using nag mv hierar cluster analysis (g03ecc). 
	A dendrogram is produced by nag mv dendrogram (g03ehc) and printed. 
	nag mv cluster indicator finds two clusters and the results are printed.
	
void test_nag_mv_cluster_indicator() 
{
	double cd[9], d[45], dord[10], s[10];
	matrix x;
	x.SetSize(10, 10);
	double dmin_;
	double dstep, ydist;
	double dlevel;
	int ic[10], ilc[9], iord[10], isx[10],iuc[9];
	int nsym;
	int i, j, k;
	int m = 3, n = 5;
	int int_method;
	int tdx = 10;
	char **c = NULL;
	char name[10] = "ABCDE";
	char char_dist;
	char char_scale;
	char char_update;
	Nag_ClusterMethod method;
	Nag_MatUpdate update;
	Nag_DistanceType dist;
	Nag_VarScaleType scale;

	int_method = 5;
	if (int_method == 1)
		method = Nag_SingleLink;
	else if (int_method == 2)
		method = Nag_CompleteLink;
	else if (int_method == 3)
		method = Nag_GroupAverage;
	else if (int_method == 4)
		method = Nag_Centroid;
	else if (int_method == 5)
		method = Nag_Median;
	else
		method = Nag_MinVariance;
	char_update = 'I';
	if (char_update == 'U')
		update = Nag_MatUp;
	else
		update = Nag_NoMatUp;
		
	char_dist = 'S';
	if (char_dist == 'A')
		dist = Nag_DistAbs;
	else if (char_dist == 'E')
		dist = Nag_DistEuclid;
	else
		dist = Nag_DistSquared;
	char_scale = 'U';
	if (char_scale == 'S')
		scale = Nag_VarScaleStd;
	else if (char_scale == 'R')
		scale = Nag_VarScaleRange;
	else if (char_scale == 'G')
		scale = Nag_VarScaleUser;
	else
		scale = Nag_NoVarScale;
			
	x[0][0] = 1;
	x[0][1] = 5.0;
	x[0][2] = 2.0;
	
	x[1][0] = 2;
	x[1][1] = 1.0;
	x[1][2] = 1.0;
	
	x[2][0] = 3;
	x[2][1] = 4.0;
	x[2][2] = 3.0;
	
	x[3][0] = 4;
	x[3][1] = 1.0;
	x[3][2] = 2.0;
	
	x[4][0] = 5;
	x[4][1] = 5.0;
	x[4][2] = 0.0;

	for (i = 0; i < m; ++i)
		isx[i] = 1;
	isx[0] = 0;
	for (i = 0; i < m; ++i)
		s[i] = 1.0;

	k = 2;
	dlevel = 0.0;
	nag_mv_distance_mat(update, dist, scale, n, m, x, tdx, isx, s, d);

	nag_mv_hierar_cluster_analysis(method, n, d, ilc, iuc, cd, iord, dord);
	
	printf("\nDistance Clusters Joined\n\n");
	for (i = 0; i < n-1; ++i)
	{
		printf("%10.3f ",cd[i]);
		printf("%3c",name[ilc[i]-1]);
		printf("%3c",name[iuc[i]-1]);
		printf("\n");
	}

	nsym = 20;
	dmin_ = 0.0;
	dstep = cd[n - 2] / (1.0* nsym);
	nag_mv_dendrogram(Nag_DendSouth, n, dord, dmin_, dstep, nsym, &c);
	printf("\n");

	printf("Dendrogram ");
	printf("\n");
	printf("\n");
	ydist = cd[n - 2];
	for (i = 0; i < nsym; ++i)
	{
		if ((i+1) % 3 == 1)
		{
			printf("%10.3f%6s",ydist,"");
			printf("%s",c[i]);
			printf("\n");
		}
		else
		{
			printf("%16s%s","", c[i]);
			printf("\n");
		}
		ydist -= dstep;
	}
	printf("\n");
	printf("%14s","");
	for (i = 0; i < n; ++i)
	{
		printf("%3c",name[iord[i]-1]);
		printf("");
	}
	printf("\n");
	
	nag_mv_cluster_indicator(n, cd, iord, dord, &k, &dlevel, ic);
	
	printf("\n%s%2i%s\n\n","Allocation to ",k," clusters");
	printf("Object Cluster\n\n");
	for (i = 0; i < n; ++i)
	{
		printf("%5s%c%5s","",name[i],"");
		printf("%ld ",ic[i]);
		printf("\n");
	}
}	

	The output is as following:
Distance Clusters Joined

     1.000   B  D
     2.000   A  C
     6.500   A  E
    14.125   A  B

Dendrogram 

    14.125            -------  
                      I     I  
                      I     I  
    12.006            I     I  
                      I     I  
                      I     I  
     9.887            I     I  
                      I     I  
                      I     I  
     7.769            I     I  
                   ---*     I  
                   I  I     I  
     5.650         I  I     I  
                   I  I     I  
                   I  I     I  
     3.531         I  I     I  
                   I  I     I  
                ---*  I     I  
     1.412      I  I  I  ---*  
                I  I  I  I  I  

                A  C  E  B  D

Allocation to  2 clusters

Object Cluster

     A     1 
     B     2 
     C     1 
     D     2 
     E     1 

Return:
	The functions returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 2: n = _value_.
	19: On entry, k = _value_ while n = _value_.  These parameters must satisfy k = n.
	118: On entry, dlevel = _value_, k = _value_.  Constraint: k = 0 and dlevel > 0.0.
	62: The sequence cd is not increasing: cd[_value_] = _value_, cd[_value_] = _value_.
	122: On entry, dlevel = _value_, cd[_value_] = _value_.  Trivial solution returned.
	120: On exit, k = 1 .  Trivial solution returned.
	121: On exit, k = _value_, n = _value_.  Trivial solution returned.
	116: Arrays cd and dord are not compatible.
	528: The precise number of clusters requested is not possible because of tied clustering distances. The actual number of clusters produced is _value_.
	74: An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
	
	successfully call of the nag_mv_cluster_indicator function.
*/
int nag_mv_cluster_indicator(
    int n, // the number of objects. 
    double cd[], // the clustering distances in increasing order as returned by nag mv hierar cluster analysis (g03ecc).
    int iord[], // the objects in the dendrogram order as returned by nag mv hierar cluster analysis
    double dord[], // the clustering distances corresponding to the order in iord.
    int *k, // Input: indicates if a speci.ed number of clusters is required.  Output: the number of clusters produced, k.
    double *dlevel, // Input: if k = 0, then dlevel must contain the distance at which clusters are produced.  Output: the distance at which the required number of clusters are found if k > 0. 
	int ic[]
);
 
/** g03fac
		performs a principal coordinate analysis also known as classical metric scaling.
		
Example:
	The data, given by Krzanowski (1990), are dissimilarities between water 
	vole populations in Europe.  The first two principal co-ordinates are 
	computed by nag mv prin coord analysis.
		
void test_nag_mv_prin_coord_analysis()
{
	double e[14];
	matrix x;
	x.SetSize(14, 14);
	int ndim;
	int tdx=14;
	int i, j, n;
	int nn;
	char char_roots;
	Nag_Eigenvalues roots;	
	double d[91] = {0.099, 0.033, 0.022, 0.183, 0.114, 0.042, 0.148, 0.224, 0.059, 0.068,
	0.198, 0.039, 0.053, 0.085, 0.051, 0.462, 0.266, 0.322, 0.435, 0.268, 0.025, 0.628, 
	0.442, 0.444, 0.406, 0.240, 0.129, 0.014, 0.113, 0.070, 0.046, 0.047, 0.034, 0.002, 
	0.106, 0.129, 0.173, 0.119, 0.162, 0.331, 0.177, 0.039, 0.089, 0.237, 0.071, 0.434, 
	0.419, 0.339, 0.505, 0.469, 0.390, 0.315, 0.349, 0.151, 0.430, 0.762, 0.633, 0.781,
	0.700, 0.758, 0.625, 0.469, 0.618, 0.440, 0.538, 0.607, 0.530, 0.389, 0.482, 0.579, 
	0.597, 0.498, 0.374, 0.562, 0.247, 0.383, 0.387, 0.084, 0.586, 0.435, 0.550, 0.530, 
	0.552, 0.509, 0.369, 0.471, 0.234, 0.346, 0.456, 0.090, 0.038};
		
	n = 14;
	ndim = 2;
	char_roots = 'L';
	nn = n * (n - 1) / 2;
	
	if (char_roots == 'L')
		roots = Nag_LargeEigVals;
	else
		roots = Nag_AllEigVals;
	nag_mv_prin_coord_analysis(roots, n, d, ndim, x, tdx, e);
	printf("\nScaled Eigenvalues\n\n");
	if (char_roots == 'L')
	{
		for (i = 0; i < ndim; ++i)
			printf("%10.4f",e[i]);
	}
	else
	{
		for (i = 0; i < n; ++i)
			printf("%10.4f",e[i]);
		printf("\n");
	}
	printf("\nCo-ordinates\n\n");
	for (i = 0; i < n; ++i)
	{
		for (j = 0; j < ndim; ++j)
			printf("%10.4f",x[i][j]);
		printf("\n");
	}
}	

	The output is as following:
Scaled Eigenvalues

    0.7871    0.2808
Co-ordinates

    0.2408    0.2337
    0.1137    0.1168
    0.2394    0.0760
    0.2129    0.0605
    0.2495   -0.0693
    0.1487   -0.0778
   -0.0514   -0.1623
    0.0115   -0.3446
   -0.0039    0.0059
    0.0386   -0.0089
   -0.0421   -0.0566
   -0.5158    0.0291
   -0.3180    0.1501
   -0.3238    0.0475

Return:		
	The function returns NAG error code or 0 if no error. 
	
	70:  On entry, parameter roots had an illegal value.
	11:  On entry, ndim must not be less than 1: ndim = _value_.
	18:  On entry, n = _value_ while ndim = _value_.  These parameters must satisfy n > ndim.
	17:  On entry, tdx = _value_ while ndim = _value_.  These parameters must satisfy tdx = ndim.
	115: On entry, d[_value_] = _value_.  Constraint: d[i - 1] = 0.0, i = 1, 2, . . . , n * (n - 1)/2.
	529: There are fewer than ndim eigenvalues greater than zero.  Try a smaller number of dimensions (ndim) or use non-metric scaling (nag mv ordinal multidimscale (g03fcc)).
	530: The computation of eigenvalues or eigenvectors has failed. 
	73:  Memory allocation failed.
	74:  An internal error has occurred in this function. Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	
	successfully call of the  nag_mv_prin_coord_analysis function.	
*/
int nag_mv_prin_coord_analysis(
    Nag_Eigenvalues roots,
    int n, // the number of objects in the distance matrix.
    double d[], // the lower triangle of the distance matrix D stored packed by rows.
    int ndim, // the number of dimensions used to represent the data.
    double x[], // the ith row contains k coordinates for the ith point, i = 1, 2, . . . , n.
    int tdx, // the last dimension of the array x.
    double eval[] // if roots = Nag AllEigVals, eval contains the n scaled eigenvalues of the matrix E.  If roots = Nag LargeEigVals, eval contains the largest k scaled eigenvalues of the matrix E.
);


/* g03fcc
	performs non-metric (ordinal) multidimensional scaling.  
		
Example
	The data, given by Krzanowski (1990), are dissimilarities between water 
	vole populations in Europe. Initial estimates are provided by the first 
	two principal co-ordinates computed by nag mv prin coord analysis (g03fac). 
	The two dimension solution is computed using nag mv ordinal multidimscale.

//note: this one requires the struct with pointer.
//	Has problem on this example code
void test_nag_mv_ordinal_multidimscale()
{
	#define NMAX 14
	#define MMAX 2
	#define NNMAX NMAX*(NMAX-1)/2
	#define X(I,J) x[(I-1)*NMAX + (J-1)]
	#define XTMP(I) xtmp[(I)-1]
	#define YTMP(I) ytmp[(I)-1]
	double dfit[364], wk[511];
	double	x[196];
	double stress;
	int ndim;
	int i, j, n;
	int nn;
	int tdx = NMAX;
	char char_type;
	Nag_ScaleCriterion type;
	double d[91] = {0.099, 0.033, 0.022, 0.183, 0.114, 0.042, 0.148, 0.224, 0.059, 0.068,
	0.198, 0.039, 0.053, 0.085, 0.051, 0.462, 0.266, 0.322, 0.435, 0.268, 0.025, 0.628, 
	0.442, 0.444, 0.406, 0.240, 0.129, 0.014, 0.113, 0.070, 0.046, 0.047, 0.034, 0.002, 
	0.106, 0.129, 0.173, 0.119, 0.162, 0.331, 0.177, 0.039, 0.089, 0.237, 0.071, 0.434, 
	0.419, 0.339, 0.505, 0.469, 0.390, 0.315, 0.349, 0.151, 0.430, 0.762, 0.633, 0.781,
	0.700, 0.758, 0.625, 0.469, 0.618, 0.440, 0.538, 0.607, 0.530, 0.389, 0.482, 0.579, 
	0.597, 0.498, 0.374, 0.562, 0.247, 0.383, 0.387, 0.084, 0.586, 0.435, 0.550, 0.530, 
	0.552, 0.509, 0.369, 0.471, 0.234, 0.346, 0.456, 0.090, 0.038};

	n = 14;
	ndim = 2;
	char_type = 'T'; 	
	
	nn = n * (n - 1) / 2;
	
	nag_mv_prin_coord_analysis(Nag_LargeEigVals, n, d, ndim, x, tdx, wk);
	if (char_type == 'T')
		type = Nag_Stress;
	else
		type = Nag_SStress;
		
	g03fcc(type, n, ndim, d, x, tdx, &stress, dfit,	E04_DEFAULT, NAGERR_DEFAULT);
	
	Vprintf("\n STRESS = %13.4e\n\n",stress);
	Vprintf("Co-ordinates\n\n");
	for (i = 1; i <= n; ++i)
	{
		for (j = 1; j <= ndim; ++j)
			Vprintf("%10.4f",X(i,j));
	3.g03fcc.4 [NP3275/5/pdf]
	g03 - Multivariate Methods g03fcc
		Vprintf("\n");
	}
	exit(EXIT_SUCCESS);
	}
	else
	{
	Vprintf("Incorrect input value of n.\n");
	exit(EXIT_FAILURE);
	}
	}
}

Return:
	The function returns NAG error code or 0 if no error. 

	70: On entry, parameter type had an illegal value.  
	11: On entry, ndim must not be less than 1: ndim = _value_.
	18: On entry, n = _value_ while ndim = _value_.  These parameters must satisfy n > ndim.
	17: On entry, tdx = _value_ while ndim = _value_.  These parameters must satisfy tdx = ndim.
	531: All elements of array d = 0.0.  Constraint: At least one element of d must be positive.
	73: Memory allocation failed.
	74: An internal error has occurred in this function.  Check the function call and array sizes. If the function call is correct, please consult NAG for assistance.  Additional error messages are output if the optimization fails to converge or if the options are set incorrectly, Details of these can be found in the nag opt conj grad (e04dgc) document.	



  int nag_mv_ordinal_multidimscale(
    Nag_ScaleCriterion type,
    int n,		// the number of objects in the distance matrix.
    int ndim,	// the number of dimensions used to represent the data.
    double d[], // the lower triangle of the distance matrix D stored packed by rows. 
    double x[],	// Input: the ith row must contain an initial estimate of the co-ordinates for the ith point, i = 1, 2, . . . , n. One method of computing these is to use nag mv prin coord analysis (g03fac).  Output: the ith row contains m co-ordinates for the ith point, i = 1, 2, . . . , n.
    int tdx,	// the last dimension of the array x as declared in the calling program.
    double *stress, // the value of STRESS or SSTRESS at the final iteration.
    double dfit[], // auxiliary outputs. 
    Nag_E04_Opt *options,
  );
*/

/** g03xzc
		frees Nag allocated memory for the dendrogram array in nag mv dendrogram (g03ehc).
	Example:
		For an example of the use of nag mv dend free, please refer 
		to the example program for nag mv dendrogram (g03ehc).
	Return:
		The function returns NAG error code or 0 if no error. 
*/
int nag_mv_dend_free(
    char ***c 	// Input: a pointer to an array of character pointers to which memory has been allocated internally by nag mv dendrogram (g03ehc).  Output: The memory allocated to each of the pointers in the pointer array is freed and the pointers are set to NULL.
);


/** g03zac
		produces standardized values (z-scores) for a data matrix.
		
Example:
	A 4 by 3 data matrix is input along with location and scaling values. 
	The first and third columns are scaled and the results printed.
	
void test_nag_mv_z_scores()
{
	double e[3] = {14.75, 0.0, 1050.0}, s[3] = {2.50, 0.0, 420.0};
	matrix z;
	z.SetSize(4,3);
	matrix x = {{15.0, 0.0, 1500.0}, {12.0, 1.0, 1000.0}, {18.0, 2.0, 1200.0}, {14.0, 3.0, 500.0}};
	int nvar = 2;
	int isx[3];
	int i, j, m = 3, n = 4;
	int tdx = 3;
	int tdz = 3;
	
	for (j = 0; j < m; ++j)
		isx[j] = 1;
	isx[1] = 0;
	
	nag_mv_z_scores(n, m, x, tdx, nvar, isx, s, e, z, tdz);
	
	printf("\nStandardized Values\n\n");
	for (i = 0; i < n; ++i)
	{
		for (j = 0; j < nvar; ++j)
			printf("%8.3f",z[i][j]);
		printf("\n");
	}	
}	

	The output is following:
	
	Standardized Values

	 0.100   1.071
  	-1.100  -0.119
  	 1.300   0.357
    -0.300  -1.310
		
Return:		
	The function returns NAG error code or 0 if no error. 
	
	11:  On entry, n must not be less than 1: n = _value_.  On entry, nvar must not be less than 1: nvar = _value_.
	17:  On entry, m = _value_ while nvar = _value_.  These parameters must satisfy m = nvar.  On entry, tdx = _value_ while m = _value_.  These parameters must satisfy tdx = m.  On entry, tdz = _value_ while nvar = _value_.  These parameters must satisfy tdz = nvar.
	501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_.  Constraint: these two numbers must be the same.
	123: On entry, isx[_value_] = _value_, s[_value_] = _value_.  Constraint: if isx[j - 1] =0 then s[j - 1] > 0.0, j=1,2,. . .,m.
	74:  An internal error has occurred in this function. Check the function call and array sizes.  If the function call is correct, please consult NAG for assistance.
	
	successfully call of the nag_mv_z_scores function.		
*/
int nag_mv_z_scores(
    int n,		// the number of observations in the data matrix.
    int m,		// the number of variables in the data array x.
    double x[], // contain the ith sample point for the jth variable x[i][j] , for i = 1, 2,. , n, j = 1, 2,  . ,m.
    int tdx,	// the last dimension of the array x.
    int nvar,	// the number of variables to be standardised.
    int isx[],	// indicates whether or not the observations on the jth variable are included in the matrix of standardized values.
    double s[], // the scaling (standard deviation), sj, for the jth variable.
    double e[], // the location shift (mean), µj, for the jth variable.
    double z[], // the matrix of standardized values (z-scores).
    int tdz	// the last dimension of the array z 
);

/* end proto */
#endif /* not NAGG03 */
