![]() ![]() | ![]() |
Passage de la branche 4.1 en branche stable.
1: /* 2: ================================================================================ 3: RPL/2 (R) version 4.1.0 4: Copyright (C) 1989-2011 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