/*------------------------------------------------------------------------------*
 * File Name: OCN_f.h															*
 * Creation: TCZ 7/20/2001														*
 * Purpose: Origin C Header for NAG functions									*
 * Copyright (c) OriginLab Corp.	2001										*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *------------------------------------------------------------------------------*/

#ifndef _O_NAG_f_H
#define _O_NAG_f_H

#importdll "ONAG" // NAG DLL prepared by OriginLab

#include <NAG\nag_types.h>

/* begin proto */

/**	f01bnc
		computes a Cholesky factorization of a complex positive-definite Hermitian matrix.
		
Example:
	This example computes the Cholesky factorization of the well-conditioned positive-definite Hermitian matrix
	.
	..
	15 1 - 2i 2 -4 + 3i
	1 + 2i 20 -2 + i 3 - 3i
	2 -2 - i 18 -1 + 2i
	-4 - 3i 3 + 3i -1 - 2i 26
	.
	..
	.

void test_nag_complex_cholesky()
{
	int i, j, n;
	matrix<complex> a;
	a.SetSize(8,8);
	double p[8];
	
	printf("f01bnc Example Program Results\n");
	a[0][0] = 15.0+0.0i;
	a[1][0] = 1.0+2.0i;
	a[1][1] = 20.0+0.0i;
	a[2][0] = 2.0+0.0i;
	a[2][1] = -2.0-1.0i;
	a[2][2] = 18.0+0.0i;
	a[3][0] = -4.0-3.0i;
	a[3][1] = 3.0+3.0i;
	a[3][2] = -1.0-2.0i;
	a[3][3] = 26.0+0.0i;
	n = 4;
	nag_complex_cholesky(n,a,8,p);
	complex gh;
	printf("\n Upper triangle of Complex matrix U by rows\n");
	for (i=0; i<n; ++i)
	{
		printf("\n");
		printf("   (%7.4f,%9.4f)\n", 1.0/p[i], 0.0);
		for (j=i+1; j<n; ++j)
		{
			gh = a[i][j];
			printf("   (%7.4f,%9.4f)\n",gh.m_re, gh.m_im);
		}
	}

}

	The output is following:
	
	Upper triangle of Complex matrix U by rows
	
	( 3.8730, 0.0000)
	( 0.2582, -0.5164)
	( 0.5164, 0.0000)
	(-1.0328, 0.7746)
	( 4.4347, 0.0000)
	(-0.4811, 0.1654)
	( 0.8268, -0.6013)
	( 4.1803, 0.0000)
	( 0.0073, 0.3463)
	( 4.8133, 0.0000)

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	370: Matrix diagonal element a [_value_][_value_] has non-zero imaginary part.
	358: The matrix is not positive-de.nite, possibly due to rounding errors.
	11: On entry, n must not be less than 1: n = _value_.  
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.

	successfully call of the nag_complex_cholesky function.


*/
int nag_complex_cholesky(
int n, // the order of the matrix A.
complex a[], // Input:the lower triangle of the n by n positive-definite Hermitian matrix A. The elements of the array above the diagonal need not be set.  Output: the off-diagonal elements of the upper triangular matrix U. The lower triangle of A is unchanged.
int tda, // the second dimension of the array a as declared in the function from which nag complex cholesky is called.
double p[] // the reciprocals of the real diagonal elements of U.
);

/**	f01mcc
		computes the Cholesky factorization of a real symmetric 
		positive-de.nite variable-bandwidth matrix.

Example:
	To obtain the Cholesky factorization of the symmetric matrix, whose lower triangle is
	.
	......
	1
	2 5
	0 3 13
	0 0 0 16
	5 14 18 8 55
	0 0 0 24 17 77
	.
	......
	.
	For this matrix, the elements of row must be set to 1, 2, 2, 1, 5, 3, and the elements
	within the envelope must be supplied in row order as 1, 2, 5, 3, 13, 16, 5, 14, 18, 8, 
	55, 24, 17, 77.

void test_nag_real_cholesky_skyline()
{

	int i, k, k1, k2, lal, n;
	double a[] = {1.0, 2.0, 5.0, 3.0, 13.0, 16.0, 5.0, 14.0,
				  18.0, 8.0, 55.0, 24.0, 17.0, 77.0};
	double al[36], d[8];
	n = 6;
	int row[] = {1, 2, 2, 1, 5, 3, 0, 0};
	k2 =14;	
	lal = k2;

	int success = nag_real_cholesky_skyline(n, a, lal, row, al, d);
	printf("\n");
	if(success == 0)
	{
		printf(" i d[i] Rowi of unit lower triangle\n\n");
		k2 = 0;
		for (i=0; i<n; ++i)
		{
			k1 = k2;
			k2 = k2+row[i];
			printf(" %3ld%8.3f", i, d[i]);
			for (k=k1; k<k2; k++)
				printf("%8.3f", al[k]);
			printf("\n");
		}
	}
	else
	printf("there is some problem with the function call");

}

	The output is following:
		
	i 	 d[i] 	Rowi of unit lower triangle
		
	0 	 1.000 	1.000
	1 	 1.000 	2.000  1.000
	2 	 4.000 	3.000  1.000
	3 	16.000 	1.000
	4 	 1.000 	5.000  4.000  1.500  0.500  1.000
	5 	16.000 	1.500  5.000  1.000

Parameters:

Return:
	This function returns NAG error code, 0 if no error.

	11: On entry, n must not be less than 1: n = _value_.  On entry, row[_value_] must not be less than 1: row[_value_] = _value_.
	13: On entry, row[_value_] = _value_ while i = _value_.  These parameters must satisfy row[i] = i + 1.
	17: On entry, lal = _value_while row[0]+ . . .+row[n-1] = _value_. These parameters must satisfy lal = row[0]+ . . .+row[n - 1].
	359: The matrix is not positive-de.nite, possibly due to rounding errors. The factorization has been completed but may be very inaccurate.
	358: The matrix is not positive-de.nite, possibly due to rounding errors.

	successfully call of the nag_real_cholesky_skyline function.

*/

int  nag_real_cholesky_skyline(
int n, // the order of the matrix A.
double a[], // the elements within the envelope of the lower triangle of the positive-de.nite symmetric matrix A, taken in row by row order.
int lal, // the smaller of the dimensions of the arrays a and al as declared in the function from which nag real cholesky skyline is called.
int row[], // contain the width of row i of the matrix A.
double al[], // : the elements within the envelope of the lower triangular matrix L, taken in row by row order.
double d[] // the diagonal elements of the diagonal matrix D. Note that the determinant of A is equal to the product of these diagonal elements.

);


/**	f01qcc
		finds the QR factorization of the real m by n matrix A, where m = n.

Example:
	To obtain the QR factorization of the 5 by 3 matrix
	A =
	.
	....
	2.0 2.5 2.5
	2.0 2.5 2.5
	1.6 -0.4 2.8
	2.0 -0.5 0.5
	1.2 -0.3 -2.9


void test_nag_real_qr()
{
	double zeta[10];
	matrix a;
	a.SetSize(20,10);
	int i, j, m, n;
	m = 5;
	n = 3; 
	a[0][0] = 2.0; 
	a[0][1] = 2.5;
	a[0][2] = 2.5;
	a[1][0] = 2.0;
	a[1][1] = 2.5;
	a[1][2] = 2.5;
	a[2][0] = 1.6;
	a[2][1] = -0.4;
	a[2][2] = 2.8;
	a[3][0] = 2.0;
	a[3][1] = -0.5;
	a[3][2] = 0.5;
	a[4][0] = 1.2;
	a[4][1] = -0.3;
	a[4][2] = -2.9;
 
	nag_real_qr(m, n, a, 10, zeta);
	printf("QR factorization of A\n\n");
	printf("Vector zeta\n");
	for (i = 0; i < n; ++i)
		printf(" %8.4f", zeta[i]);
	printf("\n\n");
	printf("Matrix A after factorization (upper triangular part is R)\n");
	for (i = 0; i < m; ++i)
	{
		for (j = 0; j < n; ++j)
			printf(" %8.4f", a[i][j]);
		printf("\n");
	}
}

	The output is following:
	
	QR factorization of A
	
	Vector zeta
	
	1.2247 1.1547 1.2649
	
	Matrix A after factorization (upper triangular part is R)
	
	-4.0000 	-2.0000 	-3.0000
	 0.4082 	-3.0000 	-2.0000
	 0.3266 	-0.4619 	-4.0000
	 0.4082 	-0.5774 	 0.0000
	 0.2449 	-0.3464 	-0.6325
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	17: On entry, m = _value_ while n = _value_. These parameters must satisfy m = n.  On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.
	11: On entry, n must not be less than 0: n = _value_.
	
	successfully call of the nag_real_qr function.


*/

int  nag_real_qr(
int m, // the number of rows of A.
int n, // the number of columns of A.
double a[], // Input: the leading m by n part of the array a must contain the matrix to be factorized.  Output: the n by n upper triangular part of a will contain the upper triangular matrix R and the m by n strictly lower triangular part of a will contain details of the factorization as described in Section 3.
int tda, // the second dimension of the array a as declared in the function from which nag real qr is called.
double zeta[] // contains the scalar zeta_k for the kth transformation. 
);


/**	f01qdc
		performs one of the transformations B := QB or B := QT B, where B 
		is an m by ncolb real matrix and Q is an m by m orthogonal matrix, 
		given as the product of Householder transformation matrices.

Example:
	To obtain the matrix QTB for the matrix B given by
	B =
	.
	....
	1.1 0 0 .00
	0.90 0.00
	0.60 1.32
	0.00 1.10
	-0.80 -0.26
	.
	....
	following the QR factorization of the 5 by 3 matrix A given by
	A =
	.
	....
	2.0 2.5 2.5
	2.0 2.5 2.5
	1.6 -0.4 2.8
	2.0 -0.5 0.5
	1.2 -0.3 -2.9
	.

void test_nag_real_apply_q()
{
	double zeta[10];
	matrix a, b;
	a.SetSize(20,10);
	b.SetSize(20,5);
	int i, j, m, n, ncolb;
	m = 5;
	n = 3;
	a[0][0] =2.0;
	a[0][1] =2.5;
	a[0][2] =2.5;
	a[1][0] =2.0;
	a[1][1] =2.5;
	a[1][2] =2.5;
	a[2][0] =1.6;
	a[2][1] =-0.4;
	a[2][2] =2.8;
	a[3][0] =2.0;
	a[3][1] =-0.5;
	a[3][2] =0.5;
	a[4][0] =1.2;
	a[4][1] =-0.3;
	a[4][2] =-2.9;
	
	b[0][0] =1.1;
	b[0][1] =0.0;

	b[1][0] =0.9;
	b[1][1] =0.0;

	b[2][0] =0.6;
	b[2][1] =1.32;

	b[3][0] =0.0;
	b[3][1] =1.1;

	b[4][0] =-0.8;
	b[4][1] =-0.26;
	ncolb = 2;
	nag_real_qr(m, n, a, 10, zeta);
	nag_real_apply_q(Transpose, Nag_ElementsSeparate, m, n, a, 10, zeta,
					ncolb, b, 5);
	printf("Matrix Q'*B\n");
	for (i = 0; i < m; ++i)
	{
		for (j = 0; j < ncolb; ++j)
			printf(" %8.4f", b[i][j]);
		printf("\n");
	}
}

	The output is following:
	
	Matrix Q'*B
	
	-1.0000 	-1.0000
	-1.0000 	 1.0000
	-1.0000 	-1.0000
	-0.1000 	 0.1000
	-0.1000	 	-0.1000

	
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	70: On entry, parameter trans had an illegal value.  On entry, parameter wheret had an illegal value.
	17: On entry, m = _value_ while n = _value_. These parameters must satisfy m = n.  On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdb = _value_ while ncolb = _value_. These parameters must satisfy tdb = ncolb.
	11: On entry, n must not be less than 0: n = _value_.  On entry, ncolb must not be less than 0: ncolb = _value_.
	73: Memory allocation failed.

	successfully call of the nag_real_apply_q function.

*/
	int  nag_real_apply_q(
	MatrixTranspose trans, 
	Nag_WhereElements wheret,
	int m, // the number of rows of A.
	int n, // the number of columns of A.
	double a[], // contain details of the matrix Q.
	int tda, // the second dimension of the array a as declared in the function from which nag real apply q is called.
	double zeta[], // contain the elements of zeta.
	int ncolb, // the number of columns of B.
	double b[], // Input: the leading m by ncolb part of the array b must contain the matrix to be transformed.  Output: b is overwritten by the transformed matrix.
	int tdb // the second dimension of the array b as declared in the function from which nag real apply q is called.
	
	);


/**	f01qec
		returns the .rst ncolq columns of the real m by m orthogonal matrix Q,
		where Q is given as the product of Householder transformation matrices.


Example:
	To obtain the 5 by 5 orthogonal matrix Q following the QR factorization 
	of the 5 by 3 matrix A given by
	A =
	.
	....
	2.0 2.5 2.5
	2.0 2.5 2.5
	1.6 -0.4 2.8
	2.0 -0.5 0.5
	1.2 -0.3 -2.9
	.
	....

void test_nag_real_form_q()
{
	double zeta[10];
	matrix a, q;
	a.SetSize(20,10);
	q.SetSize(20,20);
	int i, j, m, n, ncolq;
	m = 5;
	n = 3;
	a[0][0] =2.0;
	a[0][1] =2.5;
	a[0][2] =2.5;
	a[1][0] =2.0;
	a[1][1] =2.5;
	a[1][2] =2.5;
	a[2][0] =1.6;
	a[2][1] =-0.4;
	a[2][2] =2.8;
	a[3][0] =2.0;
	a[3][1] =-0.5;
	a[3][2] =0.5;
	a[4][0] =1.2;
	a[4][1] =-0.3;
	a[4][2] =-2.9;
	
	nag_real_qr(m, n, a, 10, zeta);

	for (j = 0; j < n; ++j)
		for (i = 0; i < m; ++i)
			q[i][j] = a[i][j];
	ncolq = m;
	nag_real_form_q(Nag_ElementsSeparate, m, n, ncolq, q, 20,
					zeta);
	printf("Matrix Q\n");
	for (i = 0; i < m; ++i)
	{
		for (j = 0; j < ncolq; ++j)
			printf(" %8.4f", q[i][j]);
		printf("\n");
	}
}

	The output is following:
	
	Matrix Q
	
	-0.5000 -0.5000  0.0000 -0.5000 -0.5000
	-0.5000 -0.5000  0.0000  0.5000  0.5000
	-0.4000  0.4000 -0.6000 -0.4000  0.4000
	-0.5000  0.5000  0.0000  0.5000 -0.5000
	-0.3000  0.3000  0.8000 -0.3000  0.3000



Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	70: On entry, parameter wheret had an illegal value.
	17: On entry, m = _value_ while n = _value_. These parameters must satisfy m = n.  On entry, tda = _value_ while max(n,ncolq) = _value_. These parameters must satisfy tda = max(n,ncolq).
	19: On entry, ncolq = _value_ while m = _value_. These parameters must satisfy ncolq = m.
	11: On entry, n must not be less than 0: n = _value_.  On entry, ncolq must not be less than 0: ncolq = _value_.
	73: Memory allocation failed.

	successfully call of the nag_real_form_q function.


*/
	int  nag_real_form_q(
	Nag_WhereElements wheret, 
	int m, // the number of rows of A.
	int n, // the number of columns of A.
	int ncolq, // the required number of columns of Q.  
	double a[], // Input: the required number of columns of Q.  When ncolq = 0 then an immediate return is effected.  Output: the first ncolq columns of the array a are overwritten by the first ncolq columns of the m by m orthogonal matrix Q. 
	int tda, // the second dimension of the array a as declared in the function from which nag real form q is called.
	double zeta[] // contain the elements of zeta.
	);


