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