/*
================================================================================
  RPL/2 (R) version 4.1.26
  Copyright (C) 1989-2017 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 réalisation la factorisation de Schur d'une matrice carrée
================================================================================
  Entrées : struct_matrice
--------------------------------------------------------------------------------
  Sorties : décomposition de Schur de la matrice d'entrée et drapeau d'erreur.
			La matrice en entrée est écrasée. La matrice de sortie est
			la forme de Schur.
			La routine renvoie aussi une matrice de complexes correspondant
			aux vecteurs de Schur. Cette matrice est allouée par
			la routine et vaut NULL sinon.
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

void
factorisation_schur(struct_processus *s_etat_processus,
		struct_matrice *s_matrice, struct_matrice **s_schur)
{
	complex16					*w;

	integer4					info;
	integer4					lwork;
	integer4					nombre_colonnes_a;
	integer4					nombre_lignes_a;
	integer4					sdim;

	real8						*rwork;
	real8						*wi;
	real8						*wr;

	unsigned char				calcul_vecteurs_schur;
	unsigned char				tri_vecteurs_schur;

	integer8					i;
	integer8					j;
	integer8					k;
	integer8					taille_matrice_f77;

	void						*matrice_a_f77;
	void						*matrice_vs_f77;
	void						*tampon;
	void						*work;

	nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
	nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
	taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;

	calcul_vecteurs_schur = 'V';
	tri_vecteurs_schur = 'N';

	switch((*s_matrice).type)
	{
		case 'I' :
		{
			/* Conversion de la matrice en matrice réelle */

			for(i = 0; i < nombre_lignes_a; i++)
			{
				tampon = (*s_matrice).tableau[i];

				if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *)
						malloc(((size_t) nombre_colonnes_a) * sizeof(real8)))
						== NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}

				for(j = 0; j < nombre_colonnes_a; j++)
				{
					((real8 **) (*s_matrice).tableau)[i][j] = (real8)
							((integer8 *) tampon)[j];
				}

				free(tampon);
			}

			(*s_matrice).type = 'R';
		}

		case 'R' :
		{
			if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((matrice_vs_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((real8 *) matrice_a_f77)[k++] = ((real8 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			if ((wr = (real8 *) malloc(((size_t) nombre_lignes_a) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((wi = (real8 *) malloc(((size_t) nombre_lignes_a) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			lwork = -1;

			if ((work = (real8 *) malloc(sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			dgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur,
					NULL, &nombre_lignes_a, matrice_a_f77,
					&nombre_colonnes_a, &sdim, wr, wi,
					matrice_vs_f77, &nombre_colonnes_a,
					work, &lwork, NULL, &info, 1, 1);

			lwork = (integer4) ((real8 *) work)[0];
			free(work);

			if ((work = (real8 *) malloc(((size_t) lwork) * sizeof(real8)))
					== NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			dgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur,
					NULL, &nombre_lignes_a, matrice_a_f77,
					&nombre_colonnes_a, &sdim, wr, wi,
					matrice_vs_f77, &nombre_colonnes_a,
					work, &lwork, NULL, &info, 1, 1);

			free(work);
			free(wr);
			free(wi);

			if (info != 0)
			{
				if (info > 0)
				{
					(*s_etat_processus).exception = d_ep_decomposition_QR;
				}
				else
				{
					(*s_etat_processus).erreur_execution =
							d_ex_routines_mathematiques;
				}

				free(matrice_a_f77);
				free(matrice_vs_f77);
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((real8 **) (*s_matrice).tableau)[j][i] =
							((real8 *) matrice_a_f77)[k++];
				}
			}

			(**s_schur).nombre_colonnes = (*s_matrice).nombre_colonnes;
			(**s_schur).nombre_lignes = (*s_matrice).nombre_lignes;
			(**s_schur).type = 'R';

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

			for(i = 0; i < (**s_schur).nombre_lignes; i++)
			{
				if ((((real8 **) (**s_schur).tableau)[i] = (real8 *)
						malloc(((size_t) (**s_schur).nombre_colonnes) *
						sizeof(real8))) == NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}
			}

			for(k = 0, i = 0; i < (**s_schur).nombre_colonnes; i++)
			{
				for(j = 0; j < (**s_schur).nombre_lignes; j++)
				{
					((real8 **) (**s_schur).tableau)[j][i] = ((real8 *)
							matrice_vs_f77)[k++];
				}
			}

			free(matrice_a_f77);
			free(matrice_vs_f77);

			break;
		}

		case 'C' :
		{
			if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((matrice_vs_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((complex16 *) matrice_a_f77)[k].partie_reelle =
							((complex16 **) (*s_matrice).tableau)[j][i]
							.partie_reelle;
					((complex16 *) matrice_a_f77)[k++].partie_imaginaire =
							((complex16 **) (*s_matrice).tableau)[j][i]
							.partie_imaginaire;
				}
			}

			if ((w = (complex16 *) malloc(((size_t) nombre_lignes_a) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			lwork = -1;

			if ((work = (complex16 *) malloc(sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((rwork = (real8 *) malloc(((size_t) nombre_lignes_a) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			zgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur,
					NULL, &nombre_lignes_a, matrice_a_f77,
					&nombre_colonnes_a, &sdim, w,
					matrice_vs_f77, &nombre_colonnes_a,
					work, &lwork, rwork, NULL, &info, 1, 1);

			lwork = (integer4) ((complex16 *) work)[0].partie_reelle;
			free(work);

			if ((work = (complex16 *) malloc(((size_t) lwork) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			zgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur,
					NULL, &nombre_lignes_a, matrice_a_f77,
					&nombre_colonnes_a, &sdim, w,
					matrice_vs_f77, &nombre_colonnes_a,
					work, &lwork, rwork, NULL, &info, 1, 1);

			free(work);
			free(rwork);
			free(w);

			if (info != 0)
			{
				if (info > 0)
				{
					(*s_etat_processus).exception = d_ep_decomposition_QR;
				}
				else
				{
					(*s_etat_processus).erreur_execution =
							d_ex_routines_mathematiques;
				}

				free(matrice_a_f77);
				free(matrice_vs_f77);
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((complex16 **) (*s_matrice).tableau)[j][i]
					   		.partie_reelle = ((complex16 *) matrice_a_f77)[k]
							.partie_reelle;
					((complex16 **) (*s_matrice).tableau)[j][i]
					   		.partie_imaginaire = ((complex16 *) matrice_a_f77)
							[k++].partie_imaginaire;
				}
			}

			(**s_schur).nombre_colonnes = (*s_matrice).nombre_colonnes;
			(**s_schur).nombre_lignes = (*s_matrice).nombre_lignes;
			(**s_schur).type = 'C';

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

			for(i = 0; i < (**s_schur).nombre_lignes; i++)
			{
				if ((((complex16 **) (**s_schur).tableau)[i] = (complex16 *)
						malloc(((size_t) (**s_schur).nombre_colonnes) *
						sizeof(complex16))) == NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}
			}

			for(k = 0, i = 0; i < (**s_schur).nombre_colonnes; i++)
			{
				for(j = 0; j < (**s_schur).nombre_lignes; j++)
				{
					((complex16 **) (**s_schur).tableau)[j][i].partie_reelle =
							((complex16 *) matrice_vs_f77)[k].partie_reelle;
					((complex16 **) (**s_schur).tableau)[j][i]
							.partie_imaginaire = ((complex16 *) matrice_vs_f77)
							[k++].partie_imaginaire;
				}
			}

			free(matrice_a_f77);
			free(matrice_vs_f77);

			break;
		}
	}

	return;
}


/*
================================================================================
  Fonction réalisation la factorisation LQ d'une matrice quelconque
================================================================================
  Entrées : struct_matrice
--------------------------------------------------------------------------------
  Sorties : décomposition de LQ de la matrice d'entrée. Le tableau tau
  est initialisé par la fonction
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

void
factorisation_lq(struct_processus *s_etat_processus, struct_matrice *s_matrice,
		void **tau)
{
	integer4					nombre_colonnes_a;
	integer4					nombre_lignes_a;
	integer4					erreur;

	integer8					i;
	integer8					j;
	integer8					k;
	integer8					taille_matrice_f77;

	void						*matrice_a_f77;
	void						*tampon;
	void						*work;

	nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
	nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
	taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;

	switch((*s_matrice).type)
	{
		case 'I' :
		{
			/* Conversion de la matrice en matrice réelle */

			for(i = 0; i < nombre_lignes_a; i++)
			{
				tampon = (*s_matrice).tableau[i];

				if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *)
						malloc(((size_t) nombre_colonnes_a) * sizeof(real8)))
						== NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}

				for(j = 0; j < nombre_colonnes_a; j++)
				{
					((real8 **) (*s_matrice).tableau)[i][j] = (real8)
							((integer8 *) tampon)[j];
				}

				free(tampon);
			}

			(*s_matrice).type = 'R';
		}

		case 'R' :
		{
			if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if (((*tau) = malloc(((size_t) ((nombre_colonnes_a <
					nombre_lignes_a) ? nombre_colonnes_a : nombre_lignes_a)) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((work = malloc(((size_t) nombre_lignes_a) * sizeof(real8)))
					== NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((real8 *) matrice_a_f77)[k++] = ((real8 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			dgelq2_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
					&nombre_lignes_a, (*((real8 **) tau)), work, &erreur);

			if (erreur != 0)
			{
				// L'erreur ne peut être que négative.

				(*s_etat_processus).erreur_execution =
						d_ex_routines_mathematiques;
				free(work);
				free(matrice_a_f77);
				return;
			}

			for(k = 0, i = 0; i < nombre_colonnes_a; i++)
			{
				for(j = 0; j < nombre_lignes_a; j++)
				{
					((real8 **) (*s_matrice).tableau)[j][i] = ((real8 *)
							matrice_a_f77)[k++];
				}
			}

			free(work);
			free(matrice_a_f77);

			break;
		}

		case 'C' :
		{
			if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if (((*tau) = malloc(((size_t) ((nombre_colonnes_a <
					nombre_lignes_a) ? nombre_colonnes_a : nombre_lignes_a)) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((work = malloc(((size_t) nombre_lignes_a) * sizeof(complex16)))
					== NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((complex16 *) matrice_a_f77)[k].partie_reelle =
							((complex16 **) (*s_matrice).tableau)[j][i]
							.partie_reelle;
					((complex16 *) matrice_a_f77)[k++].partie_imaginaire =
							((complex16 **) (*s_matrice).tableau)[j][i]
							.partie_imaginaire;
				}
			}

			zgelq2_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
					&nombre_lignes_a, (*((complex16 **) tau)), work, &erreur);

			if (erreur != 0)
			{
				// L'erreur ne peut être que négative.

				(*s_etat_processus).erreur_execution =
						d_ex_routines_mathematiques;
				free(work);
				free(matrice_a_f77);
				return;
			}

			for(k = 0, i = 0; i < nombre_colonnes_a; i++)
			{
				for(j = 0; j < nombre_lignes_a; j++)
				{
					((complex16 **) (*s_matrice).tableau)[j][i].partie_reelle =
							((complex16 *) matrice_a_f77)[k].partie_reelle;
					((complex16 **) (*s_matrice).tableau)[j][i]
							.partie_imaginaire = ((complex16 *) matrice_a_f77)
							[k++].partie_imaginaire;
				}
			}

			free(work);
			free(matrice_a_f77);

			break;
		}
	}

	return;
}


/*
================================================================================
  Fonction réalisation la factorisation QR d'une matrice quelconque
================================================================================
  Entrées : struct_matrice
--------------------------------------------------------------------------------
  Sorties : décomposition de QR de la matrice d'entrée. Le tableau tau
  est initialisé par la fonction
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

void
factorisation_qr(struct_processus *s_etat_processus, struct_matrice *s_matrice,
		void **tau)
{
	integer4					lwork;
	integer4					nombre_colonnes_a;
	integer4					nombre_lignes_a;
	integer4					erreur;
	integer4					*pivot;

	real8						*rwork;

	integer8					i;
	integer8					j;
	integer8					k;
	integer8					taille_matrice_f77;

	void						*matrice_a_f77;
	void						*registre;
	void						*tampon;
	void						*work;

	nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
	nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
	taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;

	switch((*s_matrice).type)
	{
		case 'I' :
		{
			/* Conversion de la matrice en matrice réelle */

			for(i = 0; i < nombre_lignes_a; i++)
			{
				tampon = (*s_matrice).tableau[i];

				if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *)
						malloc(((size_t) nombre_colonnes_a) * sizeof(real8)))
						== NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}

				for(j = 0; j < nombre_colonnes_a; j++)
				{
					((real8 **) (*s_matrice).tableau)[i][j] = (real8)
							((integer8 *) tampon)[j];
				}

				free(tampon);
			}

			(*s_matrice).type = 'R';
		}

		case 'R' :
		{
			if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if (((*tau) = malloc(((size_t) ((nombre_colonnes_a <
					nombre_lignes_a) ? nombre_colonnes_a : nombre_lignes_a)) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((real8 *) matrice_a_f77)[k++] = ((real8 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			if ((pivot = malloc(((size_t) nombre_colonnes_a) *
					sizeof(integer4))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(i = 0; i < nombre_colonnes_a; pivot[i++] = 0);

			lwork = -1;

			if ((work = malloc(sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			// Calcul de la taille de l'espace

			dgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
					&nombre_lignes_a, pivot, (*((real8 **) tau)),
					work, &lwork, &erreur);

			lwork = (integer4) ((real8 *) work)[0];

			free(work);

			if ((work = malloc(((size_t) lwork) * sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			// Calcul de la permutation

			if ((registre = malloc(((size_t) taille_matrice_f77) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			memcpy(registre, matrice_a_f77, ((size_t) taille_matrice_f77) *
					sizeof(real8));

			dgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, registre,
					&nombre_lignes_a, pivot, (*((real8 **) tau)),
					work, &lwork, &erreur);

			free(registre);

			if (erreur != 0)
			{
				// L'erreur ne peut être que négative.

				(*s_etat_processus).erreur_execution =
						d_ex_routines_mathematiques;
				free(work);
				free(matrice_a_f77);
				free(tau);
				return;
			}

			// La permutation doit maintenant être unitaire

			dgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
					&nombre_lignes_a, pivot, (*((real8 **) tau)),
					work, &lwork, &erreur);

			if (erreur != 0)
			{
				// L'erreur ne peut être que négative.

				(*s_etat_processus).erreur_execution =
						d_ex_routines_mathematiques;
				free(work);
				free(matrice_a_f77);
				free(tau);
				return;
			}

			for(i = 0; i < nombre_colonnes_a; i++)
			{
				if ((i + 1) != pivot[i])
				{
					(*s_etat_processus).erreur_execution =
							d_ex_routines_mathematiques;

					free(pivot);
					free(work);
					free(matrice_a_f77);
					free(tau);

					return;
				}
			}

			free(pivot);

			for(k = 0, i = 0; i < nombre_colonnes_a; i++)
			{
				for(j = 0; j < nombre_lignes_a; j++)
				{
					((real8 **) (*s_matrice).tableau)[j][i] = ((real8 *)
							matrice_a_f77)[k++];
				}
			}

			free(work);
			free(matrice_a_f77);

			break;
		}

		case 'C' :
		{
			if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if (((*tau) = malloc(((size_t) ((nombre_colonnes_a <
					nombre_lignes_a) ? nombre_colonnes_a : nombre_lignes_a)) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((complex16 *) matrice_a_f77)[k].partie_reelle =
							((complex16 **) (*s_matrice).tableau)[j][i]
							.partie_reelle;
					((complex16 *) matrice_a_f77)[k++].partie_imaginaire =
							((complex16 **) (*s_matrice).tableau)[j][i]
							.partie_imaginaire;
				}
			}

			if ((pivot = malloc(((size_t) nombre_colonnes_a) *
					sizeof(integer4))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((rwork = malloc(2 * ((size_t) nombre_colonnes_a) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(i = 0; i < nombre_colonnes_a; pivot[i++] = 0);

			lwork = -1;

			if ((work = malloc(sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			// Calcul de la taille de l'espace

			zgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
					&nombre_lignes_a, pivot, (*((complex16 **) tau)),
					work, &lwork, rwork, &erreur);

			if (erreur != 0)
			{
				// L'erreur ne peut être que négative.

				(*s_etat_processus).erreur_execution =
						d_ex_routines_mathematiques;

				free(work);
				free(rwork);
				free(pivot);
				free(matrice_a_f77);
				free(tau);
				return;
			}

			lwork = (integer4) ((complex16 *) work)[0].partie_reelle;

			free(work);

			if ((work = malloc(((size_t) lwork) * sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			// Calcul de la permutation

			if ((registre = malloc(((size_t) taille_matrice_f77) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			memcpy(registre, matrice_a_f77,
					((size_t) taille_matrice_f77) * sizeof(complex16));

			zgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, registre,
					&nombre_lignes_a, pivot, (*((complex16 **) tau)),
					work, &lwork, rwork, &erreur);

			free(registre);

			if (erreur != 0)
			{
				// L'erreur ne peut être que négative.

				(*s_etat_processus).erreur_execution =
						d_ex_routines_mathematiques;

				free(work);
				free(rwork);
				free(pivot);
				free(matrice_a_f77);
				free(tau);
				return;
			}

			// La permutation doit maintenant être unitaire

			zgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
					&nombre_lignes_a, pivot, (*((complex16 **) tau)),
					work, &lwork, rwork, &erreur);

			if (erreur != 0)
			{
				// L'erreur ne peut être que négative.

				(*s_etat_processus).erreur_execution =
						d_ex_routines_mathematiques;

				free(work);
				free(rwork);
				free(pivot);
				free(matrice_a_f77);
				free(tau);
				return;
			}

			free(rwork);

			for(i = 0; i < nombre_colonnes_a; i++)
			{
				if ((i + 1) != pivot[i])
				{
					(*s_etat_processus).erreur_execution =
							d_ex_routines_mathematiques;

					free(pivot);
					free(work);
					free(matrice_a_f77);
					free(tau);

					return;
				}
			}

			free(pivot);

			for(k = 0, i = 0; i < nombre_colonnes_a; i++)
			{
				for(j = 0; j < nombre_lignes_a; j++)
				{
					((complex16 **) (*s_matrice).tableau)[j][i].partie_reelle =
							((complex16 *) matrice_a_f77)[k].partie_reelle;
					((complex16 **) (*s_matrice).tableau)[j][i]
							.partie_imaginaire = ((complex16 *)
							matrice_a_f77)[k++].partie_imaginaire;
				}
			}

			free(work);
			free(matrice_a_f77);

			break;
		}
	}

	return;
}


/*
================================================================================
  Fonctions calculant le déterminant ou le rang d'une matrice quelconque
================================================================================
  Entrées : struct_matrice
--------------------------------------------------------------------------------
  Sorties : déterminant
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/


static integer4
calcul_rang(struct_processus *s_etat_processus, void *matrice_f77,
		integer4 nombre_lignes_a, integer4 nombre_colonnes_a,
		integer4 *pivot, integer4 dimension_vecteur_pivot, unsigned char type)
{
	integer4					erreur;
	integer4					*iwork;
	integer4					*jpvt;
	integer4					lwork;
	integer4					longueur;
	integer4					nombre_colonnes_b;
	integer4					nombre_lignes_b;
	integer4					ordre;
	integer4					rang;

	real8						anorme;
	real8						rcond;
	real8						*rwork;

	unsigned char				norme;

	integer8					i;

	void						*matrice_b;
	void						*matrice_c;
	void						*work;

#undef NORME_I
#ifdef NORME_I
	norme = 'I';

	if ((work = malloc(((size_t) nombre_lignes_a) * sizeof(real8))) == NULL)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return(-1);
	}
#else
	norme = '1';
	work = NULL;
#endif

	longueur = 1;

	if (type == 'R')
	{
		anorme = dlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a,
				matrice_f77, &nombre_lignes_a, work, longueur);

#ifndef NORME_I
		free(work);
#endif

		if ((matrice_c = malloc(((size_t) (nombre_lignes_a * nombre_colonnes_a))
				* sizeof(real8))) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		memcpy(matrice_c, matrice_f77, ((size_t) (nombre_lignes_a *
				nombre_colonnes_a)) * sizeof(real8));

		dgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77,
				&nombre_lignes_a, pivot, &erreur);

		if (erreur < 0)
		{
			(*s_etat_processus).erreur_execution =
					d_ex_routines_mathematiques;

			free(matrice_f77);
			return(-1);
		}

		if ((iwork = malloc(((size_t) nombre_colonnes_a) * sizeof(integer4)))
				== NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		if ((work = malloc(4 * ((size_t) nombre_colonnes_a) * sizeof(real8)))
				== NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		ordre = (nombre_lignes_a > nombre_colonnes_a)
				? nombre_colonnes_a : nombre_lignes_a;

		dgecon_(&norme, &ordre, matrice_f77,
				&nombre_lignes_a, &anorme, &rcond, work, iwork, &erreur,
				longueur);

		free(work);
		free(iwork);

		if (erreur < 0)
		{
			(*s_etat_processus).erreur_execution =
					d_ex_routines_mathematiques;

			free(matrice_f77);
			return(-1);
		}

		if ((jpvt = malloc(((size_t) nombre_colonnes_a) * sizeof(integer4)))
				== NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		for(i = 0; i < nombre_colonnes_a; i++)
		{
			((integer4 *) jpvt)[i] = 0;
		}

		lwork = -1;

		if ((work = malloc(sizeof(real8))) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		nombre_colonnes_b = 1;
		nombre_lignes_b = (nombre_lignes_a > nombre_colonnes_a)
				? nombre_lignes_a : nombre_colonnes_a;

		if ((matrice_b = malloc(((size_t) nombre_lignes_b) * sizeof(real8)))
				== NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		for(i = 0; i < nombre_lignes_b; i++)
		{
			((real8 *) matrice_b)[i] = 0;
		}

		dgelsy_(&nombre_lignes_a, &nombre_colonnes_a,
				&nombre_colonnes_b, matrice_c, &nombre_lignes_a,
				matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang,
				work, &lwork, &erreur);

		lwork = (integer4) ((real8 *) work)[0];
		free(work);

		if ((work = malloc(((size_t) lwork) * sizeof(real8))) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		dgelsy_(&nombre_lignes_a, &nombre_colonnes_a,
				&nombre_colonnes_b, matrice_c, &nombre_lignes_a,
				matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang,
				work, &lwork, &erreur);

		free(matrice_b);
		free(matrice_c);
		free(work);
		free(jpvt);

		if (erreur < 0)
		{
			(*s_etat_processus).erreur_execution =
					d_ex_routines_mathematiques;

			free(matrice_f77);
			return(-1);
		}
	}
	else
	{
		anorme = zlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a,
				matrice_f77, &nombre_lignes_a, work, longueur);

#ifndef NORME_I
		free(work);
#endif

		if ((matrice_c = malloc(((size_t) (nombre_lignes_a * nombre_colonnes_a))
				* sizeof(complex16))) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		memcpy(matrice_c, matrice_f77, ((size_t) (nombre_lignes_a *
				nombre_colonnes_a)) * sizeof(complex16));

		zgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77,
				&nombre_lignes_a, pivot, &erreur);

		if (erreur < 0)
		{
			(*s_etat_processus).erreur_execution =
					d_ex_routines_mathematiques;

			free(matrice_f77);
			return(-1);
		}

		if ((rwork = malloc(2 * ((size_t) nombre_colonnes_a) * sizeof(real8)))
				== NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		if ((work = malloc(2 * ((size_t) nombre_colonnes_a) *
				sizeof(complex16))) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		ordre = (nombre_lignes_a > nombre_colonnes_a)
				? nombre_colonnes_a : nombre_lignes_a;

		zgecon_(&norme, &ordre, matrice_f77,
				&nombre_lignes_a, &anorme, &rcond, work, rwork, &erreur,
				longueur);

		free(work);

		if (erreur < 0)
		{
			(*s_etat_processus).erreur_execution =
					d_ex_routines_mathematiques;

			free(matrice_f77);
			return(-1);
		}

		if ((jpvt = malloc(((size_t) nombre_colonnes_a) * sizeof(integer4)))
				== NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		for(i = 0; i < nombre_colonnes_a; i++)
		{
			((integer4 *) jpvt)[i] = 0;
		}

		lwork = -1;

		if ((work = malloc(sizeof(complex16))) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		nombre_colonnes_b = 1;
		nombre_lignes_b = (nombre_lignes_a > nombre_colonnes_a)
				? nombre_lignes_a : nombre_colonnes_a;

		if ((matrice_b = malloc(((size_t) nombre_lignes_b) *
				sizeof(complex16))) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		for(i = 0; i < nombre_lignes_b; i++)
		{
			((complex16 *) matrice_b)[i].partie_reelle = 0;
			((complex16 *) matrice_b)[i].partie_imaginaire = 0;
		}

		zgelsy_(&nombre_lignes_a, &nombre_colonnes_a,
				&nombre_colonnes_b, matrice_c, &nombre_lignes_a,
				matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang,
				work, &lwork, rwork, &erreur);

		lwork = (integer4) ((complex16 *) work)[0].partie_reelle;
		free(work);

		if ((work = malloc(((size_t) lwork) * sizeof(complex16))) == NULL)
		{
			(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
			return(-1);
		}

		zgelsy_(&nombre_lignes_a, &nombre_colonnes_a,
				&nombre_colonnes_b, matrice_c, &nombre_lignes_a,
				matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang,
				work, &lwork, rwork, &erreur);

		free(rwork);
		free(matrice_b);
		free(matrice_c);
		free(work);
		free(jpvt);

		if (erreur < 0)
		{
			(*s_etat_processus).erreur_execution =
					d_ex_routines_mathematiques;

			free(matrice_f77);
			return(-1);
		}
	}

	return(rang);
}


void
determinant(struct_processus *s_etat_processus, struct_matrice *s_matrice,
		void *valeur)
{
	complex16					*vecteur_complexe;

	integer4					dimension_vecteur_pivot;
	integer4					nombre_colonnes_a;
	integer4					nombre_lignes_a;
	integer4					*pivot;
	integer4					rang;

	integer8					i;
	integer8					j;
	integer8					k;
	integer8					signe;
	integer8					taille_matrice_f77;

	real8						*vecteur_reel;

	void						*matrice_f77;

	nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
	nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
	dimension_vecteur_pivot = (nombre_lignes_a < nombre_colonnes_a)
			? nombre_lignes_a : nombre_colonnes_a;
	taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;

	switch((*s_matrice).type)
	{
		case 'I' :
		{
			if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((real8 *) matrice_f77)[k++] = (real8) ((integer8 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			if ((pivot = (integer4 *) malloc(((size_t) dimension_vecteur_pivot)
				 	* sizeof(integer4))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((rang = calcul_rang(s_etat_processus, matrice_f77,
					nombre_lignes_a, nombre_colonnes_a, pivot,
					dimension_vecteur_pivot, 'R')) < 0)
			{
				return;
			}

			if (rang < nombre_lignes_a)
			{
				(*((real8 *) valeur)) = 0;
			}
			else
			{
				if ((vecteur_reel = malloc(((size_t) ((*s_matrice)
						.nombre_colonnes)) * sizeof(real8))) == NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}

				signe = 1;

				for(i = 0; i < (*s_matrice).nombre_colonnes; i++)
				{
					if (pivot[i] != (i + 1))
					{
						signe = (signe == 1) ? -1 : 1;
					}

					vecteur_reel[i] = ((real8 *) matrice_f77)
							[(i * nombre_colonnes_a) + i];
				}

				for(i = 1; i < (*s_matrice).nombre_colonnes; i++)
				{
					vecteur_reel[0] *= vecteur_reel[i];
				}

				(*((real8 *) valeur)) = vecteur_reel[0] * ((real8) signe);
				free(vecteur_reel);
			}

			free(matrice_f77);
			free(pivot);

			break;
		}

		case 'R' :
		{
			if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((real8 *) matrice_f77)[k++] = ((real8 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			if ((pivot = (integer4 *) malloc(((size_t) dimension_vecteur_pivot)
					* sizeof(integer4))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((rang = calcul_rang(s_etat_processus, matrice_f77,
					nombre_lignes_a, nombre_colonnes_a, pivot,
					dimension_vecteur_pivot, 'R')) < 0)
			{
				return;
			}

			if (rang < nombre_lignes_a)
			{
				(*((real8 *) valeur)) = 0;
			}
			else
			{
				if ((vecteur_reel = malloc(((size_t) (*s_matrice)
						.nombre_colonnes) * sizeof(real8))) == NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}

				signe = 1;

				for(i = 0; i < (*s_matrice).nombre_colonnes; i++)
				{
					if (pivot[i] != (i + 1))
					{
						signe = (signe == 1) ? -1 : 1;
					}

					vecteur_reel[i] = ((real8 *) matrice_f77)
							[(i * nombre_colonnes_a) + i];
				}

				for(i = 1; i < (*s_matrice).nombre_colonnes; i++)
				{
					vecteur_reel[0] *= vecteur_reel[i];
				}

				(*((real8 *) valeur)) = vecteur_reel[0] * ((real8) signe);

				free(vecteur_reel);
			}

			free(matrice_f77);
			free(pivot);

			break;
		}

		case 'C' :
		{
			if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((complex16 *) matrice_f77)[k++] = ((complex16 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			if ((pivot = (integer4 *) malloc(((size_t) dimension_vecteur_pivot)
				   	* sizeof(integer4))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			if ((rang = calcul_rang(s_etat_processus, matrice_f77,
					nombre_lignes_a, nombre_colonnes_a, pivot,
					dimension_vecteur_pivot, 'C')) < 0)
			{
				return;
			}

			if (rang < nombre_lignes_a)
			{
				(*((complex16 *) valeur)).partie_reelle = 0;
				(*((complex16 *) valeur)).partie_imaginaire = 0;
			}
			else
			{
				if ((vecteur_complexe = malloc(((size_t) (*s_matrice)
						.nombre_colonnes) * sizeof(complex16))) == NULL)
				{
					(*s_etat_processus).erreur_systeme =
							d_es_allocation_memoire;
					return;
				}

				signe = 1;

				for(i = 0; i < (*s_matrice).nombre_colonnes; i++)
				{
					if (pivot[i] != (i + 1))
					{
						signe = (signe == 1) ? -1 : 1;
					}

					vecteur_complexe[i] = ((complex16 *) matrice_f77)
							[(i * nombre_colonnes_a) + i];
				}

				for(i = 1; i < (*s_matrice).nombre_colonnes; i++)
				{
					f77multiplicationcc_(&(vecteur_complexe[0]),
							&(vecteur_complexe[i]), &(vecteur_complexe[0]));
				}

				f77multiplicationci_(&(vecteur_complexe[0]), &signe,
						((complex16 *) valeur));

				free(vecteur_complexe);
			}

			free(matrice_f77);
			free(pivot);

			break;
		}
	}

	return;
}


void
rang(struct_processus *s_etat_processus, struct_matrice *s_matrice,
		integer8 *valeur)
{
	integer4					dimension_vecteur_pivot;
	integer4					nombre_lignes_a;
	integer4					nombre_colonnes_a;
	integer4					*pivot;
	integer4					rang;
	integer4					taille_matrice_f77;

	integer8					i;
	integer8					j;
	integer8					k;

	void						*matrice_f77;

	nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
	nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
	dimension_vecteur_pivot = (nombre_lignes_a < nombre_colonnes_a)
			? nombre_lignes_a : nombre_colonnes_a;
	taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;

	if ((pivot = (integer4 *) malloc(((size_t) dimension_vecteur_pivot) *
			sizeof(integer4))) == NULL)
	{
		(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
		return;
	}

	switch((*s_matrice).type)
	{
		case 'I' :
		{
			if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((real8 *) matrice_f77)[k++] = (real8) ((integer8 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			if ((rang = calcul_rang(s_etat_processus, matrice_f77,
					nombre_lignes_a, nombre_colonnes_a, pivot,
					dimension_vecteur_pivot, 'R')) < 0)
			{
				free(pivot);
				free(matrice_f77);
				return;
			}

			free(matrice_f77);
			(*valeur) = rang;
			break;
		}

		case 'R' :
		{
			if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(real8))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((real8 *) matrice_f77)[k++] = ((real8 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			if ((rang = calcul_rang(s_etat_processus, matrice_f77,
					nombre_lignes_a, nombre_colonnes_a, pivot,
					dimension_vecteur_pivot, 'R')) < 0)
			{
				free(pivot);
				free(matrice_f77);
				return;
			}

			free(matrice_f77);
			(*valeur) = rang;
			break;
		}

		case 'C' :
		{
			if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
					sizeof(complex16))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
			{
				for(j = 0; j < (*s_matrice).nombre_lignes; j++)
				{
					((complex16 *) matrice_f77)[k++] = ((complex16 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			if ((rang = calcul_rang(s_etat_processus, matrice_f77,
					nombre_lignes_a, nombre_colonnes_a, pivot,
					dimension_vecteur_pivot, 'C')) < 0)
			{
				free(pivot);
				free(matrice_f77);
				return;
			}

			free(matrice_f77);
			(*valeur) = rang;
			break;
		}
	}

	free(pivot);

	return;
}

// vim: ts=4