/**	f01rcc
		finds the QR factorization of the complex m by n matrix A, where m = n.

Example:
	To obtain the QR factorization of the 5 by 3 matrix
	A =
	.
	....
	0.5i -0.5 + 1.5i -1.0+ 1.0i
	0.4 + 0.3i 0.9 + 1.3i 0.2+ 1.4i
	0.4 -0.4 + 0.4i 1.8
	0.3 - 0.4i 0.1 + 0.7i 0.0
	- 0.3i 0.3 + 0.3i 2.4i
	.
	....


void test_nag_complex_qr()
{	
	int i, j, m, n;

	matrix<complex> theta;
	theta.SetSize(1,10);					 //{0.0-0.3i,0.3+0.3i,0.0+2.4i}};
	matrix<complex> a;
	a.SetSize(20,10);
	a[0][0] = 0.0+0.5i;
	a[0][1] = -0.5+1.5i;
	a[0][2] = -1.0+1.0i;
	
	a[1][0] = 0.4+0.3i;
	a[1][1] = 0.9+1.3i;
	a[1][2] = 0.2+1.4i;
	
	a[2][0] = 0.4+0.0i;
	a[2][1] = -0.4+0.4i;
	a[2][2] = 1.8+0.0i;
	
	a[3][0] = 0.3-0.4i;
	a[3][1] = 0.1+0.7i;
	a[3][2] = 0.0+0.0i;
	
	a[4][0] = 0.0-0.3i;
	a[4][1] = 0.3+0.3i;
	a[4][2] = 0.0+2.4i;
	m = 5;
	n = 3;
	printf("A=\n");
	for (i=0; i<m; ++i)
	{
		for (j=0; j<n; ++j)
		{	
			complex temp = a[i][j];
			
			printf(" (%7.4f,%8.4f)", temp.m_re, temp.m_im);
			if(j%3==2 || j==n-1)
				printf("\n");
		}
	}
	nag_complex_qr(m, n, a, 10, theta);
	printf("QR factorization of A\n");
	printf("Vector THETA\n");
	for (i=0; i<n; ++i)
	{
		complex temp = theta[0][i];
		printf(" (%7.4f,%8.4f)", temp.m_re, temp.m_im);
		if(i%3==2 || i==n-1)
			printf("\n");
	}
	printf("\nMatrix A after factorization (upper triangular part is R)\n");
	for (i=0; i<m; ++i)
	{
		for (j=0; j<n; ++j)
		{	
			complex temp = a[i][j];
			
			printf(" (%7.4f,%8.4f)", temp.m_re, temp.m_im);
			if(j%3==2 || j==n-1)
				printf("\n");
		}
	}
}

	The output is following:
	A=
	 ( 0.0000,  0.5000) (-0.5000,  1.5000) (-1.0000,  1.0000)
     ( 0.4000,  0.3000) ( 0.9000,  1.3000) ( 0.2000,  1.4000)
     ( 0.4000,  0.0000) (-0.4000,  0.4000) ( 1.8000,  0.0000)
     ( 0.3000, -0.4000) ( 0.1000,  0.7000) ( 0.0000,  0.0000)
     ( 0.0000, -0.3000) ( 0.3000,  0.3000) ( 0.0000,  2.4000)

	
	QR factorization of A
	Vector THETA
	
	( 1.0000, 0.5000) ( 1.0954, -0.3333) ( 1.2649, 0.0000)
	
	Matrix A after factorization (upper triangular part is R)
	
	( 1.0000, 0.0000) 	( 1.0000, 1.0000)  ( 1.0000, 1.0000)
	(-0.2000, -0.4000)  (-2.0000, 0.0000)  (-1.0000, -1.0000)
	(-0.3200, -0.1600)  (-0.3505, 0.2629)  (-3.0000, 0.0000)
	(-0.4000, 0.2000)   ( 0.0000, 0.5477)  ( 0.0000, 0.0000)
	(-0.1200, 0.2400)   ( 0.1972, 0.2629)  ( 0.0000, 0.6325)


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	17: On entry, m = _value_ while n = _value_. These parameters must satisfy m = n.  On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.
	11: On entry, n must not be less than 0: n = _value_.
	
	successfully call of the nag_complex_qr function.

*/
	int  nag_complex_qr(
	int m, // the number of rows of A.
	int n, // the number of columns of A.
	complex a[], // Input: the leading m by n part of the array a must contain the matrix to be factorized.  Output: the n by n upper triangular part of a will contain the upper triangular matrix R, with the imaginary parts of the diagonal elements set to zero, and the m by n strictly lower triangular part of a will contain details of the factorization as described above.
	int tda, // the second dimension of the array a as declared in the function from which nag complex qr is called.
	complex theta[]  // the scalar theta_k for the kth transformation. 
	);

 
/**	f01rdc
		) performs one of the transformations B := QB or B := QHB, where B 
		is an m by ncolb complex matrix and Q is an m by m unitary matrix, 
		given as the product of Householder transformation matrices.


Example:
	To obtain the matrix QHB for the matrix B given by
	B =
	.
	....
	-0.55 + 1.05i 0.45 + 1.05i
	0.49 + 0.93i 1.09 + 0.13i
	0.56 - 0.16i 0.64 + 0.16i
	0.39 + 0.23i -0.39 - 0.23i
	1.1 3 +0.83i -1.1 3 +0.77i
	.
	....
	following the QR factorization of the 5 by 3 matrix A given by
	A =
	.
	....
	0.5i -0.5 + 1.5i -1.0+ 1.0i
	0.4 + 0.3i 0.9 + 1.3i 0.2+ 1.4i
	0.4 -0.4 + 0.4i 1.8
	0.3 - 0.4i 0.1 + 0.7i 0.0
	-0.3i 0.3 + 0.3i 2.4i
	.
	....

void test_nag_complex_apply_q()
{


	int i,j,m, n,ncolb;
	matrix<complex> a, b, theta;
	a.SetSize(20,10);
	b.SetSize(20,5);
	theta.SetSize(1,10);


	m = 5;
	n = 3;
	
	a[0][0] = 0.0+0.5i;
	a[0][1] = -0.5+1.5i;
	a[0][2] = -1.0+1.0i;
	
	a[1][0] = 0.4+0.3i;
	a[1][1] = 0.9+1.3i;
	a[1][2] = 0.2+1.4i;
	
	a[2][0] = 0.4+0.0i;
	a[2][1] = -0.4+0.4i;
	a[2][2] = 1.8+0.0i;
	
	a[3][0] = 0.3-0.4i;
	a[3][1] = 0.1+0.7i;
	a[3][2] = 0.0+0.0i;
	
	a[4][0] = 0.0-0.3i;
	a[4][1] = 0.3+0.3i;
	a[4][2] = 0.0+2.4i;
	ncolb = 2;
	b[0][0] =-0.55 +1.05i;
	b[0][1] =0.45+1.05i; 
	 
	b[1][0] = 0.49 + 0.93i;
	b[1][1] = 1.09+0.13i;
	
	b[2][0] =0.56 - 0.16i;
	b[2][1] = 0.64+0.16i;
	
	b[3][0] =0.39+0.23i;
	b[3][1] = -0.39-0.23i;
	
	b[4][0] =1.13+0.83i;
	b[4][1] = -1.13+0.77i;

	
	nag_complex_qr(m,n, a, 10, theta);
	
	nag_complex_apply_q(ConjugateTranspose, Nag_ElementsSeparate, m,n,a, 10,
						theta, ncolb, b, 5);

	printf("\nMatrix conjg( Q' )*B\n");
	for (i=0; i<m; ++i)
	{
		for (j=0; j<ncolb; ++j)
		{
			complex temp = b[i][j];
			printf(" (%7.4f,%8.4f)", temp.m_re,temp.m_im);
			if(j%2==1 || j==n-1)
				printf("\n");
			
		}
	}

}

	The output is following:
	
	Matrix conjg( Q' )*B
	
	( 1.0000,1.0000) 	( 1.0000, -1.0000)
	(-1.0000,0.0000) 	(-1.0000,0.0000)
	(-1.0000,1.0000) 	(-1.0000, -1.0000)
	(-0.0600,-0.0200) 	(-0.0400, 0.1200)
	( 0.0400,0.1200) 	(-0.0600,0.0200)


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	70: On entry, parameter trans had an illegal value.  On entry, parameter wheret had an illegal value.
	17: On entry, m = _value_ while n = _value_. These parameters must satisfy m = n.  On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdb = _value_while ncolb = _value_. These parameters must satisfy tdb = ncolb.
	11: On entry, n must not be less than 0: n = _value_.  On entry, ncolb must not be less than 0: ncolb = _value_.
	73: Memory allocation failed.

	successfully call of the nag_complex_apply_q function.
*/
	int  nag_complex_apply_q(
	MatrixTranspose trans,
	Nag_WhereElements wheret,
	int m, // the number of rows of A.
	int n, // the number of columns of A.
	complex a[], // the leading m by n strictly lower triangular part of the array a must contain details of the matrix Q.  
	int tda, // the second dimension of the array a as declared in the function from which nag complex apply q is called.
	complex theta[], // contain the elements of theta. 
	int ncolb,  // Input: ncolb, the number of columns of B. When ncolb = 0 then an immediate return is effected.
	complex b[], // Input: the leading m by ncolb part of the array b must contain the matrix to be transformed.  Output: b is overwritten by the transformed matrix.
	int tdb // the second dimension of the array b as declared in the function from which nag complex apply q is called.
	);


/**	f01rec
		returns the .rst ncolq columns of the m by m unitary matrix Q, where
		Q is given as the product of Householder transformation matrices.


Example:
	To obtain the 5 by 5 unitary matrix Q following the QR factorization 
	of the 5 by 3 matrix A given by
	A =
	.
	....
	0.5i -0.5 + 1.5i -1.0+ 1.4i
	0.4 + 0.3i 0.9 + 1.3i 0.2+ 1.4i
	0.4 -0.4 + 0.4i 1.8
	0.3 - 0.4i 0.1 + 0.7i 0.0
	-0.3i 0.3 + 0.3i 2.4i
	.
	....
	

void test_nag_complex_form_q()
{
	int i, j, m, n, ncolq;
	matrix<complex>a, q, theta;
	a.SetSize(20,10);
	q.SetSize(20,20);
	theta.SetSize(1,10);
	m = 5;
	n = 3;
	a[0][0] = 0.0+0.5i;
	a[0][1] = -0.5+1.5i;
	a[0][2] = -1.0+1.0i;
	
	a[1][0] = 0.4+0.3i;
	a[1][1] = 0.9+1.3i;
	a[1][2] = 0.2+1.4i;
	
	a[2][0] = 0.4+0.0i;
	a[2][1] = -0.4+0.4i;
	a[2][2] = 1.8+0.0i;
	
	a[3][0] = 0.3-0.4i;
	a[3][1] = 0.1+0.7i;
	a[3][2] = 0.0+0.0i;
	
	a[4][0] = 0.0-0.3i;
	a[4][1] = 0.3+0.3i;
	a[4][2] = 0.0+2.4i;
	
	
	nag_complex_qr(m, n, a, 10, theta);
	
	for (j=0; j<n; ++j)
		for (i=0; i<m; ++i)
			q[i][j] = a[i][j]; 
	

	ncolq = m;
	nag_complex_form_q(Nag_ElementsSeparate, m, n, ncolq, q, 20,theta);

	printf("\nMatrix Q\n");
	for (i=0; i<m; ++i)
	{
		for (j=0; j<ncolq; ++j)
		{
			complex temp = q[i][j];
			printf(" (%5.2f,%5.2f)", temp.m_re, temp.m_im);
			if(j%5==4 || j==ncolq-1)
				printf("\n");
		}	
	}
}


	The output is following:
	
	Matrix Q
	
	( 0.00, 0.50) ( 0.00,-0.50) ( 0.00, 0.00) ( 0.50, 0.00) ( 0.40, 0.30)
	( 0.40, 0.30) (-0.40,-0.30) ( 0.00, 0.00) (-0.30, 0.40) (-0.48, 0.14)
	( 0.40, 0.00) ( 0.40, 0.00) (-0.60, 0.00) (-0.24,-0.32) ( 0.00, 0.40)
	( 0.30,-0.40) ( 0.30,-0.40) ( 0.00, 0.00) ( 0.50, 0.00) (-0.40,-0.30)
	( 0.00,-0.30) ( 0.00,-0.30) ( 0.00,-0.80) (-0.24, 0.18) ( 0.30, 0.00)


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	70: On entry, parameter wheret had an illegal value.
	17: On entry, m = _value_ while n = _value_. These parameters must satisfy m = n.  On entry, tda = _value_ while max(n,ncolq) = _value_. These parameters must satisfy tda = n.
	11: On entry, n must not be less than 0: n = _value_.  On entry, ncolq must not be less than 0: ncolq = _value_.
	19: On entry, ncolq = _value_ while m = _value_. These parameters must satisfy ncolq = m.
	73: Memory allocation failed.

	successfully call of the nag_complex_form_q function.
*/
	int  nag_complex_form_q(
	Nag_WhereElements wheret,
	int m, // the number of rows of A.
	int n, // the number of columns of A.
	int ncolq, // the required number of columns of Q.  
	complex a[],  // Input: the leading m by n strictly lower triangular part of the array a must contain details of the matrix Q.  Output: the .rst ncolq columns of the array a are overwritten by the first ncolq columns of the m by m unitary matrix Q. When n = 0 then the .rst ncolq columns of a are overwritten by the .rst ncolq columns of the unit matrix.
	int tda,  // the second dimension of the array a as declared in the function from which nag complex form q is called.
	complex theta[] // contain the elements of theta.
	);

