/* ================================================================================ RPL/2 (R) version 4.1.32 Copyright (C) 1989-2020 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 l'inversion d'une matrice carrée ================================================================================ Entrées : struct_matrice -------------------------------------------------------------------------------- Sorties : inverse 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. -------------------------------------------------------------------------------- Effets de bord : néant ================================================================================ */ void inversion_matrice(struct_processus *s_etat_processus, struct_matrice *s_matrice) { real8 *work; integer4 dim_matrice; integer4 dim_work; integer4 erreur; integer4 *pivot; integer8 i; integer8 j; integer8 k; integer8 rang_estime; integer8 taille_matrice_f77; struct_complexe16 *c_work; void *matrice_f77; rang(s_etat_processus, s_matrice, &rang_estime); if ((*s_etat_processus).erreur_systeme != d_es) { return; } if (((*s_etat_processus).erreur_execution != d_ex) || ((*s_etat_processus).exception != d_ep)) { return; } if (rang_estime < (integer8) (*s_matrice).nombre_lignes) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; return; } taille_matrice_f77 = (*s_matrice).nombre_lignes * (*s_matrice).nombre_colonnes; dim_matrice = (integer4) (*s_matrice).nombre_lignes; 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]; } } for(i = 0; i < (*s_matrice).nombre_lignes; i++) { free(((integer8 **) (*s_matrice).tableau)[i]); } free((integer8 **) (*s_matrice).tableau); if (((*s_matrice).tableau = malloc(((size_t) (*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(((size_t) (*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(((size_t) (*s_matrice) .nombre_lignes) * sizeof(integer4))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((work = (real8 *) malloc(sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } dim_work = -1; dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work, &dim_work, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(pivot); free(work); free(matrice_f77); return; } dim_work = (integer4) work[0]; free(work); if ((work = (real8 *) malloc(((size_t) dim_work) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } dgetrf_(&dim_matrice, &dim_matrice, matrice_f77, &dim_matrice, pivot, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(pivot); free(work); free(matrice_f77); return; } dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work, &dim_work, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(pivot); free(work); free(matrice_f77); return; } free(work); free(pivot); 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); 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) (*s_matrice) .nombre_lignes) * sizeof(integer4))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((work = (real8 *) malloc(sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } dim_work = -1; dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work, &dim_work, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(pivot); free(work); free(matrice_f77); return; } dim_work = (integer4) work[0]; free(work); if ((work = (real8 *) malloc(((size_t) dim_work) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } dgetrf_(&dim_matrice, &dim_matrice, matrice_f77, &dim_matrice, pivot, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(pivot); free(work); free(matrice_f77); return; } dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work, &dim_work, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(pivot); free(work); free(matrice_f77); return; } free(work); free(pivot); for(k = 0, i = 0; i < (*s_matrice).nombre_lignes; i++) { for(j = 0; j < (*s_matrice).nombre_colonnes; j++) { ((real8 **) (*s_matrice).tableau)[j][i] = ((real8 *) matrice_f77)[k++]; } } free(matrice_f77); break; } case 'C' : { if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) * sizeof(struct_complexe16))) == 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++) { (((struct_complexe16 *) matrice_f77)[k]).partie_reelle = (((struct_complexe16 **) (*s_matrice).tableau) [j][i]).partie_reelle; (((struct_complexe16 *) matrice_f77)[k]).partie_imaginaire = (((struct_complexe16 **) (*s_matrice).tableau) [j][i]).partie_imaginaire; k++; } } if ((pivot = (integer4 *) malloc(((size_t) (*s_matrice) .nombre_lignes) * sizeof(integer4))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((c_work = (struct_complexe16 *) malloc( sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } dim_work = -1; zgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, c_work, &dim_work, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(pivot); free(c_work); free(matrice_f77); return; } dim_work = (integer4) c_work[0].partie_reelle; free(c_work); if ((c_work = (struct_complexe16 *) malloc(((size_t) dim_work) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } zgetrf_(&dim_matrice, &dim_matrice, matrice_f77, &dim_matrice, pivot, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(pivot); free(c_work); free(matrice_f77); return; } zgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, c_work, &dim_work, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_matrice_non_inversible; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(pivot); free(c_work); free(matrice_f77); return; } free(c_work); free(pivot); for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice).nombre_lignes; j++) { (((struct_complexe16 **) (*s_matrice).tableau) [j][i]).partie_reelle = (((struct_complexe16 *) matrice_f77)[k]).partie_reelle; (((struct_complexe16 **) (*s_matrice).tableau) [j][i]).partie_imaginaire = (((struct_complexe16 *) matrice_f77)[k]).partie_imaginaire; k++; } } free(matrice_f77); break; } } return; } /* ================================================================================ Fonction calculant les vecteurs propres d'une matrice carrée ================================================================================ Entrées : struct_matrice, pointeur sur la matrice d'entrée, ainsi que struct_matrice, pointeur sur le vecteur des valeurs propres (vecteur de complexes) et sur deux matrices complexes contenant les différents vecteurs propres gauches et droits. Les pointeurs sur les vecteurs propres peuvent être nuls, et seules les valeurs propres sont calculées. L'allocation des tableaux de sortie est effectuée dans la routine, mais les structures s_matrice et s_vecteurs doivent être passées en paramètre. La matrice présente en entrée n'est pas libérée. -------------------------------------------------------------------------------- Sorties : vecteur contenant les valeurs propres, matrice contenant les vecteurs propres gauches et matrice contenant les vecteurs propres droits. Toutes les sorties sont complexes. -------------------------------------------------------------------------------- Effets de bord : néant ================================================================================ */ void valeurs_propres(struct_processus *s_etat_processus, struct_matrice *s_matrice, struct_vecteur *s_valeurs_propres, struct_matrice *s_vecteurs_propres_gauches, struct_matrice *s_vecteurs_propres_droits) { real8 *rwork; integer4 dim_matrice; integer4 lwork; integer4 erreur; integer8 i; integer8 j; integer8 k; integer8 taille_matrice_f77; struct_complexe16 *matrice_f77; struct_complexe16 *vpd_f77; struct_complexe16 *vpg_f77; struct_complexe16 *work; unsigned char calcul_vp_droits; unsigned char calcul_vp_gauches; unsigned char negatif; taille_matrice_f77 = (*s_matrice).nombre_lignes * (*s_matrice).nombre_colonnes; dim_matrice = (integer4) (*s_matrice).nombre_lignes; /* * Allocation de la matrice complexe */ if ((matrice_f77 = (struct_complexe16 *) malloc(((size_t) taille_matrice_f77) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } /* * Garniture de la matrice de compatibilité Fortran */ switch((*s_matrice).type) { /* * La matrice en entrée est une matrice d'entiers. */ case 'I' : { for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice).nombre_lignes; j++) { ((struct_complexe16 *) matrice_f77)[k] .partie_reelle = (real8) ((integer8 **) (*s_matrice).tableau)[j][i]; ((struct_complexe16 *) matrice_f77)[k++] .partie_imaginaire = (real8) 0; } } break; } /* * La matrice en entrée est une matrice de réels. */ case 'R' : { for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice).nombre_lignes; j++) { ((struct_complexe16 *) matrice_f77)[k] .partie_reelle = ((real8 **) (*s_matrice).tableau)[j][i]; ((struct_complexe16 *) matrice_f77)[k++] .partie_imaginaire = (real8) 0; } } break; } /* * La matrice en entrée est une matrice de complexes. */ case 'C' : { for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice).nombre_lignes; j++) { ((struct_complexe16 *) matrice_f77)[k] .partie_reelle = ((struct_complexe16 **) (*s_matrice).tableau)[j][i].partie_reelle; ((struct_complexe16 *) matrice_f77)[k++] .partie_imaginaire = ((struct_complexe16 **) (*s_matrice).tableau)[j][i].partie_imaginaire; } } break; } } (*s_valeurs_propres).type = 'C'; (*s_valeurs_propres).taille = (*s_matrice).nombre_lignes; if (((*s_valeurs_propres).tableau = (struct_complexe16 *) malloc(((size_t) (*s_valeurs_propres).taille) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if (s_vecteurs_propres_gauches != NULL) { (*s_vecteurs_propres_gauches).type = 'C'; calcul_vp_gauches = 'V'; if ((vpg_f77 = (struct_complexe16 *) malloc(((size_t) taille_matrice_f77) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } else { vpg_f77 = NULL; calcul_vp_gauches = 'N'; } if (s_vecteurs_propres_droits != NULL) { (*s_vecteurs_propres_droits).type = 'C'; calcul_vp_droits = 'V'; if ((vpd_f77 = (struct_complexe16 *) malloc(((size_t) taille_matrice_f77) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } else { vpd_f77 = NULL; calcul_vp_droits = 'N'; } negatif = 'N'; lwork = -1; if ((rwork = (real8 *) malloc(2 * ((size_t) (*s_matrice).nombre_lignes) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((work = (struct_complexe16 *) malloc(sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } zgeev_(&negatif, &negatif, &dim_matrice, matrice_f77, &dim_matrice, (struct_complexe16 *) (*s_valeurs_propres).tableau, vpg_f77, &dim_matrice, vpd_f77, &dim_matrice, work, &lwork, rwork, &erreur, 1, 1); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_decomposition_QR; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(work); free(rwork); free(matrice_f77); if (calcul_vp_gauches == 'V') { free(vpg_f77); } if (calcul_vp_droits == 'V') { free(vpd_f77); } return; } lwork = (integer4) work[0].partie_reelle; free(work); if ((work = (struct_complexe16 *) malloc(((size_t) lwork) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } zgeev_(&calcul_vp_gauches, &calcul_vp_droits, &dim_matrice, matrice_f77, &dim_matrice, (struct_complexe16 *) (*s_valeurs_propres).tableau, vpg_f77, &dim_matrice, vpd_f77, &dim_matrice, work, &lwork, rwork, &erreur, 1, 1); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_decomposition_QR; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(work); free(rwork); free(matrice_f77); if (calcul_vp_gauches == 'V') { free(vpg_f77); } if (calcul_vp_droits == 'V') { free(vpd_f77); } return; } free(work); free(rwork); if (calcul_vp_gauches == 'V') { (*s_vecteurs_propres_gauches).type = 'C'; (*s_vecteurs_propres_gauches).nombre_lignes = (*s_matrice).nombre_lignes; (*s_vecteurs_propres_gauches).nombre_colonnes = (*s_matrice).nombre_colonnes; if (((*s_vecteurs_propres_gauches).tableau = malloc( ((size_t) (*s_vecteurs_propres_gauches).nombre_lignes) * sizeof(struct_complexe16 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } for(i = 0; i < (*s_vecteurs_propres_gauches).nombre_lignes; i++) { if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches) .tableau)[i] = (struct_complexe16 *) malloc(((size_t) (*s_vecteurs_propres_gauches).nombre_colonnes) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } for(k = 0, i = 0; i < (*s_vecteurs_propres_gauches).nombre_colonnes; i++) { for(j = 0; j < (*s_vecteurs_propres_gauches).nombre_lignes; j++) { ((struct_complexe16 **) (*s_vecteurs_propres_gauches) .tableau)[j][i].partie_reelle = vpg_f77[k].partie_reelle; ((struct_complexe16 **) (*s_vecteurs_propres_gauches) .tableau)[j][i].partie_imaginaire = vpg_f77[k++].partie_imaginaire; } } free(vpg_f77); } if (calcul_vp_droits == 'V') { (*s_vecteurs_propres_droits).type = 'C'; (*s_vecteurs_propres_droits).nombre_lignes = (*s_matrice).nombre_lignes; (*s_vecteurs_propres_droits).nombre_colonnes = (*s_matrice).nombre_colonnes; if (((*s_vecteurs_propres_droits).tableau = malloc( ((size_t) (*s_vecteurs_propres_droits).nombre_lignes) * sizeof(struct_complexe16 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } for(i = 0; i < (*s_vecteurs_propres_droits).nombre_lignes; i++) { if ((((struct_complexe16 **) (*s_vecteurs_propres_droits) .tableau)[i] = (struct_complexe16 *) malloc(((size_t) (*s_vecteurs_propres_droits).nombre_colonnes) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } for(k = 0, i = 0; i < (*s_vecteurs_propres_droits).nombre_colonnes; i++) { for(j = 0; j < (*s_vecteurs_propres_droits).nombre_lignes; j++) { ((struct_complexe16 **) (*s_vecteurs_propres_droits) .tableau)[j][i].partie_reelle = vpd_f77[k].partie_reelle; ((struct_complexe16 **) (*s_vecteurs_propres_droits) .tableau)[j][i].partie_imaginaire = vpd_f77[k++].partie_imaginaire; } } free(vpd_f77); } free(matrice_f77); return; } /* ================================================================================ Fonction calculant les vecteurs propres généralisés d'une matrice carrée dans une métrique. ================================================================================ Entrées : struct_matrice, pointeur sur la matrice d'entrée, ainsi que struct_matrice, pointeur sur la métrique et struct_matrice, pointeur sur le vecteur des valeurs propres (vecteur de complexes) et sur deux matrices complexes contenant les différents vecteurs propres gauches et droits. Les pointeurs sur les vecteurs propres peuvent être nuls, et seules les valeurs propres sont calculées. L'allocation des tableaux de sortie est effectuée dans la routine, mais les structures s_matrice et s_vecteurs doivent être passées en paramètre. La matrice présente en entrée n'est pas libérée. -------------------------------------------------------------------------------- Sorties : vecteur contenant les valeurs propres, matrice contenant les vecteurs propres gauches et matrice contenant les vecteurs propres droits. Toutes les sorties sont complexes. -------------------------------------------------------------------------------- Effets de bord : néant ================================================================================ */ void valeurs_propres_generalisees(struct_processus *s_etat_processus, struct_matrice *s_matrice, struct_matrice *s_metrique, struct_vecteur *s_valeurs_propres, struct_matrice *s_vecteurs_propres_gauches, struct_matrice *s_vecteurs_propres_droits) { real8 *rwork; integer4 dim_matrice; integer4 lwork; integer4 erreur; integer8 i; integer8 j; integer8 k; integer8 taille_matrice_f77; struct_complexe16 *alpha; struct_complexe16 *beta; struct_complexe16 *matrice_f77; struct_complexe16 *metrique_f77; struct_complexe16 *vpd_f77; struct_complexe16 *vpg_f77; struct_complexe16 *work; unsigned char calcul_vp_droits; unsigned char calcul_vp_gauches; unsigned char negatif; taille_matrice_f77 = (*s_matrice).nombre_lignes * (*s_matrice).nombre_colonnes; dim_matrice = (integer4) (*s_matrice).nombre_lignes; /* * Allocation de la matrice complexe */ if ((matrice_f77 = (struct_complexe16 *) malloc(((size_t) taille_matrice_f77) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } /* * Garniture de la matrice de compatibilité Fortran */ switch((*s_matrice).type) { /* * La matrice en entrée est une matrice d'entiers. */ case 'I' : { for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice).nombre_lignes; j++) { ((struct_complexe16 *) matrice_f77)[k] .partie_reelle = (real8) ((integer8 **) (*s_matrice).tableau)[j][i]; ((struct_complexe16 *) matrice_f77)[k++] .partie_imaginaire = (real8) 0; } } break; } /* * La matrice en entrée est une matrice de réels. */ case 'R' : { for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice).nombre_lignes; j++) { ((struct_complexe16 *) matrice_f77)[k] .partie_reelle = ((real8 **) (*s_matrice).tableau)[j][i]; ((struct_complexe16 *) matrice_f77)[k++] .partie_imaginaire = (real8) 0; } } break; } /* * La matrice en entrée est une matrice de complexes. */ case 'C' : { for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice).nombre_lignes; j++) { ((struct_complexe16 *) matrice_f77)[k] .partie_reelle = ((struct_complexe16 **) (*s_matrice).tableau)[j][i].partie_reelle; ((struct_complexe16 *) matrice_f77)[k++] .partie_imaginaire = ((struct_complexe16 **) (*s_matrice).tableau)[j][i].partie_imaginaire; } } break; } } /* * Allocation de la metrique complexe */ if ((metrique_f77 = (struct_complexe16 *) malloc(((size_t) taille_matrice_f77) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } /* * Garniture de la matrice de compatibilité Fortran */ switch((*s_metrique).type) { /* * La matrice en entrée est une matrice d'entiers. */ case 'I' : { for(k = 0, i = 0; i < (*s_metrique).nombre_colonnes; i++) { for(j = 0; j < (*s_metrique).nombre_lignes; j++) { ((struct_complexe16 *) metrique_f77)[k] .partie_reelle = (real8) ((integer8 **) (*s_metrique).tableau)[j][i]; ((struct_complexe16 *) metrique_f77)[k++] .partie_imaginaire = (real8) 0; } } break; } /* * La matrice en entrée est une matrice de réels. */ case 'R' : { for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice).nombre_lignes; j++) { ((struct_complexe16 *) metrique_f77)[k] .partie_reelle = ((real8 **) (*s_metrique).tableau)[j][i]; ((struct_complexe16 *) metrique_f77)[k++] .partie_imaginaire = (real8) 0; } } break; } /* * La matrice en entrée est une matrice de complexes. */ case 'C' : { for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice).nombre_lignes; j++) { ((struct_complexe16 *) metrique_f77)[k] .partie_reelle = ((struct_complexe16 **) (*s_metrique).tableau)[j][i].partie_reelle; ((struct_complexe16 *) metrique_f77)[k++] .partie_imaginaire = ((struct_complexe16 **) (*s_metrique).tableau)[j][i].partie_imaginaire; } } break; } } (*s_valeurs_propres).type = 'C'; (*s_valeurs_propres).taille = (*s_matrice).nombre_lignes; if (((*s_valeurs_propres).tableau = (struct_complexe16 *) malloc(((size_t) (*s_valeurs_propres).taille) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if (s_vecteurs_propres_gauches != NULL) { (*s_vecteurs_propres_gauches).type = 'C'; calcul_vp_gauches = 'V'; if ((vpg_f77 = (struct_complexe16 *) malloc(((size_t) taille_matrice_f77) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } else { vpg_f77 = NULL; calcul_vp_gauches = 'N'; } if (s_vecteurs_propres_droits != NULL) { (*s_vecteurs_propres_droits).type = 'C'; calcul_vp_droits = 'V'; if ((vpd_f77 = (struct_complexe16 *) malloc(((size_t) taille_matrice_f77) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } else { vpd_f77 = NULL; calcul_vp_droits = 'N'; } negatif = 'N'; lwork = -1; if ((rwork = (real8 *) malloc(8 * ((size_t) (*s_matrice).nombre_lignes) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((work = (struct_complexe16 *) malloc(sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((alpha = (struct_complexe16 *) malloc(((size_t) (*s_valeurs_propres) .taille) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((beta = (struct_complexe16 *) malloc(((size_t) (*s_valeurs_propres) .taille) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } zggev_(&negatif, &negatif, &dim_matrice, matrice_f77, &dim_matrice, metrique_f77, &dim_matrice, alpha, beta, vpg_f77, &dim_matrice, vpd_f77, &dim_matrice, work, &lwork, rwork, &erreur, 1, 1); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_decomposition_QZ; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(work); free(rwork); free(matrice_f77); free(metrique_f77); if (calcul_vp_gauches == 'V') { free(vpg_f77); (*s_vecteurs_propres_gauches).type = 'C'; (*s_vecteurs_propres_gauches).nombre_lignes = 1; (*s_vecteurs_propres_gauches).nombre_colonnes = 1; if (((*s_vecteurs_propres_gauches).tableau = malloc(((size_t) (*s_vecteurs_propres_gauches).nombre_lignes) * sizeof(struct_complexe16 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches) .tableau)[0] = (struct_complexe16 *) malloc(((size_t) (*s_vecteurs_propres_gauches).nombre_colonnes) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } if (calcul_vp_droits == 'V') { free(vpd_f77); (*s_vecteurs_propres_droits).type = 'C'; (*s_vecteurs_propres_droits).nombre_lignes = 1; (*s_vecteurs_propres_droits).nombre_colonnes = 1; if (((*s_vecteurs_propres_droits).tableau = malloc(((size_t) (*s_vecteurs_propres_droits).nombre_lignes) * sizeof(struct_complexe16 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((((struct_complexe16 **) (*s_vecteurs_propres_droits) .tableau)[0] = (struct_complexe16 *) malloc(((size_t) (*s_vecteurs_propres_droits).nombre_colonnes) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } return; } lwork = (integer4) work[0].partie_reelle; free(work); if ((work = (struct_complexe16 *) malloc(((size_t) lwork) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } zggev_(&calcul_vp_gauches, &calcul_vp_droits, &dim_matrice, matrice_f77, &dim_matrice, metrique_f77, &dim_matrice, alpha, beta, vpg_f77, &dim_matrice, vpd_f77, &dim_matrice, work, &lwork, rwork, &erreur, 1, 1); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_decomposition_QZ; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(work); free(rwork); free(matrice_f77); free(metrique_f77); if (calcul_vp_gauches == 'V') { free(vpg_f77); (*s_vecteurs_propres_gauches).type = 'C'; (*s_vecteurs_propres_gauches).nombre_lignes = 1; (*s_vecteurs_propres_gauches).nombre_colonnes = 1; if (((*s_vecteurs_propres_gauches).tableau = malloc(((size_t) (*s_vecteurs_propres_gauches).nombre_lignes) * sizeof(struct_complexe16 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches) .tableau)[0] = (struct_complexe16 *) malloc(((size_t) (*s_vecteurs_propres_gauches).nombre_colonnes) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } if (calcul_vp_droits == 'V') { free(vpd_f77); (*s_vecteurs_propres_droits).type = 'C'; (*s_vecteurs_propres_droits).nombre_lignes = 1; (*s_vecteurs_propres_droits).nombre_colonnes = 1; if (((*s_vecteurs_propres_droits).tableau = malloc(((size_t) (*s_vecteurs_propres_droits).nombre_lignes) * sizeof(struct_complexe16 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((((struct_complexe16 **) (*s_vecteurs_propres_droits) .tableau)[0] = (struct_complexe16 *) malloc(((size_t) (*s_vecteurs_propres_droits).nombre_colonnes) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } return; } for(i = 0; i < (*s_valeurs_propres).taille; i++) { if ((beta[i].partie_reelle != 0) || (beta[i].partie_imaginaire != 0)) { f77divisioncc_(&(alpha[i]), &(beta[i]), &(((struct_complexe16 *) (*s_valeurs_propres).tableau)[i])); } else { ((struct_complexe16 *) (*s_valeurs_propres).tableau)[i] .partie_reelle = 0; ((struct_complexe16 *) (*s_valeurs_propres).tableau)[i] .partie_imaginaire = 0; (*s_etat_processus).exception = d_ep_division_par_zero; } } free(alpha); free(beta); free(work); free(rwork); if (calcul_vp_gauches == 'V') { (*s_vecteurs_propres_gauches).type = 'C'; (*s_vecteurs_propres_gauches).nombre_lignes = (*s_matrice).nombre_lignes; (*s_vecteurs_propres_gauches).nombre_colonnes = (*s_matrice).nombre_colonnes; if (((*s_vecteurs_propres_gauches).tableau = malloc(((size_t) (*s_vecteurs_propres_gauches).nombre_lignes) * sizeof(struct_complexe16 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } for(i = 0; i < (*s_vecteurs_propres_gauches).nombre_lignes; i++) { if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches) .tableau)[i] = (struct_complexe16 *) malloc(((size_t) (*s_vecteurs_propres_gauches).nombre_colonnes) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } for(k = 0, i = 0; i < (*s_vecteurs_propres_gauches).nombre_colonnes; i++) { for(j = 0; j < (*s_vecteurs_propres_gauches).nombre_lignes; j++) { ((struct_complexe16 **) (*s_vecteurs_propres_gauches) .tableau)[j][i].partie_reelle = vpg_f77[k].partie_reelle; ((struct_complexe16 **) (*s_vecteurs_propres_gauches) .tableau)[j][i].partie_imaginaire = vpg_f77[k++].partie_imaginaire; } } free(vpg_f77); } if (calcul_vp_droits == 'V') { (*s_vecteurs_propres_droits).type = 'C'; (*s_vecteurs_propres_droits).nombre_lignes = (*s_matrice).nombre_lignes; (*s_vecteurs_propres_droits).nombre_colonnes = (*s_matrice).nombre_colonnes; if (((*s_vecteurs_propres_droits).tableau = malloc(((size_t) (*s_vecteurs_propres_droits).nombre_lignes) * sizeof(struct_complexe16 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } for(i = 0; i < (*s_vecteurs_propres_droits).nombre_lignes; i++) { if ((((struct_complexe16 **) (*s_vecteurs_propres_droits) .tableau)[i] = (struct_complexe16 *) malloc(((size_t) (*s_vecteurs_propres_droits).nombre_colonnes) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } for(k = 0, i = 0; i < (*s_vecteurs_propres_droits).nombre_colonnes; i++) { for(j = 0; j < (*s_vecteurs_propres_droits).nombre_lignes; j++) { ((struct_complexe16 **) (*s_vecteurs_propres_droits) .tableau)[j][i].partie_reelle = vpd_f77[k].partie_reelle; ((struct_complexe16 **) (*s_vecteurs_propres_droits) .tableau)[j][i].partie_imaginaire = vpd_f77[k++].partie_imaginaire; } } free(vpd_f77); } free(matrice_f77); free(metrique_f77); return; } /* ================================================================================ Fonction réalisation le calcul d'un moindres carrés (minimise [B-AX]^2) ================================================================================ Entrées : struct_matrice -------------------------------------------------------------------------------- Sorties : coefficients dans une struct_matrice allouée au vol -------------------------------------------------------------------------------- Effets de bord : néant ================================================================================ */ void moindres_carres(struct_processus *s_etat_processus, struct_matrice *s_matrice_a, struct_matrice *s_matrice_b, struct_matrice *s_matrice_x) { real8 *work; real8 rcond; real8 *rwork; real8 *vecteur_s; integer4 erreur; integer4 *iwork; integer4 lrwork; integer4 lwork; integer4 nlvl; integer4 rank; integer4 registre_1; integer4 registre_2; integer4 registre_3; integer4 registre_4; integer4 registre_5; integer4 smlsiz; complex16 *cwork; integer8 i; integer8 j; integer8 k; integer8 taille_matrice_a_f77; integer8 taille_matrice_b_f77; integer8 taille_matrice_x_f77; void *matrice_a_f77; void *matrice_b_f77; taille_matrice_a_f77 = (*s_matrice_a).nombre_lignes * (*s_matrice_a).nombre_colonnes; taille_matrice_b_f77 = (*s_matrice_b).nombre_lignes * (*s_matrice_b).nombre_colonnes; taille_matrice_x_f77 = (*s_matrice_b).nombre_colonnes * (*s_matrice_a).nombre_colonnes; /* * Résultat réel */ if (((*s_matrice_a).type != 'C') && ((*s_matrice_b).type != 'C')) { /* * Garniture de la matrice A */ if ((matrice_a_f77 = malloc(((size_t) taille_matrice_a_f77) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((*s_matrice_a).type == 'I') { for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_a).nombre_lignes; j++) { ((real8 *) matrice_a_f77)[k++] = (real8) ((integer8 **) (*s_matrice_a).tableau)[j][i]; } } } else { for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_a).nombre_lignes; j++) { ((real8 *) matrice_a_f77)[k++] = ((real8 **) (*s_matrice_a).tableau)[j][i]; } } } /* * Garniture de la matrice B */ if ((matrice_b_f77 = malloc(((size_t) ((taille_matrice_b_f77 < taille_matrice_x_f77) ? taille_matrice_x_f77 : taille_matrice_b_f77)) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((*s_matrice_b).type == 'I') { for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_b).nombre_lignes; j++) { ((real8 *) matrice_b_f77)[k++] = (real8) ((integer8 **) (*s_matrice_b).tableau)[j][i]; } for(; j < (*s_matrice_a).nombre_lignes; j++, k++); } } else { for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_b).nombre_lignes; j++) { ((real8 *) matrice_b_f77)[k++] = ((real8 **) (*s_matrice_b).tableau)[j][i]; } for(; j < (*s_matrice_a).nombre_lignes; j++, k++); } } rcond = -1; registre_1 = 9; registre_2 = 0; smlsiz = ilaenv_(®istre_1, "DGELSD", " ", ®istre_2, ®istre_2, ®istre_2, ®istre_2, 6, 1); nlvl = 1 + ((integer4) (log(((real8) (((*s_matrice_a).nombre_lignes < (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes : (*s_matrice_a).nombre_colonnes)) / (smlsiz + 1)) / log((real8) 2))); if ((iwork = (integer4 *) malloc(((size_t) ((((*s_matrice_a) .nombre_lignes < (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes : (*s_matrice_a).nombre_colonnes) * (11 + (3 * nlvl)))) * sizeof(integer4))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } registre_1 = (integer4) (*s_matrice_a).nombre_lignes; registre_2 = (integer4) (*s_matrice_a).nombre_colonnes; registre_3 = (integer4) (*s_matrice_b).nombre_colonnes; registre_4 = (registre_1 > registre_2) ? registre_1 : registre_2; registre_5 = (registre_1 > registre_2) ? registre_2 : registre_1; if ((work = malloc(sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((vecteur_s = (real8 *) malloc(((size_t) registre_5) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } lwork = -1; dgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77, ®istre_1, matrice_b_f77, ®istre_4, vecteur_s, &rcond, &rank, work, &lwork, iwork, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_decomposition_SVD; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(work); free(iwork); free(matrice_a_f77); free(matrice_b_f77); return; } lwork = (integer4) work[0]; free(work); if ((work = malloc(((size_t) lwork) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } dgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77, ®istre_1, matrice_b_f77, ®istre_4, vecteur_s, &rcond, &rank, work, &lwork, iwork, &erreur); free(vecteur_s); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_decomposition_SVD; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(work); free(iwork); free(matrice_a_f77); free(matrice_b_f77); return; } free(work); free(iwork); (*s_matrice_x).nombre_lignes = (*s_matrice_a).nombre_colonnes; (*s_matrice_x).nombre_colonnes = (*s_matrice_b).nombre_colonnes; if (((*s_matrice_x).tableau = malloc(((size_t) (*s_matrice_x) .nombre_lignes) * sizeof(real8 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } for(j = 0; j < (*s_matrice_x).nombre_lignes; j++) { if ((((real8 **) (*s_matrice_x).tableau)[j] = (real8 *) malloc(((size_t) (*s_matrice_x).nombre_colonnes) * sizeof(real8)))== NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } for(k = 0, i = 0; i < (*s_matrice_x).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_x).nombre_lignes; j++) { (((real8 **) (*s_matrice_x).tableau)[j][i]) = (((real8 *) matrice_b_f77)[k++]); } for(; j < (*s_matrice_b).nombre_lignes; j++, k++); } } /* * Résultat complexe */ else { /* * Garniture de la matrice A */ if ((matrice_a_f77 = malloc(((size_t) taille_matrice_a_f77) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((*s_matrice_a).type == 'I') { for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_a).nombre_lignes; j++) { ((struct_complexe16 *) matrice_a_f77)[k].partie_reelle = (real8) ((integer8 **) (*s_matrice_a) .tableau)[j][i]; ((struct_complexe16 *) matrice_a_f77)[k++] .partie_imaginaire = (real8) 0; } } } else if ((*s_matrice_a).type == 'R') { for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_a).nombre_lignes; j++) { ((struct_complexe16 *) matrice_a_f77)[k].partie_reelle = ((real8 **) (*s_matrice_a).tableau)[j][i]; ((struct_complexe16 *) matrice_a_f77)[k++] .partie_imaginaire = (real8) 0; } } } else { for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_a).nombre_lignes; j++) { ((struct_complexe16 *) matrice_a_f77)[k].partie_reelle = ((struct_complexe16 **) (*s_matrice_a).tableau) [j][i].partie_reelle; ((struct_complexe16 *) matrice_a_f77)[k++] .partie_imaginaire = ((struct_complexe16 **) (*s_matrice_a).tableau)[j][i].partie_imaginaire; } } } /* * Garniture de la matrice B */ if ((matrice_b_f77 = malloc(((size_t) ((taille_matrice_b_f77 < taille_matrice_x_f77) ? taille_matrice_x_f77 : taille_matrice_b_f77)) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((*s_matrice_b).type == 'I') { for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_b).nombre_lignes; j++) { ((struct_complexe16 *) matrice_b_f77)[k].partie_reelle = (real8) ((integer8 **) (*s_matrice_b) .tableau)[j][i]; ((struct_complexe16 *) matrice_b_f77)[k++] .partie_imaginaire = (real8) 0; } for(; j < (*s_matrice_a).nombre_lignes; j++, k++); } } else if ((*s_matrice_b).type == 'R') { for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_b).nombre_lignes; j++) { ((struct_complexe16 *) matrice_b_f77)[k].partie_reelle = ((real8 **) (*s_matrice_b).tableau)[j][i]; ((struct_complexe16 *) matrice_b_f77)[k++] .partie_imaginaire = (real8) 0; } for(; j < (*s_matrice_a).nombre_lignes; j++, k++); } } else { for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_b).nombre_lignes; j++) { ((struct_complexe16 *) matrice_b_f77)[k].partie_reelle = ((struct_complexe16 **) (*s_matrice_b).tableau) [j][i].partie_reelle; ((struct_complexe16 *) matrice_b_f77)[k++] .partie_imaginaire = ((struct_complexe16 **) (*s_matrice_b).tableau)[j][i].partie_imaginaire; } for(; j < (*s_matrice_a).nombre_lignes; j++, k++); } } rcond = -1; registre_1 = 9; registre_2 = 0; smlsiz = ilaenv_(®istre_1, "ZGELSD", " ", ®istre_2, ®istre_2, ®istre_2, ®istre_2, 6, 1); nlvl = 1 + ((integer4) (log(((real8) (((*s_matrice_a).nombre_lignes < (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes : (*s_matrice_a).nombre_colonnes)) / (smlsiz + 1)) / log((real8) 2))); if ((*s_matrice_a).nombre_lignes >= (*s_matrice_a).nombre_colonnes) { lrwork = (integer4) ((*s_matrice_a).nombre_colonnes * (8 + (2 * smlsiz) + (8 * nlvl) + (*s_matrice_b).nombre_colonnes)); } else { lrwork = (integer4) ((*s_matrice_a).nombre_lignes * (8 + (2 * smlsiz) + (8 * nlvl) + (*s_matrice_b).nombre_colonnes)); } if ((rwork = (real8 *) malloc(((size_t) lrwork) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((iwork = (integer4 *) malloc(((size_t) ((((*s_matrice_a) .nombre_lignes < (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes : (*s_matrice_a).nombre_colonnes) * (11 + (3 * nlvl)))) * sizeof(integer4))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } registre_1 = (integer4) (*s_matrice_a).nombre_lignes; registre_2 = (integer4) (*s_matrice_a).nombre_colonnes; registre_3 = (integer4) (*s_matrice_b).nombre_colonnes; registre_4 = (registre_1 > registre_2) ? registre_1 : registre_2; registre_5 = (registre_1 > registre_2) ? registre_2 : registre_1; if ((cwork = malloc(sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } if ((vecteur_s = (real8 *) malloc(((size_t) registre_5) * sizeof(real8))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } lwork = -1; zgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77, ®istre_1, matrice_b_f77, ®istre_4, vecteur_s, &rcond, &rank, cwork, &lwork, rwork, iwork, &erreur); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_decomposition_SVD; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(cwork); free(iwork); free(rwork); free(matrice_a_f77); free(matrice_b_f77); return; } lwork = (integer4) cwork[0].partie_reelle; free(cwork); if ((cwork = malloc(((size_t) lwork) * sizeof(struct_complexe16))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } zgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77, ®istre_1, matrice_b_f77, ®istre_4, vecteur_s, &rcond, &rank, cwork, &lwork, rwork, iwork, &erreur); free(vecteur_s); if (erreur != 0) { if (erreur > 0) { (*s_etat_processus).exception = d_ep_decomposition_SVD; } else { (*s_etat_processus).erreur_execution = d_ex_routines_mathematiques; } free(cwork); free(iwork); free(rwork); free(matrice_a_f77); free(matrice_b_f77); return; } free(cwork); free(iwork); free(rwork); (*s_matrice_x).nombre_lignes = (*s_matrice_a).nombre_colonnes; (*s_matrice_x).nombre_colonnes = (*s_matrice_b).nombre_colonnes; if (((*s_matrice_x).tableau = malloc(((size_t) (*s_matrice_x) .nombre_lignes) * sizeof(struct_complexe16 *))) == NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } for(j = 0; j < (*s_matrice_x).nombre_lignes; j++) { if ((((struct_complexe16 **) (*s_matrice_x).tableau)[j] = (struct_complexe16 *) malloc(((size_t) (*s_matrice_x).nombre_colonnes) * sizeof(struct_complexe16)))== NULL) { (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; return; } } for(k = 0, i = 0; i < (*s_matrice_x).nombre_colonnes; i++) { for(j = 0; j < (*s_matrice_x).nombre_lignes; j++) { (((struct_complexe16 **) (*s_matrice_x).tableau)[j][i]) .partie_reelle = (((struct_complexe16 *) matrice_b_f77)[k]).partie_reelle; (((struct_complexe16 **) (*s_matrice_x).tableau)[j][i]) .partie_imaginaire = (((struct_complexe16 *) matrice_b_f77)[k++]).partie_imaginaire; } for(; j < (*s_matrice_b).nombre_lignes; j++, k++); } } free(matrice_a_f77); free(matrice_b_f77); return; } // vim: ts=4