![]() ![]() | ![]() |
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 la factorisation de Schur d'une matrice carrée 29: ================================================================================ 30: Entrées : struct_matrice 31: -------------------------------------------------------------------------------- 32: Sorties : décomposition de Schur de la matrice d'entrée et drapeau d'erreur. 33: La matrice en entrée est écrasée. La matrice de sortie est 34: la forme de Schur. 35: La routine renvoie aussi une matrice de complexes correspondant 36: aux vecteurs de Schur. Cette matrice est allouée par 37: la routine et vaut NULL sinon. 38: -------------------------------------------------------------------------------- 39: Effets de bord : néant 40: ================================================================================ 41: */ 42: 43: void 44: factorisation_schur(struct_processus *s_etat_processus, 45: struct_matrice *s_matrice, struct_matrice **s_schur) 46: { 47: complex16 *w; 48: 49: integer4 info; 50: integer4 lwork; 51: integer4 nombre_colonnes_a; 52: integer4 nombre_lignes_a; 53: integer4 sdim; 54: 55: real8 *rwork; 56: real8 *wi; 57: real8 *wr; 58: 59: unsigned char calcul_vecteurs_schur; 60: unsigned char tri_vecteurs_schur; 61: 62: unsigned long i; 63: unsigned long j; 64: unsigned long k; 65: unsigned long taille_matrice_f77; 66: 67: void *matrice_a_f77; 68: void *matrice_vs_f77; 69: void *tampon; 70: void *work; 71: 72: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes; 73: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes; 74: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a; 75: 76: calcul_vecteurs_schur = 'V'; 77: tri_vecteurs_schur = 'N'; 78: 79: switch((*s_matrice).type) 80: { 81: case 'I' : 82: { 83: /* Conversion de la matrice en matrice réelle */ 84: 85: for(i = 0; i < (unsigned long) nombre_lignes_a; i++) 86: { 87: tampon = (*s_matrice).tableau[i]; 88: 89: if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *) 90: malloc(nombre_colonnes_a * sizeof(real8))) == NULL) 91: { 92: (*s_etat_processus).erreur_systeme = 93: d_es_allocation_memoire; 94: return; 95: } 96: 97: for(j = 0; j < (unsigned long) nombre_colonnes_a; j++) 98: { 99: ((real8 **) (*s_matrice).tableau)[i][j] = 100: ((integer8 *) tampon)[j]; 101: } 102: 103: free(tampon); 104: } 105: 106: (*s_matrice).type = 'R'; 107: } 108: 109: case 'R' : 110: { 111: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 * 112: sizeof(real8))) == NULL) 113: { 114: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 115: return; 116: } 117: 118: if ((matrice_vs_f77 = (void *) malloc(taille_matrice_f77 * 119: sizeof(real8))) == NULL) 120: { 121: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 122: return; 123: } 124: 125: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 126: { 127: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 128: { 129: ((real8 *) matrice_a_f77)[k++] = ((real8 **) 130: (*s_matrice).tableau)[j][i]; 131: } 132: } 133: 134: if ((wr = (real8 *) malloc(nombre_lignes_a * sizeof(real8))) 135: == NULL) 136: { 137: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 138: return; 139: } 140: 141: if ((wi = (real8 *) malloc(nombre_lignes_a * sizeof(real8))) 142: == NULL) 143: { 144: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 145: return; 146: } 147: 148: lwork = -1; 149: 150: if ((work = (real8 *) malloc(sizeof(real8))) == NULL) 151: { 152: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 153: return; 154: } 155: 156: dgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur, 157: NULL, &nombre_lignes_a, matrice_a_f77, 158: &nombre_colonnes_a, &sdim, wr, wi, 159: matrice_vs_f77, &nombre_colonnes_a, 160: work, &lwork, NULL, &info, 1, 1); 161: 162: lwork = ((real8 *) work)[0]; 163: free(work); 164: 165: if ((work = (real8 *) malloc(lwork * sizeof(real8))) == NULL) 166: { 167: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 168: return; 169: } 170: 171: dgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur, 172: NULL, &nombre_lignes_a, matrice_a_f77, 173: &nombre_colonnes_a, &sdim, wr, wi, 174: matrice_vs_f77, &nombre_colonnes_a, 175: work, &lwork, NULL, &info, 1, 1); 176: 177: free(work); 178: free(wr); 179: free(wi); 180: 181: if (info != 0) 182: { 183: if (info > 0) 184: { 185: (*s_etat_processus).exception = d_ep_decomposition_QR; 186: } 187: else 188: { 189: (*s_etat_processus).erreur_execution = 190: d_ex_routines_mathematiques; 191: } 192: 193: free(matrice_a_f77); 194: free(matrice_vs_f77); 195: return; 196: } 197: 198: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 199: { 200: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 201: { 202: ((real8 **) (*s_matrice).tableau)[j][i] = 203: ((real8 *) matrice_a_f77)[k++]; 204: } 205: } 206: 207: (**s_schur).nombre_colonnes = (*s_matrice).nombre_colonnes; 208: (**s_schur).nombre_lignes = (*s_matrice).nombre_lignes; 209: (**s_schur).type = 'R'; 210: 211: if (((**s_schur).tableau = malloc((**s_schur) 212: .nombre_lignes * sizeof(real8 *))) == NULL) 213: { 214: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 215: return; 216: } 217: 218: for(i = 0; i < (**s_schur).nombre_lignes; i++) 219: { 220: if ((((real8 **) (**s_schur).tableau)[i] = (real8 *) 221: malloc((**s_schur).nombre_colonnes * 222: sizeof(real8))) == NULL) 223: { 224: (*s_etat_processus).erreur_systeme = 225: d_es_allocation_memoire; 226: return; 227: } 228: } 229: 230: for(k = 0, i = 0; i < (**s_schur).nombre_colonnes; i++) 231: { 232: for(j = 0; j < (**s_schur).nombre_lignes; j++) 233: { 234: ((real8 **) (**s_schur).tableau)[j][i] = ((real8 *) 235: matrice_vs_f77)[k++]; 236: } 237: } 238: 239: free(matrice_a_f77); 240: free(matrice_vs_f77); 241: 242: break; 243: } 244: 245: case 'C' : 246: { 247: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 * 248: sizeof(complex16))) == NULL) 249: { 250: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 251: return; 252: } 253: 254: if ((matrice_vs_f77 = (void *) malloc(taille_matrice_f77 * 255: sizeof(complex16))) == NULL) 256: { 257: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 258: return; 259: } 260: 261: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 262: { 263: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 264: { 265: ((complex16 *) matrice_a_f77)[k].partie_reelle = 266: ((complex16 **) (*s_matrice).tableau)[j][i] 267: .partie_reelle; 268: ((complex16 *) matrice_a_f77)[k++].partie_imaginaire = 269: ((complex16 **) (*s_matrice).tableau)[j][i] 270: .partie_imaginaire; 271: } 272: } 273: 274: if ((w = (complex16 *) malloc(nombre_lignes_a * sizeof(complex16))) 275: == NULL) 276: { 277: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 278: return; 279: } 280: 281: lwork = -1; 282: 283: if ((work = (complex16 *) malloc(sizeof(complex16))) == NULL) 284: { 285: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 286: return; 287: } 288: 289: if ((rwork = (real8 *) malloc(nombre_lignes_a * sizeof(real8))) 290: == NULL) 291: { 292: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 293: return; 294: } 295: 296: zgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur, 297: NULL, &nombre_lignes_a, matrice_a_f77, 298: &nombre_colonnes_a, &sdim, w, 299: matrice_vs_f77, &nombre_colonnes_a, 300: work, &lwork, rwork, NULL, &info, 1, 1); 301: 302: lwork = ((complex16 *) work)[0].partie_reelle; 303: free(work); 304: 305: if ((work = (complex16 *) malloc(lwork * sizeof(complex16))) 306: == NULL) 307: { 308: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 309: return; 310: } 311: 312: zgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur, 313: NULL, &nombre_lignes_a, matrice_a_f77, 314: &nombre_colonnes_a, &sdim, w, 315: matrice_vs_f77, &nombre_colonnes_a, 316: work, &lwork, rwork, NULL, &info, 1, 1); 317: 318: free(work); 319: free(rwork); 320: free(w); 321: 322: if (info != 0) 323: { 324: if (info > 0) 325: { 326: (*s_etat_processus).exception = d_ep_decomposition_QR; 327: } 328: else 329: { 330: (*s_etat_processus).erreur_execution = 331: d_ex_routines_mathematiques; 332: } 333: 334: free(matrice_a_f77); 335: free(matrice_vs_f77); 336: return; 337: } 338: 339: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 340: { 341: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 342: { 343: ((complex16 **) (*s_matrice).tableau)[j][i] 344: .partie_reelle = ((complex16 *) matrice_a_f77)[k] 345: .partie_reelle; 346: ((complex16 **) (*s_matrice).tableau)[j][i] 347: .partie_imaginaire = ((complex16 *) matrice_a_f77) 348: [k++].partie_imaginaire; 349: } 350: } 351: 352: (**s_schur).nombre_colonnes = (*s_matrice).nombre_colonnes; 353: (**s_schur).nombre_lignes = (*s_matrice).nombre_lignes; 354: (**s_schur).type = 'C'; 355: 356: if (((**s_schur).tableau = malloc((**s_schur) 357: .nombre_lignes * sizeof(complex16 *))) == NULL) 358: { 359: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 360: return; 361: } 362: 363: for(i = 0; i < (**s_schur).nombre_lignes; i++) 364: { 365: if ((((complex16 **) (**s_schur).tableau)[i] = (complex16 *) 366: malloc((**s_schur).nombre_colonnes * 367: sizeof(complex16))) == NULL) 368: { 369: (*s_etat_processus).erreur_systeme = 370: d_es_allocation_memoire; 371: return; 372: } 373: } 374: 375: for(k = 0, i = 0; i < (**s_schur).nombre_colonnes; i++) 376: { 377: for(j = 0; j < (**s_schur).nombre_lignes; j++) 378: { 379: ((complex16 **) (**s_schur).tableau)[j][i].partie_reelle = 380: ((complex16 *) matrice_vs_f77)[k].partie_reelle; 381: ((complex16 **) (**s_schur).tableau)[j][i] 382: .partie_imaginaire = ((complex16 *) matrice_vs_f77) 383: [k++].partie_imaginaire; 384: } 385: } 386: 387: free(matrice_a_f77); 388: free(matrice_vs_f77); 389: 390: break; 391: } 392: } 393: 394: return; 395: } 396: 397: 398: /* 399: ================================================================================ 400: Fonction réalisation la factorisation LQ d'une matrice quelconque 401: ================================================================================ 402: Entrées : struct_matrice 403: -------------------------------------------------------------------------------- 404: Sorties : décomposition de LQ de la matrice d'entrée. Le tableau tau 405: est initialisé par la fonction 406: -------------------------------------------------------------------------------- 407: Effets de bord : néant 408: ================================================================================ 409: */ 410: 411: void 412: factorisation_lq(struct_processus *s_etat_processus, struct_matrice *s_matrice, 413: void **tau) 414: { 415: integer4 nombre_colonnes_a; 416: integer4 nombre_lignes_a; 417: integer4 erreur; 418: 419: unsigned long i; 420: unsigned long j; 421: unsigned long k; 422: unsigned long taille_matrice_f77; 423: 424: void *matrice_a_f77; 425: void *tampon; 426: void *work; 427: 428: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes; 429: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes; 430: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a; 431: 432: switch((*s_matrice).type) 433: { 434: case 'I' : 435: { 436: /* Conversion de la matrice en matrice réelle */ 437: 438: for(i = 0; i < (unsigned long) nombre_lignes_a; i++) 439: { 440: tampon = (*s_matrice).tableau[i]; 441: 442: if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *) 443: malloc(nombre_colonnes_a * sizeof(real8))) == NULL) 444: { 445: (*s_etat_processus).erreur_systeme = 446: d_es_allocation_memoire; 447: return; 448: } 449: 450: for(j = 0; j < (unsigned long) nombre_colonnes_a; j++) 451: { 452: ((real8 **) (*s_matrice).tableau)[i][j] = 453: ((integer8 *) tampon)[j]; 454: } 455: 456: free(tampon); 457: } 458: 459: (*s_matrice).type = 'R'; 460: } 461: 462: case 'R' : 463: { 464: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 * 465: sizeof(real8))) == NULL) 466: { 467: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 468: return; 469: } 470: 471: if (((*tau) = malloc(((nombre_colonnes_a < nombre_lignes_a) 472: ? nombre_colonnes_a : nombre_lignes_a) * sizeof(real8))) 473: == NULL) 474: { 475: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 476: return; 477: } 478: 479: if ((work = malloc(nombre_lignes_a * sizeof(real8))) == NULL) 480: { 481: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 482: return; 483: } 484: 485: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 486: { 487: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 488: { 489: ((real8 *) matrice_a_f77)[k++] = ((real8 **) 490: (*s_matrice).tableau)[j][i]; 491: } 492: } 493: 494: dgelq2_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77, 495: &nombre_lignes_a, (*((real8 **) tau)), work, &erreur); 496: 497: if (erreur != 0) 498: { 499: // L'erreur ne peut être que négative. 500: 501: (*s_etat_processus).erreur_execution = 502: d_ex_routines_mathematiques; 503: free(work); 504: free(matrice_a_f77); 505: return; 506: } 507: 508: for(k = 0, i = 0; i < (unsigned long) nombre_colonnes_a; i++) 509: { 510: for(j = 0; j < (unsigned long) nombre_lignes_a; j++) 511: { 512: ((real8 **) (*s_matrice).tableau)[j][i] = ((real8 *) 513: matrice_a_f77)[k++]; 514: } 515: } 516: 517: free(work); 518: free(matrice_a_f77); 519: 520: break; 521: } 522: 523: case 'C' : 524: { 525: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 * 526: sizeof(complex16))) == NULL) 527: { 528: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 529: return; 530: } 531: 532: if (((*tau) = malloc(((nombre_colonnes_a < nombre_lignes_a) 533: ? nombre_colonnes_a : nombre_lignes_a) * 534: sizeof(complex16))) == NULL) 535: { 536: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 537: return; 538: } 539: 540: if ((work = malloc(nombre_lignes_a * sizeof(complex16))) == NULL) 541: { 542: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 543: return; 544: } 545: 546: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 547: { 548: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 549: { 550: ((complex16 *) matrice_a_f77)[k].partie_reelle = 551: ((complex16 **) (*s_matrice).tableau)[j][i] 552: .partie_reelle; 553: ((complex16 *) matrice_a_f77)[k++].partie_imaginaire = 554: ((complex16 **) (*s_matrice).tableau)[j][i] 555: .partie_imaginaire; 556: } 557: } 558: 559: zgelq2_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77, 560: &nombre_lignes_a, (*((complex16 **) tau)), work, &erreur); 561: 562: if (erreur != 0) 563: { 564: // L'erreur ne peut être que négative. 565: 566: (*s_etat_processus).erreur_execution = 567: d_ex_routines_mathematiques; 568: free(work); 569: free(matrice_a_f77); 570: return; 571: } 572: 573: for(k = 0, i = 0; i < (unsigned long) nombre_colonnes_a; i++) 574: { 575: for(j = 0; j < (unsigned long) nombre_lignes_a; j++) 576: { 577: ((complex16 **) (*s_matrice).tableau)[j][i].partie_reelle = 578: ((complex16 *) matrice_a_f77)[k].partie_reelle; 579: ((complex16 **) (*s_matrice).tableau)[j][i] 580: .partie_imaginaire = ((complex16 *) matrice_a_f77) 581: [k++].partie_imaginaire; 582: } 583: } 584: 585: free(work); 586: free(matrice_a_f77); 587: 588: break; 589: } 590: } 591: 592: return; 593: } 594: 595: 596: /* 597: ================================================================================ 598: Fonction réalisation la factorisation QR d'une matrice quelconque 599: ================================================================================ 600: Entrées : struct_matrice 601: -------------------------------------------------------------------------------- 602: Sorties : décomposition de QR de la matrice d'entrée. Le tableau tau 603: est initialisé par la fonction 604: -------------------------------------------------------------------------------- 605: Effets de bord : néant 606: ================================================================================ 607: */ 608: 609: void 610: factorisation_qr(struct_processus *s_etat_processus, struct_matrice *s_matrice, 611: void **tau) 612: { 613: integer4 lwork; 614: integer4 nombre_colonnes_a; 615: integer4 nombre_lignes_a; 616: integer4 erreur; 617: integer4 *pivot; 618: 619: real8 *rwork; 620: 621: unsigned long i; 622: unsigned long j; 623: unsigned long k; 624: unsigned long taille_matrice_f77; 625: 626: void *matrice_a_f77; 627: void *registre; 628: void *tampon; 629: void *work; 630: 631: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes; 632: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes; 633: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a; 634: 635: switch((*s_matrice).type) 636: { 637: case 'I' : 638: { 639: /* Conversion de la matrice en matrice réelle */ 640: 641: for(i = 0; i < (unsigned long) nombre_lignes_a; i++) 642: { 643: tampon = (*s_matrice).tableau[i]; 644: 645: if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *) 646: malloc(nombre_colonnes_a * sizeof(real8))) == NULL) 647: { 648: (*s_etat_processus).erreur_systeme = 649: d_es_allocation_memoire; 650: return; 651: } 652: 653: for(j = 0; j < (unsigned long) nombre_colonnes_a; j++) 654: { 655: ((real8 **) (*s_matrice).tableau)[i][j] = 656: ((integer8 *) tampon)[j]; 657: } 658: 659: free(tampon); 660: } 661: 662: (*s_matrice).type = 'R'; 663: } 664: 665: case 'R' : 666: { 667: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 * 668: sizeof(real8))) == NULL) 669: { 670: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 671: return; 672: } 673: 674: if (((*tau) = malloc(((nombre_colonnes_a < nombre_lignes_a) 675: ? nombre_colonnes_a : nombre_lignes_a) * sizeof(real8))) 676: == NULL) 677: { 678: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 679: return; 680: } 681: 682: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 683: { 684: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 685: { 686: ((real8 *) matrice_a_f77)[k++] = ((real8 **) 687: (*s_matrice).tableau)[j][i]; 688: } 689: } 690: 691: if ((pivot = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL) 692: { 693: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 694: return; 695: } 696: 697: for(i = 0; i < (unsigned long) nombre_colonnes_a; pivot[i++] = 0); 698: 699: lwork = -1; 700: 701: if ((work = malloc(sizeof(real8))) == NULL) 702: { 703: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 704: return; 705: } 706: 707: // Calcul de la taille de l'espace 708: 709: dgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77, 710: &nombre_lignes_a, pivot, (*((real8 **) tau)), 711: work, &lwork, &erreur); 712: 713: lwork = ((real8 *) work)[0]; 714: 715: free(work); 716: 717: if ((work = malloc(lwork * sizeof(real8))) == NULL) 718: { 719: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 720: return; 721: } 722: 723: // Calcul de la permutation 724: 725: if ((registre = (void *) malloc(taille_matrice_f77 * 726: sizeof(real8))) == NULL) 727: { 728: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 729: return; 730: } 731: 732: memcpy(registre, matrice_a_f77, taille_matrice_f77 * sizeof(real8)); 733: 734: dgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, registre, 735: &nombre_lignes_a, pivot, (*((real8 **) tau)), 736: work, &lwork, &erreur); 737: 738: free(registre); 739: 740: if (erreur != 0) 741: { 742: // L'erreur ne peut être que négative. 743: 744: (*s_etat_processus).erreur_execution = 745: d_ex_routines_mathematiques; 746: free(work); 747: free(matrice_a_f77); 748: free(tau); 749: return; 750: } 751: 752: // La permutation doit maintenant être unitaire 753: 754: dgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77, 755: &nombre_lignes_a, pivot, (*((real8 **) tau)), 756: work, &lwork, &erreur); 757: 758: if (erreur != 0) 759: { 760: // L'erreur ne peut être que négative. 761: 762: (*s_etat_processus).erreur_execution = 763: d_ex_routines_mathematiques; 764: free(work); 765: free(matrice_a_f77); 766: free(tau); 767: return; 768: } 769: 770: for(i = 0; i < (unsigned long) nombre_colonnes_a; i++) 771: { 772: if ((i + 1) != (unsigned long) pivot[i]) 773: { 774: (*s_etat_processus).erreur_execution = 775: d_ex_routines_mathematiques; 776: 777: free(pivot); 778: free(work); 779: free(matrice_a_f77); 780: free(tau); 781: 782: return; 783: } 784: } 785: 786: free(pivot); 787: 788: for(k = 0, i = 0; i < (unsigned long) nombre_colonnes_a; i++) 789: { 790: for(j = 0; j < (unsigned long) nombre_lignes_a; j++) 791: { 792: ((real8 **) (*s_matrice).tableau)[j][i] = ((real8 *) 793: matrice_a_f77)[k++]; 794: } 795: } 796: 797: free(work); 798: free(matrice_a_f77); 799: 800: break; 801: } 802: 803: case 'C' : 804: { 805: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 * 806: sizeof(complex16))) == NULL) 807: { 808: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 809: return; 810: } 811: 812: if (((*tau) = malloc(((nombre_colonnes_a < nombre_lignes_a) 813: ? nombre_colonnes_a : nombre_lignes_a) * sizeof(complex16))) 814: == NULL) 815: { 816: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 817: return; 818: } 819: 820: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 821: { 822: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 823: { 824: ((complex16 *) matrice_a_f77)[k].partie_reelle = 825: ((complex16 **) (*s_matrice).tableau)[j][i] 826: .partie_reelle; 827: ((complex16 *) matrice_a_f77)[k++].partie_imaginaire = 828: ((complex16 **) (*s_matrice).tableau)[j][i] 829: .partie_imaginaire; 830: } 831: } 832: 833: if ((pivot = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL) 834: { 835: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 836: return; 837: } 838: 839: if ((rwork = malloc(2 * nombre_colonnes_a * sizeof(real8))) == NULL) 840: { 841: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 842: return; 843: } 844: 845: for(i = 0; i < (unsigned long) nombre_colonnes_a; pivot[i++] = 0); 846: 847: lwork = -1; 848: 849: if ((work = malloc(sizeof(complex16))) == NULL) 850: { 851: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 852: return; 853: } 854: 855: // Calcul de la taille de l'espace 856: 857: zgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77, 858: &nombre_lignes_a, pivot, (*((complex16 **) tau)), 859: work, &lwork, rwork, &erreur); 860: 861: if (erreur != 0) 862: { 863: // L'erreur ne peut être que négative. 864: 865: (*s_etat_processus).erreur_execution = 866: d_ex_routines_mathematiques; 867: 868: free(work); 869: free(rwork); 870: free(pivot); 871: free(matrice_a_f77); 872: free(tau); 873: return; 874: } 875: 876: lwork = ((complex16 *) work)[0].partie_reelle; 877: 878: free(work); 879: 880: if ((work = malloc(lwork * sizeof(complex16))) == NULL) 881: { 882: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 883: return; 884: } 885: 886: // Calcul de la permutation 887: 888: if ((registre = (void *) malloc(taille_matrice_f77 * 889: sizeof(complex16))) == NULL) 890: { 891: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 892: return; 893: } 894: 895: memcpy(registre, matrice_a_f77, 896: taille_matrice_f77 * sizeof(complex16)); 897: 898: zgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, registre, 899: &nombre_lignes_a, pivot, (*((complex16 **) tau)), 900: work, &lwork, rwork, &erreur); 901: 902: free(registre); 903: 904: if (erreur != 0) 905: { 906: // L'erreur ne peut être que négative. 907: 908: (*s_etat_processus).erreur_execution = 909: d_ex_routines_mathematiques; 910: 911: free(work); 912: free(rwork); 913: free(pivot); 914: free(matrice_a_f77); 915: free(tau); 916: return; 917: } 918: 919: // La permutation doit maintenant être unitaire 920: 921: zgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77, 922: &nombre_lignes_a, pivot, (*((complex16 **) tau)), 923: work, &lwork, rwork, &erreur); 924: 925: if (erreur != 0) 926: { 927: // L'erreur ne peut être que négative. 928: 929: (*s_etat_processus).erreur_execution = 930: d_ex_routines_mathematiques; 931: 932: free(work); 933: free(rwork); 934: free(pivot); 935: free(matrice_a_f77); 936: free(tau); 937: return; 938: } 939: 940: free(rwork); 941: 942: for(i = 0; i < (unsigned long) nombre_colonnes_a; i++) 943: { 944: if ((i + 1) != (unsigned long) pivot[i]) 945: { 946: (*s_etat_processus).erreur_execution = 947: d_ex_routines_mathematiques; 948: 949: free(pivot); 950: free(work); 951: free(matrice_a_f77); 952: free(tau); 953: 954: return; 955: } 956: } 957: 958: free(pivot); 959: 960: for(k = 0, i = 0; i < (unsigned long) nombre_colonnes_a; i++) 961: { 962: for(j = 0; j < (unsigned long) nombre_lignes_a; j++) 963: { 964: ((complex16 **) (*s_matrice).tableau)[j][i].partie_reelle = 965: ((complex16 *) matrice_a_f77)[k].partie_reelle; 966: ((complex16 **) (*s_matrice).tableau)[j][i] 967: .partie_imaginaire = ((complex16 *) 968: matrice_a_f77)[k++].partie_imaginaire; 969: } 970: } 971: 972: free(work); 973: free(matrice_a_f77); 974: 975: break; 976: } 977: } 978: 979: return; 980: } 981: 982: 983: /* 984: ================================================================================ 985: Fonctions calculant le déterminant ou le rang d'une matrice quelconque 986: ================================================================================ 987: Entrées : struct_matrice 988: -------------------------------------------------------------------------------- 989: Sorties : déterminant 990: -------------------------------------------------------------------------------- 991: Effets de bord : néant 992: ================================================================================ 993: */ 994: 995: 996: static integer4 997: calcul_rang(struct_processus *s_etat_processus, void *matrice_f77, 998: integer4 nombre_lignes_a, integer4 nombre_colonnes_a, 999: integer4 *pivot, integer4 dimension_vecteur_pivot, unsigned char type) 1000: { 1001: integer4 erreur; 1002: integer4 *iwork; 1003: integer4 *jpvt; 1004: integer4 lwork; 1005: integer4 longueur; 1006: integer4 nombre_colonnes_b; 1007: integer4 nombre_lignes_b; 1008: integer4 ordre; 1009: integer4 rang; 1010: 1011: real8 anorme; 1012: real8 rcond; 1013: real8 *rwork; 1014: 1015: unsigned char norme; 1016: 1017: unsigned long i; 1018: 1019: void *matrice_b; 1020: void *matrice_c; 1021: void *work; 1022: 1023: #undef NORME_I 1024: #ifdef NORME_I 1025: norme = 'I'; 1026: 1027: if ((work = malloc(nombre_lignes_a * sizeof(real8))) == NULL) 1028: { 1029: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1030: return(-1); 1031: } 1032: #else 1033: norme = '1'; 1034: work = NULL; 1035: #endif 1036: 1037: longueur = 1; 1038: 1039: if (type == 'R') 1040: { 1041: anorme = dlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a, 1042: matrice_f77, &nombre_lignes_a, work, longueur); 1043: 1044: #ifndef NORME_I 1045: free(work); 1046: #endif 1047: 1048: if ((matrice_c = malloc(nombre_lignes_a * nombre_colonnes_a * 1049: sizeof(real8))) == NULL) 1050: { 1051: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1052: return(-1); 1053: } 1054: 1055: memcpy(matrice_c, matrice_f77, nombre_lignes_a * nombre_colonnes_a * 1056: sizeof(real8)); 1057: 1058: dgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77, 1059: &nombre_lignes_a, pivot, &erreur); 1060: 1061: if (erreur < 0) 1062: { 1063: (*s_etat_processus).erreur_execution = 1064: d_ex_routines_mathematiques; 1065: 1066: free(matrice_f77); 1067: return(-1); 1068: } 1069: 1070: if ((iwork = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL) 1071: { 1072: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1073: return(-1); 1074: } 1075: 1076: if ((work = malloc(4 * nombre_colonnes_a * sizeof(real8))) == NULL) 1077: { 1078: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1079: return(-1); 1080: } 1081: 1082: ordre = (nombre_lignes_a > nombre_colonnes_a) 1083: ? nombre_colonnes_a : nombre_lignes_a; 1084: 1085: dgecon_(&norme, &ordre, matrice_f77, 1086: &nombre_lignes_a, &anorme, &rcond, work, iwork, &erreur, 1087: longueur); 1088: 1089: free(work); 1090: free(iwork); 1091: 1092: if (erreur < 0) 1093: { 1094: (*s_etat_processus).erreur_execution = 1095: d_ex_routines_mathematiques; 1096: 1097: free(matrice_f77); 1098: return(-1); 1099: } 1100: 1101: if ((jpvt = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL) 1102: { 1103: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1104: return(-1); 1105: } 1106: 1107: for(i = 0; i < (unsigned long) nombre_colonnes_a; i++) 1108: { 1109: ((integer4 *) jpvt)[i] = 0; 1110: } 1111: 1112: lwork = -1; 1113: 1114: if ((work = malloc(sizeof(real8))) == NULL) 1115: { 1116: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1117: return(-1); 1118: } 1119: 1120: nombre_colonnes_b = 1; 1121: nombre_lignes_b = (nombre_lignes_a > nombre_colonnes_a) 1122: ? nombre_lignes_a : nombre_colonnes_a; 1123: 1124: if ((matrice_b = malloc(nombre_lignes_b * sizeof(real8))) == NULL) 1125: { 1126: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1127: return(-1); 1128: } 1129: 1130: for(i = 0; i < (unsigned long) nombre_lignes_b; i++) 1131: { 1132: ((real8 *) matrice_b)[i] = 0; 1133: } 1134: 1135: dgelsy_(&nombre_lignes_a, &nombre_colonnes_a, 1136: &nombre_colonnes_b, matrice_c, &nombre_lignes_a, 1137: matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang, 1138: work, &lwork, &erreur); 1139: 1140: lwork = ((real8 *) work)[0]; 1141: free(work); 1142: 1143: if ((work = malloc(lwork * sizeof(real8))) == NULL) 1144: { 1145: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1146: return(-1); 1147: } 1148: 1149: dgelsy_(&nombre_lignes_a, &nombre_colonnes_a, 1150: &nombre_colonnes_b, matrice_c, &nombre_lignes_a, 1151: matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang, 1152: work, &lwork, &erreur); 1153: 1154: free(matrice_b); 1155: free(matrice_c); 1156: free(work); 1157: free(jpvt); 1158: 1159: if (erreur < 0) 1160: { 1161: (*s_etat_processus).erreur_execution = 1162: d_ex_routines_mathematiques; 1163: 1164: free(matrice_f77); 1165: return(-1); 1166: } 1167: } 1168: else 1169: { 1170: anorme = zlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a, 1171: matrice_f77, &nombre_lignes_a, work, longueur); 1172: 1173: #ifndef NORME_I 1174: free(work); 1175: #endif 1176: 1177: if ((matrice_c = malloc(nombre_lignes_a * nombre_colonnes_a * 1178: sizeof(complex16))) == NULL) 1179: { 1180: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1181: return(-1); 1182: } 1183: 1184: memcpy(matrice_c, matrice_f77, nombre_lignes_a * nombre_colonnes_a * 1185: sizeof(complex16)); 1186: 1187: zgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77, 1188: &nombre_lignes_a, pivot, &erreur); 1189: 1190: if (erreur < 0) 1191: { 1192: (*s_etat_processus).erreur_execution = 1193: d_ex_routines_mathematiques; 1194: 1195: free(matrice_f77); 1196: return(-1); 1197: } 1198: 1199: if ((rwork = malloc(2 * nombre_colonnes_a * sizeof(real8))) == NULL) 1200: { 1201: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1202: return(-1); 1203: } 1204: 1205: if ((work = malloc(2 * nombre_colonnes_a * sizeof(complex16))) == NULL) 1206: { 1207: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1208: return(-1); 1209: } 1210: 1211: ordre = (nombre_lignes_a > nombre_colonnes_a) 1212: ? nombre_colonnes_a : nombre_lignes_a; 1213: 1214: zgecon_(&norme, &ordre, matrice_f77, 1215: &nombre_lignes_a, &anorme, &rcond, work, rwork, &erreur, 1216: longueur); 1217: 1218: free(work); 1219: 1220: if (erreur < 0) 1221: { 1222: (*s_etat_processus).erreur_execution = 1223: d_ex_routines_mathematiques; 1224: 1225: free(matrice_f77); 1226: return(-1); 1227: } 1228: 1229: if ((jpvt = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL) 1230: { 1231: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1232: return(-1); 1233: } 1234: 1235: for(i = 0; i < (unsigned long) nombre_colonnes_a; i++) 1236: { 1237: ((integer4 *) jpvt)[i] = 0; 1238: } 1239: 1240: lwork = -1; 1241: 1242: if ((work = malloc(sizeof(complex16))) == NULL) 1243: { 1244: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1245: return(-1); 1246: } 1247: 1248: nombre_colonnes_b = 1; 1249: nombre_lignes_b = (nombre_lignes_a > nombre_colonnes_a) 1250: ? nombre_lignes_a : nombre_colonnes_a; 1251: 1252: if ((matrice_b = malloc(nombre_lignes_b * sizeof(complex16))) == NULL) 1253: { 1254: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1255: return(-1); 1256: } 1257: 1258: for(i = 0; i < (unsigned long) nombre_lignes_b; i++) 1259: { 1260: ((complex16 *) matrice_b)[i].partie_reelle = 0; 1261: ((complex16 *) matrice_b)[i].partie_imaginaire = 0; 1262: } 1263: 1264: zgelsy_(&nombre_lignes_a, &nombre_colonnes_a, 1265: &nombre_colonnes_b, matrice_c, &nombre_lignes_a, 1266: matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang, 1267: work, &lwork, rwork, &erreur); 1268: 1269: lwork = ((complex16 *) work)[0].partie_reelle; 1270: free(work); 1271: 1272: if ((work = malloc(lwork * sizeof(complex16))) == NULL) 1273: { 1274: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1275: return(-1); 1276: } 1277: 1278: zgelsy_(&nombre_lignes_a, &nombre_colonnes_a, 1279: &nombre_colonnes_b, matrice_c, &nombre_lignes_a, 1280: matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang, 1281: work, &lwork, rwork, &erreur); 1282: 1283: free(rwork); 1284: free(matrice_b); 1285: free(matrice_c); 1286: free(work); 1287: free(jpvt); 1288: 1289: if (erreur < 0) 1290: { 1291: (*s_etat_processus).erreur_execution = 1292: d_ex_routines_mathematiques; 1293: 1294: free(matrice_f77); 1295: return(-1); 1296: } 1297: } 1298: 1299: return(rang); 1300: } 1301: 1302: 1303: void 1304: determinant(struct_processus *s_etat_processus, struct_matrice *s_matrice, 1305: void *valeur) 1306: { 1307: complex16 *vecteur_complexe; 1308: 1309: integer4 dimension_vecteur_pivot; 1310: integer4 nombre_colonnes_a; 1311: integer4 nombre_lignes_a; 1312: integer4 *pivot; 1313: integer4 rang; 1314: 1315: integer8 signe; 1316: 1317: real8 *vecteur_reel; 1318: 1319: unsigned long i; 1320: unsigned long j; 1321: unsigned long k; 1322: unsigned long taille_matrice_f77; 1323: 1324: void *matrice_f77; 1325: 1326: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes; 1327: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes; 1328: dimension_vecteur_pivot = (nombre_lignes_a < nombre_colonnes_a) 1329: ? nombre_lignes_a : nombre_colonnes_a; 1330: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a; 1331: 1332: switch((*s_matrice).type) 1333: { 1334: case 'I' : 1335: { 1336: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 1337: sizeof(real8))) == NULL) 1338: { 1339: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1340: return; 1341: } 1342: 1343: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 1344: { 1345: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 1346: { 1347: ((real8 *) matrice_f77)[k++] = ((integer8 **) 1348: (*s_matrice).tableau)[j][i]; 1349: } 1350: } 1351: 1352: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot * 1353: sizeof(integer4))) == NULL) 1354: { 1355: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1356: return; 1357: } 1358: 1359: if ((rang = calcul_rang(s_etat_processus, matrice_f77, 1360: nombre_lignes_a, nombre_colonnes_a, pivot, 1361: dimension_vecteur_pivot, 'R')) < 0) 1362: { 1363: return; 1364: } 1365: 1366: if (rang < nombre_lignes_a) 1367: { 1368: (*((real8 *) valeur)) = 0; 1369: } 1370: else 1371: { 1372: if ((vecteur_reel = malloc((*s_matrice).nombre_colonnes * 1373: sizeof(real8))) == NULL) 1374: { 1375: (*s_etat_processus).erreur_systeme = 1376: d_es_allocation_memoire; 1377: return; 1378: } 1379: 1380: signe = 1; 1381: 1382: for(i = 0; i < (*s_matrice).nombre_colonnes; i++) 1383: { 1384: if ((unsigned long) pivot[i] != (i + 1)) 1385: { 1386: signe = (signe == 1) ? -1 : 1; 1387: } 1388: 1389: vecteur_reel[i] = ((real8 *) matrice_f77) 1390: [(i * nombre_colonnes_a) + i]; 1391: } 1392: 1393: for(i = 1; i < (*s_matrice).nombre_colonnes; i++) 1394: { 1395: vecteur_reel[0] *= vecteur_reel[i]; 1396: } 1397: 1398: (*((real8 *) valeur)) = vecteur_reel[0] * signe; 1399: 1400: free(vecteur_reel); 1401: } 1402: 1403: free(matrice_f77); 1404: free(pivot); 1405: 1406: break; 1407: } 1408: 1409: case 'R' : 1410: { 1411: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 1412: sizeof(real8))) == NULL) 1413: { 1414: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1415: return; 1416: } 1417: 1418: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 1419: { 1420: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 1421: { 1422: ((real8 *) matrice_f77)[k++] = ((real8 **) 1423: (*s_matrice).tableau)[j][i]; 1424: } 1425: } 1426: 1427: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot * 1428: sizeof(integer4))) == NULL) 1429: { 1430: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1431: return; 1432: } 1433: 1434: if ((rang = calcul_rang(s_etat_processus, matrice_f77, 1435: nombre_lignes_a, nombre_colonnes_a, pivot, 1436: dimension_vecteur_pivot, 'R')) < 0) 1437: { 1438: return; 1439: } 1440: 1441: if (rang < nombre_lignes_a) 1442: { 1443: (*((real8 *) valeur)) = 0; 1444: } 1445: else 1446: { 1447: if ((vecteur_reel = malloc((*s_matrice).nombre_colonnes * 1448: sizeof(real8))) == NULL) 1449: { 1450: (*s_etat_processus).erreur_systeme = 1451: d_es_allocation_memoire; 1452: return; 1453: } 1454: 1455: signe = 1; 1456: 1457: for(i = 0; i < (*s_matrice).nombre_colonnes; i++) 1458: { 1459: if ((unsigned long) pivot[i] != (i + 1)) 1460: { 1461: signe = (signe == 1) ? -1 : 1; 1462: } 1463: 1464: vecteur_reel[i] = ((real8 *) matrice_f77) 1465: [(i * nombre_colonnes_a) + i]; 1466: } 1467: 1468: for(i = 1; i < (*s_matrice).nombre_colonnes; i++) 1469: { 1470: vecteur_reel[0] *= vecteur_reel[i]; 1471: } 1472: 1473: (*((real8 *) valeur)) = vecteur_reel[0] * signe; 1474: 1475: free(vecteur_reel); 1476: } 1477: 1478: free(matrice_f77); 1479: free(pivot); 1480: 1481: break; 1482: } 1483: 1484: case 'C' : 1485: { 1486: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 1487: sizeof(complex16))) == NULL) 1488: { 1489: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1490: return; 1491: } 1492: 1493: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 1494: { 1495: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 1496: { 1497: ((complex16 *) matrice_f77)[k++] = ((complex16 **) 1498: (*s_matrice).tableau)[j][i]; 1499: } 1500: } 1501: 1502: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot * 1503: sizeof(integer4))) == NULL) 1504: { 1505: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1506: return; 1507: } 1508: 1509: if ((rang = calcul_rang(s_etat_processus, matrice_f77, 1510: nombre_lignes_a, nombre_colonnes_a, pivot, 1511: dimension_vecteur_pivot, 'C')) < 0) 1512: { 1513: return; 1514: } 1515: 1516: if (rang < nombre_lignes_a) 1517: { 1518: (*((complex16 *) valeur)).partie_reelle = 0; 1519: (*((complex16 *) valeur)).partie_imaginaire = 0; 1520: } 1521: else 1522: { 1523: if ((vecteur_complexe = malloc((*s_matrice).nombre_colonnes * 1524: sizeof(complex16))) == NULL) 1525: { 1526: (*s_etat_processus).erreur_systeme = 1527: d_es_allocation_memoire; 1528: return; 1529: } 1530: 1531: signe = 1; 1532: 1533: for(i = 0; i < (*s_matrice).nombre_colonnes; i++) 1534: { 1535: if ((unsigned long) pivot[i] != (i + 1)) 1536: { 1537: signe = (signe == 1) ? -1 : 1; 1538: } 1539: 1540: vecteur_complexe[i] = ((complex16 *) matrice_f77) 1541: [(i * nombre_colonnes_a) + i]; 1542: } 1543: 1544: for(i = 1; i < (*s_matrice).nombre_colonnes; i++) 1545: { 1546: f77multiplicationcc_(&(vecteur_complexe[0]), 1547: &(vecteur_complexe[i]), &(vecteur_complexe[0])); 1548: } 1549: 1550: f77multiplicationci_(&(vecteur_complexe[0]), &signe, 1551: ((complex16 *) valeur)); 1552: 1553: free(vecteur_complexe); 1554: } 1555: 1556: free(matrice_f77); 1557: free(pivot); 1558: 1559: break; 1560: } 1561: } 1562: 1563: return; 1564: } 1565: 1566: 1567: void 1568: rang(struct_processus *s_etat_processus, struct_matrice *s_matrice, 1569: integer8 *valeur) 1570: { 1571: integer4 dimension_vecteur_pivot; 1572: integer4 nombre_lignes_a; 1573: integer4 nombre_colonnes_a; 1574: integer4 *pivot; 1575: integer4 rang; 1576: integer4 taille_matrice_f77; 1577: 1578: unsigned long i; 1579: unsigned long j; 1580: unsigned long k; 1581: 1582: void *matrice_f77; 1583: 1584: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes; 1585: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes; 1586: dimension_vecteur_pivot = (nombre_lignes_a < nombre_colonnes_a) 1587: ? nombre_lignes_a : nombre_colonnes_a; 1588: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a; 1589: 1590: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot * 1591: sizeof(integer4))) == NULL) 1592: { 1593: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1594: return; 1595: } 1596: 1597: switch((*s_matrice).type) 1598: { 1599: case 'I' : 1600: { 1601: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 1602: sizeof(real8))) == NULL) 1603: { 1604: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1605: return; 1606: } 1607: 1608: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 1609: { 1610: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 1611: { 1612: ((real8 *) matrice_f77)[k++] = ((integer8 **) 1613: (*s_matrice).tableau)[j][i]; 1614: } 1615: } 1616: 1617: if ((rang = calcul_rang(s_etat_processus, matrice_f77, 1618: nombre_lignes_a, nombre_colonnes_a, pivot, 1619: dimension_vecteur_pivot, 'R')) < 0) 1620: { 1621: free(pivot); 1622: free(matrice_f77); 1623: return; 1624: } 1625: 1626: free(matrice_f77); 1627: (*valeur) = rang; 1628: break; 1629: } 1630: 1631: case 'R' : 1632: { 1633: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 1634: sizeof(real8))) == NULL) 1635: { 1636: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1637: return; 1638: } 1639: 1640: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 1641: { 1642: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 1643: { 1644: ((real8 *) matrice_f77)[k++] = ((real8 **) 1645: (*s_matrice).tableau)[j][i]; 1646: } 1647: } 1648: 1649: if ((rang = calcul_rang(s_etat_processus, matrice_f77, 1650: nombre_lignes_a, nombre_colonnes_a, pivot, 1651: dimension_vecteur_pivot, 'R')) < 0) 1652: { 1653: free(pivot); 1654: free(matrice_f77); 1655: return; 1656: } 1657: 1658: free(matrice_f77); 1659: (*valeur) = rang; 1660: break; 1661: } 1662: 1663: case 'C' : 1664: { 1665: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 1666: sizeof(complex16))) == NULL) 1667: { 1668: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 1669: return; 1670: } 1671: 1672: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 1673: { 1674: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 1675: { 1676: ((complex16 *) matrice_f77)[k++] = ((complex16 **) 1677: (*s_matrice).tableau)[j][i]; 1678: } 1679: } 1680: 1681: if ((rang = calcul_rang(s_etat_processus, matrice_f77, 1682: nombre_lignes_a, nombre_colonnes_a, pivot, 1683: dimension_vecteur_pivot, 'C')) < 0) 1684: { 1685: free(pivot); 1686: free(matrice_f77); 1687: return; 1688: } 1689: 1690: free(matrice_f77); 1691: (*valeur) = rang; 1692: break; 1693: } 1694: } 1695: 1696: free(pivot); 1697: 1698: return; 1699: } 1700: 1701: // vim: ts=4