/**	f02aac
		calculates all the eigenvalues of a real symmetric matrix.
		
Example:
	To calculate all the eigenvalues of the real symmetric matrix
	.
	..
	0.5 0.0 2.3 -2.6
	0.0 0.5 -1.4 -0.7
	2.3 -1.4 0.5 0.0
	-2.6 -0.7 0.0 0.5
	.
	..
	.

void test_nag_real_symm_eigenvalues()
{	
	int i, j, n;
	matrix a;
	a.SetSize(8,8);
	double r[8];

	n = 4;
	a[0][0] = 0.5;
	a[0][1] = 0.0;
	a[0][2] = 2.3;
	a[0][3] = -2.6;
	
	a[1][0] = 0.0;
	a[1][1] = 0.5;
	a[1][2] = -1.4;
	a[1][3] = -0.7;
	
	a[2][0] = 2.3;
	a[2][1] = -1.4;
	a[2][2] = 0.5;
	a[2][3] = 0.0;
	
	a[3][0] = -2.6;
	a[3][1] = -0.7;
	a[3][2] = 0.0;
	a[3][3] = 0.5;
		
	nag_real_symm_eigenvalues(n, a, 8, r);
	printf("Eigenvalues\n");
	for (i=0; i<n; i++)
	{
		printf("%9.4f",r[i]);
		if(i%8==7 || i==n-1)
			printf("\n");
	}	
}

	The output is following:
	
	Eigenvalues
	
	-3.0000 -1.0000 2.0000 4.0000

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	362: More than _value_ iterations are required to isolate all the eigenvalues.
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.
	73: Memory allocation failed.

	successfully call of the nag_real_symm_eigenvalues function.



*/

	int  nag_real_symm_eigenvalues(
	int n,  // the order of the matrix A.
	double a[], // Input: the lower triangle of the n by n symmetric matrix A. The elements of the array above the diagonal need not be set.  Output: the elements of A below the diagonal are overwritten, and the rest of the array is unchanged.
	int tda, // the second dimension of the array a as declared in the function from which nag real symm eigenvalues is called.
	double r[] // the eigenvalues in ascending order.
	);



/**	f02abc
		calculates all the eigenvalues and eigenvectors of a real
		symmetric matrix.

Example:
	To calculate all the eigenvalues and eigenvectors of the real symmetric matrix
	.
	..
	0.5 0.0 2.3 -2.6
	0.0 0.5 -1.4 -0.7
	2.3 -1.4 0.5 0.0
	-2.6 -0.7 0.0 0.5
	.
	..
	.

void test_nag_real_symm_eigensystem()
{	
	
	int i, j, n;
	matrix a, v;
	a.SetSize(8,8);
	double  r[8];
	v.SetSize(8,8);

	n = 4;
	a[0][0] = 0.5;
	a[0][1] = 0.0;
	a[0][2] = 2.3;
	a[0][3] = -2.6;
	
	a[1][0] = 0.0;
	a[1][1] = 0.5;
	a[1][2] = -1.4;
	a[1][3] = -0.7;
	
	a[2][0] = 2.3;
	a[2][1] = -1.4;
	a[2][2] = 0.5;
	a[2][3] = 0.0;
	
	a[3][0] = -2.6;
	a[3][1] = -0.7;
	a[3][2] = 0.0;
	a[3][3] = 0.5;

	
	nag_real_symm_eigensystem(n, a, 8, r, v, 8);
	printf("Eigenvalues\n");
	for (i=0; i<n; i++)
	{
		printf("%9.4f", r[i]);
		if(i%8==7 || i==n-1)
			printf("\n");
	}
		
	printf("Eigenvectors\n");
	for ( i=0; i<n; i++)
		for (j=0; j<n; j++)
		{
			printf("%9.4f", v[i][j]);
			if(j%8==7 || j==n-1)
				printf("\n");
		}
		
}

	The output is following:
	
	Eigenvalues
	
	-3.0000 -1.0000 2.0000 4.0000
	
	Eigenvectors
	
     0.7000   0.1000   0.1000  -0.7000
    -0.1000   0.7000   0.7000   0.1000
    -0.5000   0.5000  -0.5000  -0.5000
     0.5000   0.5000  -0.5000   0.5000


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	362: More than _value_ iterations are required to isolate all the eigenvalues.
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdv = _value_ while n = _value_. These parameters must satisfy tdv = n.
	73: Memory allocation failed.

	successfully call of the  nag_real_symm_eigensystem function.

*/

	int  nag_real_symm_eigensystem(
	int n, // the order of the matrix A.
	double a[], // the lower triangle of the n by n symmetric matrix A. The elements of the array above the diagonal need not be set. 
	int tda, // the second dimension of the array a as declared in the function from which nag real symm eigensystem is called.
	double r[], // the eigenvalues in ascending order.
	double v[], // the normalised eigenvectors, stored by columns; the ith column corresponds to the ith eigenvalue. The eigenvectors are normalised so that the sum of squares of the elements is equal to 1.
	int tdv // the second dimension of the array v as declared in the function from which nag real symm eigensystem is called.
	
	);


/**	f02adc
		calculates all the eigenvalues of Ax = lamdaBx, where A is a real 
		symmetric matrix and B is a real symmetric positive-definite matrix.


Example:
	To calculate all the eigenvalues of the general symmetric 
	eigenproblem Ax = lamdaBx where A is the symmetric matrix
	.
	..
	0.5 1.5 6.6 4.8
	1.5 6.5 16.2 8.6
	6.6 16.2 37.6 9.8
	4.8 8.6 9.8 -17.1
	.
	..
	and B is the symmetric positive-de.nite matrix
	.
	..
	1 3 4 1
	3 13 16 11
	4 16 24 18
	1 11 18 27
	.
	..
	.
	

void test_nag_real_symm_general_eigenvalues()
{		
	int i, j, n;
	matrix a, b;
	a.SetSize(8,8);
	b.SetSize(8,8);
	double r[8];
	n = 4;
	a[0][0] = 0.5;
	a[0][1] = 1.5;
	a[0][2] = 6.6; 
	a[0][3] = 4.8;
	
	
	a[1][0] = 1.5;
	a[1][1] = 6.5;
	a[1][2] = 16.2;
	a[1][3] = 8.6;
	
	a[2][0] = 6.6;
	a[2][1] = 16.2;
	a[2][2] = 37.6;
	a[2][3] = 9.8; 
	
	a[3][0] = 4.8;
	a[3][1] = 8.6;
	a[3][2] = 9.8; 
	a[3][3] = -17.1;
	
	b[0][0] = 1.0;
	b[0][1] = 3.0;
	b[0][2] = 4.0;
	b[0][3] = 1.0;
	
	
	b[1][0] = 3.0;
	b[1][1] = 13.0; 
	b[1][2] = 16.0;
	b[1][3] = 11.0;
	
	b[2][0] = 4.0;
	b[2][1] = 16.0; 
	b[2][2] = 24.0;
	b[2][3] = 18.0;
  	
	b[3][0] = 1.0;
	b[3][1] = 11.0;
	b[3][2] = 18.0;
	b[3][3] = 27.0;
	

	nag_real_symm_general_eigenvalues(n, a, 8, b, 8, r);
	printf("Eigenvalues\n");
	for (i=0; i<n; i++)
	{
		printf("%9.4f",r[i]);
		if(i%8==7 || i==n-1)
			printf("\n");
	}	
}

	The output is following:
	
	Eigenvalues
	
	-3.0000 -1.0000 2.0000 4.0000


Parameters:

Return:
	This function returns NAG error code, 0 if no error.

	358: The matrix B is not positive-de.nite, possibly due to rounding errors.
	362: More than _value_ iterations are required to isolate all the eigenvalues.
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdb = _value_ while n = _value_. These parameters must satisfy tdb = n.
	73: Memory allocation failed.

	successfully call of the nag_real_symm_general_eigenvalues function.


*/
	
	int  nag_real_symm_general_eigenvalues(
	int n, // the order of the matrices A and B.
	double a[], // Input: the upper triangle of the n by n symmetric matrix A. The elements of the array below the diagonal need not be set.  Output: the lower triangle of the array is overwritten. The rest of the array is unchanged.
	int tda, // the second dimension of the array a as declared in the function from which nag real symm general eigenvalues is called.
	double b[], // Input: the upper triangle of the n by n symmetric positive-de.nite matrix B. The elements of the array below the diagonal need not be set.  Output: the elements below the diagonal are overwritten. The rest of the array is unchanged.
	int tdb, // the second dimension of the array b as declared in the function from which nag real symm general eigenvalues is called.
	double r[] // the eigenvalues in ascending order.
	);


/**	f02aec
		calculates all the eigenvalues and eigenvectors of Ax = lamdaBx, where A 
		is a real symmetric matrix and B is a real symmetric positive-definite matrix.


Example:
	To calculate all the eigenvalues and eigenvectors of the general 
	symmetric eigenproblem Ax = lamdaBx where A is the symmetric matrix
	.
	..
	0.5 1.5 6.6 4.8
	1.5 6.5 16.2 8.6
	6.6 16.2 37.6 9.8
	4.8 8.6 9.8 -17.1
	.
	..
	and B is the symmetric positive-de.nite matrix
	.
	..
	1 3 4 1
	3 13 16 11
	4 16 24 18
	1 11 18 27
	.
	..
	.

void test_nag_real_symm_general_eigensystem()
{
	int i, j, n;
	matrix a, b, v;
	a.SetSize(8,8);
	b.SetSize(8,8);
	v.SetSize(8,8);
	double  r[8];
	n = 4;
	a[0][0] = 0.5;
	a[0][1] = 1.5;
	a[0][2] = 6.6; 
	a[0][3] = 4.8;
	
	
	a[1][0] = 1.5;
	a[1][1] = 6.5;
	a[1][2] = 16.2;
	a[1][3] = 8.6;
	
	a[2][0] = 6.6;
	a[2][1] = 16.2;
	a[2][2] = 37.6;
	a[2][3] = 9.8; 
	
	a[3][0] = 4.8;
	a[3][1] = 8.6;
	a[3][2] = 9.8; 
	a[3][3] = -17.1;
	
	b[0][0] = 1.0;
	b[0][1] = 3.0;
	b[0][2] = 4.0;
	b[0][3] = 1.0;	
	
	b[1][0] = 3.0;
	b[1][1] = 13.0; 
	b[1][2] = 16.0;
	b[1][3] = 11.0;
	
	b[2][0] = 4.0;
	b[2][1] = 16.0; 
	b[2][2] = 24.0;
	b[2][3] = 18.0;
  	
	b[3][0] = 1.0;
	b[3][1] = 11.0;
	b[3][2] = 18.0;
	b[3][3] = 27.0;
	nag_real_symm_general_eigensystem(n, a, 8, b, 8, r, v, 8);
	printf("Eigenvalues\n");
	for (i=0; i<n; i++)
	{
		printf("%9.4f",r[i]);
		if(i%8==7 || i==n-1)
			printf("\n");
	}
	printf("Eigenvectors\n");
	for (i=0; i<n; i++)
	{
		for (j=0; j<n; j++)
		{
			printf("%9.4f",v[i][j]);
			if(j%8==7 || j==n-1)
				printf("\n");
		}
	}
}

	The output is following:
	
	Eigenvalues
	
	-3.0000 -1.0000 2.0000 4.0000
	
	Eigenvectors
	
	-4.3500 -2.0500 -3.9500  2.6500
	 0.0500  0.1500  0.8500  0.0500
	 1.0000  0.5000  0.5000 -1.0000
	-0.5000 -0.5000 -0.5000  0.5000


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	358:  The matrix B is not positive-de.nite, possibly due to rounding errors.
	362: More than _value_ iterations are required to isolate all the eigenvalues.
	11: On entry, n must not be less than 1: n = _value_.  
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdb = _value_ while n = _value_. These parameters must satisfy tdb = n.  On entry, tdv = _value_ while n = _value_. These parameters must satisfy tdv = n.
	73: Memory allocation failed.

	successfully call of the nag_real_symm_general_eigensystem function.

*/

	int  nag_real_symm_general_eigensystem(
	int n, // the order of the matrices A and B.
	double a[],  // Input: the upper triangle of the n by n symmetric matrix A. The elements of the array below the diagonal need not be set.  Output: the lower triangle of the array is overwritten. The rest of the array is unchanged.
	int tda, // the second dimension of the array a as declared in the function from which nag real symm general eigensystem is called.
	double b[], // Input: the upper triangle of the n by n symmetric positive-de.nite matrix B. The elements of the array below the diagonal need not be set.  Output: the elements below the diagonal are overwritten. The rest of the array is unchanged.
	int tdb, // the second dimension of the array b as declared in the function from which nag real symm general eigensystem is called.
	double r[], // the eigenvalues in ascending order.
	double v[], // the normalised eigenvectors, stored by columns; the ith column corresponds to the ith eigenvalue. The eigenvectors x are normalised so that xTBx = 1. 
	int tdv // the second dimension of the array v as declared in the function from which nag real symm general eigensystem is called.

	);
/**	f02afc
		calculates all the eigenvalues of a real unsymmetric matrix.

Example:

	To calculate all the eigenvalues of the real matrix.
	.
	..
	1.5 0.1 4.5 -1.5
	-22.5 3.5 12.5 -2.5
	-2.5 0.3 4.5 -2.5
	-2.5 0.1 4.5 2.5
	.
	..
	.


void test_nag_real_eigenvalues()
{
	int i, j, n;
	matrix a;
	a.SetSize(4,4);
	matrix<complex> r;
	r.SetSize(1,4);
	int iter[4];
	n = 4;
	a[0][0] = 1.5;
	a[0][1] = 0.1;
	a[0][2] = 4.5;
	a[0][3] = -1.5;
		
	a[1][0] = -22.5;
	a[1][1] = 3.5;
	a[1][2] = 12.5;
	a[1][3] = -2.5;
		
	a[2][0] = -2.5;
	a[2][1] = 0.3;
	a[2][2] = 4.5;
	a[2][3] = -2.5;
		
	a[3][0] = -2.5;
	a[3][1] = 0.1; 
	a[3][2] = 4.5;
	a[3][3] = 2.5;
	nag_real_eigenvalues(n, a, 4, r, iter);
	printf("Eigenvalues\n");
	for (i=0; i<n; i++)
	{
		complex temp =  r[0][i];
		printf("( %7.3f , %7.3f ) \n", temp.m_re, temp.m_im);
	}
}

	The output is following:
	
	Eigenvalues
	
	( 3.000 ,  4.000 )
	( 3.000 , -4.000 )
	( 4.000 ,  0.000 )
	( 2.000 ,  0.000 )
	

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	362: More than _value_ iterations are required to isolate all the eigenvalues.
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry, tda = _value_ while n = value. These parameters must satisfy tda = n.
	73: Memory allocation failed.

	successfully call of the nag_real_eigenvalues function.

*/
int  nag_real_eigenvalues(
int n, // the order of the matrix A.
double a[], // Input: the n by n matrix A.  Output: the array is overwritten.
int tda, // the second dimension of the array a as declared in the function from which nag real eigenvalues is called.
complex r[], // the eigenvalues.
int iter[] // contains the number of iterations used to find the ith eigenvalue. If iter[i - 1] is negative, the ith eigenvalue is the second of a pair found simultaneously.
);
/**	f02agc
		calculates all the eigenvalues and eigenvectors of 
		a real unsymmetric matrix.

Example:
	To calculate all the eigenvalues and eigenvectors of the real matrix.
	.
	..
	1.5 0.1 4.5 -1.5
	-22.5 3.5 1 2 .5 -2.5
	-2.5 0.3 4.5 -2.5
	-2.5 0.1 4.5 2.5
	.
	..
	.

void test_nag_real_eigensystem()
{

	int i, j, n;
	matrix a;
	a.SetSize(4,4);
	matrix<complex>r, v;
	r.SetSize(1,4);
	v.SetSize(4,4);
	int iter[4];
  
	a[0][0] = 1.5;	
	a[0][1] = 0.1;
	a[0][2] = 4.5;
	a[0][3] = -1.5;
		
	a[1][0] = -22.5;
	a[1][1] = 3.5;
	a[1][2] = 12.5;
	a[1][3] = -2.5;
		
	a[2][0] = -2.5;
	a[2][1] = 0.3;
	a[2][2] = 4.5;
	a[2][3] = -2.5;
		
	a[3][0] = -2.5;
	a[3][1] = 0.1; 
	a[3][2] = 4.5;
	a[3][3] = 2.5;
	n = 4;
	nag_real_eigensystem(n, a, 4, r, v, 4, iter);
	printf("Eigenvalues\n");
	for (i=0; i<n; i++)
	{
		complex temp = r[0][i];
		printf("(%7.3f, %7.3f) \n", temp.m_re, temp.m_im);
	
	}
	printf("\nEigenvectors\n");
	for (i=0; i<n; i++)
		for (j=0; j<n; j++)
		{
			complex temp = v[i][j];
			printf("(%7.3f, %7.3f)", temp.m_re, temp.m_im);
			if(j%4==3 || j==n-1)
				printf("\n");
		}

}


	The output is following:
	
	Eigenvalues
	
	( 3.000,  4.000)
	( 3.000, -4.000)
	( 4.000,  0.000)
	( 2.000,  0.000)
	
	Eigenvectors
	
	( 0.113, -0.151) ( 0.113, 0.151) ( -0.033, 0.000) ( 0.063, 0.000)
	( 0.945,  0.000) ( 0.945, 0.000) (  0.988, 0.000) ( 0.996, 0.000)
	( 0.189,  0.000) ( 0.189, 0.000) (  0.011, 0.000) ( 0.006, 0.000)
	( 0.113, -0.151) ( 0.113, 0.151) (  0.154, 0.000) ( 0.063, 0.000)
		

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	362: More than _value_ iterations are required to isolate all the eigenvalues.
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdv = _value_ while n = _value_. These parameters must satisfy tdv = n.
	73: Memory allocation failed.

	successfully call of the nag_real_eigensystem function.


*/
	
	int  nag_real_eigensystem(
	int n, // the order of the matrix A.
	double a[], // Input: the n by n matrix A.  Output: the array is overwritten.
	int tda, // the second dimension of the array a as declared in the function from which nag real eigensystem is called.
	complex r[], // the eigenvalues.
	complex v[], // the eigenvectors, stored by columns. The ith column corresponds to the ith eigenvalue. The eigenvectors are normalised so that the sum of the squares of the moduli of the elements is equal to 1an d the element of largest modulus is real. This ensures that real eigenvalues have real eigenvectors.
	int tdv, // the second dimension of the array v as declared in the function from which nag real eigensystem is called.
	int iter[] // contains the number of iterations used to find the ith eigenvalue. If iter[i - 1] is negative, the ith eigenvalue is the second of a pair found simultaneously.
	);


/**	f02awc
		calculates all the eigenvalues of a complex Hermitian matrix.

Example:
	To calculate all the eigenvalues of the complex Hermitian matrix:
	.
	..
	0.50 0.00 1.84 + 1.38i 2.08 - 1.56i
	0.00 0.50 1.1 2 +0.84i -0.56 + 0.42i
	1.84 - 1.38i 1.12 - 0.84i 0.50 0.00
	2.08 + 1.56i -0.56 - 0.42i 0.00 0.50
	.
	..
	.

void test_nag_hermitian_eigenvalues()
{	
	int i, j, n;
	matrix<complex>a;
	a.SetSize(4,4);
	double r[4];
	n = 4;

	a[0][0] = 0.50 + 0.00i;
	a[0][1] =  0.00 + 0.00i;
	a[0][2] = 1.84 + 1.38i;
	a[0][3] = 2.08 - 1.56i;
	
	a[1][0] = 0.00 + 0.00i;
	a[1][1] = 0.50 + 0.00i;
	a[1][2] = 1.12 + 0.84i;
	a[1][3] = -0.56 + 0.42i;
	
	a[2][0] = 1.84 - 1.38i;
	a[2][1] = 1.12 - 0.84i;
	a[2][2] = 0.50 + 0.00i;
	a[2][3] = 0.00 + 0.00i;
	
	a[3][0] = 2.08 + 1.56i;
	a[3][1] = -0.56 - 0.42i;
	a[3][2] = 0.00 + 0.00i;
	a[3][3] = 0.50 + 0.00i;	

	nag_hermitian_eigenvalues(n, a, 4, r);
	printf("Eigenvalues\n");
	for (i=0; i<n; i++)
		printf("%9.4f", r[i]);
	printf("\n");
}

	The output is following:
	
	Eigenvalues
	
	-3.0000 -1.0000 2.0000 4.0000
	

Parameters:

Return:
	This function returns NAG error code, 0 if no error.

	11: On entry, n must not be less than 1: n = _value_.
	17: On entry tda = _value_ while n = _value_. These parameters must satisfy tda = n.
	73: Memory allocation failed.
	362: More than _value_ iterations are required to isolate all the eigenvalues.

	successfully call of the nag_hermitian_eigenvalues function.

*/

	int  nag_hermitian_eigenvalues(
	int n, // the order of the matrix A.
	complex a[], // Input: the elements of the lower triangle of the n by n complex Hermitian matrix A. Elements of the array above the diagonal need not be set.  Output: the array is overwritten.
	int tda, // the last dimension of the array a as declared in the function from which nag hermitian eigenvalues is called.
	double r[] // the eigenvalues in ascending order.
	);

/**	f02axc
		calculates all the eigenvalues and eigenvectors of 
		a complex Hermitian matrix.


Example:
	To calculate the eigenvalues and eigenvectors of the complex Hermitian matrix:
	.
	..
	0.50 0.00 1.84 + 1.38i 2.08 - 1.56i
	0.00 0.50 1.1 2 +0.84i -0.56 + 0.42i
	1.84 - 1.38i 1.12 - 0.84i 0.50 0.00
	2.08 + 1.56i -0.56 - 0.42i 0.00 0.50
	.
	..
	.

void test_nag_hermitian_eigensystem()
{	
	int i, j, n;
	matrix<complex> a, v;
	a.SetSize(4,4);
	v.SetSize(4,4);
	double r[4];
	n = 4;

	a[0][0] = 0.50 + 0.00i;
	a[0][1] =  0.00 + 0.00i;
	a[0][2] = 1.84 + 1.38i;
	a[0][3] = 2.08 - 1.56i;
	
	a[1][0] = 0.00 + 0.00i;
	a[1][1] = 0.50 + 0.00i;
	a[1][2] = 1.12 + 0.84i;
	a[1][3] = -0.56 + 0.42i;
	
	a[2][0] = 1.84 - 1.38i;
	a[2][1] = 1.12 - 0.84i;
	a[2][2] = 0.50 + 0.00i;
	a[2][3] = 0.00 + 0.00i;
	
	a[3][0] = 2.08 + 1.56i;
	a[3][1] = -0.56 - 0.42i;
	a[3][2] = 0.00 + 0.00i;
	a[3][3] = 0.50 + 0.00i;
	nag_hermitian_eigensystem(n, a, 4, r, v, 4);
	printf("Eigenvalues\n");
	for (i=0; i<n; i++)
		printf("%9.4f", r[i]);
	printf("\nEigenvectors\n");
	for (i=0; i<n; i++)
		for (j=0; j<n; j++)
		{
			complex temp = v[i][j];
			printf("(%7.3f %7.3f )",temp.m_re, temp.m_im);
			if(j%4==3 || j==n-1)
				printf("\n");
		}
}

	The output is following:
		
	Eigenvalues
	
	-3.0000 -1.0000 2.0000 4.0000
	
	Eigenvectors
	
	(  0.700  0.000 ) ( -0.100 0.000 ) ( -0.100  0.000 ) ( 0.700  0.000 )
	(  0.100  0.000 ) (  0.700 0.000 ) (  0.700  0.000 ) ( 0.100  0.000 )
	( -0.400  0.300 ) ( -0.400 0.300 ) (  0.400 -0.300 ) ( 0.400 -0.300 )
	( -0.400 -0.300 ) (  0.400 0.300 ) ( -0.400 -0.300 ) ( 0.400  0.300 )

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry tdv = _value_ while n = _value_. These parameters must satisfy tdv = n.
	73: Memory allocation failed.

	successfully call of the nag_hermitian_eigensystem function.

*/

	int  nag_hermitian_eigensystem(
	int n, // the order of the matrix A.
	complex a[], // the elements of the lower triangle of the n by n complex Hermitian matrix A. Elements of the array above the diagonal need not be set. 
	int tda, // the last dimension of the array a as declared in the function from which nag hermitian eigensystem is called.
	double r[], // the eigenvalues in ascending order.
	complex v[], // the eigenvectors, stored by columns. The ith column corresponds to the ith eigenvector. The eigenvectors are normalised so that the sum of the squares of the moduli of the elements is equal to 1an d the element of largest modulus is real.
	int tdv // the last dimension of the array v as declared in the function from which nag hermitian eigensystem is called.
	);

/**	f02bjc
		calculates all the eigenvalues and, if required, all the 
		eigenvectors of the generalized eigenproblem Ax = lamdaBx 
		where A and B are real, square matrices, using the QZ algorithm.

Example:
	To find all the eigenvalues and eigenvectors of Ax = ?Bx where
	A =
	.
	..
	3.9 4.3 4.3 4.4
	12.5 21.5 21.5 26.0
	-34.5 -47.5 -43.5 -46.0
	-0.5 7.5 3.5 6.0
	.
	..
	and
	B =
	.
	..
	1 1 1 1
	2 3 3 3
	-3 -5 -4 -4
	1 4 3 4
	.
	..
	.

void test_nag_real_general_eigensystem()
{
	int i, j, k, n, ip, iter[8];
	matrix<complex> alfa;
	alfa.SetSize(1,8);
	double beta[8], eps1;
	matrix a, b, z;
	a.SetSize(8,8);
	b.SetSize(8,8);
	z.SetSize(8,8);
	BOOL matz;	
	n = 4;
	
	a[0][0] = 3.9;
	a[0][1] = 12.5;
	a[0][2] = -34.5;
	a[0][3] = -0.5;
	
	a[1][0] = 4.3;  
	a[1][1] = 21.5;
	a[1][2] = -47.5;
	a[1][3] = 7.5;
	
	a[2][0] = 4.3;
	a[2][1] = 21.5;
	a[2][2] = -43.5;
	a[2][3] = 3.5;
	
	a[3][0] = 4.4;
	a[3][1] = 26.0;
	a[3][2] = -46.0;
	a[3][3] = 6.0;
	
	b[0][0] = 1.0;
	b[0][1] = 2.0; 
	b[0][2] = -3.0; 
	b[0][3] = 1.0;
	
	b[1][0] = 1.0; 
	b[1][1] = 3.0; 
	b[1][2] = -5.0;
	b[1][3] = 4.0;
	
	b[2][0] = 1.0;
	b[2][1] = 3.0;
	b[2][2] = -4.0;
	b[2][3] = 3.0;
	
	b[3][0] = 1.0;
	b[3][1] = 3.0;
	b[3][2] = -4.0;
	b[3][3] = 4.0;
	  	
	matz = TRUE;
	eps1 = 1e-8;
	nag_real_general_eigensystem(n, a, 8, b, 8, eps1, alfa, beta, matz, z, 8, iter);
	ip = 0;
	for (i=0; i<n; ++i)
	{
		complex tmp = alfa[0][i];
		printf("Eigensolution %4ld\n",i+1);
		printf("alfa[%ld].re %7.3f",i,tmp.m_re);
		printf(" alfa[%ld].im %7.3f",i,tmp.m_im);
		printf(" beta[%ld] %7.3f\n",i,beta[i]);
		if (beta[i] == 0.0)
			printf("lambda is infinite");
		else if (tmp.m_im == 0.0)
		{
			printf("lambda %7.3f\n",tmp.m_re/beta[i]);
			printf("Eigenvector\n");
			for (j=0; j<n; ++j)
				printf("%7.3f\n", z[j][i]);
		}
		else
		{
			complex temp = alfa[0][i];
			printf("lambda %7.3f %7.3f\n",temp.m_re/beta[i], temp.m_im/beta[i] );
			printf("Eigenvector\n");
			k = pow(-1, ip+2);
			for (j=0; j<n; ++j)
			{
				printf("%7.3f",z[j][i-ip]);
				printf("%7.3f\n",k*z[j][i-ip+1]);
			}
			ip = 1-ip;
		}
	}
	printf("Number of iterations (machine-dependent)\n");
	for (i=0; i<n; ++i)
		printf("%2ld",iter[i]);
	printf("\n");	
}

	The output is following:
	
	Eigensolution 1
	
	alfa[0].re 3.801 alfa[0].im 0.000 beta[0] 1.900
	
	lambda 2.000
	
	Eigenvector
	
	0.996
	0.006
	0.063
	0.063
	
	Eigensolution 2
	
	alfa[1].re 1.563 alfa[1].im 2.084 beta[1] 0.521
	
	lambda 3.000 4.000
	
	Eigenvector
	
	0.945 0.000
	0.189 0.000
	0.113 -0.151
	0.113 -0.151
	
	Eigensolution 3
	
	alfa[2].re 3.030 alfa[2].im -4.040 beta[2] 1.010
	
	lambda 3.000 -4.000
	
	Eigenvector
	
	0.945 0.000
	0.189 0.000
	0.113 0.151
	0.113 0.151
	
	Eigensolution 4
	
	alfa[3].re 4.000 alfa[3].im 0.000 beta[3] 1.000
	
	lambda 4.000
	
	Eigenvector
	
	0.988
	0.011
	-0.033
	0.154
	
	Number of iterations (machine-dependent)
	0 0 5 0


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry tdb = _value_ while n = _value_. These parameters must satisfy tdb = n.  On entry tdv = _value_ while n = _value_. These parameters must satisfy tdv = n.
	374: More than n × 30 iterations are required to determine all the diagonal 1 by 1 or 2 by 2 blocks of the quasi-triangular form in the second step of the QZ algorithm. This failure occurs at the ith eigenvalue, i = _value_. aj and ßj are correct for j = i+1, i+2, . . ., n but v does not contain any correct eigenvectors.  The value of i will be returned in member errnum of the NAG error structure provided NAGERR DEFAULT is not used as the error parameter.
	
	successfully call of the nag_real_general_eigensystem function.

*/

	int  nag_real_general_eigensystem(
	int n, // the order of the matrices A and B.
	double a[], // Input: the n by n matrix A.  Output: the array is overwritten.
	int	tda, // the second dimension of the array a as declared in the function from which nag real general eigensystem is called.
	double b[], // Input: the n by n matrix B.  Output: the array is overwritten.
	int tdb, // the second dimension of the array b as declared in the function from which nag real general eigensystem is called.
	double tol, // the tolerance used to determine negligible elements. If tol > 0.0, an element will be considered negligible if it is less than tol times the norm of its matrix. If tol = 0.0, machine precision is used in place of tol. A value of tol greater than machine precision may result in faster execution but less accurate results.
	complex alfa[], // Output: alpha_j, for j = 1, 2, . . . , n.
	double beta[], // Output: beta_j, for j = 1, 2, . . . , n.
	bool wantv, // wantv must be set to TRUE if the eigenvectors are required. If wantv is set to FALSE then the array v is not referenced.
	double v[], // Output: if wantv = TRUE, then (i) if the jth eigenvalue is real, the jth column of v contains its eigenvector; (ii) if the jth and (j + 1)th eigenvalues form a complex pair, the jth and (j + 1)th columns of v contain the real and imaginary parts of the eigenvector associated with the first eigenvalue of the pair. The conjugate of this vector is the eigenvector for the conjugate eigenvalue.  Each eigenvector is normalised so that the component of largest modulus is real and the sum of squares of the moduli equal one.  If wantv = FALSE, v is not referenced and may be set to the null pointer
	int tdv, // the second dimension of the array v as declared in the function from which nag real general eigensystem is called.
	int iter[] // contains the number of iterations needed to obtain the jth eigenvalue.
	);

