/*
================================================================================
  RPL/2 (R) version 4.0.12
  Copyright (C) 1989-2010 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 LU d'une matrice quelconque
================================================================================
  Entrées : struct_matrice
--------------------------------------------------------------------------------
  Sorties : décomposition A=PLU de la matrice d'entrée et drapeau d'erreur.
  			La matrice en entrée est écrasée. La matrice de sortie est réelle
			si la matrice d'entrée est entière ou réelle, et complexe
			si la matrice d'entrée est complexe.
			La routine renvoie aussi une matrice d'entiers correspondant
			à la matrice de permutation. Cette matrice est allouée par
			la routine et vaut NULL sinon.
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

void
factorisation_lu(struct_processus *s_etat_processus,
		struct_matrice *s_matrice, struct_matrice **s_permutation)
{
	integer4					dimension_vecteur_pivot;
	integer4					erreur;
	integer4					nombre_colonnes_a;
	integer4					nombre_lignes_a;
	integer4					*pivot;

	long						n;

	unsigned long				i;
	unsigned long				j;
	unsigned long				k;
	unsigned long				taille_matrice_f77;

	void						*matrice_f77;
	void						*tampon;

	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 = (void *) malloc(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++] = ((integer8 **)
							(*s_matrice).tableau)[j][i];
				}
			}

			for(i = 0; i < (*s_matrice).nombre_lignes; i++)
			{
				free(((integer8 **) (*s_matrice).tableau)[i]);
			}

			free((integer8 **) (*s_matrice).tableau);

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

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

			(*s_matrice).type = 'R';

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

			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(pivot);
				free(matrice_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_f77)[k++];
				}
			}

			free(matrice_f77);

			(**s_permutation).type = 'I';
			(**s_permutation).nombre_lignes = dimension_vecteur_pivot;
			(**s_permutation).nombre_colonnes = dimension_vecteur_pivot;

			if (((**s_permutation).tableau = malloc((**s_permutation)
					.nombre_lignes * sizeof(integer8 *))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

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

				for(j = 0; j < (**s_permutation).nombre_colonnes; j++)
				{
					((integer8 **) (**s_permutation).tableau)[i][j] =
							(j == i) ? 1 : 0;
				}
			}

			for(n = dimension_vecteur_pivot - 1; n >= 0; n--)
			{
				if ((pivot[n] - 1) != n)
				{
					tampon = ((integer8 **) (**s_permutation).tableau)[n];
					((integer8 **) (**s_permutation).tableau)[n] =
							((integer8 **) (**s_permutation).tableau)
							[pivot[n] - 1];
					((integer8 **) (**s_permutation).tableau)[pivot[n] - 1] =
							tampon;
				}
			}

			free(pivot);

			break;
		}

		case 'R' :
		{
			if ((matrice_f77 = (void *) malloc(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(dimension_vecteur_pivot *
					sizeof(integer4))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			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(pivot);
				free(matrice_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_f77)[k++];
				}
			}

			free(matrice_f77);

			(**s_permutation).type = 'I';
			(**s_permutation).nombre_lignes = dimension_vecteur_pivot;
			(**s_permutation).nombre_colonnes = dimension_vecteur_pivot;

			if (((**s_permutation).tableau = malloc((**s_permutation)
					.nombre_lignes * sizeof(integer8 *))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

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

				for(j = 0; j < (**s_permutation).nombre_colonnes; j++)
				{
					((integer8 **) (**s_permutation).tableau)[i][j] =
							(j == i) ? 1 : 0;
				}
			}

			for(n = dimension_vecteur_pivot - 1; n >= 0; n--)
			{
				if ((pivot[n] - 1) != n)
				{
					tampon = ((integer8 **) (**s_permutation).tableau)[n];
					((integer8 **) (**s_permutation).tableau)[n] =
							((integer8 **) (**s_permutation).tableau)
							[pivot[n] - 1];
					((integer8 **) (**s_permutation).tableau)[pivot[n] - 1] =
							tampon;
				}
			}

			free(pivot);

			break;
		}

		case 'C' :
		{
			if ((matrice_f77 = (void *) malloc(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(dimension_vecteur_pivot *
					sizeof(integer4))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

			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(pivot);
				free(matrice_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] =
							((complex16 *) matrice_f77)[k++];
				}
			}

			free(matrice_f77);

			(**s_permutation).type = 'I';
			(**s_permutation).nombre_lignes = dimension_vecteur_pivot;
			(**s_permutation).nombre_colonnes = dimension_vecteur_pivot;

			if (((**s_permutation).tableau = malloc((**s_permutation)
					.nombre_lignes * sizeof(integer8 *))) == NULL)
			{
				(*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
				return;
			}

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

				for(j = 0; j < (**s_permutation).nombre_colonnes; j++)
				{
					((integer8 **) (**s_permutation).tableau)[i][j] =
							(j == i) ? 1 : 0;
				}
			}

			for(n = dimension_vecteur_pivot - 1; n >= 0; n--)
			{
				if ((pivot[n] - 1) != n)
				{
					tampon = ((integer8 **) (**s_permutation).tableau)[n];
					((integer8 **) (**s_permutation).tableau)[n] =
							((integer8 **) (**s_permutation).tableau)
							[pivot[n] - 1];
					((integer8 **) (**s_permutation).tableau)[pivot[n] - 1] =
							tampon;
				}
			}

			free(pivot);

			break;
		}
	}

	return;
}


/*
================================================================================
  Fonction réalisation la factorisation de Cholesky d'une matrice symétrique
  définie positive
================================================================================
  Entrées : struct_matrice, mode (U ou L selon la décomposition demandée)
--------------------------------------------------------------------------------
  Sorties : décomposition de Cholesky de la matrice d'entrée et drapeau
  			d'erreur.  La matrice en entrée est écrasée.
--------------------------------------------------------------------------------
  Effets de bord : néant
================================================================================
*/

void
factorisation_cholesky(struct_processus *s_etat_processus,
		struct_matrice *s_matrice, unsigned char mode)
{
	integer4					erreur;
	integer4					i;
	integer4					j;
	integer4					nombre_colonnes_a;
	integer4					nombre_lignes_a;

	unsigned long				taille_matrice_f77;

	void						*matrice_f77;

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

	if (nombre_lignes_a != nombre_colonnes_a)
	{
		(*s_etat_processus).erreur_execution = d_ex_dimensions_invalides;
		return;
	}

	taille_matrice_f77 = (nombre_lignes_a * (nombre_colonnes_a + 1)) / 2;

	switch((*s_matrice).type)
	{
		case 'I' :
		{
			/*
			 * Allocation du vecteur représentant la matrice triangulaire
			 */

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

			if (mode == 'U')
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = 1; i <= j; i++)
					{
						((real8 *) matrice_f77)[i + (((j - 1) * j) / 2) - 1] =
								(real8) ((integer8 **) (*s_matrice).tableau)
								[i - 1][j - 1];
					}
				}
			}
			else
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = j; i <= nombre_colonnes_a; i++)
					{
						((real8 *) matrice_f77)[(i - 1) + (((j - 1) *
								((nombre_colonnes_a * 2) - j)) / 2)] = (real8)
								((integer8 **) (*s_matrice).tableau)
								[i - 1][j - 1];
					}
				}
			}

			dpptrf_(&mode, &nombre_colonnes_a, matrice_f77, &erreur, 1);

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

				free(matrice_f77);
				return;
			}

			for(i = 0; i < (integer4) (*s_matrice).nombre_lignes; i++)
			{
				free(((integer8 **) (*s_matrice).tableau)[i]);
			}

			free((integer8 **) (*s_matrice).tableau);

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

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

			(*s_matrice).type = 'R';

			if (mode == 'U')
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = 1; i <= j; i++)
					{
						((real8 **) (*s_matrice).tableau)[i - 1][j - 1] =
								((real8 *) matrice_f77)
								[i + (((j - 1) * j) / 2) - 1];
					}
				}

				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = j + 1; i <= nombre_colonnes_a; i++)
					{
						((real8 **) (*s_matrice).tableau)[i - 1][j - 1] = 0;
					}
				}
			}
			else
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = j; i <= nombre_colonnes_a; i++)
					{
						((real8 **) (*s_matrice).tableau)[i - 1][j - 1] =
								((real8 *) matrice_f77)[(i - 1) + (((j - 1) *
								((nombre_colonnes_a * 2) - j)) / 2)];
					}
				}

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

			free(matrice_f77);
			break;
		}

		case 'R' :
		{
			/*
			 * Allocation du vecteur représentant la matrice triangulaire
			 */

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

			if (mode == 'U')
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = 1; i <= j; i++)
					{
						((real8 *) matrice_f77)[i + (((j - 1) * j) / 2) - 1] =
								((real8 **) (*s_matrice).tableau)[i - 1][j - 1];
					}
				}
			}
			else
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = j; i <= nombre_colonnes_a; i++)
					{
						((real8 *) matrice_f77)[(i - 1) + (((j - 1) *
								((nombre_colonnes_a * 2) - j)) / 2)] = (real8)
								((real8 **) (*s_matrice).tableau)[i - 1][j - 1];
					}
				}
			}

			dpptrf_(&mode, &nombre_colonnes_a, matrice_f77, &erreur, 1);

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

				free(matrice_f77);
				return;
			}

			if (mode == 'U')
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = 1; i <= j; i++)
					{
						((real8 **) (*s_matrice).tableau)[i - 1][j - 1] =
								((real8 *) matrice_f77)
								[i + (((j - 1) * j) / 2) - 1];
					}
				}

				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = j + 1; i <= nombre_colonnes_a; i++)
					{
						((real8 **) (*s_matrice).tableau)[i - 1][j - 1] = 0;
					}
				}
			}
			else
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = j; i <= nombre_colonnes_a; i++)
					{
						((real8 **) (*s_matrice).tableau)[i - 1][j - 1] =
								((real8 *) matrice_f77)[(i - 1) + (((j - 1) *
								((nombre_colonnes_a * 2) - j)) / 2)];
					}
				}

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

			free(matrice_f77);
			break;
		}

		case 'C' :
		{
			/*
			 * Allocation du vecteur représentant la matrice triangulaire
			 */

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

			if (mode == 'U')
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = 1; i <= j; i++)
					{
						((complex16 *) matrice_f77)[i + (((j - 1) * j) / 2)
						   		- 1].partie_reelle = ((complex16 **)
								(*s_matrice).tableau)[i - 1][j - 1]
								.partie_reelle;
						((complex16 *) matrice_f77)[i + (((j - 1) * j) / 2)
						   		- 1].partie_imaginaire = ((complex16 **)
								(*s_matrice).tableau)[i - 1][j - 1]
								.partie_imaginaire;
					}
				}
			}
			else
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = j; i <= nombre_colonnes_a; i++)
					{
						((complex16 *) matrice_f77)[(i - 1) + (((j - 1) *
								((nombre_colonnes_a * 2) - j)) / 2)]
						   		.partie_reelle = ((complex16 **)
								(*s_matrice).tableau)[i - 1][j - 1]
								.partie_reelle;
						((complex16 *) matrice_f77)[(i - 1) + (((j - 1) *
								((nombre_colonnes_a * 2) - j)) / 2)]
						   		.partie_imaginaire = ((complex16 **)
								(*s_matrice).tableau)[i - 1][j - 1]
								.partie_imaginaire;
					}
				}
			}

			zpptrf_(&mode, &nombre_colonnes_a, matrice_f77, &erreur, 1);

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

				free(matrice_f77);
				return;
			}

			if (mode == 'U')
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = 1; i <= j; i++)
					{
						((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
								.partie_reelle = ((complex16 *) matrice_f77)
								[i + (((j - 1) * j) / 2) - 1].partie_reelle;
						((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
								.partie_imaginaire = ((complex16 *) matrice_f77)
								[i + (((j - 1) * j) / 2) - 1].partie_imaginaire;
					}
				}

				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = j + 1; i <= nombre_colonnes_a; i++)
					{
						((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
								.partie_reelle = 0;
						((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
								.partie_imaginaire = 0;
					}
				}
			}
			else
			{
				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = j; i <= nombre_colonnes_a; i++)
					{
						((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
								.partie_reelle = ((complex16 *)
								matrice_f77)[(i - 1) + (((j - 1) *
								((nombre_colonnes_a * 2) - j)) / 2)]
								.partie_reelle;
						((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
								.partie_imaginaire = ((complex16 *)
								matrice_f77)[(i - 1) + (((j - 1) *
								((nombre_colonnes_a * 2) - j)) / 2)]
								.partie_imaginaire;
					}
				}

				for(j = 1; j <= nombre_colonnes_a; j++)
				{
					for(i = 1; i < j; i++)
					{
						((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
								.partie_reelle = 0;
						((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
								.partie_imaginaire = 0;
					}
				}
			}

			free(matrice_f77);
			break;
		}
	}

	return;
}

// vim: ts=4