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