/**	f02ecc
		computes selected eigenvalues and eigenvectors of 
		a real general matrix.


Example:
	To compute those eigenvalues of the matrix A whose moduli lie in 
	the range [0.2,0.5], and their corresponding eigenvectors, where
	A =
	.
	..
	0.35 0.45 -0.14 -0.17
	0.09 0.07 -0.54 0.35
	-0.44 -0.33 -0.03 0.17
	0.25 -0.32 -0.13 0.11
	.


void test_nag_real_eigensystem_sel()
{
	matrix a;
	a.SetSize(8,8);
	double wl, wu;
	int i, j, m, n;
	int mest = 3;
	
	matrix<complex>w, v;
	w.SetSize(1,8);
	v.SetSize(8,3);
	n = 4;
	wl = 0.2;
	wu = 0.5;
	a[0][0] =0.35;
	a[0][1] =0.45;
	a[0][2] =-0.14;
	a[0][3] = -0.17;
		
	a[1][0] =0.09;
	a[1][1] =0.07;
	a[1][2] = -0.54;
	a[1][3] =0.35;
	
	a[2][0] =-0.44;
	a[2][1] =-0.33;
	a[2][2] =-0.03;
	a[2][3] =0.17;
	
	a[3][0] =0.25;
	a[3][1] =-0.32;
	a[3][2] =-0.13;
	a[3][3] =0.11;
	
	nag_real_eigensystem_sel(Nag_Select_Modulus, n, a, 8, wl, wu, mest, &m, w,v, 3);
	printf("\n\nEigenvalues\n");
	for (i = 0; i < m; ++i)
	{
		complex temp;
		temp = w[0][i];
		printf("(%7.4f, %7.4f) ", temp.m_re, temp.m_im);
	}
	printf("\n");
	printf(" Eigenvectors\n ");
	for (i = 1; i <= m; i++)
	{
		printf("%15ld", i);
		if(i%m == 0 )
			printf("\n");
	}
	for (i = 0; i < n; i++)
	{
		printf("%ld ", i + 1);
		for (j = 0; j < m; j++)
		{
			complex temp;
			temp = v[i][j];
			printf("(%8.4f, %8.4f)", temp.m_re, temp.m_im);
			if((j + 1) % m == 0)
				printf("\n");
		}
	}
}

	The output is following:

	Eigenvalues

	(-0.0994, 0.4008) (-0.0994, -0.4008)

	Eigenvectors

			1 					2
	
	1 ( -0.1933,  0.2546) ( -0.1933, -0.2546)
	2 (  0.2519, -0.5224) (  0.2519,  0.5224)
	3 (  0.0972, -0.3084) (  0.0972,  0.3084)
    4 (  0.6760,  0.0000) (  0.6760,  0.0000)
 

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	70: On entry, parameter crit had an illegal value.
	11: On entry, n must not be less than 0: n = _value_.  On entry, mest must not be less than 1: mest = _value_.
	91: On entry, tda = _value_ while n = _value_.  Constraint: tda = max(1,n).  On entry, tdv = _value_ while mest = _value_.  Constraint: tdv = max(1,mest).
	24: On entry, wu = _value_ while wl = _value_.  These parameters must satisfy wu > wl.
	378: The QR algorithm failed to compute all the eigenvalues. No eigenvectors have been computed.
	379: There are more than mest eigenvalues in the speci.ed range. The actual number of eigenvalues in the range is returned in m. No eigenvectors have been computed.  Rerun with the second dimension of v = mest = m.
	
	successfully call of the nag_real_eigensystem_sel function.

*/
	
	int  nag_real_eigensystem_sel(
	Nag_Select_Eigenvalues crit,
	int n, // the order of the matrix A.
	double a[], // Input: the n by n general matrix A.  Output: a contains the Hessenberg form of the balanced input matrix A_
	int tda,  // the last dimension of the array a as declared in the calling program.
	double wl, // the lower bound on the criterion for the selected eigenvalues.
	double wu, // the upper bound on the criterion for the selected eigenvalues.
	int mest, // an upper bound on m, the number of eigenvalues and eigenvectors selected. No eigenvectors are computed if mest < m.
	int *m, // the number of eigenvalues actually selected
	complex w[], // the first m elements of w hold the values of the selected eigenvalues; elements from the index m to n-1 contain the other eigenvalues. Complex conjugate pairs of eigenvalues are stored in consecutive elements of the array, with the eigenvalue having the positive imaginary part first.
	complex v[], // the selected eigenvectors, with the ith column holding the eigenvector associated with the eigenvalue lamda.
	int tdv // the second dimension of the array v as declared in the calling program.
	);

/**	f02gcc
		computes selected eigenvalues and eigenvectors of 
		a complex general matrix.


Example:
	To compute those eigenvalues of the matrix A whose moduli lie 
	in the range [-5.5,+5.5], and their corresponding eigenvectors, where
	A =
	.
	..
	-3.97 - 5.04i -4.11 + 3.70i -0.34 + 1.01i 1.29 - 0.86i
	0.34 - 1.50i 1.52 - 0.43i 1.88 - 5.38i 3.36 + 0.65i
	3.31 - 3.85i 2.50 + 3.45i 0.88 - 1.08i 0.64 - 1.48i
	-1.10 + 0.82i 1.81 - 1.59i 3.25 + 1.33i 1.57 - 3.44i
	.
	..
	
void test_nag_complex_eigensystem_sel()
{	
	double wl,wu;
	int mest = 3;
	int i,j,m, n;
	matrix<complex> a, v, w;
	a.SetSize(8,8);
	v.SetSize(8,3);
	w.SetSize(1,8);	
	n = 4;
	wl = -5.5;
	wu = 5.5;		
	
	a[0][0] = -3.97 - 5.04i;
	a[0][1] = -4.11 + 3.70i;
	a[0][2] = -0.34 + 1.01i;
	a[0][3] = 1.29 - 0.86i;
		
	a[1][0] = 0.34 - 1.50i;
	a[1][1] = 1.52 - 0.43i;
	a[1][2] = 1.88 - 5.38i;
	a[1][3] = 3.36 + 0.65i;
	
	a[2][0] = 3.31 - 3.85i;
	a[2][1] = 2.50 + 3.45i;
	a[2][2] = 0.88 - 1.08i;
	a[2][3] = 0.64 - 1.48i;
		
	a[3][0] = -1.10 + 0.82i;
	a[3][1] = 1.81 - 1.59i;
	a[3][2] = 3.25 + 1.33i;
	a[3][3] = 1.57 - 3.44i;
		
	nag_complex_eigensystem_sel(Nag_Select_Modulus, n, a, 8, wl, wu, mest, &m, w, v,3);
	printf("\n\nEigenvalues\n\n");
	for (i = 0; i < m; ++i)
	{
		complex temp = w[0][i];
		printf("(%7.4f,%7.4f)\t", temp.m_re, temp.m_im);
		if((i+1)%2 == 0)
			printf("\n");
	}
	printf("\nEigenvectors\n");
	for (i=1; i<=m; i++)
	{
		printf("%15ld",i);
		if(i%m == 0)
			printf("\n");
	}
	for (i = 0; i < n; i++)
	{
		printf("% ld ",i + 1);
		for (j = 0; j < m; j++)
		{
			complex temp = v[i][j];
			printf("(%8.4f,%8.4f)", temp.m_re, temp.m_im);
			if((j + 1)%m == 0)
				printf("\n");
		}		
	}
}
	

	The output is following:
	
	Eigenvalues
	
	(-5.0000,2.0060) ( 3.0023, -3.9998)
	
	Eigenvectors
	
				1 			2
	
	1 ( -0.3865, 0.1732) ( -0.0356, -0.1782)
	2 ( -0.3539, 0.4529) (  0.1264,  0.2666)
	3 (  0.6124, 0.0000) (  0.0129, -0.2966)
	4 ( -0.0859,-0.3284) (  0.8898,  0.0000)
		

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	70: On entry, parameter crit had an illegal value.
	11: On entry, n must not be less than 0: n = _value_.  On entry, mest must not be less than 1: mest = _value_.
	91: On entry, tda = _value_ while n = _value_. Constraint: tda = max (1,n).
	24: On entry, wu = _value_ while wl = _value_.  These parameters must satisfy wu > wl.
	17: On entry, tdv = _value_ while mest = _value_.  These parameters must satisfy tdv = mest.
	378: The QR algorithm failed to compute all the eigenvalues. No eigenvectors have been computed.
	379: There are more than mest eigenvalues in the speci.ed range. The actual number of eigenvalues in the range is returned in m. No eigenvectors have been computed.  Rerun with the second dimension of v = mest = m.
	380: Inverse iteration failed to compute all the specified eigenvectors. If an eigenvector failed to converge, the corresponding column of v is set to zero.

	successfully call of the nag_complex_eigensystem_sel function.

	
*/
	
	int  nag_complex_eigensystem_sel(
	Nag_Select_Eigenvalues crit,
	int n, // the order of the matrix A.
	complex a[], // Input: the n by n general matrix A.  Output: a contains the Hessenberg form of the balanced input matrix A
	int tda,  // the last dimension of the array a as declared in the calling program.
	double wl, // the lower bound on the criterion for the selected eigenvalues.
	double wu, // the upper bound on the criterion for the selected eigenvalues.
	int mest, // an upper bound on m, the number of eigenvalues and eigenvectors selected. No eigenvectors are computed if mest < m.
	int *m, // the number of eigenvalues actually selected.
	complex w[], // the first m elements of w hold the selected eigenvalues; elements from the index m to n-1 contain the other eigenvalues.
	complex v[], // the selected eigenvectors, with the ith column holding the eigenvector associated with the eigenvalue lamda.
	int tdv // the second dimension of the array v as declared in the calling program.
	);
	
