Annotation of rpl/src/algebre_lineaire1.c, revision 1.1
1.1 ! bertrand 1: /*
! 2: ================================================================================
! 3: RPL/2 (R) version 4.0.9
! 4: Copyright (C) 1989-2010 Dr. BERTRAND Joël
! 5:
! 6: This file is part of RPL/2.
! 7:
! 8: RPL/2 is free software; you can redistribute it and/or modify it
! 9: under the terms of the CeCILL V2 License as published by the french
! 10: CEA, CNRS and INRIA.
! 11:
! 12: RPL/2 is distributed in the hope that it will be useful, but WITHOUT
! 13: ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! 14: FITNESS FOR A PARTICULAR PURPOSE. See the CeCILL V2 License
! 15: for more details.
! 16:
! 17: You should have received a copy of the CeCILL License
! 18: along with RPL/2. If not, write to info@cecill.info.
! 19: ================================================================================
! 20: */
! 21:
! 22:
! 23: #include "rpl.conv.h"
! 24:
! 25:
! 26: /*
! 27: ================================================================================
! 28: Fonction réalisation l'inversion d'une matrice carrée
! 29: ================================================================================
! 30: Entrées : struct_matrice
! 31: --------------------------------------------------------------------------------
! 32: Sorties : inverse de la matrice d'entrée et drapeau d'erreur. La matrice
! 33: en entrée est écrasée. La matrice de sortie est réelle si
! 34: la matrice d'entrée est entière ou réelle, et complexe
! 35: si la matrice d'entrée est complexe.
! 36: --------------------------------------------------------------------------------
! 37: Effets de bord : néant
! 38: ================================================================================
! 39: */
! 40:
! 41: void
! 42: inversion_matrice(struct_processus *s_etat_processus,
! 43: struct_matrice *s_matrice)
! 44: {
! 45: real8 *work;
! 46:
! 47: integer4 dim_matrice;
! 48: integer4 dim_work;
! 49: integer4 erreur;
! 50: integer4 *pivot;
! 51:
! 52: integer8 rang_estime;
! 53:
! 54: struct_complexe16 *c_work;
! 55:
! 56: unsigned long i;
! 57: unsigned long j;
! 58: unsigned long k;
! 59: unsigned long taille_matrice_f77;
! 60:
! 61: void *matrice_f77;
! 62:
! 63: rang(s_etat_processus, s_matrice, &rang_estime);
! 64:
! 65: if ((*s_etat_processus).erreur_systeme != d_es)
! 66: {
! 67: return;
! 68: }
! 69:
! 70: if (((*s_etat_processus).erreur_execution != d_ex) ||
! 71: ((*s_etat_processus).exception != d_ep))
! 72: {
! 73: return;
! 74: }
! 75:
! 76: if (rang_estime < (integer8) (*s_matrice).nombre_lignes)
! 77: {
! 78: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 79: return;
! 80: }
! 81:
! 82: taille_matrice_f77 = (*s_matrice).nombre_lignes *
! 83: (*s_matrice).nombre_colonnes;
! 84: dim_matrice = (integer4) (*s_matrice).nombre_lignes;
! 85:
! 86: switch((*s_matrice).type)
! 87: {
! 88: case 'I' :
! 89: {
! 90: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 91: sizeof(real8))) == NULL)
! 92: {
! 93: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 94: return;
! 95: }
! 96:
! 97: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 98: {
! 99: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 100: {
! 101: ((real8 *) matrice_f77)[k++] = ((integer8 **)
! 102: (*s_matrice).tableau)[j][i];
! 103: }
! 104: }
! 105:
! 106: for(i = 0; i < (*s_matrice).nombre_lignes; i++)
! 107: {
! 108: free(((integer8 **) (*s_matrice).tableau)[i]);
! 109: }
! 110:
! 111: free((integer8 **) (*s_matrice).tableau);
! 112:
! 113: if (((*s_matrice).tableau = (void **) malloc((*s_matrice)
! 114: .nombre_lignes * sizeof(real8 *))) == NULL)
! 115: {
! 116: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 117: return;
! 118: }
! 119:
! 120: for(i = 0; i < (*s_matrice).nombre_lignes; i++)
! 121: {
! 122: if ((((*s_matrice).tableau)[i] =
! 123: (real8 *) malloc((*s_matrice)
! 124: .nombre_colonnes * sizeof(real8))) == NULL)
! 125: {
! 126: (*s_etat_processus).erreur_systeme =
! 127: d_es_allocation_memoire;
! 128: return;
! 129: }
! 130: }
! 131:
! 132: (*s_matrice).type = 'R';
! 133:
! 134: if ((pivot = (integer4 *) malloc((*s_matrice).nombre_lignes *
! 135: sizeof(integer4))) == NULL)
! 136: {
! 137: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 138: return;
! 139: }
! 140:
! 141: if ((work = (real8 *) malloc(sizeof(real8))) == NULL)
! 142: {
! 143: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 144: return;
! 145: }
! 146:
! 147: dim_work = -1;
! 148:
! 149: dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work,
! 150: &dim_work, &erreur);
! 151:
! 152: if (erreur != 0)
! 153: {
! 154: if (erreur > 0)
! 155: {
! 156: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 157: }
! 158: else
! 159: {
! 160: (*s_etat_processus).erreur_execution =
! 161: d_ex_routines_mathematiques;
! 162: }
! 163:
! 164: free(pivot);
! 165: free(work);
! 166: free(matrice_f77);
! 167: return;
! 168: }
! 169:
! 170: dim_work = (integer4) work[0];
! 171:
! 172: free(work);
! 173:
! 174: if ((work = (real8 *) malloc(dim_work * sizeof(real8))) == NULL)
! 175: {
! 176: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 177: return;
! 178: }
! 179:
! 180: dgetrf_(&dim_matrice, &dim_matrice, matrice_f77,
! 181: &dim_matrice, pivot, &erreur);
! 182:
! 183: if (erreur != 0)
! 184: {
! 185: if (erreur > 0)
! 186: {
! 187: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 188: }
! 189: else
! 190: {
! 191: (*s_etat_processus).erreur_execution =
! 192: d_ex_routines_mathematiques;
! 193: }
! 194:
! 195: free(pivot);
! 196: free(work);
! 197: free(matrice_f77);
! 198: return;
! 199: }
! 200:
! 201: dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work,
! 202: &dim_work, &erreur);
! 203:
! 204: if (erreur != 0)
! 205: {
! 206: if (erreur > 0)
! 207: {
! 208: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 209: }
! 210: else
! 211: {
! 212: (*s_etat_processus).erreur_execution =
! 213: d_ex_routines_mathematiques;
! 214: }
! 215:
! 216: free(pivot);
! 217: free(work);
! 218: free(matrice_f77);
! 219: return;
! 220: }
! 221:
! 222: free(work);
! 223: free(pivot);
! 224:
! 225: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 226: {
! 227: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 228: {
! 229: ((real8 **) (*s_matrice).tableau)[j][i] =
! 230: ((real8 *) matrice_f77)[k++];
! 231: }
! 232: }
! 233:
! 234: free(matrice_f77);
! 235:
! 236: break;
! 237: }
! 238:
! 239: case 'R' :
! 240: {
! 241: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 242: sizeof(real8))) == NULL)
! 243: {
! 244: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 245: return;
! 246: }
! 247:
! 248: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 249: {
! 250: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 251: {
! 252: ((real8 *) matrice_f77)[k++] = ((real8 **)
! 253: (*s_matrice).tableau)[j][i];
! 254: }
! 255: }
! 256:
! 257: if ((pivot = (integer4 *) malloc((*s_matrice).nombre_lignes *
! 258: sizeof(integer4))) == NULL)
! 259: {
! 260: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 261: return;
! 262: }
! 263:
! 264: if ((work = (real8 *) malloc(sizeof(real8))) == NULL)
! 265: {
! 266: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 267: return;
! 268: }
! 269:
! 270: dim_work = -1;
! 271:
! 272: dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work,
! 273: &dim_work, &erreur);
! 274:
! 275: if (erreur != 0)
! 276: {
! 277: if (erreur > 0)
! 278: {
! 279: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 280: }
! 281: else
! 282: {
! 283: (*s_etat_processus).erreur_execution =
! 284: d_ex_routines_mathematiques;
! 285: }
! 286:
! 287: free(pivot);
! 288: free(work);
! 289: free(matrice_f77);
! 290: return;
! 291: }
! 292:
! 293: dim_work = (integer4) work[0];
! 294:
! 295: free(work);
! 296:
! 297: if ((work = (real8 *) malloc(dim_work * sizeof(real8))) == NULL)
! 298: {
! 299: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 300: return;
! 301: }
! 302:
! 303: dgetrf_(&dim_matrice, &dim_matrice, matrice_f77,
! 304: &dim_matrice, pivot, &erreur);
! 305:
! 306: if (erreur != 0)
! 307: {
! 308: if (erreur > 0)
! 309: {
! 310: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 311: }
! 312: else
! 313: {
! 314: (*s_etat_processus).erreur_execution =
! 315: d_ex_routines_mathematiques;
! 316: }
! 317:
! 318: free(pivot);
! 319: free(work);
! 320: free(matrice_f77);
! 321: return;
! 322: }
! 323:
! 324: dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work,
! 325: &dim_work, &erreur);
! 326:
! 327: if (erreur != 0)
! 328: {
! 329: if (erreur > 0)
! 330: {
! 331: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 332: }
! 333: else
! 334: {
! 335: (*s_etat_processus).erreur_execution =
! 336: d_ex_routines_mathematiques;
! 337: }
! 338:
! 339: free(pivot);
! 340: free(work);
! 341: free(matrice_f77);
! 342: return;
! 343: }
! 344:
! 345: free(work);
! 346: free(pivot);
! 347:
! 348: for(k = 0, i = 0; i < (*s_matrice).nombre_lignes; i++)
! 349: {
! 350: for(j = 0; j < (*s_matrice).nombre_colonnes; j++)
! 351: {
! 352: ((real8 **) (*s_matrice).tableau)[j][i] =
! 353: ((real8 *) matrice_f77)[k++];
! 354: }
! 355: }
! 356:
! 357: free(matrice_f77);
! 358:
! 359: break;
! 360: }
! 361:
! 362: case 'C' :
! 363: {
! 364: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 365: sizeof(struct_complexe16))) == NULL)
! 366: {
! 367: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 368: return;
! 369: }
! 370:
! 371: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 372: {
! 373: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 374: {
! 375: (((struct_complexe16 *) matrice_f77)[k]).partie_reelle =
! 376: (((struct_complexe16 **) (*s_matrice).tableau)
! 377: [j][i]).partie_reelle;
! 378: (((struct_complexe16 *) matrice_f77)[k]).partie_imaginaire =
! 379: (((struct_complexe16 **) (*s_matrice).tableau)
! 380: [j][i]).partie_imaginaire;
! 381: k++;
! 382: }
! 383: }
! 384:
! 385: if ((pivot = (integer4 *) malloc((*s_matrice).nombre_lignes *
! 386: sizeof(integer4))) == NULL)
! 387: {
! 388: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 389: return;
! 390: }
! 391:
! 392: if ((c_work = (struct_complexe16 *) malloc(
! 393: sizeof(struct_complexe16))) == NULL)
! 394: {
! 395: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 396: return;
! 397: }
! 398:
! 399: dim_work = -1;
! 400:
! 401: zgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, c_work,
! 402: &dim_work, &erreur);
! 403:
! 404: if (erreur != 0)
! 405: {
! 406: if (erreur > 0)
! 407: {
! 408: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 409: }
! 410: else
! 411: {
! 412: (*s_etat_processus).erreur_execution =
! 413: d_ex_routines_mathematiques;
! 414: }
! 415:
! 416: free(pivot);
! 417: free(c_work);
! 418: free(matrice_f77);
! 419: return;
! 420: }
! 421:
! 422: dim_work = (integer4) c_work[0].partie_reelle;
! 423:
! 424: free(c_work);
! 425:
! 426: if ((c_work = (struct_complexe16 *) malloc(dim_work *
! 427: sizeof(struct_complexe16))) == NULL)
! 428: {
! 429: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 430: return;
! 431: }
! 432:
! 433: zgetrf_(&dim_matrice, &dim_matrice, matrice_f77,
! 434: &dim_matrice, pivot, &erreur);
! 435:
! 436: if (erreur != 0)
! 437: {
! 438: if (erreur > 0)
! 439: {
! 440: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 441: }
! 442: else
! 443: {
! 444: (*s_etat_processus).erreur_execution =
! 445: d_ex_routines_mathematiques;
! 446: }
! 447:
! 448: free(pivot);
! 449: free(c_work);
! 450: free(matrice_f77);
! 451: return;
! 452: }
! 453:
! 454: zgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, c_work,
! 455: &dim_work, &erreur);
! 456:
! 457: if (erreur != 0)
! 458: {
! 459: if (erreur > 0)
! 460: {
! 461: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
! 462: }
! 463: else
! 464: {
! 465: (*s_etat_processus).erreur_execution =
! 466: d_ex_routines_mathematiques;
! 467: }
! 468:
! 469: free(pivot);
! 470: free(c_work);
! 471: free(matrice_f77);
! 472: return;
! 473: }
! 474:
! 475: free(c_work);
! 476: free(pivot);
! 477:
! 478: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 479: {
! 480: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 481: {
! 482: (((struct_complexe16 **) (*s_matrice).tableau)
! 483: [j][i]).partie_reelle = (((struct_complexe16 *)
! 484: matrice_f77)[k]).partie_reelle;
! 485: (((struct_complexe16 **) (*s_matrice).tableau)
! 486: [j][i]).partie_imaginaire = (((struct_complexe16 *)
! 487: matrice_f77)[k]).partie_imaginaire;
! 488: k++;
! 489: }
! 490: }
! 491:
! 492: free(matrice_f77);
! 493:
! 494: break;
! 495: }
! 496: }
! 497:
! 498: return;
! 499: }
! 500:
! 501:
! 502: /*
! 503: ================================================================================
! 504: Fonction calculant les vecteurs propres d'une matrice carrée
! 505: ================================================================================
! 506: Entrées : struct_matrice, pointeur sur la matrice d'entrée, ainsi que
! 507: struct_matrice, pointeur sur le vecteur des valeurs propres
! 508: (vecteur de complexes) et sur deux matrices complexes
! 509: contenant les différents vecteurs propres gauches et droits.
! 510: Les pointeurs sur les vecteurs propres peuvent être nuls,
! 511: et seules les valeurs propres sont calculées.
! 512: L'allocation des tableaux de sortie est effectuée dans la
! 513: routine, mais les structures s_matrice et s_vecteurs
! 514: doivent être passées en paramètre.
! 515: La matrice présente en entrée n'est pas libérée.
! 516: --------------------------------------------------------------------------------
! 517: Sorties : vecteur contenant les valeurs propres, matrice contenant
! 518: les vecteurs propres gauches et matrice contenant les vecteurs
! 519: propres droits. Toutes les sorties sont complexes.
! 520: --------------------------------------------------------------------------------
! 521: Effets de bord : néant
! 522: ================================================================================
! 523: */
! 524:
! 525: void
! 526: valeurs_propres(struct_processus *s_etat_processus,
! 527: struct_matrice *s_matrice,
! 528: struct_vecteur *s_valeurs_propres,
! 529: struct_matrice *s_vecteurs_propres_gauches,
! 530: struct_matrice *s_vecteurs_propres_droits)
! 531: {
! 532: real8 *rwork;
! 533:
! 534: integer4 dim_matrice;
! 535: integer4 lwork;
! 536: integer4 erreur;
! 537:
! 538: struct_complexe16 *matrice_f77;
! 539: struct_complexe16 *vpd_f77;
! 540: struct_complexe16 *vpg_f77;
! 541: struct_complexe16 *work;
! 542:
! 543: unsigned char calcul_vp_droits;
! 544: unsigned char calcul_vp_gauches;
! 545: unsigned char negatif;
! 546:
! 547: unsigned long i;
! 548: unsigned long j;
! 549: unsigned long k;
! 550: unsigned long taille_matrice_f77;
! 551:
! 552: taille_matrice_f77 = (*s_matrice).nombre_lignes *
! 553: (*s_matrice).nombre_colonnes;
! 554: dim_matrice = (integer4) (*s_matrice).nombre_lignes;
! 555:
! 556: /*
! 557: * Allocation de la matrice complexe
! 558: */
! 559:
! 560: if ((matrice_f77 = (struct_complexe16 *) malloc(taille_matrice_f77 *
! 561: sizeof(struct_complexe16))) == NULL)
! 562: {
! 563: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 564: return;
! 565: }
! 566:
! 567: /*
! 568: * Garniture de la matrice de compatibilité Fortran
! 569: */
! 570:
! 571: switch((*s_matrice).type)
! 572: {
! 573: /*
! 574: * La matrice en entrée est une matrice d'entiers.
! 575: */
! 576:
! 577: case 'I' :
! 578: {
! 579: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 580: {
! 581: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 582: {
! 583: ((struct_complexe16 *) matrice_f77)[k]
! 584: .partie_reelle = (real8) ((integer8 **)
! 585: (*s_matrice).tableau)[j][i];
! 586: ((struct_complexe16 *) matrice_f77)[k++]
! 587: .partie_imaginaire = (real8) 0;
! 588: }
! 589: }
! 590:
! 591: break;
! 592: }
! 593:
! 594: /*
! 595: * La matrice en entrée est une matrice de réels.
! 596: */
! 597:
! 598: case 'R' :
! 599: {
! 600: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 601: {
! 602: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 603: {
! 604: ((struct_complexe16 *) matrice_f77)[k]
! 605: .partie_reelle = ((real8 **)
! 606: (*s_matrice).tableau)[j][i];
! 607: ((struct_complexe16 *) matrice_f77)[k++]
! 608: .partie_imaginaire = (real8) 0;
! 609: }
! 610: }
! 611:
! 612: break;
! 613: }
! 614:
! 615: /*
! 616: * La matrice en entrée est une matrice de complexes.
! 617: */
! 618:
! 619: case 'C' :
! 620: {
! 621: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 622: {
! 623: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 624: {
! 625: ((struct_complexe16 *) matrice_f77)[k]
! 626: .partie_reelle = ((struct_complexe16 **)
! 627: (*s_matrice).tableau)[j][i].partie_reelle;
! 628: ((struct_complexe16 *) matrice_f77)[k++]
! 629: .partie_imaginaire = ((struct_complexe16 **)
! 630: (*s_matrice).tableau)[j][i].partie_imaginaire;
! 631: }
! 632: }
! 633:
! 634: break;
! 635: }
! 636: }
! 637:
! 638: (*s_valeurs_propres).type = 'C';
! 639: (*s_valeurs_propres).taille = (*s_matrice).nombre_lignes;
! 640:
! 641: if (((*s_valeurs_propres).tableau = (struct_complexe16 *)
! 642: malloc((*s_valeurs_propres).taille * sizeof(struct_complexe16)))
! 643: == NULL)
! 644: {
! 645: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 646: return;
! 647: }
! 648:
! 649: if (s_vecteurs_propres_gauches != NULL)
! 650: {
! 651: (*s_vecteurs_propres_gauches).type = 'C';
! 652: calcul_vp_gauches = 'V';
! 653:
! 654: if ((vpg_f77 = (struct_complexe16 *) malloc(taille_matrice_f77 *
! 655: sizeof(struct_complexe16))) == NULL)
! 656: {
! 657: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 658: return;
! 659: }
! 660: }
! 661: else
! 662: {
! 663: vpg_f77 = NULL;
! 664: calcul_vp_gauches = 'N';
! 665: }
! 666:
! 667: if (s_vecteurs_propres_droits != NULL)
! 668: {
! 669: (*s_vecteurs_propres_droits).type = 'C';
! 670: calcul_vp_droits = 'V';
! 671:
! 672: if ((vpd_f77 = (struct_complexe16 *) malloc(taille_matrice_f77 *
! 673: sizeof(struct_complexe16))) == NULL)
! 674: {
! 675: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 676: return;
! 677: }
! 678: }
! 679: else
! 680: {
! 681: vpd_f77 = NULL;
! 682: calcul_vp_droits = 'N';
! 683: }
! 684:
! 685: negatif = 'N';
! 686: lwork = -1;
! 687:
! 688: if ((rwork = (real8 *) malloc(2 * (*s_matrice).nombre_lignes *
! 689: sizeof(real8))) == NULL)
! 690: {
! 691: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 692: return;
! 693: }
! 694:
! 695: if ((work = (struct_complexe16 *) malloc(sizeof(struct_complexe16)))
! 696: == NULL)
! 697: {
! 698: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 699: return;
! 700: }
! 701:
! 702: zgeev_(&negatif, &negatif, &dim_matrice, matrice_f77, &dim_matrice,
! 703: (struct_complexe16 *) (*s_valeurs_propres).tableau,
! 704: vpg_f77, &dim_matrice, vpd_f77, &dim_matrice,
! 705: work, &lwork, rwork, &erreur, 1, 1);
! 706:
! 707: if (erreur != 0)
! 708: {
! 709: if (erreur > 0)
! 710: {
! 711: (*s_etat_processus).exception = d_ep_decomposition_QR;
! 712: }
! 713: else
! 714: {
! 715: (*s_etat_processus).erreur_execution =
! 716: d_ex_routines_mathematiques;
! 717: }
! 718:
! 719: free(work);
! 720: free(rwork);
! 721: free(matrice_f77);
! 722:
! 723: if (calcul_vp_gauches == 'V')
! 724: {
! 725: free(vpg_f77);
! 726: }
! 727:
! 728: if (calcul_vp_droits == 'V')
! 729: {
! 730: free(vpd_f77);
! 731: }
! 732:
! 733: return;
! 734: }
! 735:
! 736: lwork = (integer4) work[0].partie_reelle;
! 737: free(work);
! 738:
! 739: if ((work = (struct_complexe16 *) malloc(lwork * sizeof(struct_complexe16)))
! 740: == NULL)
! 741: {
! 742: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 743: return;
! 744: }
! 745:
! 746: zgeev_(&calcul_vp_gauches, &calcul_vp_droits, &dim_matrice,
! 747: matrice_f77, &dim_matrice,
! 748: (struct_complexe16 *) (*s_valeurs_propres).tableau,
! 749: vpg_f77, &dim_matrice, vpd_f77, &dim_matrice,
! 750: work, &lwork, rwork, &erreur, 1, 1);
! 751:
! 752: if (erreur != 0)
! 753: {
! 754: if (erreur > 0)
! 755: {
! 756: (*s_etat_processus).exception = d_ep_decomposition_QR;
! 757: }
! 758: else
! 759: {
! 760: (*s_etat_processus).erreur_execution =
! 761: d_ex_routines_mathematiques;
! 762: }
! 763:
! 764: free(work);
! 765: free(rwork);
! 766: free(matrice_f77);
! 767:
! 768: if (calcul_vp_gauches == 'V')
! 769: {
! 770: free(vpg_f77);
! 771: }
! 772:
! 773: if (calcul_vp_droits == 'V')
! 774: {
! 775: free(vpd_f77);
! 776: }
! 777:
! 778: return;
! 779: }
! 780:
! 781: free(work);
! 782: free(rwork);
! 783:
! 784: if (calcul_vp_gauches == 'V')
! 785: {
! 786: (*s_vecteurs_propres_gauches).type = 'C';
! 787: (*s_vecteurs_propres_gauches).nombre_lignes =
! 788: (*s_matrice).nombre_lignes;
! 789: (*s_vecteurs_propres_gauches).nombre_colonnes =
! 790: (*s_matrice).nombre_colonnes;
! 791:
! 792: if (((*s_vecteurs_propres_gauches).tableau = malloc(
! 793: (*s_vecteurs_propres_gauches).nombre_lignes *
! 794: sizeof(struct_complexe16 *))) == NULL)
! 795: {
! 796: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 797: return;
! 798: }
! 799:
! 800: for(i = 0; i < (*s_vecteurs_propres_gauches).nombre_lignes; i++)
! 801: {
! 802: if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches)
! 803: .tableau)[i] = (struct_complexe16 *) malloc(
! 804: (*s_vecteurs_propres_gauches).nombre_colonnes *
! 805: sizeof(struct_complexe16))) == NULL)
! 806: {
! 807: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 808: return;
! 809: }
! 810: }
! 811:
! 812: for(k = 0, i = 0; i < (*s_vecteurs_propres_gauches).nombre_colonnes;
! 813: i++)
! 814: {
! 815: for(j = 0; j < (*s_vecteurs_propres_gauches).nombre_lignes; j++)
! 816: {
! 817: ((struct_complexe16 **) (*s_vecteurs_propres_gauches)
! 818: .tableau)[j][i].partie_reelle =
! 819: vpg_f77[k].partie_reelle;
! 820: ((struct_complexe16 **) (*s_vecteurs_propres_gauches)
! 821: .tableau)[j][i].partie_imaginaire =
! 822: vpg_f77[k++].partie_imaginaire;
! 823: }
! 824: }
! 825:
! 826: free(vpg_f77);
! 827: }
! 828:
! 829: if (calcul_vp_droits == 'V')
! 830: {
! 831: (*s_vecteurs_propres_droits).type = 'C';
! 832: (*s_vecteurs_propres_droits).nombre_lignes =
! 833: (*s_matrice).nombre_lignes;
! 834: (*s_vecteurs_propres_droits).nombre_colonnes =
! 835: (*s_matrice).nombre_colonnes;
! 836:
! 837: if (((*s_vecteurs_propres_droits).tableau = malloc(
! 838: (*s_vecteurs_propres_droits).nombre_lignes *
! 839: sizeof(struct_complexe16 *))) == NULL)
! 840: {
! 841: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 842: return;
! 843: }
! 844:
! 845: for(i = 0; i < (*s_vecteurs_propres_droits).nombre_lignes; i++)
! 846: {
! 847: if ((((struct_complexe16 **) (*s_vecteurs_propres_droits)
! 848: .tableau)[i] = (struct_complexe16 *) malloc(
! 849: (*s_vecteurs_propres_droits).nombre_colonnes *
! 850: sizeof(struct_complexe16))) == NULL)
! 851: {
! 852: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 853: return;
! 854: }
! 855: }
! 856:
! 857: for(k = 0, i = 0; i < (*s_vecteurs_propres_droits).nombre_colonnes; i++)
! 858: {
! 859: for(j = 0; j < (*s_vecteurs_propres_droits).nombre_lignes; j++)
! 860: {
! 861: ((struct_complexe16 **) (*s_vecteurs_propres_droits)
! 862: .tableau)[j][i].partie_reelle =
! 863: vpd_f77[k].partie_reelle;
! 864: ((struct_complexe16 **) (*s_vecteurs_propres_droits)
! 865: .tableau)[j][i].partie_imaginaire =
! 866: vpd_f77[k++].partie_imaginaire;
! 867: }
! 868: }
! 869:
! 870: free(vpd_f77);
! 871: }
! 872:
! 873: free(matrice_f77);
! 874:
! 875: return;
! 876: }
! 877:
! 878:
! 879: /*
! 880: ================================================================================
! 881: Fonction calculant les vecteurs propres généralisés d'une matrice carrée
! 882: dans une métrique.
! 883: ================================================================================
! 884: Entrées : struct_matrice, pointeur sur la matrice d'entrée, ainsi que
! 885: struct_matrice, pointeur sur la métrique et
! 886: struct_matrice, pointeur sur le vecteur des valeurs propres
! 887: (vecteur de complexes) et sur deux matrices complexes
! 888: contenant les différents vecteurs propres gauches et droits.
! 889: Les pointeurs sur les vecteurs propres peuvent être nuls,
! 890: et seules les valeurs propres sont calculées.
! 891: L'allocation des tableaux de sortie est effectuée dans la
! 892: routine, mais les structures s_matrice et s_vecteurs
! 893: doivent être passées en paramètre.
! 894: La matrice présente en entrée n'est pas libérée.
! 895: --------------------------------------------------------------------------------
! 896: Sorties : vecteur contenant les valeurs propres, matrice contenant
! 897: les vecteurs propres gauches et matrice contenant les vecteurs
! 898: propres droits. Toutes les sorties sont complexes.
! 899: --------------------------------------------------------------------------------
! 900: Effets de bord : néant
! 901: ================================================================================
! 902: */
! 903:
! 904: void
! 905: valeurs_propres_generalisees(struct_processus *s_etat_processus,
! 906: struct_matrice *s_matrice,
! 907: struct_matrice *s_metrique,
! 908: struct_vecteur *s_valeurs_propres,
! 909: struct_matrice *s_vecteurs_propres_gauches,
! 910: struct_matrice *s_vecteurs_propres_droits)
! 911: {
! 912: real8 *rwork;
! 913:
! 914: integer4 dim_matrice;
! 915: integer4 lwork;
! 916: integer4 erreur;
! 917:
! 918: struct_complexe16 *alpha;
! 919: struct_complexe16 *beta;
! 920: struct_complexe16 *matrice_f77;
! 921: struct_complexe16 *metrique_f77;
! 922: struct_complexe16 *vpd_f77;
! 923: struct_complexe16 *vpg_f77;
! 924: struct_complexe16 *work;
! 925:
! 926: unsigned char calcul_vp_droits;
! 927: unsigned char calcul_vp_gauches;
! 928: unsigned char negatif;
! 929:
! 930: unsigned long i;
! 931: unsigned long j;
! 932: unsigned long k;
! 933: unsigned long taille_matrice_f77;
! 934:
! 935: taille_matrice_f77 = (*s_matrice).nombre_lignes *
! 936: (*s_matrice).nombre_colonnes;
! 937: dim_matrice = (integer4) (*s_matrice).nombre_lignes;
! 938:
! 939: /*
! 940: * Allocation de la matrice complexe
! 941: */
! 942:
! 943: if ((matrice_f77 = (struct_complexe16 *) malloc(taille_matrice_f77 *
! 944: sizeof(struct_complexe16))) == NULL)
! 945: {
! 946: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 947: return;
! 948: }
! 949:
! 950: /*
! 951: * Garniture de la matrice de compatibilité Fortran
! 952: */
! 953:
! 954: switch((*s_matrice).type)
! 955: {
! 956: /*
! 957: * La matrice en entrée est une matrice d'entiers.
! 958: */
! 959:
! 960: case 'I' :
! 961: {
! 962: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 963: {
! 964: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 965: {
! 966: ((struct_complexe16 *) matrice_f77)[k]
! 967: .partie_reelle = (real8) ((integer8 **)
! 968: (*s_matrice).tableau)[j][i];
! 969: ((struct_complexe16 *) matrice_f77)[k++]
! 970: .partie_imaginaire = (real8) 0;
! 971: }
! 972: }
! 973:
! 974: break;
! 975: }
! 976:
! 977: /*
! 978: * La matrice en entrée est une matrice de réels.
! 979: */
! 980:
! 981: case 'R' :
! 982: {
! 983: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 984: {
! 985: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 986: {
! 987: ((struct_complexe16 *) matrice_f77)[k]
! 988: .partie_reelle = ((real8 **)
! 989: (*s_matrice).tableau)[j][i];
! 990: ((struct_complexe16 *) matrice_f77)[k++]
! 991: .partie_imaginaire = (real8) 0;
! 992: }
! 993: }
! 994:
! 995: break;
! 996: }
! 997:
! 998: /*
! 999: * La matrice en entrée est une matrice de complexes.
! 1000: */
! 1001:
! 1002: case 'C' :
! 1003: {
! 1004: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1005: {
! 1006: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 1007: {
! 1008: ((struct_complexe16 *) matrice_f77)[k]
! 1009: .partie_reelle = ((struct_complexe16 **)
! 1010: (*s_matrice).tableau)[j][i].partie_reelle;
! 1011: ((struct_complexe16 *) matrice_f77)[k++]
! 1012: .partie_imaginaire = ((struct_complexe16 **)
! 1013: (*s_matrice).tableau)[j][i].partie_imaginaire;
! 1014: }
! 1015: }
! 1016:
! 1017: break;
! 1018: }
! 1019: }
! 1020:
! 1021: /*
! 1022: * Allocation de la metrique complexe
! 1023: */
! 1024:
! 1025: if ((metrique_f77 = (struct_complexe16 *) malloc(taille_matrice_f77 *
! 1026: sizeof(struct_complexe16))) == NULL)
! 1027: {
! 1028: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1029: return;
! 1030: }
! 1031:
! 1032: /*
! 1033: * Garniture de la matrice de compatibilité Fortran
! 1034: */
! 1035:
! 1036: switch((*s_metrique).type)
! 1037: {
! 1038: /*
! 1039: * La matrice en entrée est une matrice d'entiers.
! 1040: */
! 1041:
! 1042: case 'I' :
! 1043: {
! 1044: for(k = 0, i = 0; i < (*s_metrique).nombre_colonnes; i++)
! 1045: {
! 1046: for(j = 0; j < (*s_metrique).nombre_lignes; j++)
! 1047: {
! 1048: ((struct_complexe16 *) metrique_f77)[k]
! 1049: .partie_reelle = (real8) ((integer8 **)
! 1050: (*s_metrique).tableau)[j][i];
! 1051: ((struct_complexe16 *) metrique_f77)[k++]
! 1052: .partie_imaginaire = (real8) 0;
! 1053: }
! 1054: }
! 1055:
! 1056: break;
! 1057: }
! 1058:
! 1059: /*
! 1060: * La matrice en entrée est une matrice de réels.
! 1061: */
! 1062:
! 1063: case 'R' :
! 1064: {
! 1065: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1066: {
! 1067: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 1068: {
! 1069: ((struct_complexe16 *) metrique_f77)[k]
! 1070: .partie_reelle = ((real8 **)
! 1071: (*s_metrique).tableau)[j][i];
! 1072: ((struct_complexe16 *) metrique_f77)[k++]
! 1073: .partie_imaginaire = (real8) 0;
! 1074: }
! 1075: }
! 1076:
! 1077: break;
! 1078: }
! 1079:
! 1080: /*
! 1081: * La matrice en entrée est une matrice de complexes.
! 1082: */
! 1083:
! 1084: case 'C' :
! 1085: {
! 1086: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1087: {
! 1088: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 1089: {
! 1090: ((struct_complexe16 *) metrique_f77)[k]
! 1091: .partie_reelle = ((struct_complexe16 **)
! 1092: (*s_metrique).tableau)[j][i].partie_reelle;
! 1093: ((struct_complexe16 *) metrique_f77)[k++]
! 1094: .partie_imaginaire = ((struct_complexe16 **)
! 1095: (*s_metrique).tableau)[j][i].partie_imaginaire;
! 1096: }
! 1097: }
! 1098:
! 1099: break;
! 1100: }
! 1101: }
! 1102:
! 1103: (*s_valeurs_propres).type = 'C';
! 1104: (*s_valeurs_propres).taille = (*s_matrice).nombre_lignes;
! 1105:
! 1106: if (((*s_valeurs_propres).tableau = (struct_complexe16 *)
! 1107: malloc((*s_valeurs_propres).taille * sizeof(struct_complexe16)))
! 1108: == NULL)
! 1109: {
! 1110: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1111: return;
! 1112: }
! 1113:
! 1114: if (s_vecteurs_propres_gauches != NULL)
! 1115: {
! 1116: (*s_vecteurs_propres_gauches).type = 'C';
! 1117: calcul_vp_gauches = 'V';
! 1118:
! 1119: if ((vpg_f77 = (struct_complexe16 *) malloc(taille_matrice_f77 *
! 1120: sizeof(struct_complexe16))) == NULL)
! 1121: {
! 1122: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1123: return;
! 1124: }
! 1125: }
! 1126: else
! 1127: {
! 1128: vpg_f77 = NULL;
! 1129: calcul_vp_gauches = 'N';
! 1130: }
! 1131:
! 1132: if (s_vecteurs_propres_droits != NULL)
! 1133: {
! 1134: (*s_vecteurs_propres_droits).type = 'C';
! 1135: calcul_vp_droits = 'V';
! 1136:
! 1137: if ((vpd_f77 = (struct_complexe16 *) malloc(taille_matrice_f77 *
! 1138: sizeof(struct_complexe16))) == NULL)
! 1139: {
! 1140: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1141: return;
! 1142: }
! 1143: }
! 1144: else
! 1145: {
! 1146: vpd_f77 = NULL;
! 1147: calcul_vp_droits = 'N';
! 1148: }
! 1149:
! 1150: negatif = 'N';
! 1151: lwork = -1;
! 1152:
! 1153: if ((rwork = (real8 *) malloc(8 * (*s_matrice).nombre_lignes *
! 1154: sizeof(real8))) == NULL)
! 1155: {
! 1156: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1157: return;
! 1158: }
! 1159:
! 1160: if ((work = (struct_complexe16 *) malloc(sizeof(struct_complexe16)))
! 1161: == NULL)
! 1162: {
! 1163: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1164: return;
! 1165: }
! 1166:
! 1167: if ((alpha = (struct_complexe16 *) malloc((*s_valeurs_propres).taille *
! 1168: sizeof(struct_complexe16))) == NULL)
! 1169: {
! 1170: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1171: return;
! 1172: }
! 1173:
! 1174: if ((beta = (struct_complexe16 *) malloc((*s_valeurs_propres).taille *
! 1175: sizeof(struct_complexe16))) == NULL)
! 1176: {
! 1177: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1178: return;
! 1179: }
! 1180:
! 1181: zggev_(&negatif, &negatif, &dim_matrice, matrice_f77, &dim_matrice,
! 1182: metrique_f77, &dim_matrice, alpha, beta,
! 1183: vpg_f77, &dim_matrice, vpd_f77, &dim_matrice,
! 1184: work, &lwork, rwork, &erreur, 1, 1);
! 1185:
! 1186: if (erreur != 0)
! 1187: {
! 1188: if (erreur > 0)
! 1189: {
! 1190: (*s_etat_processus).exception = d_ep_decomposition_QZ;
! 1191: }
! 1192: else
! 1193: {
! 1194: (*s_etat_processus).erreur_execution =
! 1195: d_ex_routines_mathematiques;
! 1196: }
! 1197:
! 1198: free(work);
! 1199: free(rwork);
! 1200: free(matrice_f77);
! 1201: free(metrique_f77);
! 1202:
! 1203: if (calcul_vp_gauches == 'V')
! 1204: {
! 1205: free(vpg_f77);
! 1206:
! 1207: (*s_vecteurs_propres_gauches).type = 'C';
! 1208: (*s_vecteurs_propres_gauches).nombre_lignes = 1;
! 1209: (*s_vecteurs_propres_gauches).nombre_colonnes = 1;
! 1210:
! 1211: if (((*s_vecteurs_propres_gauches).tableau = malloc(
! 1212: (*s_vecteurs_propres_gauches).nombre_lignes *
! 1213: sizeof(struct_complexe16 *))) == NULL)
! 1214: {
! 1215: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1216: return;
! 1217: }
! 1218:
! 1219: if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches)
! 1220: .tableau)[0] = (struct_complexe16 *) malloc(
! 1221: (*s_vecteurs_propres_gauches).nombre_colonnes *
! 1222: sizeof(struct_complexe16))) == NULL)
! 1223: {
! 1224: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1225: return;
! 1226: }
! 1227: }
! 1228:
! 1229: if (calcul_vp_droits == 'V')
! 1230: {
! 1231: free(vpd_f77);
! 1232:
! 1233: (*s_vecteurs_propres_droits).type = 'C';
! 1234: (*s_vecteurs_propres_droits).nombre_lignes = 1;
! 1235: (*s_vecteurs_propres_droits).nombre_colonnes = 1;
! 1236:
! 1237: if (((*s_vecteurs_propres_droits).tableau = malloc(
! 1238: (*s_vecteurs_propres_droits).nombre_lignes *
! 1239: sizeof(struct_complexe16 *))) == NULL)
! 1240: {
! 1241: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1242: return;
! 1243: }
! 1244:
! 1245: if ((((struct_complexe16 **) (*s_vecteurs_propres_droits)
! 1246: .tableau)[0] = (struct_complexe16 *) malloc(
! 1247: (*s_vecteurs_propres_droits).nombre_colonnes *
! 1248: sizeof(struct_complexe16))) == NULL)
! 1249: {
! 1250: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1251: return;
! 1252: }
! 1253: }
! 1254:
! 1255: return;
! 1256: }
! 1257:
! 1258: lwork = (integer4) work[0].partie_reelle;
! 1259: free(work);
! 1260:
! 1261: if ((work = (struct_complexe16 *) malloc(lwork * sizeof(struct_complexe16)))
! 1262: == NULL)
! 1263: {
! 1264: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1265: return;
! 1266: }
! 1267:
! 1268: zggev_(&calcul_vp_gauches, &calcul_vp_droits, &dim_matrice,
! 1269: matrice_f77, &dim_matrice,
! 1270: metrique_f77, &dim_matrice, alpha, beta,
! 1271: vpg_f77, &dim_matrice, vpd_f77, &dim_matrice,
! 1272: work, &lwork, rwork, &erreur, 1, 1);
! 1273:
! 1274: if (erreur != 0)
! 1275: {
! 1276: if (erreur > 0)
! 1277: {
! 1278: (*s_etat_processus).exception = d_ep_decomposition_QZ;
! 1279: }
! 1280: else
! 1281: {
! 1282: (*s_etat_processus).erreur_execution =
! 1283: d_ex_routines_mathematiques;
! 1284: }
! 1285:
! 1286: free(work);
! 1287: free(rwork);
! 1288: free(matrice_f77);
! 1289: free(metrique_f77);
! 1290:
! 1291: if (calcul_vp_gauches == 'V')
! 1292: {
! 1293: free(vpg_f77);
! 1294:
! 1295: (*s_vecteurs_propres_gauches).type = 'C';
! 1296: (*s_vecteurs_propres_gauches).nombre_lignes = 1;
! 1297: (*s_vecteurs_propres_gauches).nombre_colonnes = 1;
! 1298:
! 1299: if (((*s_vecteurs_propres_gauches).tableau = malloc(
! 1300: (*s_vecteurs_propres_gauches).nombre_lignes *
! 1301: sizeof(struct_complexe16 *))) == NULL)
! 1302: {
! 1303: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1304: return;
! 1305: }
! 1306:
! 1307: if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches)
! 1308: .tableau)[0] = (struct_complexe16 *) malloc(
! 1309: (*s_vecteurs_propres_gauches).nombre_colonnes *
! 1310: sizeof(struct_complexe16))) == NULL)
! 1311: {
! 1312: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1313: return;
! 1314: }
! 1315: }
! 1316:
! 1317: if (calcul_vp_droits == 'V')
! 1318: {
! 1319: free(vpd_f77);
! 1320:
! 1321: (*s_vecteurs_propres_droits).type = 'C';
! 1322: (*s_vecteurs_propres_droits).nombre_lignes = 1;
! 1323: (*s_vecteurs_propres_droits).nombre_colonnes = 1;
! 1324:
! 1325: if (((*s_vecteurs_propres_droits).tableau = malloc(
! 1326: (*s_vecteurs_propres_droits).nombre_lignes *
! 1327: sizeof(struct_complexe16 *))) == NULL)
! 1328: {
! 1329: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1330: return;
! 1331: }
! 1332:
! 1333: if ((((struct_complexe16 **) (*s_vecteurs_propres_droits)
! 1334: .tableau)[0] = (struct_complexe16 *) malloc(
! 1335: (*s_vecteurs_propres_droits).nombre_colonnes *
! 1336: sizeof(struct_complexe16))) == NULL)
! 1337: {
! 1338: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1339: return;
! 1340: }
! 1341: }
! 1342:
! 1343: return;
! 1344: }
! 1345:
! 1346: for(i = 0; i < (*s_valeurs_propres).taille; i++)
! 1347: {
! 1348: if ((beta[i].partie_reelle != 0) || (beta[i].partie_imaginaire != 0))
! 1349: {
! 1350: f77divisioncc_(&(alpha[i]), &(beta[i]), &(((struct_complexe16 *)
! 1351: (*s_valeurs_propres).tableau)[i]));
! 1352: }
! 1353: else
! 1354: {
! 1355: ((struct_complexe16 *) (*s_valeurs_propres).tableau)[i]
! 1356: .partie_reelle = 0;
! 1357: ((struct_complexe16 *) (*s_valeurs_propres).tableau)[i]
! 1358: .partie_imaginaire = 0;
! 1359:
! 1360: (*s_etat_processus).exception = d_ep_division_par_zero;
! 1361: }
! 1362: }
! 1363:
! 1364: free(alpha);
! 1365: free(beta);
! 1366:
! 1367: free(work);
! 1368: free(rwork);
! 1369:
! 1370: if (calcul_vp_gauches == 'V')
! 1371: {
! 1372: (*s_vecteurs_propres_gauches).type = 'C';
! 1373: (*s_vecteurs_propres_gauches).nombre_lignes =
! 1374: (*s_matrice).nombre_lignes;
! 1375: (*s_vecteurs_propres_gauches).nombre_colonnes =
! 1376: (*s_matrice).nombre_colonnes;
! 1377:
! 1378: if (((*s_vecteurs_propres_gauches).tableau = malloc(
! 1379: (*s_vecteurs_propres_gauches).nombre_lignes *
! 1380: sizeof(struct_complexe16 *))) == NULL)
! 1381: {
! 1382: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1383: return;
! 1384: }
! 1385:
! 1386: for(i = 0; i < (*s_vecteurs_propres_gauches).nombre_lignes; i++)
! 1387: {
! 1388: if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches)
! 1389: .tableau)[i] = (struct_complexe16 *) malloc(
! 1390: (*s_vecteurs_propres_gauches).nombre_colonnes *
! 1391: sizeof(struct_complexe16))) == NULL)
! 1392: {
! 1393: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1394: return;
! 1395: }
! 1396: }
! 1397:
! 1398: for(k = 0, i = 0; i < (*s_vecteurs_propres_gauches).nombre_colonnes;
! 1399: i++)
! 1400: {
! 1401: for(j = 0; j < (*s_vecteurs_propres_gauches).nombre_lignes; j++)
! 1402: {
! 1403: ((struct_complexe16 **) (*s_vecteurs_propres_gauches)
! 1404: .tableau)[j][i].partie_reelle =
! 1405: vpg_f77[k].partie_reelle;
! 1406: ((struct_complexe16 **) (*s_vecteurs_propres_gauches)
! 1407: .tableau)[j][i].partie_imaginaire =
! 1408: vpg_f77[k++].partie_imaginaire;
! 1409: }
! 1410: }
! 1411:
! 1412: free(vpg_f77);
! 1413: }
! 1414:
! 1415: if (calcul_vp_droits == 'V')
! 1416: {
! 1417: (*s_vecteurs_propres_droits).type = 'C';
! 1418: (*s_vecteurs_propres_droits).nombre_lignes =
! 1419: (*s_matrice).nombre_lignes;
! 1420: (*s_vecteurs_propres_droits).nombre_colonnes =
! 1421: (*s_matrice).nombre_colonnes;
! 1422:
! 1423: if (((*s_vecteurs_propres_droits).tableau = malloc(
! 1424: (*s_vecteurs_propres_droits).nombre_lignes *
! 1425: sizeof(struct_complexe16 *))) == NULL)
! 1426: {
! 1427: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1428: return;
! 1429: }
! 1430:
! 1431: for(i = 0; i < (*s_vecteurs_propres_droits).nombre_lignes; i++)
! 1432: {
! 1433: if ((((struct_complexe16 **) (*s_vecteurs_propres_droits)
! 1434: .tableau)[i] = (struct_complexe16 *) malloc(
! 1435: (*s_vecteurs_propres_droits).nombre_colonnes *
! 1436: sizeof(struct_complexe16))) == NULL)
! 1437: {
! 1438: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1439: return;
! 1440: }
! 1441: }
! 1442:
! 1443: for(k = 0, i = 0; i < (*s_vecteurs_propres_droits).nombre_colonnes; i++)
! 1444: {
! 1445: for(j = 0; j < (*s_vecteurs_propres_droits).nombre_lignes; j++)
! 1446: {
! 1447: ((struct_complexe16 **) (*s_vecteurs_propres_droits)
! 1448: .tableau)[j][i].partie_reelle =
! 1449: vpd_f77[k].partie_reelle;
! 1450: ((struct_complexe16 **) (*s_vecteurs_propres_droits)
! 1451: .tableau)[j][i].partie_imaginaire =
! 1452: vpd_f77[k++].partie_imaginaire;
! 1453: }
! 1454: }
! 1455:
! 1456: free(vpd_f77);
! 1457: }
! 1458:
! 1459: free(matrice_f77);
! 1460: free(metrique_f77);
! 1461:
! 1462: return;
! 1463: }
! 1464:
! 1465:
! 1466: /*
! 1467: ================================================================================
! 1468: Fonction réalisation le calcul d'un moindres carrés
! 1469: (minimise [B-AX]^2)
! 1470: ================================================================================
! 1471: Entrées : struct_matrice
! 1472: --------------------------------------------------------------------------------
! 1473: Sorties : coefficients dans une struct_matrice allouée au vol
! 1474: --------------------------------------------------------------------------------
! 1475: Effets de bord : néant
! 1476: ================================================================================
! 1477: */
! 1478:
! 1479: void
! 1480: moindres_carres(struct_processus *s_etat_processus,
! 1481: struct_matrice *s_matrice_a, struct_matrice *s_matrice_b,
! 1482: struct_matrice *s_matrice_x)
! 1483: {
! 1484: real8 *work;
! 1485: real8 rcond;
! 1486: real8 *rwork;
! 1487: real8 *vecteur_s;
! 1488:
! 1489: integer4 erreur;
! 1490: integer4 *iwork;
! 1491: integer4 lrwork;
! 1492: integer4 lwork;
! 1493: integer4 nlvl;
! 1494: integer4 rank;
! 1495: integer4 registre_1;
! 1496: integer4 registre_2;
! 1497: integer4 registre_3;
! 1498: integer4 registre_4;
! 1499: integer4 registre_5;
! 1500: integer4 smlsiz;
! 1501:
! 1502: complex16 *cwork;
! 1503:
! 1504: unsigned long i;
! 1505: unsigned long j;
! 1506: unsigned long k;
! 1507: unsigned long taille_matrice_a_f77;
! 1508: unsigned long taille_matrice_b_f77;
! 1509: unsigned long taille_matrice_x_f77;
! 1510:
! 1511: void *matrice_a_f77;
! 1512: void *matrice_b_f77;
! 1513:
! 1514: taille_matrice_a_f77 = (*s_matrice_a).nombre_lignes *
! 1515: (*s_matrice_a).nombre_colonnes;
! 1516: taille_matrice_b_f77 = (*s_matrice_b).nombre_lignes *
! 1517: (*s_matrice_b).nombre_colonnes;
! 1518: taille_matrice_x_f77 = (*s_matrice_b).nombre_colonnes *
! 1519: (*s_matrice_a).nombre_colonnes;
! 1520:
! 1521: /*
! 1522: * Résultat réel
! 1523: */
! 1524:
! 1525: if (((*s_matrice_a).type != 'C') && ((*s_matrice_b).type != 'C'))
! 1526: {
! 1527:
! 1528: /*
! 1529: * Garniture de la matrice A
! 1530: */
! 1531:
! 1532: if ((matrice_a_f77 = (void *) malloc(taille_matrice_a_f77 *
! 1533: sizeof(real8))) == NULL)
! 1534: {
! 1535: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1536: return;
! 1537: }
! 1538:
! 1539: if ((*s_matrice_a).type == 'I')
! 1540: {
! 1541: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
! 1542: {
! 1543: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
! 1544: {
! 1545: ((real8 *) matrice_a_f77)[k++] = ((integer8 **)
! 1546: (*s_matrice_a).tableau)[j][i];
! 1547: }
! 1548: }
! 1549: }
! 1550: else
! 1551: {
! 1552: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
! 1553: {
! 1554: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
! 1555: {
! 1556: ((real8 *) matrice_a_f77)[k++] = ((real8 **)
! 1557: (*s_matrice_a).tableau)[j][i];
! 1558: }
! 1559: }
! 1560: }
! 1561:
! 1562: /*
! 1563: * Garniture de la matrice B
! 1564: */
! 1565:
! 1566: if ((matrice_b_f77 = (void *) malloc(((taille_matrice_b_f77
! 1567: < taille_matrice_x_f77) ? taille_matrice_x_f77
! 1568: : taille_matrice_b_f77) * sizeof(real8))) == NULL)
! 1569: {
! 1570: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1571: return;
! 1572: }
! 1573:
! 1574: if ((*s_matrice_b).type == 'I')
! 1575: {
! 1576: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
! 1577: {
! 1578: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
! 1579: {
! 1580: ((real8 *) matrice_b_f77)[k++] = ((integer8 **)
! 1581: (*s_matrice_b).tableau)[j][i];
! 1582: }
! 1583:
! 1584: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
! 1585: }
! 1586: }
! 1587: else
! 1588: {
! 1589: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
! 1590: {
! 1591: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
! 1592: {
! 1593: ((real8 *) matrice_b_f77)[k++] = ((real8 **)
! 1594: (*s_matrice_b).tableau)[j][i];
! 1595: }
! 1596:
! 1597: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
! 1598: }
! 1599: }
! 1600:
! 1601: rcond = -1;
! 1602:
! 1603: registre_1 = 9;
! 1604: registre_2 = 0;
! 1605:
! 1606: smlsiz = ilaenv_(®istre_1, "DGELSD", " ", ®istre_2,
! 1607: ®istre_2, ®istre_2, ®istre_2, 6, 1);
! 1608:
! 1609: nlvl = 1 + ((integer4) log(((real8) (((*s_matrice_a).nombre_lignes <
! 1610: (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes
! 1611: : (*s_matrice_a).nombre_colonnes)) / (smlsiz + 1)) /
! 1612: log((real8) 2));
! 1613:
! 1614: if ((iwork = (integer4 *) malloc(((((*s_matrice_a).nombre_lignes <
! 1615: (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes
! 1616: : (*s_matrice_a).nombre_colonnes) * (11 + (3 * nlvl))) *
! 1617: sizeof(integer4))) == NULL)
! 1618: {
! 1619: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1620: return;
! 1621: }
! 1622:
! 1623: registre_1 = (integer4) (*s_matrice_a).nombre_lignes;
! 1624: registre_2 = (integer4) (*s_matrice_a).nombre_colonnes;
! 1625: registre_3 = (integer4) (*s_matrice_b).nombre_colonnes;
! 1626: registre_4 = (registre_1 > registre_2) ? registre_1 : registre_2;
! 1627: registre_5 = (registre_1 > registre_2) ? registre_2 : registre_1;
! 1628:
! 1629: if ((work = malloc(sizeof(real8))) == NULL)
! 1630: {
! 1631: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1632: return;
! 1633: }
! 1634:
! 1635: if ((vecteur_s = (real8 *) malloc(registre_5 * sizeof(real8)))
! 1636: == NULL)
! 1637: {
! 1638: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1639: return;
! 1640: }
! 1641:
! 1642: lwork = -1;
! 1643:
! 1644: dgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77,
! 1645: ®istre_1, matrice_b_f77, ®istre_4, vecteur_s,
! 1646: &rcond, &rank, work, &lwork, iwork, &erreur);
! 1647:
! 1648: if (erreur != 0)
! 1649: {
! 1650: if (erreur > 0)
! 1651: {
! 1652: (*s_etat_processus).exception = d_ep_decomposition_SVD;
! 1653: }
! 1654: else
! 1655: {
! 1656: (*s_etat_processus).erreur_execution =
! 1657: d_ex_routines_mathematiques;
! 1658: }
! 1659:
! 1660: free(work);
! 1661: free(iwork);
! 1662: free(matrice_a_f77);
! 1663: free(matrice_b_f77);
! 1664: return;
! 1665: }
! 1666:
! 1667: lwork = (integer4) work[0];
! 1668: free(work);
! 1669:
! 1670: if ((work = malloc(lwork * sizeof(real8))) == NULL)
! 1671: {
! 1672: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1673: return;
! 1674: }
! 1675:
! 1676: dgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77,
! 1677: ®istre_1, matrice_b_f77, ®istre_4, vecteur_s,
! 1678: &rcond, &rank, work, &lwork, iwork, &erreur);
! 1679:
! 1680: free(vecteur_s);
! 1681:
! 1682: if (erreur != 0)
! 1683: {
! 1684: if (erreur > 0)
! 1685: {
! 1686: (*s_etat_processus).exception = d_ep_decomposition_SVD;
! 1687: }
! 1688: else
! 1689: {
! 1690: (*s_etat_processus).erreur_execution =
! 1691: d_ex_routines_mathematiques;
! 1692: }
! 1693:
! 1694: free(work);
! 1695: free(iwork);
! 1696: free(matrice_a_f77);
! 1697: free(matrice_b_f77);
! 1698: return;
! 1699: }
! 1700:
! 1701: free(work);
! 1702: free(iwork);
! 1703:
! 1704: (*s_matrice_x).nombre_lignes = (*s_matrice_a).nombre_colonnes;
! 1705: (*s_matrice_x).nombre_colonnes = (*s_matrice_b).nombre_colonnes;
! 1706:
! 1707: if (((*s_matrice_x).tableau = malloc((*s_matrice_x).nombre_lignes *
! 1708: sizeof(real8 *))) == NULL)
! 1709: {
! 1710: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1711: return;
! 1712: }
! 1713:
! 1714: for(j = 0; j < (*s_matrice_x).nombre_lignes; j++)
! 1715: {
! 1716: if ((((real8 **) (*s_matrice_x).tableau)[j] = (real8 *)
! 1717: malloc((*s_matrice_x).nombre_colonnes *
! 1718: sizeof(real8)))== NULL)
! 1719: {
! 1720: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1721: return;
! 1722: }
! 1723: }
! 1724:
! 1725: for(k = 0, i = 0; i < (*s_matrice_x).nombre_colonnes; i++)
! 1726: {
! 1727: for(j = 0; j < (*s_matrice_x).nombre_lignes; j++)
! 1728: {
! 1729: (((real8 **) (*s_matrice_x).tableau)[j][i]) = (((real8 *)
! 1730: matrice_b_f77)[k++]);
! 1731: }
! 1732:
! 1733: for(; j < (*s_matrice_b).nombre_lignes; j++, k++);
! 1734: }
! 1735: }
! 1736:
! 1737: /*
! 1738: * Résultat complexe
! 1739: */
! 1740:
! 1741: else
! 1742: {
! 1743:
! 1744: /*
! 1745: * Garniture de la matrice A
! 1746: */
! 1747:
! 1748: if ((matrice_a_f77 = (void *) malloc(taille_matrice_a_f77 *
! 1749: sizeof(struct_complexe16))) == NULL)
! 1750: {
! 1751: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1752: return;
! 1753: }
! 1754:
! 1755: if ((*s_matrice_a).type == 'I')
! 1756: {
! 1757: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
! 1758: {
! 1759: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
! 1760: {
! 1761: ((struct_complexe16 *) matrice_a_f77)[k].partie_reelle =
! 1762: (real8) ((integer8 **) (*s_matrice_a)
! 1763: .tableau)[j][i];
! 1764: ((struct_complexe16 *) matrice_a_f77)[k++]
! 1765: .partie_imaginaire = (real8) 0;
! 1766: }
! 1767: }
! 1768: }
! 1769: else if ((*s_matrice_a).type == 'R')
! 1770: {
! 1771: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
! 1772: {
! 1773: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
! 1774: {
! 1775: ((struct_complexe16 *) matrice_a_f77)[k].partie_reelle =
! 1776: ((real8 **) (*s_matrice_a).tableau)[j][i];
! 1777: ((struct_complexe16 *) matrice_a_f77)[k++]
! 1778: .partie_imaginaire = (real8) 0;
! 1779: }
! 1780: }
! 1781: }
! 1782: else
! 1783: {
! 1784: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
! 1785: {
! 1786: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
! 1787: {
! 1788: ((struct_complexe16 *) matrice_a_f77)[k].partie_reelle =
! 1789: ((struct_complexe16 **) (*s_matrice_a).tableau)
! 1790: [j][i].partie_reelle;
! 1791: ((struct_complexe16 *) matrice_a_f77)[k++]
! 1792: .partie_imaginaire = ((struct_complexe16 **)
! 1793: (*s_matrice_a).tableau)[j][i].partie_imaginaire;
! 1794: }
! 1795: }
! 1796: }
! 1797:
! 1798: /*
! 1799: * Garniture de la matrice B
! 1800: */
! 1801:
! 1802: if ((matrice_b_f77 = (void *) malloc(((taille_matrice_b_f77
! 1803: < taille_matrice_x_f77) ? taille_matrice_x_f77
! 1804: : taille_matrice_b_f77) * sizeof(struct_complexe16))) == NULL)
! 1805: {
! 1806: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1807: return;
! 1808: }
! 1809:
! 1810: if ((*s_matrice_b).type == 'I')
! 1811: {
! 1812: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
! 1813: {
! 1814: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
! 1815: {
! 1816: ((struct_complexe16 *) matrice_b_f77)[k].partie_reelle =
! 1817: (real8) ((integer8 **) (*s_matrice_b)
! 1818: .tableau)[j][i];
! 1819: ((struct_complexe16 *) matrice_b_f77)[k++]
! 1820: .partie_imaginaire = (real8) 0;
! 1821: }
! 1822:
! 1823: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
! 1824: }
! 1825: }
! 1826: else if ((*s_matrice_b).type == 'R')
! 1827: {
! 1828: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
! 1829: {
! 1830: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
! 1831: {
! 1832: ((struct_complexe16 *) matrice_b_f77)[k].partie_reelle =
! 1833: ((real8 **) (*s_matrice_b).tableau)[j][i];
! 1834: ((struct_complexe16 *) matrice_b_f77)[k++]
! 1835: .partie_imaginaire = (real8) 0;
! 1836: }
! 1837:
! 1838: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
! 1839: }
! 1840: }
! 1841: else
! 1842: {
! 1843: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
! 1844: {
! 1845: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
! 1846: {
! 1847: ((struct_complexe16 *) matrice_b_f77)[k].partie_reelle =
! 1848: ((struct_complexe16 **) (*s_matrice_b).tableau)
! 1849: [j][i].partie_reelle;
! 1850: ((struct_complexe16 *) matrice_b_f77)[k++]
! 1851: .partie_imaginaire = ((struct_complexe16 **)
! 1852: (*s_matrice_b).tableau)[j][i].partie_imaginaire;
! 1853: }
! 1854:
! 1855: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
! 1856: }
! 1857: }
! 1858:
! 1859: rcond = -1;
! 1860:
! 1861: registre_1 = 9;
! 1862: registre_2 = 0;
! 1863:
! 1864: smlsiz = ilaenv_(®istre_1, "ZGELSD", " ", ®istre_2,
! 1865: ®istre_2, ®istre_2, ®istre_2, 6, 1);
! 1866:
! 1867: nlvl = 1 + ((integer4) log(((real8) (((*s_matrice_a).nombre_lignes <
! 1868: (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes
! 1869: : (*s_matrice_a).nombre_colonnes)) / (smlsiz + 1))
! 1870: / log((real8) 2));
! 1871:
! 1872: if ((*s_matrice_a).nombre_lignes >= (*s_matrice_a).nombre_colonnes)
! 1873: {
! 1874: lrwork = (integer4) ((*s_matrice_a).nombre_colonnes * (8 + (2 *
! 1875: smlsiz) + (8 * nlvl) + (*s_matrice_b).nombre_colonnes));
! 1876: }
! 1877: else
! 1878: {
! 1879: lrwork = (integer4) ((*s_matrice_a).nombre_lignes * (8 + (2 *
! 1880: smlsiz) + (8 * nlvl) + (*s_matrice_b).nombre_colonnes));
! 1881: }
! 1882:
! 1883: if ((rwork = (real8 *) malloc(lrwork * sizeof(real8))) == NULL)
! 1884: {
! 1885: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1886: return;
! 1887: }
! 1888:
! 1889: if ((iwork = (integer4 *) malloc(((((*s_matrice_a).nombre_lignes <
! 1890: (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes
! 1891: : (*s_matrice_a).nombre_colonnes) * (11 + (3 * nlvl))) *
! 1892: sizeof(integer4))) == NULL)
! 1893: {
! 1894: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1895: return;
! 1896: }
! 1897:
! 1898: registre_1 = (integer4) (*s_matrice_a).nombre_lignes;
! 1899: registre_2 = (integer4) (*s_matrice_a).nombre_colonnes;
! 1900: registre_3 = (integer4) (*s_matrice_b).nombre_colonnes;
! 1901: registre_4 = (registre_1 > registre_2) ? registre_1 : registre_2;
! 1902: registre_5 = (registre_1 > registre_2) ? registre_2 : registre_1;
! 1903:
! 1904: if ((cwork = malloc(sizeof(struct_complexe16))) == NULL)
! 1905: {
! 1906: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1907: return;
! 1908: }
! 1909:
! 1910: if ((vecteur_s = (real8 *) malloc(registre_5 * sizeof(real8)))
! 1911: == NULL)
! 1912: {
! 1913: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1914: return;
! 1915: }
! 1916:
! 1917: lwork = -1;
! 1918:
! 1919: zgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77,
! 1920: ®istre_1, matrice_b_f77, ®istre_4, vecteur_s,
! 1921: &rcond, &rank, cwork, &lwork, rwork, iwork, &erreur);
! 1922:
! 1923: if (erreur != 0)
! 1924: {
! 1925: if (erreur > 0)
! 1926: {
! 1927: (*s_etat_processus).exception = d_ep_decomposition_SVD;
! 1928: }
! 1929: else
! 1930: {
! 1931: (*s_etat_processus).erreur_execution =
! 1932: d_ex_routines_mathematiques;
! 1933: }
! 1934:
! 1935: free(cwork);
! 1936: free(iwork);
! 1937: free(rwork);
! 1938: free(matrice_a_f77);
! 1939: free(matrice_b_f77);
! 1940: return;
! 1941: }
! 1942:
! 1943: lwork = (integer4) cwork[0].partie_reelle;
! 1944: free(cwork);
! 1945:
! 1946: if ((cwork = malloc(lwork * sizeof(struct_complexe16))) == NULL)
! 1947: {
! 1948: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1949: return;
! 1950: }
! 1951:
! 1952: zgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77,
! 1953: ®istre_1, matrice_b_f77, ®istre_4, vecteur_s,
! 1954: &rcond, &rank, cwork, &lwork, rwork, iwork, &erreur);
! 1955:
! 1956: free(vecteur_s);
! 1957:
! 1958: if (erreur != 0)
! 1959: {
! 1960: if (erreur > 0)
! 1961: {
! 1962: (*s_etat_processus).exception = d_ep_decomposition_SVD;
! 1963: }
! 1964: else
! 1965: {
! 1966: (*s_etat_processus).erreur_execution =
! 1967: d_ex_routines_mathematiques;
! 1968: }
! 1969:
! 1970: free(cwork);
! 1971: free(iwork);
! 1972: free(rwork);
! 1973: free(matrice_a_f77);
! 1974: free(matrice_b_f77);
! 1975: return;
! 1976: }
! 1977:
! 1978: free(cwork);
! 1979: free(iwork);
! 1980: free(rwork);
! 1981:
! 1982: (*s_matrice_x).nombre_lignes = (*s_matrice_a).nombre_colonnes;
! 1983: (*s_matrice_x).nombre_colonnes = (*s_matrice_b).nombre_colonnes;
! 1984:
! 1985: if (((*s_matrice_x).tableau = malloc((*s_matrice_x).nombre_lignes *
! 1986: sizeof(struct_complexe16 *))) == NULL)
! 1987: {
! 1988: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1989: return;
! 1990: }
! 1991:
! 1992: for(j = 0; j < (*s_matrice_x).nombre_lignes; j++)
! 1993: {
! 1994: if ((((struct_complexe16 **) (*s_matrice_x).tableau)[j] =
! 1995: (struct_complexe16 *) malloc((*s_matrice_x).nombre_colonnes
! 1996: * sizeof(struct_complexe16)))== NULL)
! 1997: {
! 1998: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1999: return;
! 2000: }
! 2001: }
! 2002:
! 2003: for(k = 0, i = 0; i < (*s_matrice_x).nombre_colonnes; i++)
! 2004: {
! 2005: for(j = 0; j < (*s_matrice_x).nombre_lignes; j++)
! 2006: {
! 2007: (((struct_complexe16 **) (*s_matrice_x).tableau)[j][i])
! 2008: .partie_reelle = (((struct_complexe16 *)
! 2009: matrice_b_f77)[k]).partie_reelle;
! 2010: (((struct_complexe16 **) (*s_matrice_x).tableau)[j][i])
! 2011: .partie_imaginaire = (((struct_complexe16 *)
! 2012: matrice_b_f77)[k++]).partie_imaginaire;
! 2013: }
! 2014:
! 2015: for(; j < (*s_matrice_b).nombre_lignes; j++, k++);
! 2016: }
! 2017: }
! 2018:
! 2019: free(matrice_a_f77);
! 2020: free(matrice_b_f77);
! 2021:
! 2022: return;
! 2023: }
! 2024:
! 2025: // vim: ts=4
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to algebre_lineaire1.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / src