//  matrix_m03.cpp test the matrix inversion in c03.

//  this file has a triple purpose:  1) my testing.  2) demonstration.  3) your testing.



#include <stdio.h>

#include "matrix_dll.h"
#include "StdAfx.h"
//#include "matrix_h0s.h"
//#include "matrix_h01.h"
//#include "matrix_h02.h"
//#include "matrix_h03.h"
//#include "matrix_h04.h"
//#include "matrix_h05.h"
//#include "matrix_h06.h"
//#include "matrix_h07.h"
//#include "matrix_h08.h"
#include "matrix_h09.h"


// get yourself enough arrays.  any unused ones consume neglegible storage.
#define amt_of_arrays 15


matrix_preferences mat_pref;		// declare the global preferences.




int main() {
	FILE *outfile;
	const double ZERO=0.0;
	MAT array[amt_of_arrays], scratch[amt_of_arrays];
	CHN chan[amt_of_arrays];
	double determinant, norm_left, norm_right, trace, relative, residue_bra, residue_ket,
		relative_array, relative_bra, relative_diag, relative_ket, sum, R, Theta, rho;
	unsigned long i, n, m0, m1, m2, mlicr, mlicr_mx, step, sort=0, style, complex, imp;
	long nn;
	int y=0;
	char *string0[]={"without pivoting", "with pivoting"},
		*string1[]={"system", "Gauss-Jordan", "LDU", "partition & conquer"},
		*string2[]={"system", "once", "twice"},
		*string3[]={"system", "Jacobi", "squaring"},
		*string4[]={"original", "by rows", "by diagonals"},
		*string5[]={"positive symmetric", "symmetric", "positive symmetric tested as general",
			"symmetric tested as real", "general"},
		*string6[]={"real", "complex"},
		*string7[]={"always +  +", "always -  +", "small  +", "large  +",
			"always +  -", "always -  -", "small  -", "large  -"};
	extern matrix_preferences mat_pref;		// so that you may access them.

	outfile = fopen("test.txt", "w");

	for(i=0; i<amt_of_arrays; i++) {	// required for the automatic storage allocation to work.
		array[i].mx = 0; scratch[i].mx = 0;
		matrix_initialize_channel(&chan[i]); };



	// *********** the preferences ***************

	matrix_set_preferences();			// required to set the defaults.
//	for safety, do this even if you intend to over-ride all of them.

	// over-ride any (or all) of them.
//	mat_pref.eigenpairs=1;		// 0=system default 1=Jacobi 2=squaring.
	mat_pref.Jacobi_style=0;	// 0=original 1=simplified - by rows  2=simplified - by super-diagonals.
//	mat_pref.test_Jac = 1.0e-18;	// test criterion for Jacobi.
//	mat_pref.Jacobi_angle=0;	// 0=always+ 1=always- 2=small 3=large.  add 4 to have cosine negative.
	mat_pref.random=1;			// 0=system default 1=199017 2=352345383.
	mat_pref.seed=54321;		// the seed for the pseudo-random number generator.

//	mat_pref.test_Jac = 1.0e-10;	// temporarily. <<<<<<<<<<<<<<<<<
//	mat_pref.singular = 1.0e-10;	// temp 


	// ********** ask for the test situation to be specified.  ******************
	
	matrix_about();		// displays the copyright notice.
capo:
	printf("enter the size of the matrix -- negative to test inversion -- (zero to quit) "); scanf("%d", &nn);
	if(nn == 0) goto tie_point;		// user requested to quit.
	n = matrix_abs(nn);
	printf("size %d   enter the mlicr_mx ", n); scanf("%d", &mlicr_mx);
	if(mlicr_mx == 0) {
		printf("mlicr_mx %d   enter the style 0=pos sym 1=sym 2=pos sym as general 3=sym as general 4=general ",
			mlicr_mx);
		step = 1; }
	else {
		printf("mlicr_mx %d   enter its step ", mlicr_mx); scanf("%d", &step); if(step == 0) step = 1;
		printf("step %d   enter the style 0=pos sym 1=sym 2=pos sym as general 3=sym as general 4=general ", step); };
			scanf("%d", &style);
	printf("style %d %s   enter kill-style 0=largest 1=rows 2=super-diagonals +4=principal ", style, string5[style]);
		scanf("%d", &mat_pref.Jacobi_style);
	printf("kill-style %d %s   enter complex 0=real 1=complex ",
		mat_pref.Jacobi_style, string4[mat_pref.Jacobi_style%4]); scanf("%d", &complex);
	printf("complex %d %s   enter Jacobi angle 0=always+ 1=always- 2=small 3=large ", complex, string6[complex]);
		scanf("%d", &mat_pref.Jacobi_angle);
	printf("Jacobi angle %d %s   do you want to sort the eigenvalues?  0=no 1=yes ",
		mat_pref.Jacobi_angle, string7[mat_pref.Jacobi_angle]);
		scanf("%d", &sort);
	printf("sort %d \n", sort);



	// ********* generate the requested test matrix. *****************************************

	for(mlicr=0; mlicr<=mlicr_mx; mlicr+=step) {
		printf("\nthis is a test of a matrix size %d with mlicr of %d style is %s %s\n\n",
				n, mlicr, string5[style], string6[complex]);
		if(mlicr == 0) {
			switch(style) {
				case 0:		// positive symmetric.
				case 2:		// but, test as general.
					if(complex == 0) {
						y = matrix_random_positive_symmetric(n, &array[0], &scratch[0]); }
					else {
						y = matrix_random_positive_Hermitian(n, &array[0], &scratch[0],
								&scratch[1], &scratch[2], &scratch[3]); };
					if(style == 2) (array[0]).sym = 0;
					break;
				case 1:		// symmetric.
				case 3:		// but, test as general.
					if(complex == 0) {
						y = matrix_random_symmetric_and_skew(n, &array[0], &scratch[0], &scratch[1],
								&scratch[2], &scratch[3]); }
					else {
						y = matrix_random_Hermitian_and_skew(n, &array[0], &scratch[0], &scratch[1],
							&scratch[2], &scratch[3], &scratch[4]); };
					if(style == 3) (array[0]).sym = 0;
					break;
				case 4:		// general.
				default:
					array[0].imx = n, array[0].jmx = n;
					if(complex == 0) {
						y = matrix_uniformly_random(&array[0]); }
					else {
						y = matrix_random_complex(n, n, &array[0], &scratch[0], &scratch[1], &scratch[2]); };
					break; };
			if(y != 0) goto tie_point; }
		else {
			switch(style) {
				case 0:		// positive symmetric.
				case 2:		// but, test as general.
					if(complex == 0) {
						y = matrix_random_positive_symmetric(n, &array[0], &scratch[0]); }
					else {
						y = matrix_random_positive_Hermitian(n, &array[0], &scratch[0],
								&scratch[1], &scratch[2], &scratch[3]); };
					if(style == 2) (array[0]).sym = 0;
					break;
				case 1:		// symmetric (and factors).
				case 3:		// but, test as general.
					if(complex == 0) {
						//												  array		 bra	    diag	   ket
						y = matrix_random_symmetic_and_factors(n, mlicr, &array[0], &array[1], &array[2], &array[3],
								&scratch[0], &scratch[1], &scratch[2]); }
					else {
						//												  array		  bra	     diag		ket
						y = matrix_random_Hermitian_and_factors(n, mlicr, &array[0], &array[1], &array[2], &array[3],
								&scratch[0], &scratch[1], &scratch[2], &scratch[3]); };
					if(style == 2) (array[0]).sym = 0;
					break;
				case 4:		// general (and factors).
				default:
					if(complex == 0) {
						//											  array      bra        diag       ket
						y = matrix_random_real_and_factors(n, mlicr, &array[0], &array[1], &array[2], &array[3],
								&scratch[0], &scratch[1], &scratch[2], &scratch[3]); }
					else {
						//											     array      bra        diag       ket
						y = matrix_random_complex_and_factors(n, mlicr, &array[0], &array[1], &array[2], &array[3],
								&scratch[0], &scratch[1], &scratch[2], &scratch[3], &scratch[4]); };
						break; };
			if(y != 0) goto tie_point;
			if(style == 1 || style == 3 || style == 4) {
				y = matrix_isit_orthonormal(&array[1], &residue_bra, 
						&scratch[0], &scratch[1], &scratch[2], &scratch[3]);
				if(y != 0) goto tie_point;
				y = matrix_isit_orthonormal(&array[3], &residue_ket, 
						&scratch[0], &scratch[1], &scratch[2], &scratch[3]);
				if(y != 0) goto tie_point;
				printf("are the bra and ket ortho-normal?  bra=%e  ket=%e\n\n", residue_bra, residue_ket);
				y = matrix_mat_diag_matt_multiply(&array[1], &array[2], &array[1], &scratch[0],
						&scratch[1], &scratch[2]);
				if(y != 0) goto tie_point;
				y = matrix_isit_equal(&array[0], &scratch[0], &relative, &scratch[1], &scratch[2]);
				if(y != 0) goto tie_point;
				printf("reconstituted array, relative=%e \n", relative); }; };
		y = matrix_trace(&array[0], &trace);
			if(y != 0) goto tie_point;




	// *************** run the tests.  *************************************
		


		// ********** Gramm-Schmidt ortho-normalization and random permutation. ******
		
	y = matrix_copy(&array[0], &array[1]);
		if(y != 0) goto tie_point;
//	if(n == 3) {
//		y = matrix_zero(&array[1]);
//			if(y != 0) goto tie_point;
//		array[1].a[0]=2.0; array[1].a[2]=-3.0;
//		array[1].a[3]=3.0; array[1].a[5]=2.0;
//		matrix_display(&array[1], "input to GSch"); };
	y = matrix_insitu_GrammSchmidt(&array[1], 1);
		if(y != 0) goto tie_point;
	y = matrix_isit_orthonormal(&array[1], &relative, &scratch[0], &scratch[1],
			&scratch[2], &scratch[3]);
		if(y != 0) goto tie_point;
	if(array->cmp == 0) {
		relative_bra = ZERO; }
	else {
		y = matrix_isit_complex(&array[1], &relative_bra, &scratch[0], &scratch[1], &scratch[2],
				&scratch[3], &scratch[4], &scratch[5]);
		if(y != 0) goto tie_point; };
	printf("test of Gramm-Schmidt ortho-normalization relative=%e complex=%e\n",
		relative, relative_bra);
	if(n <= 4) {
		matrix_display(&array[1], "ortho-normalized"); };
	y = matrix_random_permutation(n, complex, &array[1]);
		if(y != 0) goto tie_point;
	y = matrix_isit_orthonormal(&array[1], &relative, &scratch[0], &scratch[1],
			&scratch[2], &scratch[3]);
		if(y != 0) goto tie_point;
	if(array->cmp == 0) {
		relative_bra = ZERO; }
	else {
		y = matrix_isit_complex(&array[1], &relative_bra, &scratch[0], &scratch[1], &scratch[2],
				&scratch[3], &scratch[4], &scratch[5]);
		if(y != 0) goto tie_point; };
	printf("test of random permutation relative=%e complex=%e\n", relative, relative_bra);
	if(n <= 4) {
		matrix_display(&array[1], "permutation"); };
	printf("\n");



	// *********** LDU and Gauss-Jordan matrix inversion. ***************


	if(nn < 0) {
		for(m1=1; m1<=2; m1++) {			// algorithm
			for(m0=0; m0<=1; m0++) {		// pivot.
				for(m2=1; m2<=2; m2++) {	// twice.
					mat_pref.inversion = m1;	// set the algorithm.
					mat_pref.pivot = m0;		// set the pivoting.
					mat_pref.twice = m2;		// set the twice.
					printf("inverting by the %s %s algorithm, %s. trace is actual %f claimed %f \n",
						string1[m1], string2[m2], string0[m0], trace, array[0].inh.trace);
					y = matrix_copy(&array[0], &array[1]);
					if(y != 0) goto tie_point;
					y = matrix_insitu_inverse(&array[1], &determinant, &scratch[0], &scratch[1]);
						if(y != 0) goto tie_point;
					y = matrix_isit_inverse(&array[0], &array[1], &norm_left, &norm_right,
							&scratch[0], &scratch[1], &scratch[2], &scratch[3]);
					if(y != 0) goto tie_point;
					printf("test of matrix inversion:  determinant=%e \n  norm left=%e right=%e \n\n",
						determinant, norm_left, norm_right); }; }; }; };


		// ************ Jacobi eigenvalue-eigenvector (eigen-pair) factorization. ************

						//			array	   bra		  diag		 ket			  inverse
	y = matrix_real_factor_inverse(&array[0], &array[2], &array[1], &array[3], sort, &array[5],
			&scratch[1], &scratch[2], &scratch[3], &scratch[4],
			&scratch[5], &scratch[6], &scratch[7], &scratch[8]);	// tests the Jacobi algorithm.
		if(y != 0) goto tie_point;
	y = matrix_isit_orthonormal(&array[2], &relative_bra, &scratch[0], &scratch[1], &scratch[2], &scratch[3]);
		if(y != 0) goto tie_point;
	y = matrix_isit_diagonal(&array[1], &relative_diag, &scratch[0], &scratch[1]);
		if(y != 0) goto tie_point;
	y = matrix_isit_orthonormal(&array[3], &relative_ket, &scratch[0], &scratch[1], &scratch[2], &scratch[3]);
		if(y != 0) goto tie_point;
	y = matrix_mat_diag_mat_multiply(&array[2], &array[1], &array[3], &array[4], &scratch[1]);
		if(y != 0) goto tie_point;
	y = matrix_isit_equal(&array[0], &array[4], &relative, &scratch[1], &scratch[2]);
		if(y != 0) goto tie_point;
	sum = relative_bra + relative_diag + relative_ket;
	printf("Jacobi:  bra=%e diag=%e ket=%e sum=%e relative=%e \n",
				relative_bra, relative_diag, relative_ket, sum, relative);
	y = matrix_isit_inverse(&array[0], &array[5], &norm_left, &norm_right,
			&scratch[1], &scratch[2], &scratch[3], &scratch[4]);
		if(y != 0) goto tie_point;
	printf("    test of matrix inversion by Jacobi:  norm left=%e right=%e \n",
				norm_left, norm_right);
	if(complex != 0) {
		y = matrix_isit_complex(&array[0], &relative_array, &scratch[0], &scratch[1],
				&scratch[2], &scratch[3], &scratch[4], &scratch[5]);
			if(y != 0) goto tie_point;
		y = matrix_isit_complex(&array[2], &relative_bra, &scratch[0], &scratch[1],
				&scratch[2], &scratch[3], &scratch[4], &scratch[5]);
			if(y != 0) goto tie_point;
		y = matrix_isit_complex(&array[1], &relative_diag, &scratch[0], &scratch[1],
				&scratch[2], &scratch[3], &scratch[4], &scratch[5]);
			if(y != 0) goto tie_point;
		y = matrix_isit_complex(&array[3], &relative_ket, &scratch[0], &scratch[1],
				&scratch[2], &scratch[3], &scratch[4], &scratch[5]);
			if(y != 0) goto tie_point;
		sum = relative_array + relative_bra + relative_diag + relative_ket;
		printf("rel complex?  array=%e bra=%e diag=%e ket=%e their sum=%e \n",
				relative_array, relative_bra, relative_diag, relative_ket, sum); };
	if(n <= 4) {
		matrix_display(&array[0], "array");
		matrix_display(&array[2], "bra  ");
		matrix_display(&array[1], "diag ");
		matrix_display(&array[3], "ket  ");
		matrix_display(&array[4], "reconstructed"); };
	printf("\n"); };


	// ********** channel. ***************


	if((&array[0])->cmp == 0 && (&array[0])->sym == 1 && (&array[0])->positive != 0 && n >= 2) {
//		if((n %2) == 0) {
//			y = matrix_copy(&array[0], &array[1]); }
//		else {
//			y = matrix_augment(&array[0], 1, n + 3, 1, &array[1]); };
		y = matrix_copy(&array[0], &array[8]);
		if(y != 0) goto tie_point;
		imp = array[8].imx - array[8].jmx / 2; array[8].part = 1; array[8].imp = imp; array[8].jmp = imp;
		y = matrix_augment_channel(&array[8], 1, n+3, &array[1]);
		if(y != 0) goto tie_point;
		y = matrix_channel_rate(&array[1], &R, &Theta, &rho, &scratch[0], &scratch[1], &scratch[2],
				&scratch[3], &scratch[4], &scratch[5]);
			if(y != 0) goto tie_point;
		printf("predicted channel rate R=%f Theta=%f rho=%f \n", R, Theta, rho);
		array[2].imx = 1; array[2].jmx = array[1].jmx;
		y = matrix_uniformly_random(&array[2]);
			if(y != 0) goto tie_point;
		array[2].jmp = array[1].jmp; array[2].part = array[1].part;
		y = matrix_array_bar_channel(&array[1], &array[2], &chan[0],
			//	 bra		diag	   bra		  diag	  reconstituted		these are large-style.
				&array[3], &array[4], &array[5], &array[6], &array[7], &relative,
				&scratch[0], &scratch[1], &scratch[2], &scratch[3],
				&scratch[4], &scratch[5], &scratch[6], &scratch[7],
				&scratch[8], &scratch[9], &scratch[10], &scratch[11]);
			if(y != 0) goto tie_point;
		y = matrix_channel_rate(&array[7], &R, &Theta, &rho, &scratch[4], &scratch[5], &scratch[6],
				&scratch[7], &scratch[8], &scratch[9]);
			if(y != 0) goto tie_point;
		printf("channel (small style)  R=%f Theta=%f rho=%f \n",
			chan[0].R, chan[0].Theta, chan[0].rho);
		printf("reconstituted relative=%e  R=%f Theta=%f rho=%f \n\n", relative, R, Theta, rho);
		if(n < 6) {
			matrix_display(&array[1], "the augmented array");
			matrix_display(&array[3], "exterior bra");
			matrix_display(&array[4], "exterior diag");
			matrix_display(&array[5], "interior bra");
			matrix_display(&array[6], "interior diag");
			matrix_display(&array[7], "reconstructed");
			matrix_display(&scratch[3], "the d small array");
			matrix_display(&array[2], "the given average");
			matrix_display(&(chan[0].eleven[AVE]), "the large average");
			matrix_display(&(chan[0].eleven[LAVE]), "the left average");
			matrix_display(&(chan[0].eleven[RAVE]), "the right average");

		
			matrix_print(outfile, &array[1], "the augmented array");
			matrix_print(outfile, &array[3], "exterior bra");
			matrix_print(outfile, &array[4], "exterior diag");
			matrix_print(outfile, &array[5], "interior bra");
			matrix_print(outfile, &array[6], "interior diag");
			matrix_print(outfile, &array[7], "reconstructed");
			matrix_print(outfile, &scratch[3], "the d small array"); }; };





goto capo;		// del capo al fine.



	// ********* finish.  ************************************8

tie_point:		// fine.
	for(i=0; i<amt_of_arrays; i++) {	// clean-up the storage.
		y += matrix_dealloc(&array[i]) + matrix_dealloc(&scratch[i]);
		y += matrix_channel_dealloc(&chan[i]); };
	fclose(outfile);
	return(0); };





// copyright (c) 2003,4 by R.I. 'Scibor-Marchocki.  last modified Tuesday 20-th April 2004.

// eof