/**	f02wec
			returns all, or part, of the singular value 
			decomposition of a general real matrix.

Example:
	To find the singular value decomposition of the 5 by 3 matrix
	A =
	.
	....
	2.0 2.5 2.5
	2.0 2.5 2.5
	1.6 -0.4 2.8
	2.0 -0.5 0.5
	1.2 -0.3 -2.9
	.
	....
	together with the vector QT b for the vector
	b =
	.
	....
	1.1
	0.9
	0.6
	0.0
	-0.8
	.
	....
	.

void test_nag_real_svd()
{
	int tda = 10;
	int tdpt = 10;
	matrix a;
	double b[20],e[9];
	a.SetSize(20,10);
	
	matrix pt;
	pt.SetSize(10,10);
	double sv[10], dummy[1];
	int i, j, m, n, iter, failinfo;
	BOOL wantp, wantq;
	m = 5;
	n = 3;
	a[0][0] = 2.0;
	a[0][1] = 2.5;
	a[0][2] = 2.5;
	
	a[1][0] = 2.0;
	a[1][1] = 2.5;
	a[1][2] = 2.5;
	
	a[2][0] = 1.6;
	a[2][1] = -0.4;
	a[2][2] = 2.8;
	
	a[3][0] = 2.0;
	a[3][1] = -0.5;
	a[3][2] = 0.5;
	
	a[4][0] = 1.2;
	a[4][1] = -0.3;
	a[4][2] = -2.9;

	b[0] = 1.1;
	b[1] = 0.9;
	b[2] = 0.6;
	b[3] = 0.0;
	b[4] = -0.8;

	wantq = TRUE;
	wantp = TRUE;
	nag_real_svd(m, n, a, tda, 1, b, 1, wantq, dummy, 1, sv, wantp, pt, tdpt, &iter, e, &failinfo);
	
	printf("Singular value decomposition of A\n\n");
	printf("Singular values\n");
	for (i = 0; i < n; ++i)
		printf(" %8.4f", sv[i]);
	printf("\n\n");
	printf("Left-hand singular vectors, by column\n");
	for (i = 0; i < m; ++i)
	{
		for (j = 0; j < n; ++j)
			printf(" %8.4f", a[i][j]);
		printf("\n");
	}
	printf("\n");
	printf("Right-hand singular vectors, by column\n");
	for (i = 0; i < n; ++i)
	{
		for (j = 0; j < n; ++j)
			printf(" %8.4f", pt[j][i]);
		printf("\n");
	}
	printf("\n");
	printf("Vector Q'*B\n");
	for (i = 0; i < m; ++i)
		printf(" %8.4f", b[i]);
	printf("\n\n");
}

	The output is following:
		
	Example 2
	Singular value decomposition of A
	
	Singular values
    6.5616   3.0000   2.4384

	Left-hand singular vectors, by column
	   0.6011  -0.1961  -0.3165
	   0.6011  -0.1961  -0.3165
	   0.4166   0.1569   0.6941
	   0.1688  -0.3922   0.5636
	  -0.2742  -0.8629   0.0139
	
	Right-hand singular vectors, by column
	   0.4694  -0.7845   0.4054
	   0.4324  -0.1961  -0.8801
	   0.7699   0.5883   0.2471
	
	Vector Q'*B
	   1.6716   0.3922  -0.2276  -0.1000  -0.1000


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, m must not be less than 0: m = _value_.  On entry, n must not be less than 0: n = _value_.  On entry, ncolb must not be less than 0: ncolb = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n. On entry, tdb = _value_ while ncolb = _value_. These parameters must satisfy tdb = ncolb.
	363: On entry, tdq = _value_ while m = _value_. Whenwantq is TRUE and m < n then relationship tdq = m must be satis.ed.
	364: On entry, tdpt = _value_ while n = _value_. When wantq and wantp are TRUE and m = n then relationship tdpt = n must be satis.ed.
	366: The QR algorithm has failed to converge in _value_ iterations. Singular values 1,2,. . .,failinfo may not have been found correctly and the remaining singular values may not be the smallest.  The matrix A will nevertheless have been factorized as A = QEPT , where the leading min(m, n) by min(m, n) part of E is a bidiagonal matrix with sv[0], sv[1], . . ., sv[min(m,n-1)] as the diagonal elements and e[0], e[1], . . ., e[min(m,n-2)] as the superdiagonal elements.  This failure is not likely to occur.
	73: Memory allocation failed.

	successfully call of the nag_real_svd function.

*/
	
	int  nag_real_svd(
	int m, // the number of rows, m, of the matrix A.
	int n, // the number of columns, n, of the matrix A.
	double a[], // Input: the leading m by n part of the array a must contain the matrix A whose singular value decomposition is required.  Output: if m = n and wantq = TRUE, then the leading m by n part of a will contain the first n columns of the orthogonal matrix Q.  If m < n and wantp = TRUE, then the leading m by n part of a will contain the .rst m rows of the orthogonal matrix PT .  If m = n and wantq = FALSE and wantp = TRUE, then the leading n by n part of a will contain the first n rows of the orthogonal matrix PT .  Otherwise the contents of the leading m by n part of a are indeterminate.
	int tda, // the second dimension of the array a as declared in the function from which nag real svd is called.
	int ncolb, // the number of columns of the matrix B. When ncolb = 0 the array b is not referenced.
	double b[], // Input: if ncolb > 0, the leading m by ncolb part of the array b must contain the matrix to be transformed. If ncolb = 0 the array b is not referenced and may be set to the null pointer, i.e., (double *)0.  Output: b is overwritten by the m by ncolb matrix QTB.
	int tdb, // the second dimension of the array b as declared in the function from which nag real svd is called.
	bool wantq, // wantq must be TRUE, if the left-hand singular vectors are required. If wantq = FALSE, then the array q is not referenced.
	double q[], // if m < n and wantq = TRUE, the leading m by m part of the array q will contain the orthogonal matrix Q. Otherwise the array q is not referenced and may be set to the null pointer, i.e., (double *)0.
	int tdq, // the second dimension of the array q as declared in the function from which nag real svd is called.
	double sv[], // the min(m,n) diagonal elements of the matrix S.
	bool wantp, // wantp must be TRUE if the right-hand singular vectors are required. If wantp = FALSE, then the array pt is not referenced.
	double pt[], // if m = n and wantq and wantp are TRUE, the leading n by n part of the array pt will contain the orthogonal matrix PT . Otherwise the array pt is not referenced and may be set to the null pointer, i.e., (double *)0.
	int tdpt, // the second dimension of the array pt as declared in the function from which nag real svd is called.
	int *iter, // the total number of iterations taken by the QR algorithm.
	double e[], // if the error NE QR NOT CONV occurs the array e contains the super diagonal elements of matrix E in the factorisation of A according to A = QEPT . See Section 5 for further details.
	int *failinfo // if the error NE QR NOT CONV occurs failinfo contains the number of singular values which may not have been found correctly. See Section 5 for details.
	);

/**	f02xec
		returns all, or part, of the singular value 
		decomposition of a general complex matrix.


Example:
	To find the singular value decomposition of the 5 by 3 matrix
	A =
	.
	....
	0.5i -0.5 + 1.5i -1.0+ 1.0i
	0.4 + 0.3i 0.9 + 1.3i 0.2+ 1.4i
	0.4 -0.4 + 0.4i 1.8
	0.3 - 0.4i 0.1 + 0.7i 0.0
	-0.3i 0.3 + 0.3i 2.4i
	.
	....
	together with the vector QHb for the vector
	b =
	.
	....
	-0.55 + 1.05i
	0.49 + 0.93i
	0.56 - 0.16i
	0.39 + 0.23i
	1.13 + 0.83i
	.
	....
	.

void test_nag_complex_svd()
{
	int tda = 10;
	int tdph = 10;
	matrix<complex>a, b, ph, dummy;
	a.SetSize(20,10);
	b.SetSize(1,20);
	ph.SetSize(10,10);
	dummy.SetSize(1,1);
	double e[9], sv[10];
	int i, j, m, n, iter, failinfo;
	BOOL wantp, wantq;
	m = 5;
	n = 3;
	
	a[0][0] =0.00 + 0.50i;
	a[0][1] = -0.50 + 1.50i;
	a[0][2] = -1.00 + 1.00i;
		
	a[1][0] =0.40 + 0.30i;
	a[1][1] = 0.90 + 1.30i;
	a[1][2] = 0.20 + 1.40i;
	
	a[2][0] =0.40 + 0.00i;
	a[2][1] = -0.40 + 0.40i;
	a[2][2] = 1.80 + 0.00i;
	
	a[3][0] =0.30 - 0.40i;
	a[3][1] = 0.10 + 0.70i;
	a[3][2] = 0.00 + 0.00i;
	
	a[4][0] = 0.00 - 0.30i;
	a[4][1] = 0.30 + 0.30i;
	a[4][2] = 0.00 + 2.40i;
	
	b[0][0] = -0.55 + 1.05i;
	b[0][1] = 0.49 + 0.93i;
	b[0][2] = 0.56 - 0.16i;
	b[0][3] = 0.39 + 0.23i;
	b[0][4] = 1.13 + 0.83i;
	
	wantq = TRUE;
	wantp = TRUE;
	nag_complex_svd(m, n, a, tda, 1, b, 1, wantq, dummy, 1, sv, wantp, ph, tdph, &iter, e, &failinfo);
	printf("Singular value decomposition of A\n\nSingular values\n");
	for (i=0; i<n; ++i)
	{	
		printf("%9.4f", sv[i]);
		if(i%5==4 || i==n-1)
			printf("\n");
	}
	printf("\nLeft-hand singular vectors, by column\n");
	for (i=0; i<m; ++i)
		for (j=0; j<n; ++j)
		{
			complex temp = a[i][j];
			printf("(%7.4f %7.4f)  ", temp.m_re, temp.m_im);
			if(j%3==2 || j==n-1)
				printf("\n")
		}
		
	printf("\nRight-hand singular vectors, by column\n");
	for (i=0; i<n; ++i)
		for (j=0; j<n; ++j)
		{
			complex temp = ph[i][j];
			printf("(%7.4f %7.4f)  ", temp.m_re, temp.m_im);
			if(j%3==2 || j==n-1)
				printf("\n");
		}
		
	printf("\nVector conjg(Q')*B\n");
	for (i=0; i<m; ++i)
	{
		complex temp = b[0][i];	
		printf("(%7.4f %7.4f)  ", temp.m_re, temp.m_im);
		if(i%3==2 || i==m-1)
			printf("\n");			
	}
}

	The output is following:
		
	Example 2
	Singular value decomposition of A

	Singular values
   	   3.9263   2.0000   0.7641

	Left-hand singular vectors, by column
	(-0.0757 -0.5079)  (-0.2831 -0.2831)  (-0.2251  0.1594)  
	(-0.4517 -0.2441)  (-0.3963  0.0566)  (-0.0075  0.2757)  
	(-0.2366  0.2669)  (-0.1359 -0.6341)  ( 0.2983 -0.2082)  
	(-0.0561 -0.0513)  (-0.3284 -0.0340)  ( 0.1670 -0.5978)  
	(-0.4820 -0.3277)  ( 0.3737  0.1019)  (-0.0976 -0.5664)  

	Right-hand singular vectors, by column
	(-0.1275  0.0000)  (-0.3899 -0.2046)  (-0.5289 -0.7142)  
	(-0.2265  0.0000)  (-0.3397 -0.7926)  ( 0.0000  0.4529)  
	( 0.9656  0.0000)  (-0.1311 -0.2129)  (-0.0698  0.0119)  

	Vector conjg(Q')*B
	(-1.9656 -0.7935)  ( 0.1132 -0.3397)  ( 0.0915  0.6086)  
	(-0.0600 -0.0200)  ( 0.0400  0.1200)  



Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, m must not be less than 0: m = _value_.  On entry, n must not be less than 0: n = _value_.  On entry, ncolb must not be less than 0: ncolb = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdb = _value_ while ncolb = _value_. These parameters must satisfy tdb = ncolb.
	363: On entry, tdq = _value_ while m = _value_. Whenwantq is TRUE and m < n then relationship tdq = m must be satis.ed.
	364: On entry, tdph = _value_ while n = _value_. When wantq and wantp are TRUE and m = n then relationship tdph = n must be satis.ed.
	366: The QR algorithm has failed to converge in _value_ iterations. Singular values 1, 2, . . . ,fail info may not have been found correctly and the remaining singular values may not be the smallest.  The matrix A will nevertheless have been factorized as A = QEPT , where the leading min(m, n) by min(m, n) part of E is a bidiagonal matrix with sv[0], sv[1], . . ., sv[min(m,n-1)] as the diagonal elements and e[0], e[1], . . ., e[min(m,n-2)] as the super-diagonal elements.  This failure is not likely to occur.  
	73: Memory allocation failed.

	successfully call of the nag_complex_svd function.

*/
	
	int  nag_complex_svd(
	int m, // the number of rows, m, of the matrix A.
	int n, // the number of columns, n, of the matrix A.
	complex a[], // Input: the leading m by n part of the array a must contain the matrix A whose singular value decomposition is required.  Output: if m = n and wantq = TRUE, then the leading m by n part of a will contain the first n columns of the unitary matrix Q.  If m < n and wantp = TRUE, then the leading m by n part of a will contain the .rst m rows of the unitary matrix PH.  If m = n and wantq = FALSE and wantp = TRUE, then the leading n by n part of a will contain the .rst n rows of the unitary matrix PH.  Otherwise the contents of the leading m by n part of a are indeterminate.
	int tda, // the second dimension of the array a as declared in the function from which nag complex svd is called.
	int ncolb, // the number of columns of the matrix B. When ncolb = 0 the array b is not referenced.
	complex b[], // Input: if ncolb > 0, the leading m by ncolb part of the array b must contain the matrix to be transformed. If ncolb = 0 the array b is not referenced and may be set to the null pointer, i.e., (Complex *)0.  Output: b is overwritten by the m by ncolb matrix QHB.
	int tdb, // the second dimension of the array b as declared in the function from which nag complex svd is called.
	bool wantq, // Input: wantq must be TRUE if the left-hand singular vectors are required. If wantq = FALSE then the array q is not referenced.  
	complex q[], // Output: if m < n and wantq = TRUE, the leading m by m part of the array q will contain the unitary matrix Q. Otherwise the array q is not referenced and may be set to the null pointer, i.e., (Complex *)0.
	int tdq, // the second dimension of the array q as declared in the function from which nag complex svd is called.
	double sv[], // the min(m,n) diagonal elements of the matrix S.
	bool wantp, // wantp must be TRUE if the right-hand singular vectors are required. If wantp = FALSE then the array ph is not referenced
	complex ph[], // if m = n and wantq and wantp are TRUE, the leading n by n part of the array ph will contain the unitary matrix PH. Otherwise the array ph is not referenced and may be set to the null pointer, i.e., (Complex *)0.
	int tdph, // the second dimension of the array ph as declared in the function from which nag complex svd is called.
	int *iter, // the total number of iterations taken by the QR algorithm.
	double e[], // if the error NE QR NOT CONV occurs the array e contains the super-diagonal elements of matrix E in the factorisation of A according to A = QEPH. See Section 5 for further details.
	int *failinfo // if the error NE QR NOT CONV occurs failinfo contains the number of singular values which may not have been found correctly. See Section 5 for details.
	);

/** f03aec
*/

int  nag_real_cholesky(
int n, 
double a[], 
int tda, 
double p[],
double *detf, 
int *dete
);

/** f03afc
*/
int  nag_real_lu(
int n, 
double a[], 
int tda, 
int pivot[],
double *detf, 
int *dete
);

/**	f03ahc
*/
int  nag_complex_lu(
int n, 
complex a[], 
int tda, 
int pivot[],
complex *detf, 
int *dete
);

/** f04adc
*/
int  nag_complex_lin_eqn_mult_rhs(
int n, 
int nrhs, 
complex a[],
int tda, 
complex b[], 
int tdb, 
complex x[], 
int tdx
);
/** f04agc
*/

int  nag_real_cholesky_solve_mult_rhs(
int n, 
int nrhs, 
double a[],
int tda, 
double p[], 
double b[], 
int tdb, 
double x[], 
int tdx
);


/** f04ajc
*/
int  nag_real_lu_solve_mult_rhs(
int n, 
int nrhs,
double a[],
int tda, 
int pivot[], 
double b[], 
int tdb
);


/** f04akc
*/
int  nag_complex_lu_solve_mult_rhs(
int n, 
int nrhs,
complex a[],
int tda, 
int pivot[], 
complex b[], 
int tdb
);


/** f04arc
*/
int  nag_real_lin_eqn(
int n, 
double a[], 
int tda, 
double b[],
double x[]
);
/** f04awc
*/
int  nag_hermitian_lin_eqn_mult_rhs(
int n, 
int nrhs, 
complex a[],
int tda, 
double p[], 
complex b[], 
int tdb, 
complex x[],
int tdx
);

/** f04mcc
*/
int  nag_real_cholesky_skyline_solve(
Nag_SolveSystem selct, 
int n,
int nrhs, 
double al[], 
int lal, 
double d[], 
int row[],
double b[], 
int tdb, 
double x[], 
int tdx
);

