/*
================================================================================
  RPL/2 (R) version 4.1.20
  Copyright (C) 1989-2015 Dr. BERTRAND Joël

  This file is part of RPL/2.

  RPL/2 is free software; you can redistribute it and/or modify it
  under the terms of the CeCILL V2 License as published by the french
  CEA, CNRS and INRIA.
 
  RPL/2 is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  FITNESS FOR A PARTICULAR PURPOSE.  See the CeCILL V2 License
  for more details.
 
  You should have received a copy of the CeCILL License
  along with RPL/2. If not, write to info@cecill.info.
================================================================================
*/


#include "rpl-conv.h"


/*
================================================================================
  Fonction 'qr'
================================================================================
  Entrées : pointeur sur une structure struct_processus
--------------------------------------------------------------------------------
  Sorties :
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

void
instruction_qr(struct_processus *s_etat_processus)
{
	complex16					registre;
	complex16					*tau_complexe;
	complex16					*vecteur_complexe;

	real8						*tau_reel;
	real8						*vecteur_reel;

	struct_liste_chainee		*registre_pile_last;

	struct_objet				*s_copie_argument;
	struct_objet				*s_matrice_identite;
	struct_objet				*s_objet;
	struct_objet				*s_objet_argument;
	struct_objet				*s_objet_resultat;

	integer8					i;
	integer8					j;
	integer8					k;
	integer8					nombre_reflecteurs_elementaires;

	void						*tau;

	(*s_etat_processus).erreur_execution = d_ex;

	if ((*s_etat_processus).affichage_arguments == 'Y')
	{
		printf("\n  QR ");
		
		if ((*s_etat_processus).langue == 'F')
		{
			printf("(décomposition QR)\n\n");
		}
		else
		{
			printf("(QR décomposition)\n\n");
		}

		printf("    1: %s, %s\n", d_MIN, d_MRL);
		printf("->  2: %s\n", d_MRL);
		printf("    1: %s\n\n", d_MRL);

		printf("    1: %s\n", d_MCX);
		printf("->  2: %s\n", d_MCX);
		printf("    1: %s\n", d_MCX);

		return;
	}
	else if ((*s_etat_processus).test_instruction == 'Y')
	{
		(*s_etat_processus).nombre_arguments = -1;
		return;
	}

	if (test_cfsf(s_etat_processus, 31) == d_vrai)
	{
		if (empilement_pile_last(s_etat_processus, 1) == d_erreur)
		{
			return;
		}
	}

	if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
			&s_objet_argument) == d_erreur)
	{
		(*s_etat_processus).erreur_execution = d_ex_manque_argument;
		return;
	}

	if (((*s_objet_argument).type == MIN) ||
			((*s_objet_argument).type == MRL))
	{
		/*
		 * Matrice entière ou réelle
		 */

		if ((s_copie_argument = copie_objet(s_etat_processus,
				s_objet_argument, 'Q')) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return;
		}

		factorisation_qr(s_etat_processus, (*s_copie_argument).objet, &tau);
		(*s_copie_argument).type = MRL;

		tau_reel = (real8 *) tau;

		if ((*s_etat_processus).erreur_systeme != d_es)
		{
			return;
		}

		if (((*s_etat_processus).exception != d_ep) ||
				((*s_etat_processus).erreur_execution != d_ex))
		{
			free(tau);
			liberation(s_etat_processus, s_objet_argument);
			liberation(s_etat_processus, s_copie_argument);
			return;
		}

		if ((s_objet_resultat = copie_objet(s_etat_processus,
				s_copie_argument, 'O')) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return;
		}

		// Matrice Q

		nombre_reflecteurs_elementaires = ((*((struct_matrice *)
				(*s_copie_argument).objet)).nombre_colonnes <
				(*((struct_matrice *) (*s_copie_argument).objet))
				.nombre_lignes) ? (*((struct_matrice *)
				(*s_copie_argument).objet)).nombre_colonnes
				: (*((struct_matrice *) (*s_copie_argument).objet))
				.nombre_lignes;

		registre_pile_last = NULL;

		if (test_cfsf(s_etat_processus, 31) == d_vrai)
		{
			registre_pile_last = (*s_etat_processus).l_base_pile_last;
			(*s_etat_processus).l_base_pile_last = NULL;
		}

		if ((s_objet = allocation(s_etat_processus, INT)) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return;
		}

		(*((integer8 *) (*s_objet).objet)) = (*((struct_matrice *)
				(*s_copie_argument).objet)).nombre_lignes;

		if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
				s_objet) == d_erreur)
		{
			return;
		}

		instruction_idn(s_etat_processus);

		if (((*s_etat_processus).erreur_systeme != d_es) ||
				((*s_etat_processus).erreur_execution != d_ex) ||
				((*s_etat_processus).exception != d_ep))
		{
			liberation(s_etat_processus, s_copie_argument);
			free(tau);

			if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
			{
				return;
			}

			(*s_etat_processus).l_base_pile_last = registre_pile_last;
			return;
		}

		if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
				&s_matrice_identite) == d_erreur)
		{
			(*s_etat_processus).erreur_execution = d_ex_manque_argument;
			return;
		}

		for(i = 0; i < nombre_reflecteurs_elementaires; i++)
		{
			// Calcul de H(i) = I - tau * v * v'

			if ((s_objet = copie_objet(s_etat_processus, s_matrice_identite,
					'P')) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
					s_objet) == d_erreur)
			{
				return;
			}

			if ((s_objet = allocation(s_etat_processus, REL)) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			(*((real8 *) (*s_objet).objet)) = tau_reel[i];

			if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
					s_objet) == d_erreur)
			{
				return;
			}

			if ((s_objet = allocation(s_etat_processus, MRL)) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			(*((struct_matrice *) (*s_objet).objet)).nombre_lignes =
					(*((struct_matrice *) (*s_copie_argument).objet))
					.nombre_lignes;
			(*((struct_matrice *) (*s_objet).objet)).nombre_colonnes =
					(*((struct_matrice *) (*s_copie_argument).objet))
					.nombre_lignes;

			if ((vecteur_reel = malloc(((size_t) (*((struct_matrice *)
					(*s_objet).objet)).nombre_lignes) * sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
					.nombre_lignes; j++)
			{
				if (j < i)
				{
					vecteur_reel[j] = 0;
				}
				else if (j == i)
				{
					vecteur_reel[j] = 1;
				}
				else
				{
					vecteur_reel[j] = ((real8 **) (*((struct_matrice *)
							(*s_copie_argument).objet)).tableau)[j][i];
				}
			}

			if (((*((struct_matrice *) (*s_objet).objet)).tableau =
					malloc(((size_t) (*((struct_matrice *) (*s_objet).objet))
					.nombre_lignes) * sizeof(real8 *))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
					.nombre_lignes; j++)
			{
				if ((((real8 **) (*((struct_matrice *) (*s_objet).objet))
						.tableau)[j] = malloc(((size_t) (*((struct_matrice *)
						(*s_objet).objet)).nombre_lignes) * sizeof(real8)))
						== NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}

				for(k = 0; k < (*((struct_matrice *) (*s_objet).objet))
						.nombre_colonnes; k++)
				{
					((real8 **) (*((struct_matrice *) (*s_objet).objet))
							.tableau)[j][k] = vecteur_reel[j] * vecteur_reel[k];
				}
			}

			free(vecteur_reel);

			if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
					s_objet) == d_erreur)
			{
				return;
			}

			instruction_multiplication(s_etat_processus);

			if (((*s_etat_processus).erreur_systeme != d_es) ||
					((*s_etat_processus).erreur_execution != d_ex) ||
					((*s_etat_processus).exception != d_ep))
			{
				liberation(s_etat_processus, s_copie_argument);
				liberation(s_etat_processus, s_matrice_identite);
				free(tau);

				if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
				{
					return;
				}

				(*s_etat_processus).l_base_pile_last = registre_pile_last;
				return;
			}

			instruction_moins(s_etat_processus);

			if (((*s_etat_processus).erreur_systeme != d_es) ||
					((*s_etat_processus).erreur_execution != d_ex) ||
					((*s_etat_processus).exception != d_ep))
			{
				liberation(s_etat_processus, s_copie_argument);
				liberation(s_etat_processus, s_matrice_identite);
				free(tau);

				if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
				{
					return;
				}

				(*s_etat_processus).l_base_pile_last = registre_pile_last;
				return;
			}

			if (i > 0)
			{
				instruction_multiplication(s_etat_processus);

				if (((*s_etat_processus).erreur_systeme != d_es) ||
						((*s_etat_processus).erreur_execution != d_ex) ||
						((*s_etat_processus).exception != d_ep))
				{
					liberation(s_etat_processus, s_copie_argument);
					liberation(s_etat_processus, s_matrice_identite);
					free(tau);

					if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
					{
						return;
					}

					(*s_etat_processus).l_base_pile_last = registre_pile_last;
					return;
				}
			}
		}

		if (test_cfsf(s_etat_processus, 31) == d_vrai)
		{
			if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
			{
				return;
			}

			(*s_etat_processus).l_base_pile_last = registre_pile_last;
		}

		// Matrice R

		for(i = 0; i < (*((struct_matrice *) (*s_objet_resultat).objet))
				.nombre_lignes; i++)
		{
			for(j = 0; j < i; j++)
			{
				((real8 **) (*((struct_matrice *) (*s_objet_resultat).objet))
						.tableau)[i][j] = 0;
			}
		}

		if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
				s_objet_resultat) == d_erreur)
		{
			return;
		}

		liberation(s_etat_processus, s_matrice_identite);
		liberation(s_etat_processus, s_copie_argument);
		free(tau);
	}
	else if ((*s_objet_argument).type == MCX)
	{
		/*
		 * Matrice complexe
		 */

		if ((s_copie_argument = copie_objet(s_etat_processus,
				s_objet_argument, 'Q')) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return;
		}

		factorisation_qr(s_etat_processus, (*s_copie_argument).objet, &tau);

		tau_complexe = (complex16 *) tau;

		if ((*s_etat_processus).erreur_systeme != d_es)
		{
			return;
		}

		if (((*s_etat_processus).exception != d_ep) ||
				((*s_etat_processus).erreur_execution != d_ex))
		{
			free(tau);
			liberation(s_etat_processus, s_objet_argument);
			liberation(s_etat_processus, s_copie_argument);
			return;
		}

		if ((s_objet_resultat = copie_objet(s_etat_processus,
				s_copie_argument, 'O')) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return;
		}

		// Matrice Q

		nombre_reflecteurs_elementaires = ((*((struct_matrice *)
				(*s_copie_argument).objet)).nombre_colonnes <
				(*((struct_matrice *) (*s_copie_argument).objet))
				.nombre_lignes) ? (*((struct_matrice *)
				(*s_copie_argument).objet)).nombre_colonnes
				: (*((struct_matrice *) (*s_copie_argument).objet))
				.nombre_lignes;

		registre_pile_last = NULL;

		if (test_cfsf(s_etat_processus, 31) == d_vrai)
		{
			registre_pile_last = (*s_etat_processus).l_base_pile_last;
			(*s_etat_processus).l_base_pile_last = NULL;
		}

		if ((s_objet = allocation(s_etat_processus, INT)) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return;
		}

		(*((integer8 *) (*s_objet).objet)) = (*((struct_matrice *)
				(*s_copie_argument).objet)).nombre_lignes;

		if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
				s_objet) == d_erreur)
		{
			return;
		}

		instruction_idn(s_etat_processus);

		if (((*s_etat_processus).erreur_systeme != d_es) ||
				((*s_etat_processus).erreur_execution != d_ex) ||
				((*s_etat_processus).exception != d_ep))
		{
			liberation(s_etat_processus, s_copie_argument);
			free(tau);

			if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
			{
				return;
			}

			(*s_etat_processus).l_base_pile_last = registre_pile_last;
			return;
		}

		if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
				&s_matrice_identite) == d_erreur)
		{
			(*s_etat_processus).erreur_execution = d_ex_manque_argument;
			return;
		}

		for(i = 0; i < nombre_reflecteurs_elementaires; i++)
		{
			// Calcul de H(i) = I - tau * v * v'

			if ((s_objet = copie_objet(s_etat_processus, s_matrice_identite,
					'P')) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
					s_objet) == d_erreur)
			{
				return;
			}

			if ((s_objet = allocation(s_etat_processus, CPL)) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			(*((complex16 *) (*s_objet).objet)) = tau_complexe[i];

			if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
					s_objet) == d_erreur)
			{
				return;
			}

			if ((s_objet = allocation(s_etat_processus, MCX)) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			(*((struct_matrice *) (*s_objet).objet)).nombre_lignes =
					(*((struct_matrice *) (*s_copie_argument).objet))
					.nombre_lignes;
			(*((struct_matrice *) (*s_objet).objet)).nombre_colonnes =
					(*((struct_matrice *) (*s_copie_argument).objet))
					.nombre_lignes;

			if ((vecteur_complexe = malloc(((size_t) (*((struct_matrice *)
					(*s_objet).objet)).nombre_lignes) * sizeof(complex16)))
					== NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
					.nombre_lignes; j++)
			{
				if (j < i)
				{
					vecteur_complexe[j].partie_reelle = 0;
					vecteur_complexe[j].partie_imaginaire = 0;
				}
				else if (j == i)
				{
					vecteur_complexe[j].partie_reelle = 1;
					vecteur_complexe[j].partie_imaginaire = 0;
				}
				else
				{
					vecteur_complexe[j].partie_reelle = ((complex16 **)
							(*((struct_matrice *) (*s_copie_argument).objet))
							.tableau)[j][i].partie_reelle;
					vecteur_complexe[j].partie_imaginaire = ((complex16 **)
							(*((struct_matrice *) (*s_copie_argument).objet))
							.tableau)[j][i].partie_imaginaire;
				}
			}

			if (((*((struct_matrice *) (*s_objet).objet)).tableau =
					malloc(((size_t) (*((struct_matrice *) (*s_objet).objet))
					.nombre_lignes) * sizeof(complex16 *))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
					.nombre_lignes; j++)
			{
				if ((((complex16 **) (*((struct_matrice *) (*s_objet).objet))
						.tableau)[j] = malloc(((size_t) (*((struct_matrice *)
						(*s_objet).objet)).nombre_lignes) * sizeof(complex16)))
						== NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}

				for(k = 0; k < (*((struct_matrice *) (*s_objet).objet))
						.nombre_colonnes; k++)
				{
					registre = vecteur_complexe[k];
					registre.partie_imaginaire =
							-vecteur_complexe[k].partie_imaginaire;

					f77multiplicationcc_(&(vecteur_complexe[j]),
							&registre, &(((complex16 **)
							(*((struct_matrice *) (*s_objet).objet)).tableau)
							[j][k]));
				}
			}

			free(vecteur_complexe);

			if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
					s_objet) == d_erreur)
			{
				return;
			}

			instruction_multiplication(s_etat_processus);

			if (((*s_etat_processus).erreur_systeme != d_es) ||
					((*s_etat_processus).erreur_execution != d_ex) ||
					((*s_etat_processus).exception != d_ep))
			{
				liberation(s_etat_processus, s_copie_argument);
				liberation(s_etat_processus, s_matrice_identite);
				free(tau);

				if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
				{
					return;
				}

				(*s_etat_processus).l_base_pile_last = registre_pile_last;
				return;
			}

			instruction_moins(s_etat_processus);

			if (((*s_etat_processus).erreur_systeme != d_es) ||
					((*s_etat_processus).erreur_execution != d_ex) ||
					((*s_etat_processus).exception != d_ep))
			{
				liberation(s_etat_processus, s_copie_argument);
				liberation(s_etat_processus, s_matrice_identite);
				free(tau);

				if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
				{
					return;
				}

				(*s_etat_processus).l_base_pile_last = registre_pile_last;
				return;
			}

			if (i > 0)
			{
				instruction_multiplication(s_etat_processus);

				if (((*s_etat_processus).erreur_systeme != d_es) ||
						((*s_etat_processus).erreur_execution != d_ex) ||
						((*s_etat_processus).exception != d_ep))
				{
					liberation(s_etat_processus, s_copie_argument);
					liberation(s_etat_processus, s_matrice_identite);
					free(tau);

					if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
					{
						return;
					}

					(*s_etat_processus).l_base_pile_last = registre_pile_last;
					return;
				}
			}
		}

		if (test_cfsf(s_etat_processus, 31) == d_vrai)
		{
			if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
			{
				return;
			}

			(*s_etat_processus).l_base_pile_last = registre_pile_last;
		}

		// Matrice R

		for(i = 0; i < (*((struct_matrice *) (*s_objet_resultat).objet))
				.nombre_lignes; i++)
		{
			for(j = 0; j < i; j++)
			{
				((complex16 **) (*((struct_matrice *) (*s_objet_resultat)
						.objet)).tableau)[i][j].partie_reelle = 0;
				((complex16 **) (*((struct_matrice *) (*s_objet_resultat)
						.objet)).tableau)[i][j].partie_imaginaire = 0;
			}
		}

		if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
				s_objet_resultat) == d_erreur)
		{
			return;
		}

		liberation(s_etat_processus, s_matrice_identite);
		liberation(s_etat_processus, s_copie_argument);
		free(tau);
	}

	/*
	 * Type d'argument invalide
	 */

	else
	{
		liberation(s_etat_processus, s_objet_argument);

		(*s_etat_processus).erreur_execution = d_ex_erreur_type_argument;
		return;
	}

	liberation(s_etat_processus, s_objet_argument);

	return;
}

// vim: ts=4