![]() ![]() | ![]() |
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 calculant le nombre de condition d'une matrice 29: ================================================================================ 30: Entrées : struct_matrice 31: -------------------------------------------------------------------------------- 32: Sorties : nombre de condition de la matrice 33: -------------------------------------------------------------------------------- 34: Effets de bord : néant 35: ================================================================================ 36: */ 37: 38: static integer4 39: calcul_cond(struct_processus *s_etat_processus, void *matrice_f77, 40: integer4 nombre_lignes_a, integer4 nombre_colonnes_a, 41: integer4 *pivot, unsigned char type, real8 *cond) 42: { 43: integer4 erreur; 44: integer4 *iwork; 45: integer4 longueur; 46: integer4 ordre; 47: 48: real8 anorme; 49: real8 rcond; 50: real8 *rwork; 51: 52: unsigned char norme; 53: 54: void *work; 55: 56: norme = '1'; 57: longueur = 1; 58: 59: if (type == 'R') 60: { 61: // work est NULL dans le cas de la norme '1' 62: anorme = dlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a, 63: matrice_f77, &nombre_lignes_a, NULL, longueur); 64: 65: dgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77, 66: &nombre_lignes_a, pivot, &erreur); 67: 68: if (erreur < 0) 69: { 70: (*s_etat_processus).erreur_execution = 71: d_ex_routines_mathematiques; 72: 73: free(matrice_f77); 74: return(-1); 75: } 76: 77: if ((iwork = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL) 78: { 79: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 80: return(-1); 81: } 82: 83: if ((work = malloc(4 * nombre_colonnes_a * sizeof(real8))) == NULL) 84: { 85: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 86: return(-1); 87: } 88: 89: ordre = (nombre_lignes_a > nombre_colonnes_a) 90: ? nombre_colonnes_a : nombre_lignes_a; 91: 92: dgecon_(&norme, &ordre, matrice_f77, 93: &nombre_lignes_a, &anorme, &rcond, work, iwork, &erreur, 94: longueur); 95: 96: free(work); 97: free(iwork); 98: 99: if (erreur < 0) 100: { 101: (*s_etat_processus).erreur_execution = 102: d_ex_routines_mathematiques; 103: 104: free(matrice_f77); 105: return(-1); 106: } 107: } 108: else 109: { 110: // work est NULL dans le cas de la norme '1' 111: anorme = zlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a, 112: matrice_f77, &nombre_lignes_a, NULL, longueur); 113: 114: zgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77, 115: &nombre_lignes_a, pivot, &erreur); 116: 117: if (erreur < 0) 118: { 119: (*s_etat_processus).erreur_execution = 120: d_ex_routines_mathematiques; 121: 122: free(matrice_f77); 123: return(-1); 124: } 125: 126: if ((rwork = malloc(2 * nombre_colonnes_a * sizeof(real8))) == NULL) 127: { 128: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 129: return(-1); 130: } 131: 132: if ((work = malloc(2 * nombre_colonnes_a * sizeof(complex16))) == NULL) 133: { 134: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 135: return(-1); 136: } 137: 138: ordre = (nombre_lignes_a > nombre_colonnes_a) 139: ? nombre_colonnes_a : nombre_lignes_a; 140: 141: zgecon_(&norme, &ordre, matrice_f77, 142: &nombre_lignes_a, &anorme, &rcond, work, rwork, &erreur, 143: longueur); 144: 145: free(work); 146: free(rwork); 147: 148: if (erreur < 0) 149: { 150: (*s_etat_processus).erreur_execution = 151: d_ex_routines_mathematiques; 152: 153: free(matrice_f77); 154: return(-1); 155: } 156: } 157: 158: (*cond) = ((real8) 1 / rcond); 159: return(0); 160: } 161: 162: 163: void 164: cond(struct_processus *s_etat_processus, 165: struct_matrice *s_matrice, real8 *condition) 166: { 167: integer4 dimension_vecteur_pivot; 168: integer4 nombre_lignes_a; 169: integer4 nombre_colonnes_a; 170: integer4 *pivot; 171: integer4 rang; 172: integer4 taille_matrice_f77; 173: 174: real8 cond; 175: 176: unsigned long i; 177: unsigned long j; 178: unsigned long k; 179: 180: void *matrice_f77; 181: 182: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes; 183: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes; 184: dimension_vecteur_pivot = (nombre_lignes_a < nombre_colonnes_a) 185: ? nombre_lignes_a : nombre_colonnes_a; 186: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a; 187: 188: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot * 189: sizeof(integer4))) == NULL) 190: { 191: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 192: return; 193: } 194: 195: switch((*s_matrice).type) 196: { 197: case 'I' : 198: { 199: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 200: sizeof(real8))) == NULL) 201: { 202: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 203: return; 204: } 205: 206: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 207: { 208: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 209: { 210: ((real8 *) matrice_f77)[k++] = ((integer8 **) 211: (*s_matrice).tableau)[j][i]; 212: } 213: } 214: 215: if ((rang = calcul_cond(s_etat_processus, matrice_f77, 216: nombre_lignes_a, nombre_colonnes_a, pivot, 217: 'R', &cond)) < 0) 218: { 219: free(pivot); 220: free(matrice_f77); 221: return; 222: } 223: 224: free(matrice_f77); 225: (*condition) = cond; 226: break; 227: } 228: 229: case 'R' : 230: { 231: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 232: sizeof(real8))) == NULL) 233: { 234: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 235: return; 236: } 237: 238: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 239: { 240: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 241: { 242: ((real8 *) matrice_f77)[k++] = ((real8 **) 243: (*s_matrice).tableau)[j][i]; 244: } 245: } 246: 247: if ((rang = calcul_cond(s_etat_processus, matrice_f77, 248: nombre_lignes_a, nombre_colonnes_a, pivot, 249: 'R', &cond)) < 0) 250: { 251: free(pivot); 252: free(matrice_f77); 253: return; 254: } 255: 256: free(matrice_f77); 257: (*condition) = cond; 258: break; 259: } 260: 261: case 'C' : 262: { 263: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 264: sizeof(complex16))) == NULL) 265: { 266: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 267: return; 268: } 269: 270: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 271: { 272: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 273: { 274: ((complex16 *) matrice_f77)[k++] = ((complex16 **) 275: (*s_matrice).tableau)[j][i]; 276: } 277: } 278: 279: if ((rang = calcul_cond(s_etat_processus, matrice_f77, 280: nombre_lignes_a, nombre_colonnes_a, pivot, 281: 'C', &cond)) < 0) 282: { 283: free(pivot); 284: free(matrice_f77); 285: return; 286: } 287: 288: free(matrice_f77); 289: (*condition) = cond; 290: break; 291: } 292: } 293: 294: free(pivot); 295: 296: return; 297: } 298: 299: 300: /* 301: ================================================================================ 302: Fonction effectuant une décomposition en valeurs singulières 303: ================================================================================ 304: Entrées : struct_matrice 305: -------------------------------------------------------------------------------- 306: Sorties : valeurs singulières dans le vecteur S. Si les pointeurs sur U 307: et VH ne sont pas nul, les matrices U et VH sont aussi calculées. 308: -------------------------------------------------------------------------------- 309: Effets de bord : néant 310: ================================================================================ 311: */ 312: 313: void valeurs_singulieres(struct_processus *s_etat_processus, 314: struct_matrice *s_matrice, struct_matrice *matrice_u, 315: struct_vecteur *vecteur_s, struct_matrice *matrice_vh) 316: { 317: integer4 erreur; 318: integer4 longueur; 319: integer4 lwork; 320: integer4 nombre_colonnes_a; 321: integer4 nombre_lignes_a; 322: integer4 nombre_valeurs_singulieres; 323: integer4 taille_matrice_f77; 324: 325: real8 *rwork; 326: 327: unsigned char jobu; 328: unsigned char jobvh; 329: 330: unsigned long i; 331: unsigned long j; 332: unsigned long k; 333: 334: void *matrice_f77; 335: void *matrice_f77_u; 336: void *matrice_f77_vh; 337: void *vecteur_f77_s; 338: void *work; 339: 340: longueur = 1; 341: 342: if (matrice_u != NULL) 343: { 344: jobu = 'A'; 345: } 346: else 347: { 348: jobu = 'N'; 349: } 350: 351: if (matrice_vh != NULL) 352: { 353: jobvh = 'A'; 354: } 355: else 356: { 357: jobvh = 'N'; 358: } 359: 360: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes; 361: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes; 362: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a; 363: nombre_valeurs_singulieres = (nombre_lignes_a < nombre_colonnes_a) 364: ? nombre_lignes_a : nombre_colonnes_a; 365: 366: switch((*s_matrice).type) 367: { 368: case 'I' : 369: { 370: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 371: sizeof(real8))) == NULL) 372: { 373: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 374: return; 375: } 376: 377: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 378: { 379: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 380: { 381: ((real8 *) matrice_f77)[k++] = ((integer8 **) 382: (*s_matrice).tableau)[j][i]; 383: } 384: } 385: 386: lwork = -1; 387: 388: if ((work = malloc(sizeof(real8))) == NULL) 389: { 390: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 391: return; 392: } 393: 394: if (matrice_u != NULL) 395: { 396: if ((matrice_f77_u = malloc(nombre_lignes_a * nombre_lignes_a * 397: sizeof(real8))) == NULL) 398: { 399: (*s_etat_processus).erreur_systeme = 400: d_es_allocation_memoire; 401: return; 402: } 403: } 404: else 405: { 406: matrice_f77_u = NULL; 407: } 408: 409: if ((vecteur_f77_s = malloc(nombre_valeurs_singulieres * 410: sizeof(real8))) == NULL) 411: { 412: (*s_etat_processus).erreur_systeme = 413: d_es_allocation_memoire; 414: return; 415: } 416: 417: if (matrice_vh != NULL) 418: { 419: if ((matrice_f77_vh = malloc(nombre_colonnes_a 420: * nombre_colonnes_a * sizeof(real8))) == NULL) 421: { 422: (*s_etat_processus).erreur_systeme = 423: d_es_allocation_memoire; 424: return; 425: } 426: } 427: else 428: { 429: matrice_f77_vh = NULL; 430: } 431: 432: dgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a, 433: matrice_f77, &nombre_lignes_a, vecteur_f77_s, 434: matrice_f77_u, &nombre_lignes_a, 435: matrice_f77_vh, &nombre_colonnes_a, 436: work, &lwork, &erreur, longueur, longueur); 437: 438: lwork = ((real8 *) work)[0]; 439: free(work); 440: 441: if ((work = malloc(lwork * sizeof(real8))) == NULL) 442: { 443: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 444: return; 445: } 446: 447: dgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a, 448: matrice_f77, &nombre_lignes_a, vecteur_f77_s, 449: matrice_f77_u, &nombre_lignes_a, 450: matrice_f77_vh, &nombre_colonnes_a, 451: work, &lwork, &erreur, longueur, longueur); 452: 453: free(work); 454: free(matrice_f77); 455: 456: if (erreur != 0) 457: { 458: if (erreur > 0) 459: { 460: (*s_etat_processus).exception = d_ep_decomposition_SVD; 461: } 462: else 463: { 464: (*s_etat_processus).erreur_execution = 465: d_ex_routines_mathematiques; 466: } 467: 468: free(matrice_f77_u); 469: free(matrice_f77_vh); 470: free(vecteur_f77_s); 471: return; 472: } 473: 474: if (matrice_u != NULL) 475: { 476: (*matrice_u).nombre_lignes = nombre_lignes_a; 477: (*matrice_u).nombre_colonnes = nombre_lignes_a; 478: 479: if (((*matrice_u).tableau = malloc((*matrice_u).nombre_lignes * 480: sizeof(real8 *))) == NULL) 481: { 482: (*s_etat_processus).erreur_systeme = 483: d_es_allocation_memoire; 484: return; 485: } 486: 487: for(i = 0; i < (*matrice_u).nombre_lignes; i++) 488: { 489: if ((((real8 **) (*matrice_u).tableau)[i] = 490: malloc((*matrice_u).nombre_colonnes * 491: sizeof(real8))) == NULL) 492: { 493: (*s_etat_processus).erreur_systeme = 494: d_es_allocation_memoire; 495: return; 496: } 497: } 498: 499: for(k = 0, i = 0; i < (*matrice_u).nombre_colonnes; i++) 500: { 501: for(j = 0; j < (*matrice_u).nombre_lignes; j++) 502: { 503: ((real8 **) (*matrice_u).tableau)[j][i] = 504: ((real8 *) matrice_f77_u)[k++]; 505: } 506: } 507: 508: free(matrice_f77_u); 509: } 510: 511: if (matrice_vh != NULL) 512: { 513: (*matrice_vh).nombre_lignes = nombre_colonnes_a; 514: (*matrice_vh).nombre_colonnes = nombre_colonnes_a; 515: 516: if (((*matrice_vh).tableau = malloc((*matrice_vh) 517: .nombre_lignes * sizeof(real8 *))) == NULL) 518: { 519: (*s_etat_processus).erreur_systeme = 520: d_es_allocation_memoire; 521: return; 522: } 523: 524: for(i = 0; i < (*matrice_vh).nombre_lignes; i++) 525: { 526: if ((((real8 **) (*matrice_vh).tableau)[i] = 527: malloc((*matrice_vh).nombre_colonnes * 528: sizeof(real8))) == NULL) 529: { 530: (*s_etat_processus).erreur_systeme = 531: d_es_allocation_memoire; 532: return; 533: } 534: } 535: 536: for(k = 0, i = 0; i < (*matrice_vh).nombre_colonnes; i++) 537: { 538: for(j = 0; j < (*matrice_vh).nombre_lignes; j++) 539: { 540: ((real8 **) (*matrice_vh).tableau)[j][i] = 541: ((real8 *) matrice_f77_vh)[k++]; 542: } 543: } 544: 545: free(matrice_f77_vh); 546: } 547: 548: (*vecteur_s).taille = nombre_valeurs_singulieres; 549: (*vecteur_s).type = 'R'; 550: (*vecteur_s).tableau = vecteur_f77_s; 551: 552: break; 553: } 554: 555: case 'R' : 556: { 557: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 558: sizeof(real8))) == NULL) 559: { 560: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 561: return; 562: } 563: 564: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 565: { 566: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 567: { 568: ((real8 *) matrice_f77)[k++] = ((real8 **) 569: (*s_matrice).tableau)[j][i]; 570: } 571: } 572: 573: lwork = -1; 574: 575: if ((work = malloc(sizeof(real8))) == NULL) 576: { 577: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 578: return; 579: } 580: 581: if (matrice_u != NULL) 582: { 583: if ((matrice_f77_u = malloc(nombre_lignes_a * nombre_lignes_a * 584: sizeof(real8))) == NULL) 585: { 586: (*s_etat_processus).erreur_systeme = 587: d_es_allocation_memoire; 588: return; 589: } 590: } 591: else 592: { 593: matrice_f77_u = NULL; 594: } 595: 596: if ((vecteur_f77_s = malloc(nombre_valeurs_singulieres * 597: sizeof(real8))) == NULL) 598: { 599: (*s_etat_processus).erreur_systeme = 600: d_es_allocation_memoire; 601: return; 602: } 603: 604: if (matrice_vh != NULL) 605: { 606: if ((matrice_f77_vh = malloc(nombre_colonnes_a 607: * nombre_colonnes_a * sizeof(real8))) == NULL) 608: { 609: (*s_etat_processus).erreur_systeme = 610: d_es_allocation_memoire; 611: return; 612: } 613: } 614: else 615: { 616: matrice_f77_vh = NULL; 617: } 618: 619: dgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a, 620: matrice_f77, &nombre_lignes_a, vecteur_f77_s, 621: matrice_f77_u, &nombre_lignes_a, 622: matrice_f77_vh, &nombre_colonnes_a, 623: work, &lwork, &erreur, longueur, longueur); 624: 625: lwork = ((real8 *) work)[0]; 626: free(work); 627: 628: if ((work = malloc(lwork * sizeof(real8))) == NULL) 629: { 630: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 631: return; 632: } 633: 634: dgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a, 635: matrice_f77, &nombre_lignes_a, vecteur_f77_s, 636: matrice_f77_u, &nombre_lignes_a, 637: matrice_f77_vh, &nombre_colonnes_a, 638: work, &lwork, &erreur, longueur, longueur); 639: 640: free(work); 641: free(matrice_f77); 642: 643: if (erreur != 0) 644: { 645: if (erreur > 0) 646: { 647: (*s_etat_processus).exception = d_ep_decomposition_SVD; 648: } 649: else 650: { 651: (*s_etat_processus).erreur_execution = 652: d_ex_routines_mathematiques; 653: } 654: 655: free(matrice_f77_u); 656: free(matrice_f77_vh); 657: free(vecteur_f77_s); 658: return; 659: } 660: 661: if (matrice_u != NULL) 662: { 663: (*matrice_u).nombre_lignes = nombre_lignes_a; 664: (*matrice_u).nombre_colonnes = nombre_lignes_a; 665: 666: if (((*matrice_u).tableau = malloc((*matrice_u).nombre_lignes * 667: sizeof(real8 *))) == NULL) 668: { 669: (*s_etat_processus).erreur_systeme = 670: d_es_allocation_memoire; 671: return; 672: } 673: 674: for(i = 0; i < (*matrice_u).nombre_lignes; i++) 675: { 676: if ((((real8 **) (*matrice_u).tableau)[i] = 677: malloc((*matrice_u).nombre_colonnes * 678: sizeof(real8))) == NULL) 679: { 680: (*s_etat_processus).erreur_systeme = 681: d_es_allocation_memoire; 682: return; 683: } 684: } 685: 686: for(k = 0, i = 0; i < (*matrice_u).nombre_colonnes; i++) 687: { 688: for(j = 0; j < (*matrice_u).nombre_lignes; j++) 689: { 690: ((real8 **) (*matrice_u).tableau)[j][i] = 691: ((real8 *) matrice_f77_u)[k++]; 692: } 693: } 694: 695: free(matrice_f77_u); 696: } 697: 698: if (matrice_vh != NULL) 699: { 700: (*matrice_vh).nombre_lignes = nombre_colonnes_a; 701: (*matrice_vh).nombre_colonnes = nombre_colonnes_a; 702: 703: if (((*matrice_vh).tableau = malloc((*matrice_vh) 704: .nombre_lignes * sizeof(real8 *))) == NULL) 705: { 706: (*s_etat_processus).erreur_systeme = 707: d_es_allocation_memoire; 708: return; 709: } 710: 711: for(i = 0; i < (*matrice_vh).nombre_lignes; i++) 712: { 713: if ((((real8 **) (*matrice_vh).tableau)[i] = 714: malloc((*matrice_vh).nombre_colonnes * 715: sizeof(real8))) == NULL) 716: { 717: (*s_etat_processus).erreur_systeme = 718: d_es_allocation_memoire; 719: return; 720: } 721: } 722: 723: for(k = 0, i = 0; i < (*matrice_vh).nombre_colonnes; i++) 724: { 725: for(j = 0; j < (*matrice_vh).nombre_lignes; j++) 726: { 727: ((real8 **) (*matrice_vh).tableau)[j][i] = 728: ((real8 *) matrice_f77_vh)[k++]; 729: } 730: } 731: 732: free(matrice_f77_vh); 733: } 734: 735: (*vecteur_s).taille = nombre_valeurs_singulieres; 736: (*vecteur_s).type = 'R'; 737: (*vecteur_s).tableau = vecteur_f77_s; 738: 739: break; 740: } 741: 742: case 'C' : 743: { 744: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 * 745: sizeof(complex16))) == NULL) 746: { 747: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 748: return; 749: } 750: 751: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++) 752: { 753: for(j = 0; j < (*s_matrice).nombre_lignes; j++) 754: { 755: ((complex16 *) matrice_f77)[k++] = ((complex16 **) 756: (*s_matrice).tableau)[j][i]; 757: } 758: } 759: 760: lwork = -1; 761: 762: if ((work = malloc(sizeof(complex16))) == NULL) 763: { 764: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 765: return; 766: } 767: 768: if (matrice_u != NULL) 769: { 770: if ((matrice_f77_u = malloc(nombre_lignes_a * nombre_lignes_a * 771: sizeof(complex16))) == NULL) 772: { 773: (*s_etat_processus).erreur_systeme = 774: d_es_allocation_memoire; 775: return; 776: } 777: } 778: else 779: { 780: matrice_f77_u = NULL; 781: } 782: 783: if ((vecteur_f77_s = malloc(nombre_valeurs_singulieres * 784: sizeof(real8))) == NULL) 785: { 786: (*s_etat_processus).erreur_systeme = 787: d_es_allocation_memoire; 788: return; 789: } 790: 791: if (matrice_vh != NULL) 792: { 793: if ((matrice_f77_vh = malloc(nombre_colonnes_a 794: * nombre_colonnes_a * sizeof(complex16))) == NULL) 795: { 796: (*s_etat_processus).erreur_systeme = 797: d_es_allocation_memoire; 798: return; 799: } 800: } 801: else 802: { 803: matrice_f77_vh = NULL; 804: } 805: 806: if ((rwork = malloc(5 * nombre_valeurs_singulieres * sizeof(real8))) 807: == NULL) 808: { 809: (*s_etat_processus).erreur_systeme = 810: d_es_allocation_memoire; 811: return; 812: } 813: 814: zgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a, 815: matrice_f77, &nombre_lignes_a, vecteur_f77_s, 816: matrice_f77_u, &nombre_lignes_a, 817: matrice_f77_vh, &nombre_colonnes_a, 818: work, &lwork, rwork, &erreur, longueur, longueur); 819: 820: lwork = ((real8 *) work)[0]; 821: free(work); 822: 823: if ((work = malloc(lwork * sizeof(real8))) == NULL) 824: { 825: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire; 826: return; 827: } 828: 829: zgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a, 830: matrice_f77, &nombre_lignes_a, vecteur_f77_s, 831: matrice_f77_u, &nombre_lignes_a, 832: matrice_f77_vh, &nombre_colonnes_a, 833: work, &lwork, rwork, &erreur, longueur, longueur); 834: 835: free(work); 836: free(rwork); 837: free(matrice_f77); 838: 839: if (erreur != 0) 840: { 841: if (erreur > 0) 842: { 843: (*s_etat_processus).exception = d_ep_decomposition_SVD; 844: } 845: else 846: { 847: (*s_etat_processus).erreur_execution = 848: d_ex_routines_mathematiques; 849: } 850: 851: free(matrice_f77_u); 852: free(matrice_f77_vh); 853: free(vecteur_f77_s); 854: return; 855: } 856: 857: if (matrice_u != NULL) 858: { 859: (*matrice_u).nombre_lignes = nombre_lignes_a; 860: (*matrice_u).nombre_colonnes = nombre_lignes_a; 861: 862: if (((*matrice_u).tableau = malloc((*matrice_u).nombre_lignes * 863: sizeof(complex16 *))) == NULL) 864: { 865: (*s_etat_processus).erreur_systeme = 866: d_es_allocation_memoire; 867: return; 868: } 869: 870: for(i = 0; i < (*matrice_u).nombre_lignes; i++) 871: { 872: if ((((complex16 **) (*matrice_u).tableau)[i] = 873: malloc((*matrice_u).nombre_colonnes * 874: sizeof(complex16))) == NULL) 875: { 876: (*s_etat_processus).erreur_systeme = 877: d_es_allocation_memoire; 878: return; 879: } 880: } 881: 882: for(k = 0, i = 0; i < (*matrice_u).nombre_colonnes; i++) 883: { 884: for(j = 0; j < (*matrice_u).nombre_lignes; j++) 885: { 886: ((complex16 **) (*matrice_u).tableau)[j][i] = 887: ((complex16 *) matrice_f77_u)[k++]; 888: } 889: } 890: 891: free(matrice_f77_u); 892: } 893: 894: if (matrice_vh != NULL) 895: { 896: (*matrice_vh).nombre_lignes = nombre_colonnes_a; 897: (*matrice_vh).nombre_colonnes = nombre_colonnes_a; 898: 899: if (((*matrice_vh).tableau = malloc((*matrice_vh) 900: .nombre_lignes * sizeof(complex16 *))) == NULL) 901: { 902: (*s_etat_processus).erreur_systeme = 903: d_es_allocation_memoire; 904: return; 905: } 906: 907: for(i = 0; i < (*matrice_vh).nombre_lignes; i++) 908: { 909: if ((((complex16 **) (*matrice_vh).tableau)[i] = 910: malloc((*matrice_vh).nombre_colonnes * 911: sizeof(complex16))) == NULL) 912: { 913: (*s_etat_processus).erreur_systeme = 914: d_es_allocation_memoire; 915: return; 916: } 917: } 918: 919: for(k = 0, i = 0; i < (*matrice_vh).nombre_colonnes; i++) 920: { 921: for(j = 0; j < (*matrice_vh).nombre_lignes; j++) 922: { 923: ((complex16 **) (*matrice_vh).tableau)[j][i] = 924: ((complex16 *) matrice_f77_vh)[k++]; 925: } 926: } 927: 928: free(matrice_f77_vh); 929: } 930: 931: (*vecteur_s).taille = nombre_valeurs_singulieres; 932: (*vecteur_s).type = 'R'; 933: (*vecteur_s).tableau = vecteur_f77_s; 934: 935: break; 936: } 937: } 938: 939: return; 940: } 941: 942: // vim: ts=4