/** f03aec
		computes a Cholesky factorization of a real symmetric 
		positive-definite matrix, and evaluates the determinant.

Example:
	To compute a Cholesky factorization and calculate the determinant of the 
	real symmetric positive-definite matrix
	.
	..
	6 7 6 5
	7 11 8 7
	6 8 11 9
	5 7 9 11
	.
	..
	.

void test_nag_real_cholesky()
{	
	double detf, determ;
	matrix a;
	a.SetSize(8,8);
	double p[8];
	int i, dete, j, n;
	n = 4;
	a[0][0] = 6;
	a[0][1] = 7;
	a[0][2] = 6;
	a[0][3] = 5;
	
	a[1][0] = 7;
	a[1][1] = 11;
	a[1][2] = 8;
	a[1][3] = 7;
	
	a[2][0] = 6;
	a[2][1] = 8;
	a[2][2] = 11;
	a[2][3] = 9;
	
	a[3][0] = 5;
	a[3][1] = 7;
	a[3][2] = 9;
	a[3][3] = 11;
	nag_real_cholesky(n, a, 8, p, &detf, &dete);
	printf("Array A after factorization\n");
	for (i=0; i<n; i++)
		for (j=0; j<n; j++)
		{
			printf("%9.4f", a[i][j]);
			if(j%8==7 || j==n-1)
				printf("\n");
		}
	printf("\nArray p\n");
	for (i=0; i<n; i++)
	{
		printf("%9.4f", p[i]);
		if(i%8==7 || i==n-1)
			printf("\n");
	}
	printf("\ndetf = %9.4f dete = %2ld\n\n", detf, dete);
	determ = detf*pow(2.0,1.0*dete);
	printf("Value of determinant = %9.4f\n", determ);

}

	The output is following:
	
	Array A after factorization
	
	6.0000  7.0000  6.0000  5.0000
	2.8577 11.0000  8.0000  7.0000
	2.4495  0.5941 11.0000  9.0000
	2.0412  0.6931  1.6645 11.0000
	
	Array p
	
	0.4082 0.5941 0.4639 0.5283
	
	detf = 0.0691 dete = 12
	
	Value of determinant = 283.0000


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	358: The matrix is not positive-de.nite, possibly due to rounding errors. The factorization could not be completed. detf and dete are set to zero.
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.
	
	successfully call of the nag_real_cholesky function.


*/

	int  nag_real_cholesky(
	int n, // the order of the matrix A.
	double a[], // Input: the upper triangle of the n by n positive-de.nite symmetric matrix A. The elements of the array below the diagonal need not be set.  Output: the sub-diagonal elements of the lower triangular matrix L. The upper triangle of A is unchanged.
	int tda, // the second dimension of the array a as declared in the function from which nag real cholesky is called.
	double p[], // the reciprocals of the diagonal elements of L.
	double *detf, // 
	int *dete // the determinant of A is given by detf × 2.0dete. It is given in this form to avoid overflow or underflow.
	);

/** f03afc
		computes an LU factorization of a real matrix, with 
		partial pivoting, and evaluates the determinant.

Example:
	To compute the LU factorization with partial pivoting, and 
	calculate the determinant, of the real matrix
	.
	.
	33 16 72
   -24 -10 -57
   -8 -4 -17
    .
    ..
    
void test_nag_real_lu()
{
	double detf, two=2.0;
	int i, dete, j, n;
	matrix a;
	a.SetSize(8,8);
	int pivot[8];
	
	n = 3;
	a[0][0] = 33;
	a[0][1] = 16;
	a[0][2] = 72;
	
	a[1][0] = -24;
	a[1][1] = -10;
	a[1][2] = -57;
	
	a[2][0] = -8;
	a[2][1] = -4;
	a[2][2] = -17;
		

	nag_real_lu(n,a,8,pivot,&detf, &dete);

	printf("Array A after factorization\n");
	for (i=0; i<n; i++)
		for (j=0; j<n; j++)
		{
			printf("%9.4f",a[i][j]);
			if(j%8==7 || j==n-1)
				printf("\n");
		}
	printf("\nArray P\n");
	for (i=0; i<n; i++)
	{
		printf("%3ld",pivot[i]);
		if(i%8==7 || i==n-1) 
			printf("\n");
	}		
	printf("\ndetf = %9.4f dete = %2ld\n", detf, dete);
	detf = detf * pow(two, 1.0*dete);
	printf("\nValue of determinant = %9.4f\n", detf);	
}

	The output is following:
	
	Array A after factorization
	
	 -8.0000  0.5000  2.1250
	-24.0000  2.0000 -3.0000
	 33.0000 -0.5000  0.3750
	 
	Array P
	
	3 2 3
	
	detf = 0.3750 dete = 4
	
	Value of determinant = 6.0000
	

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	357: The matrix A is singular, possibly due to rounding errors. The factorization could not be completed. detf and dete are set to zero.
	11: On entry, n must not be less than 1: n = _value_.
	17: O n entry, tda = _value_ while n = _value_. The parameters must satisfy tda = n.
	73: Memory allocation failed.

	successfully call of the nag_real_lu function.


*/
	int  nag_real_lu(
	int n, // the order of the matrix A.
	double a[], // Input: the n by n matrix A.  Output: A is overwritten by the lower triangular matrix L and the o.-diagonal elements of the upper triangular matrix U. The unit diagonal elements of U are not stored.
	int tda, // the second dimension of the array a as declared in the function from which nag real lu is called.
	int pivot[], // pivot[i - 1] gives the row index of the ith pivot.
	double *detf, //
	int *dete // the determinant of A is given by detf × 2.0dete. It is given in this form to avoid overflow or underflow.
	);

/**	f03ahc
		computes an LU factorization of a complex matrix, 
		with partial pivoting, and evaluates the determinant.

Example:
	To compute an LU factorization, with partial pivoting, and 
	calculate the determinant, of the complex matrix
	.
	.
	2 1+2i 2 + 10i
	1 + i 1 + 3i -5+ 14i
	1 + i 5i -7+ 20i
	.
	..

void test_nag_complex_lu()
{

	int pivot[5];
	matrix<complex> a;
	a.SetSize(5,5);
	complex det;
	int i, j, n, dete;
	n = 3;
	a[0][0] = 2.0 + 0.0i;
	a[0][1] = 1.0 + 2.0i;
	a[0][2] = 2.0 + 10.0i;
	
	a[1][0] = 1.0 + 1.0i;
	a[1][1] = 1.0 + 3.0i;
	a[1][2] = -5.0 +14.0i;
	
	a[2][0] = 1.0 + 1.0i;
	a[2][1] = 0.0 + 5.0i;
	a[2][2] = -7.0 + 20.0i;
		
	nag_complex_lu(n, a, 5, pivot, &det, &dete);
	
	printf("Array a after factorization\n");
	for (i=0; i<n; i++)
	{
		for (j=0; j<n; j++)
		{
			complex temp  =  a[i][j];
			printf("(%7.3f, %7.3f) ", temp.m_re, temp.m_im);
		}
		printf("\n");
	}
	printf("\nArray pivot\n");
	for (i=0; i<n; i++)
		printf("%5ld",pivot[i]);
	printf("\n\ndet.re = %7.4f, det.im = %7.4f, dete = %2ld.\n",det.m_re, det.m_im, dete);
 	//the determinent given in this form to avoid overflow and underflow
 	det.m_re=det.m_re * pow(2,dete);
	det.m_im=det.m_im * pow(2,dete);	
	
	printf("\nValue of determinant = (%7.4f, %7.4f)\n", det.m_re, det.m_im);
	
}

	The output is following:
	
	Array a after factorization
	
	( 2.000, 0.000) ( 0.500, 1.000) ( 1.000, 5.000)
	( 1.000, 1.000) ( 0.500, 3.500) ( 3.800, 1.400)
	( 1.000, 1.000) ( 1.500, 1.500) ( -4.600, 0.200)
	
	Array pivot
	
	1 3 3
	
	det.re = 0.0234, det.im = 0.1250, dete = 8.
	Value of determinant = ( 6.0000, 32.0000)
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	357: The matrix A is singular, possibly due to rounding errors. The factorization could not be completed. detf and dete are set to zero.
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry, tda = _value_ while n = _value_. The parameters must satisfy tda = n.
	73: Memory allocation failed.

	successfully call of the nag_complex_lu function.
	

*/

int  nag_complex_lu(
	int n, // the order of the matrix A. 
	complex a[], // Input: the n by n matrix A.  Output: A is overwritten by the lower triangular matrix L and the o.-diagonal elements of the upper triangular matrix U. The unit diagonal elements of U are not stored.
	int tda, // the second dimension of the array a as declared in the function from which nag complex lu is called.
	int pivot[], // pivot[i - 1] gives the row index of the ith pivot.
	complex *detf, //
	int *dete // the determinant of A is given by (det.re + idet.im) × 2.0dete. It is given in this form to avoid overflow and underflow.  
);
	
/** f04adc
		calculates the solution of a set of complex linear equations
		with multiple right-hand sides.

Example:
	To solve the set of linear equations AX = B where
	A = .
	.
	1     1+2i    2+ 10i
	1 + i 3i     -5+ 14i
	1 + i 5i     -8+ 20i
	.
	.
	and
	B = .
	.
	1
	0
	0
	.
	..

void test_nag_complex_lin_eqn_mult_rhs()
{
	matrix<complex> a, b, x;
	a.SetSize(5,5);
	b.SetSize(5,1);
	x.SetSize(5,1);
	int i, j, n, l = 1;

	n = 3;
	a[0][0] = 1.0 + 0.0i; 
	a[0][1] = 1.0 + 2.0i;
	a[0][2] = 2.0 +10.0i;
	
	a[1][0] = 1.0 + 1.0i;
	a[1][1] = 0.0 + 3.0i;
	a[1][2] = -5.0 +14.0i;
	
	a[2][0] = 1.0 + 1.0i;
	a[2][1] = 0.0 + 5.0i;
	a[2][2] = -8.0 + 20.0i;
	
	
	b[0][0] = 1.0 + 0.0i;
	b[1][0] = 0.0 + 0.0i;
	b[2][0] = 0.0 + 0.0i;
	
	nag_complex_lin_eqn_mult_rhs(n, 1, a, 5, b, 1, x, 1);
	
	printf("Solution\n");
	for (i=0; i<n; i++)
	{
		complex temp = x[i][0];
		printf("(%7.4f, %7.4f)\n", temp.m_re, temp.m_im);	
	}
}

The output is following:

Solution
(10.0000, 1.0000)
( 9.0000, -3.0000)
(-2.0000, 2.0000)

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	357: The matrix A is singular, possibly due to rounding errors.
	11: On entry, n must not be less than 1: n = _value_.  On entry, nrhs must not be less than 1: nrhs = _value_.
	17: On entry, tda = _value_ while n = _value_. The parameters must satisfy tda = n.  On entry, tdb = _value_ while nrhs = _value_. These parameters must satisfy tdb = nrhs.  On entry, tdx = _value_ while nphs = _value_. These parameters must satisfy tdx = nrhs.
	73: Memory allocation failed.
	
	successfully call of the nag_complex_lin_eqn_mult_rhs function.
*/

	int  nag_complex_lin_eqn_mult_rhs(
	int n, // the order of the matrix A.
	int nrhs, // the number of right-hand sides.
	complex a[], // Input: the n by n matrix A.  Output: A is overwritten by the lower triangular matrix L and the o.-diagonal elements of the upper triangular matrix U. The unit diagonal elements of U are not stored.
	int tda, // the last dimension of the array a as declared in the function from which nag complex lin eqn mult rhs is called.
	complex b[], // the n by r right-hand side matrix B.
	int tdb,  // the last dimension of the array b as declared in the function from which nag complex lin eqn mult rhs is called.
	complex x[], // the n by r solution matrix X. 
	int tdx // the last dimension of the array x as declared in the function from which nag complex lin eqn mult rhs is called.
	);
	
/** f04agc
		calculates the approximate solution of a set of real symmetric 
		positive-definite linear equations with multiple right-hand sides, 
		AX = B, where A has been factorized by nag real cholesky (f03aec).

Example:
	To solve the set of linear equations AX = B where
	A =
	.
	..
	5 7 6 5
	7 10 8 7
	6 8 10 9
	5 7 9 10
	.
	..
	and
	B =
	.
	..
	23
	32
	33
	31
	.
	..
	.


void test_nag_real_cholesky_solve_mult_rhs()
{	
	double d1;
	int i, id, j, n;
	matrix a, b, x;
	a.SetSize(8,8);
	b.SetSize(8,1);
	double  p[8];
	x.SetSize(8,1);
	n = 4;
	a[0][0] = 5;
	a[0][1] = 7; 
	a[0][2] = 6;
	a[0][3] = 5;
	
	a[1][0] = 7;
	a[1][1] = 10;
	a[1][2] = 8;
	a[1][3] = 7;
	
	a[2][0] = 6;
	a[2][1] = 8;
	a[2][2] = 10;
	a[2][3] = 9;
	
	a[3][0] = 5;
	a[3][1] = 7;
	a[3][2] = 9;
	a[3][3] = 10;
	
		
	b[0][0] = 23;
	b[1][0] = 32;
	b[2][0] = 33;
	b[3][0] = 31;
	
	
	nag_real_cholesky(n, a, 8, p, &d1, &id);
	nag_real_cholesky_solve_mult_rhs(n, 1, a, 8, p, b, 1, x, 1);

	printf("\n Solution\n");
	for (i=0; i<n; ++i)
	{
		for (j=0; j <1 ; ++j)
			printf("%9.4f", x[i][j]);
		printf("\n");
	}
}

The output is following:

Solution
1.0000
1.0000
1.0000
1.0000

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 1: n = _value_.  On entry, nrhs must not be less than 1: nrhs = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdb = _value_while nrhs = _value_. These parameters must satisfy tdb = nrhs.  On entry, tdx = _value_while nrhs = _value_. These parameters must satisfy tdx = nrhs.

	successfully call of the nag_real_cholesky_solve_mult_rhs function.
*/

	int  nag_real_cholesky_solve_mult_rhs(
	int n, // the order of the matrix A.
	int nrhs, // the number of right-hand sides.
	double a[], // the upper triangle of the n by n positive-de.nite symmetric matrix A, and the subdiagonal elements of its Cholesky factor L, as returned by nag real cholesky (f03aec).
	int tda, // : the second dimension of the array a as declared in the function from which nag real cholesky solve mult rhs is called.
	double p[], // the reciprocals of the diagonal elements of L, as returned by nag real cholesky (f03aec).
	double b[], // the n by r right-hand side matrix B.
	int tdb, // the second dimension of the array b as declared in the function from which nag real cholesky solve mult rhs is called.
	double x[], // the n by r solution matrix X. 
	int tdx // the second dimension of the array x as declared in the function from which nag real cholesky solve mult rhs is called.
	);


