Annotation of rpl/src/algebre_lineaire4.c, revision 1.12
1.12 ! bertrand 1: /*
! 2: ================================================================================
! 3: RPL/2 (R) version 4.0.18
! 4: Copyright (C) 1989-2010 Dr. BERTRAND Joël
! 5:
! 6: This file is part of RPL/2.
! 7:
! 8: RPL/2 is free software; you can redistribute it and/or modify it
! 9: under the terms of the CeCILL V2 License as published by the french
! 10: CEA, CNRS and INRIA.
! 11:
! 12: RPL/2 is distributed in the hope that it will be useful, but WITHOUT
! 13: ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! 14: FITNESS FOR A PARTICULAR PURPOSE. See the CeCILL V2 License
! 15: for more details.
! 16:
! 17: You should have received a copy of the CeCILL License
! 18: along with RPL/2. If not, write to info@cecill.info.
! 19: ================================================================================
! 20: */
! 21:
! 22:
! 23: #include "rpl-conv.h"
! 24:
! 25:
! 26: /*
! 27: ================================================================================
! 28: Fonction 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
CVSweb interface <joel.bertrand@systella.fr>
![[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