/** f04ajc
		calculates the approximate solution of a set of real linear
		equations with multiple right-hand sides, AX = B, 
		where A has been factorized by nag real lu (f03afc).

Example:
	To solve the set of linear equations AX = B where
	A = .
	.
	 33 16 72
	-24 -10 -57
	-8 -4 -17
	.
	.
	and
	B = .
	.
	-359
	 281
	  85
	 .
	..

void test_nag_real_lu_solve_mult_rhs()
{
	int TDB =  3;
	int NMAX = 8;
	int TDA = 8;
	double detf;
	int i, dete, j, n, nrhs = 1;
	matrix a, b;
	a.SetSize(NMAX,TDA);
	b.SetSize(NMAX, TDB);
	int pivot[8];

	n = 3;
	a[0][0] = 33;
	a[0][1] = 16;
	a[0][2] = 72;
	
	a[1][0] = -24;
	a[1][1] = -10;
	a[1][2] = -57;
	
	a[2][0] = -8;
	a[2][1] = -4;
	a[2][2] = -17;
	
	b[0][0] = -359;
	b[1][0] = 281;
	b[2][0] = 85;
	
	nag_real_lu(n, a, TDA, pivot, &detf,&dete);
	
	nag_real_lu_solve_mult_rhs(n, nrhs, a, TDA, pivot, b, TDB);
	
	printf("Solution\n");
	for (i=0; i<n; i++)
		for (j=0; j<nrhs; j++)
			printf("%9.4f\n",b[i][j]);
}

	The output is following:
	
	Solution
	1.0000
	-2.0000
	-5.0000

Parameters:

Return:
	This function returns NAG error code, 0 if no error.

	11: On entry, n must not be less than 1: n = _value_.  On entry, nrhs must not be less than 1: nrhs = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdb = _value_ while nrhs = _value_. These parameters must satisfy tdb = nrhs.

	successfully call of the nag_real_lu_solve_mult_rhs function.

*/

	int  nag_real_lu_solve_mult_rhs(
	int n, // the order of the matrix A.
	int nrhs, // the number of right-hand sides.
	double a[], // details of the LU factorization, as returned by nag real lu (f03afc).
	int tda, // the second dimension of the array a as declared in the function from which nag real lu solve mult rhs is called.
	int pivot[], // details of the row interchanges as returned by nag real lu (f03afc).
	double b[], // Input: the n by r right-hand side matrix B.  Output: B is overwritten by the solution matrix X.
	int tdb // the second dimension of the array b as declared in the function from which nag real lu solve mult rhs is called.
	);


/** f04akc
		calculates the approximate solution of a set of complex
		linear equations with multiple right-hand sides AX = B, 
		where A has been factorized by nag complex lu (f03ahc).

Example:
	To solve the set of linear equations AX = B where
	A = .
	.
	1 1+2i 2+ 10i
	1 + i 3i -5+ 14i
	1 + i 5i -8+ 20i
	.
	.
	and
	B = .
	.
	1
	0
	0
	.
	..


void test_nag_complex_lu_solve_mult_rhs()
{
	int NMAX = 5;
	int TDA =  NMAX;
	int TDB  = NMAX;
	
	complex det;
	matrix<complex>a, b;
	a.SetSize(NMAX,TDA);
	b.SetSize(NMAX,TDB);
	int i, j, n, dete, nrhs = 1;
	int pivot[5];


	n = 3;
	a[0][0] = 1.0 + 0.0i;
	a[0][1] = 1.0 + 2.0i;
	a[0][2] = 2.0 +10.0i;
	
	a[1][0] = 1.0 + 1.0i;
	a[1][1] = 0.0 + 3.0i;
	a[1][2] = -5.0 + 14.0i;
	
	a[2][0] = 1.0 + 1.0i;
	a[2][1] = 0.0 + 5.0i;
	a[2][2] = -8.0 + 20.0i;
	
	
	b[0][0] = 1.0 + 0.0i;
	b[1][0] = 0.0 + 0.0i;
	b[2][0] = 0.0 + 0.0i;
	
	nag_complex_lu(n, a, TDA, pivot, &det, &dete);
	
	nag_complex_lu_solve_mult_rhs(n, nrhs, a, TDA, pivot, b, TDB);
	printf("Solution\n");
	for (i=0; i<n; i++)
	{
		for (j=0; j<nrhs; j++)
		{
			complex temp = b[i][j];
			printf("(%7.3f, %7.3f) ", temp.m_re, temp.m_im);
		}
		printf("\n");
	}

}

	The output is following:
	
	Solution
	( 10.000, 1.000)
	( 9.000, -3.000)
	( -2.000, 2.000)


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 1: n = _value_.  On entry, nrhs must not be less than 1: nrhs = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.  On entry, tdb = _value_ while nrhs = _value_. These parameters must satisfy tdb = nrhs.
	
	successfully call of the  nag_complex_lu_solve_mult_rhs function.

*/

	int  nag_complex_lu_solve_mult_rhs(
	int n, // the order of the matrix A.
	int nrhs, // the number of right-hand sides.
	complex a[], // details of the LU factorization, as returned by nag complex lu (f03ahc).
	int tda, // the second dimension of the array a as declared in the function from which nag complex lu solve mult rhs is called.
	int pivot[], // details of the row interchanges as returned by nag complex lu (f03ahc).
	complex b[], // Input: the n by r right-hand side matrix B.  Output: B is overwritten by the solution matrix X.
	int tdb // the second dimension of the array b as declared in the function from which nag complex lu solve mult rhs is called.
	);


/** f04arc
		calculates the approximate solution of a set of real linear equations with
		a single right-hand side, using an LU factorization with partial pivoting.

Example:
	To solve the set of linear equations Ax = b where
	A = .
	.
	 33 16 72
	-24 -10 -57
	-8 -4 -17
	.
	.
	and
	B = .
	.
	-359
	 281
	  85
	 .
	 ..

void test_nag_real_lin_eqn()
{
	int NMAX = 8;
	int TDA = NMAX;
	
	matrix a;
	a.SetSize(NMAX,TDA);
	double  b[8], x[8];
	int i, j, n;
	n = 3;
	a[0][0] = 33;
	a[0][1] = 16;
	a[0][2] = 72;
	
	
	a[1][0] = -24;
	a[1][1] = -10;
	a[1][2] = -57;
	
	a[2][0] = -8;
	a[2][1] = -4;
	a[2][2] = -17;
	
	b[0] = -359;
	b[1] = 281;
	b[2] = 85;
	

	nag_real_lin_eqn(n, a, TDA, b, x);
	printf("Solution\n");
	for (i=0; i<n; i++)
		printf("%9.4f\n",x[i]);
	

}

	The output is following:
	
	Solution
	1.0000
	-2.0000
	-5.0000

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 1: n = _value_.
	17: On entry, tda = _value_ while n = _value_. These parameters must satisfy tda = n.
	357: The matrix A is singular, possibly due to rounding errors.
	73: Memory allocation failed.

	successfully call of the nag_real_lin_eqn function.



*/

	int  nag_real_lin_eqn(
	int n, // the order of the matrix A.
	double a[], // Input: the n by n matrix A.  Output: A is overwritten by the lower triangular matrix L and the o.-diagonal elements of the upper triangular matrix U. The unit diagonal elements of U are not stored.
	int tda, // the second dimension of the array a as declared in the function from which nag real lin eqn is called.
	double b[], // the right-hand side vector b.
	double x[] // the solution vector x.
	);

/** f04awc
		calculates the approximate solution of a set of complex Hermitian 
		positive-definite linear equations with multiple right-hand sides, AX = B, 
		where A has been factorized by nag complex cholesky (f01bnc).

Example:
	To solve the set of linear equations AX = B where A is the 
	positive-definite Hermitian matrix:
	.
	..
	15 1 - 2i 2 -4 + 3i
	1 + 2i 20 -2 + i 3 - 3i
	2 -2 - i 18 -1 + 2i
	-4 - 3i 3 + 3i -1 - 2i 26
	.
	..
	and B is the single column vector:
	.
	..
	25 + 34i
	21 + 19i
	12 - 21i
	21 - 27i
	.
	..
	.

void test_nag_hermitian_lin_eqn_mult_rhs()
{
	int NMAX =  8;
	int TDA =  NMAX;
	int TDB =  2;
	int TDX = 2;
	int i, j, n;
	matrix<complex> a, b, x;
	a.SetSize(NMAX,TDA);
	b.SetSize(NMAX,TDB);
	x.SetSize(NMAX,TDX);
	double p[8];
	int nrhs = 1;
	
	n = 4;
	a[0][0] = 15.0 + 0.0i;
	
	a[1][0] = 1.0 + 2.0i;
	a[1][1]	= 20.0 + 0.0i;
	
	a[2][0] = 2.0 + 0.0i;
	a[2][1] = -2.0 -1.0i;
	a[2][2] = 18.0 + 0.0i;
	
	a[3][0] = -4.0 - 3.0i;
	a[3][1] = 3.0 + 3.0i;
	a[3][2] = -1.0 - 2.0i;
	a[3][3] = 26.0 + 0.0i;
	
	b[0][0] = 25.0 + 34.0i;
	b[1][0] = 21.0 + 19.0i;
	b[2][0] = 12.0 - 21.0i;
	b[3][0] = 21.0 - 27.0i; 
	
	
	
	nag_complex_cholesky(n, a, TDA, p);
	
	nag_hermitian_lin_eqn_mult_rhs(n, nrhs, a, TDA, p, b, TDB, x, TDX);
	
	printf("\nSolution\n\n");
	for (i=0; i<n; ++i)
	{
		for (j=0; j<nrhs; ++j)
		{
			complex temp = x[i][j];
	 		printf(" (%7.4f,%7.4f)", temp.m_re, temp.m_im);
	 	}
		printf("\n");
	}

}

	The output is following:
	
	Solution
	( 1.4917, 2.1788)
	( 1.1629, 0.7698)
	( 0.5506,-1.3977)
	( 0.8691,-0.7655)

Parameters:

Return:
	This function returns NAG error code, 0 if no error.

	11: On entry, nrhs must not be less than 1: nrhs = _value_.  On entry, n must not be less than 1: n = _value_.
	17: On entry, tda = _value_ while n = _value_. These paremters must satisfy tda = n.  On entry, tdb = _value_ while nrhs = _value_. These parameters must satisfy tdb = nrhs.  On entry, tdx = _value_ while nrhs = _value_. These parameters must satisfy tdx = nrhs.

	successfully call of the nag_hermitian_lin_eqn_mult_rhs function.
*/

int  nag_hermitian_lin_eqn_mult_rhs(
	int n, // the order of the matrix A.
	int nrhs, // the number of right-hand sides.
	complex a[], // the o.-diagonal elements of the upper triangular matrix U as returned by nag complex cholesky (f01bnc). The lower triangle of the array is not used.
	int tda,  // the second dimension of the array a as declared in the function from which nag hermitian lin eqn mult rhs is called.
	double p[], // the reciprocals of the diagonal elements of U, as returned by nag complex cholesky (f01bnc).
	complex b[], // the n by r right-hand side matrix B. 
	int tdb, // the second dimension of the array b as declared in the function from which nag hermitian lin eqn mult rhs is called.
	complex x[], // the n by r solution matrix X. 
	int tdx // the second dimension of the array x as declared in the function from which nag hermitian lin eqn mult rhs is called.
);

/** f04mcc
		computes the approximate solution of a system of real linear equations 
		with multiple right-hand sides, AX = B, where A is a symmetric 
		positive-definite variable-bandwidth matrix, which has previously 
		been factorized by nag real cholesky skyline (f01mcc). Related 
		systems may also be solved.

Example:
	To solve the system of equations AX = B, where
	A =
	.
	......
	1 2 0 0 5 0
	2 5 3 0 1 40
	0 3 1 3 01 8 0
	0 0 0 1 68 24
	5 14 18 8 55 17
	0 0 0 24 1 77
	.
	......
	and
	B =
	.
	......
	6 -10
	15 -21
	11 -3
	0 24
	51 -39
	46 67
	.
	......
	.
	Here A is symmetric and positive-definite and must first 
	be factorized by nag real cholesky skyline


void test_nag_real_cholesky_skyline_solve()
{
	int NMAX = 6;
	int NRHSMAX  = 2;
	int TDB  = NRHSMAX;
	int TDX =  NRHSMAX;
	int LALMAX  = 14;
	
	int i, nrhs, k, k1, k2, lal, n;
	matrix b, x;
	b.SetSize(NMAX,TDB);
	x.SetSize(NMAX, TDX);
	double a[14] = {1.0, 2.0, 5.0, 3.0, 13.0, 16.0, 5.0, 14.0, 18.0, 8.0, 55.0, 24.0, 17.0, 77.0};
	double al[14], d[6];
	int row[6] = {1, 2, 2, 1, 5, 3};
	Nag_SolveSystem select;
	n = 6;
	k2 = 0;
	nrhs = 2;	
	lal = 14;
	
	b[0][0] = 6.0;
	b[0][1] = -10.0;
	
	b[1][0] = 15.0;
	b[1][1] = -21.0;
	
	b[2][0] = 11.0;
	b[2][1] = -3.0;
	
	b[3][0] = 0.0;
	b[3][1] = 24.0;
	
	b[4][0] = 51.0;
	b[4][1] = -39.0;
	
	b[5][0] = 46.0;
	b[5][1] = 67.0;		
	
	nag_real_cholesky_skyline(n, a, lal, row, al, d);
	select = Nag_LDLTX;
	nag_real_cholesky_skyline_solve(select, n, nrhs, al, lal, d, row, b, TDB, x, TDX);
	printf("\n Solution\n");
	for (i=0; i<n; ++i)
	{
		for (k=0; k<nrhs; ++k)
			printf("%9.3f",x[i][k]);
		printf("\n");
	}
}

The output is following:

Solution
   -3.000    4.000
    2.000   -2.000
   -1.000    3.000
   -2.000    1.000
    1.000   -2.000
    1.000    1.000
 

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 1: n = _value_.  On entry, row[_value_] must not be less than 1: row[_value_] = _value_.  On entry, nrhs must not be less than 1: nrhs = _value_.
	19: On entry, row[i] = _value_ while i = _value_. These parameters must satisfy row[i] = i + 1 .
	17: On entry, lal = _value_ while row[0] + . . .+row[n - 1 ] =_value_. These parameters must satisfy lal = row[0] + . . .+row[n - 1].  On entry, tdb = _value_ while nrhs = _value_. These parameters must satisfy tdb = nrhs.  On entry, tdx = _value_ while nrhs = _value_. These parameters must satisfy tdx = nrhs.
	70: On entry, parameter selct had an illegal value.
	370: The diagonal matrix D is singular as it has at least one zero element. The .rst zero element has been located in the array d[_value_]
	371: The lower triangular matrix L has at least one diagonal element which is not equal to unity.  The first non-unit element has been located in the array al[_value_]

	successfully call of the nag_real_cholesky_skyline_solve function.
*/

int  nag_real_cholesky_skyline_solve(
	Nag_SolveSystem selct, 
	int n, // the order of the matrix L.
	int nrhs, // the number of right-hand sides.
	double al[], // the elements within the envelope of the lower triangular matrix L, taken in row by row order, as returned by nag real cholesky skyline (f01mcc). The unit diagonal elements of L must be stored explicitly.
	int lal, // the dimension of the array al as declared in the function from which nag real cholesky skyline solve is called.
	double d[], // the diagonal elements of the diagonal matrix D. d is not referenced if selct = Nag LLTX, Nag LX or Nag LTX
	int row[], // row[i] must contain the width of row i of L, i.e., the number of elements between the first (left-most) non-zero element and the element on the diagonal, inclusive.
	double b[], // the n by r right-hand side matrix B. 
	int tdb, // the second dimension of the array b as declared in the function from which nag real cholesky skyline solve is called.
	double x[], // the n by r solution matrix X. 
	int tdx // the second dimension of the array x as declared in the function from which nag real cholesky skyline solve is called.
);

/* end proto */
#endif //!_O_NAG_f_H