Annotation of rpl/src/algebre_lineaire2.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 réalisation la factorisation LU d'une matrice quelconque
! 29: ================================================================================
! 30: Entrées : struct_matrice
! 31: --------------------------------------------------------------------------------
! 32: Sorties : décomposition A=PLU de la matrice d'entrée et drapeau d'erreur.
! 33: La matrice en entrée est écrasée. La matrice de sortie est réelle
! 34: si la matrice d'entrée est entière ou réelle, et complexe
! 35: si la matrice d'entrée est complexe.
! 36: La routine renvoie aussi une matrice d'entiers correspondant
! 37: à la matrice de permutation. Cette matrice est allouée par
! 38: la routine et vaut NULL sinon.
! 39: --------------------------------------------------------------------------------
! 40: Effets de bord : néant
! 41: ================================================================================
! 42: */
! 43:
! 44: void
! 45: factorisation_lu(struct_processus *s_etat_processus,
! 46: struct_matrice *s_matrice, struct_matrice **s_permutation)
! 47: {
! 48: integer4 dimension_vecteur_pivot;
! 49: integer4 erreur;
! 50: integer4 nombre_colonnes_a;
! 51: integer4 nombre_lignes_a;
! 52: integer4 *pivot;
! 53:
! 54: long n;
! 55:
! 56: unsigned long i;
! 57: unsigned long j;
! 58: unsigned long k;
! 59: unsigned long taille_matrice_f77;
! 60:
! 61: void *matrice_f77;
! 62: void *tampon;
! 63:
! 64: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
! 65: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
! 66: dimension_vecteur_pivot = (nombre_lignes_a < nombre_colonnes_a)
! 67: ? nombre_lignes_a : nombre_colonnes_a;
! 68: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;
! 69:
! 70: switch((*s_matrice).type)
! 71: {
! 72: case 'I' :
! 73: {
! 74: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 75: sizeof(real8))) == NULL)
! 76: {
! 77: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 78: return;
! 79: }
! 80:
! 81: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 82: {
! 83: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 84: {
! 85: ((real8 *) matrice_f77)[k++] = ((integer8 **)
! 86: (*s_matrice).tableau)[j][i];
! 87: }
! 88: }
! 89:
! 90: for(i = 0; i < (*s_matrice).nombre_lignes; i++)
! 91: {
! 92: free(((integer8 **) (*s_matrice).tableau)[i]);
! 93: }
! 94:
! 95: free((integer8 **) (*s_matrice).tableau);
! 96:
! 97: if (((*s_matrice).tableau = (void **) malloc((*s_matrice)
! 98: .nombre_lignes * sizeof(real8 *))) == NULL)
! 99: {
! 100: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 101: return;
! 102: }
! 103:
! 104: for(i = 0; i < (*s_matrice).nombre_lignes; i++)
! 105: {
! 106: if ((((*s_matrice).tableau)[i] =
! 107: (real8 *) malloc((*s_matrice)
! 108: .nombre_colonnes * sizeof(real8))) == NULL)
! 109: {
! 110: (*s_etat_processus).erreur_systeme =
! 111: d_es_allocation_memoire;
! 112: return;
! 113: }
! 114: }
! 115:
! 116: (*s_matrice).type = 'R';
! 117:
! 118: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot *
! 119: sizeof(integer4))) == NULL)
! 120: {
! 121: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 122: return;
! 123: }
! 124:
! 125: dgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77,
! 126: &nombre_lignes_a, pivot, &erreur);
! 127:
! 128: if (erreur < 0)
! 129: {
! 130: (*s_etat_processus).erreur_execution =
! 131: d_ex_routines_mathematiques;
! 132:
! 133: free(pivot);
! 134: free(matrice_f77);
! 135: return;
! 136: }
! 137:
! 138: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 139: {
! 140: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 141: {
! 142: ((real8 **) (*s_matrice).tableau)[j][i] =
! 143: ((real8 *) matrice_f77)[k++];
! 144: }
! 145: }
! 146:
! 147: free(matrice_f77);
! 148:
! 149: (**s_permutation).type = 'I';
! 150: (**s_permutation).nombre_lignes = dimension_vecteur_pivot;
! 151: (**s_permutation).nombre_colonnes = dimension_vecteur_pivot;
! 152:
! 153: if (((**s_permutation).tableau = malloc((**s_permutation)
! 154: .nombre_lignes * sizeof(integer8 *))) == NULL)
! 155: {
! 156: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 157: return;
! 158: }
! 159:
! 160: for(i = 0; i < (**s_permutation).nombre_lignes; i++)
! 161: {
! 162: if ((((integer8 **) (**s_permutation).tableau)[i] =
! 163: (integer8 *) malloc((**s_permutation).nombre_colonnes *
! 164: sizeof(integer8))) == NULL)
! 165: {
! 166: (*s_etat_processus).erreur_systeme =
! 167: d_es_allocation_memoire;
! 168: return;
! 169: }
! 170:
! 171: for(j = 0; j < (**s_permutation).nombre_colonnes; j++)
! 172: {
! 173: ((integer8 **) (**s_permutation).tableau)[i][j] =
! 174: (j == i) ? 1 : 0;
! 175: }
! 176: }
! 177:
! 178: for(n = dimension_vecteur_pivot - 1; n >= 0; n--)
! 179: {
! 180: if ((pivot[n] - 1) != n)
! 181: {
! 182: tampon = ((integer8 **) (**s_permutation).tableau)[n];
! 183: ((integer8 **) (**s_permutation).tableau)[n] =
! 184: ((integer8 **) (**s_permutation).tableau)
! 185: [pivot[n] - 1];
! 186: ((integer8 **) (**s_permutation).tableau)[pivot[n] - 1] =
! 187: tampon;
! 188: }
! 189: }
! 190:
! 191: free(pivot);
! 192:
! 193: break;
! 194: }
! 195:
! 196: case 'R' :
! 197: {
! 198: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 199: sizeof(real8))) == NULL)
! 200: {
! 201: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 202: return;
! 203: }
! 204:
! 205: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 206: {
! 207: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 208: {
! 209: ((real8 *) matrice_f77)[k++] = ((real8 **)
! 210: (*s_matrice).tableau)[j][i];
! 211: }
! 212: }
! 213:
! 214: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot *
! 215: sizeof(integer4))) == NULL)
! 216: {
! 217: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 218: return;
! 219: }
! 220:
! 221: dgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77,
! 222: &nombre_lignes_a, pivot, &erreur);
! 223:
! 224: if (erreur < 0)
! 225: {
! 226: (*s_etat_processus).erreur_execution =
! 227: d_ex_routines_mathematiques;
! 228:
! 229: free(pivot);
! 230: free(matrice_f77);
! 231: return;
! 232: }
! 233:
! 234: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 235: {
! 236: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 237: {
! 238: ((real8 **) (*s_matrice).tableau)[j][i] =
! 239: ((real8 *) matrice_f77)[k++];
! 240: }
! 241: }
! 242:
! 243: free(matrice_f77);
! 244:
! 245: (**s_permutation).type = 'I';
! 246: (**s_permutation).nombre_lignes = dimension_vecteur_pivot;
! 247: (**s_permutation).nombre_colonnes = dimension_vecteur_pivot;
! 248:
! 249: if (((**s_permutation).tableau = malloc((**s_permutation)
! 250: .nombre_lignes * sizeof(integer8 *))) == NULL)
! 251: {
! 252: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 253: return;
! 254: }
! 255:
! 256: for(i = 0; i < (**s_permutation).nombre_lignes; i++)
! 257: {
! 258: if ((((integer8 **) (**s_permutation).tableau)[i] =
! 259: (integer8 *) malloc((**s_permutation).nombre_colonnes *
! 260: sizeof(integer8))) == NULL)
! 261: {
! 262: (*s_etat_processus).erreur_systeme =
! 263: d_es_allocation_memoire;
! 264: return;
! 265: }
! 266:
! 267: for(j = 0; j < (**s_permutation).nombre_colonnes; j++)
! 268: {
! 269: ((integer8 **) (**s_permutation).tableau)[i][j] =
! 270: (j == i) ? 1 : 0;
! 271: }
! 272: }
! 273:
! 274: for(n = dimension_vecteur_pivot - 1; n >= 0; n--)
! 275: {
! 276: if ((pivot[n] - 1) != n)
! 277: {
! 278: tampon = ((integer8 **) (**s_permutation).tableau)[n];
! 279: ((integer8 **) (**s_permutation).tableau)[n] =
! 280: ((integer8 **) (**s_permutation).tableau)
! 281: [pivot[n] - 1];
! 282: ((integer8 **) (**s_permutation).tableau)[pivot[n] - 1] =
! 283: tampon;
! 284: }
! 285: }
! 286:
! 287: free(pivot);
! 288:
! 289: break;
! 290: }
! 291:
! 292: case 'C' :
! 293: {
! 294: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 295: sizeof(complex16))) == NULL)
! 296: {
! 297: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 298: return;
! 299: }
! 300:
! 301: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 302: {
! 303: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 304: {
! 305: ((complex16 *) matrice_f77)[k++] = ((complex16 **)
! 306: (*s_matrice).tableau)[j][i];
! 307: }
! 308: }
! 309:
! 310: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot *
! 311: sizeof(integer4))) == NULL)
! 312: {
! 313: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 314: return;
! 315: }
! 316:
! 317: zgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77,
! 318: &nombre_lignes_a, pivot, &erreur);
! 319:
! 320: if (erreur < 0)
! 321: {
! 322: (*s_etat_processus).erreur_execution =
! 323: d_ex_routines_mathematiques;
! 324: free(pivot);
! 325: free(matrice_f77);
! 326: return;
! 327: }
! 328:
! 329: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 330: {
! 331: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 332: {
! 333: ((complex16 **) (*s_matrice).tableau)[j][i] =
! 334: ((complex16 *) matrice_f77)[k++];
! 335: }
! 336: }
! 337:
! 338: free(matrice_f77);
! 339:
! 340: (**s_permutation).type = 'I';
! 341: (**s_permutation).nombre_lignes = dimension_vecteur_pivot;
! 342: (**s_permutation).nombre_colonnes = dimension_vecteur_pivot;
! 343:
! 344: if (((**s_permutation).tableau = malloc((**s_permutation)
! 345: .nombre_lignes * sizeof(integer8 *))) == NULL)
! 346: {
! 347: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 348: return;
! 349: }
! 350:
! 351: for(i = 0; i < (**s_permutation).nombre_lignes; i++)
! 352: {
! 353: if ((((integer8 **) (**s_permutation).tableau)[i] =
! 354: (integer8 *) malloc((**s_permutation).nombre_colonnes *
! 355: sizeof(integer8))) == NULL)
! 356: {
! 357: (*s_etat_processus).erreur_systeme =
! 358: d_es_allocation_memoire;
! 359: return;
! 360: }
! 361:
! 362: for(j = 0; j < (**s_permutation).nombre_colonnes; j++)
! 363: {
! 364: ((integer8 **) (**s_permutation).tableau)[i][j] =
! 365: (j == i) ? 1 : 0;
! 366: }
! 367: }
! 368:
! 369: for(n = dimension_vecteur_pivot - 1; n >= 0; n--)
! 370: {
! 371: if ((pivot[n] - 1) != n)
! 372: {
! 373: tampon = ((integer8 **) (**s_permutation).tableau)[n];
! 374: ((integer8 **) (**s_permutation).tableau)[n] =
! 375: ((integer8 **) (**s_permutation).tableau)
! 376: [pivot[n] - 1];
! 377: ((integer8 **) (**s_permutation).tableau)[pivot[n] - 1] =
! 378: tampon;
! 379: }
! 380: }
! 381:
! 382: free(pivot);
! 383:
! 384: break;
! 385: }
! 386: }
! 387:
! 388: return;
! 389: }
! 390:
! 391:
! 392: /*
! 393: ================================================================================
! 394: Fonction réalisation la factorisation de Cholesky d'une matrice symétrique
! 395: définie positive
! 396: ================================================================================
! 397: Entrées : struct_matrice, mode (U ou L selon la décomposition demandée)
! 398: --------------------------------------------------------------------------------
! 399: Sorties : décomposition de Cholesky de la matrice d'entrée et drapeau
! 400: d'erreur. La matrice en entrée est écrasée.
! 401: --------------------------------------------------------------------------------
! 402: Effets de bord : néant
! 403: ================================================================================
! 404: */
! 405:
! 406: void
! 407: factorisation_cholesky(struct_processus *s_etat_processus,
! 408: struct_matrice *s_matrice, unsigned char mode)
! 409: {
! 410: integer4 erreur;
! 411: integer4 i;
! 412: integer4 j;
! 413: integer4 nombre_colonnes_a;
! 414: integer4 nombre_lignes_a;
! 415:
! 416: unsigned long taille_matrice_f77;
! 417:
! 418: void *matrice_f77;
! 419:
! 420: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
! 421: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
! 422:
! 423: if (nombre_lignes_a != nombre_colonnes_a)
! 424: {
! 425: (*s_etat_processus).erreur_execution = d_ex_dimensions_invalides;
! 426: return;
! 427: }
! 428:
! 429: taille_matrice_f77 = (nombre_lignes_a * (nombre_colonnes_a + 1)) / 2;
! 430:
! 431: switch((*s_matrice).type)
! 432: {
! 433: case 'I' :
! 434: {
! 435: /*
! 436: * Allocation du vecteur représentant la matrice triangulaire
! 437: */
! 438:
! 439: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 440: sizeof(real8))) == NULL)
! 441: {
! 442: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 443: return;
! 444: }
! 445:
! 446: if (mode == 'U')
! 447: {
! 448: for(j = 1; j <= nombre_colonnes_a; j++)
! 449: {
! 450: for(i = 1; i <= j; i++)
! 451: {
! 452: ((real8 *) matrice_f77)[i + (((j - 1) * j) / 2) - 1] =
! 453: (real8) ((integer8 **) (*s_matrice).tableau)
! 454: [i - 1][j - 1];
! 455: }
! 456: }
! 457: }
! 458: else
! 459: {
! 460: for(j = 1; j <= nombre_colonnes_a; j++)
! 461: {
! 462: for(i = j; i <= nombre_colonnes_a; i++)
! 463: {
! 464: ((real8 *) matrice_f77)[(i - 1) + (((j - 1) *
! 465: ((nombre_colonnes_a * 2) - j)) / 2)] = (real8)
! 466: ((integer8 **) (*s_matrice).tableau)
! 467: [i - 1][j - 1];
! 468: }
! 469: }
! 470: }
! 471:
! 472: dpptrf_(&mode, &nombre_colonnes_a, matrice_f77, &erreur, 1);
! 473:
! 474: if (erreur != 0)
! 475: {
! 476: if (erreur > 0)
! 477: {
! 478: (*s_etat_processus).exception =
! 479: d_ep_matrice_non_definie_positive;
! 480: }
! 481: else
! 482: {
! 483: (*s_etat_processus).erreur_execution =
! 484: d_ex_routines_mathematiques;
! 485: }
! 486:
! 487: free(matrice_f77);
! 488: return;
! 489: }
! 490:
! 491: for(i = 0; i < (integer4) (*s_matrice).nombre_lignes; i++)
! 492: {
! 493: free(((integer8 **) (*s_matrice).tableau)[i]);
! 494: }
! 495:
! 496: free((integer8 **) (*s_matrice).tableau);
! 497:
! 498: if (((*s_matrice).tableau = (void **) malloc((*s_matrice)
! 499: .nombre_lignes * sizeof(real8 *))) == NULL)
! 500: {
! 501: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 502: return;
! 503: }
! 504:
! 505: for(i = 0; i < (integer4) (*s_matrice).nombre_lignes; i++)
! 506: {
! 507: if ((((*s_matrice).tableau)[i] =
! 508: (real8 *) malloc((*s_matrice)
! 509: .nombre_colonnes * sizeof(real8))) == NULL)
! 510: {
! 511: (*s_etat_processus).erreur_systeme =
! 512: d_es_allocation_memoire;
! 513: return;
! 514: }
! 515: }
! 516:
! 517: (*s_matrice).type = 'R';
! 518:
! 519: if (mode == 'U')
! 520: {
! 521: for(j = 1; j <= nombre_colonnes_a; j++)
! 522: {
! 523: for(i = 1; i <= j; i++)
! 524: {
! 525: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1] =
! 526: ((real8 *) matrice_f77)
! 527: [i + (((j - 1) * j) / 2) - 1];
! 528: }
! 529: }
! 530:
! 531: for(j = 1; j <= nombre_colonnes_a; j++)
! 532: {
! 533: for(i = j + 1; i <= nombre_colonnes_a; i++)
! 534: {
! 535: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1] = 0;
! 536: }
! 537: }
! 538: }
! 539: else
! 540: {
! 541: for(j = 1; j <= nombre_colonnes_a; j++)
! 542: {
! 543: for(i = j; i <= nombre_colonnes_a; i++)
! 544: {
! 545: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1] =
! 546: ((real8 *) matrice_f77)[(i - 1) + (((j - 1) *
! 547: ((nombre_colonnes_a * 2) - j)) / 2)];
! 548: }
! 549: }
! 550:
! 551: for(j = 1; j <= nombre_colonnes_a; j++)
! 552: {
! 553: for(i = 1; i < j; i++)
! 554: {
! 555: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1] = 0;
! 556: }
! 557: }
! 558: }
! 559:
! 560: free(matrice_f77);
! 561: break;
! 562: }
! 563:
! 564: case 'R' :
! 565: {
! 566: /*
! 567: * Allocation du vecteur représentant la matrice triangulaire
! 568: */
! 569:
! 570: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 571: sizeof(real8))) == NULL)
! 572: {
! 573: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 574: return;
! 575: }
! 576:
! 577: if (mode == 'U')
! 578: {
! 579: for(j = 1; j <= nombre_colonnes_a; j++)
! 580: {
! 581: for(i = 1; i <= j; i++)
! 582: {
! 583: ((real8 *) matrice_f77)[i + (((j - 1) * j) / 2) - 1] =
! 584: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1];
! 585: }
! 586: }
! 587: }
! 588: else
! 589: {
! 590: for(j = 1; j <= nombre_colonnes_a; j++)
! 591: {
! 592: for(i = j; i <= nombre_colonnes_a; i++)
! 593: {
! 594: ((real8 *) matrice_f77)[(i - 1) + (((j - 1) *
! 595: ((nombre_colonnes_a * 2) - j)) / 2)] = (real8)
! 596: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1];
! 597: }
! 598: }
! 599: }
! 600:
! 601: dpptrf_(&mode, &nombre_colonnes_a, matrice_f77, &erreur, 1);
! 602:
! 603: if (erreur != 0)
! 604: {
! 605: if (erreur > 0)
! 606: {
! 607: (*s_etat_processus).exception =
! 608: d_ep_matrice_non_definie_positive;
! 609: }
! 610: else
! 611: {
! 612: (*s_etat_processus).erreur_execution =
! 613: d_ex_routines_mathematiques;
! 614: }
! 615:
! 616: free(matrice_f77);
! 617: return;
! 618: }
! 619:
! 620: if (mode == 'U')
! 621: {
! 622: for(j = 1; j <= nombre_colonnes_a; j++)
! 623: {
! 624: for(i = 1; i <= j; i++)
! 625: {
! 626: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1] =
! 627: ((real8 *) matrice_f77)
! 628: [i + (((j - 1) * j) / 2) - 1];
! 629: }
! 630: }
! 631:
! 632: for(j = 1; j <= nombre_colonnes_a; j++)
! 633: {
! 634: for(i = j + 1; i <= nombre_colonnes_a; i++)
! 635: {
! 636: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1] = 0;
! 637: }
! 638: }
! 639: }
! 640: else
! 641: {
! 642: for(j = 1; j <= nombre_colonnes_a; j++)
! 643: {
! 644: for(i = j; i <= nombre_colonnes_a; i++)
! 645: {
! 646: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1] =
! 647: ((real8 *) matrice_f77)[(i - 1) + (((j - 1) *
! 648: ((nombre_colonnes_a * 2) - j)) / 2)];
! 649: }
! 650: }
! 651:
! 652: for(j = 1; j <= nombre_colonnes_a; j++)
! 653: {
! 654: for(i = 1; i < j; i++)
! 655: {
! 656: ((real8 **) (*s_matrice).tableau)[i - 1][j - 1] = 0;
! 657: }
! 658: }
! 659: }
! 660:
! 661: free(matrice_f77);
! 662: break;
! 663: }
! 664:
! 665: case 'C' :
! 666: {
! 667: /*
! 668: * Allocation du vecteur représentant la matrice triangulaire
! 669: */
! 670:
! 671: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 672: sizeof(complex16))) == NULL)
! 673: {
! 674: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 675: return;
! 676: }
! 677:
! 678: if (mode == 'U')
! 679: {
! 680: for(j = 1; j <= nombre_colonnes_a; j++)
! 681: {
! 682: for(i = 1; i <= j; i++)
! 683: {
! 684: ((complex16 *) matrice_f77)[i + (((j - 1) * j) / 2)
! 685: - 1].partie_reelle = ((complex16 **)
! 686: (*s_matrice).tableau)[i - 1][j - 1]
! 687: .partie_reelle;
! 688: ((complex16 *) matrice_f77)[i + (((j - 1) * j) / 2)
! 689: - 1].partie_imaginaire = ((complex16 **)
! 690: (*s_matrice).tableau)[i - 1][j - 1]
! 691: .partie_imaginaire;
! 692: }
! 693: }
! 694: }
! 695: else
! 696: {
! 697: for(j = 1; j <= nombre_colonnes_a; j++)
! 698: {
! 699: for(i = j; i <= nombre_colonnes_a; i++)
! 700: {
! 701: ((complex16 *) matrice_f77)[(i - 1) + (((j - 1) *
! 702: ((nombre_colonnes_a * 2) - j)) / 2)]
! 703: .partie_reelle = ((complex16 **)
! 704: (*s_matrice).tableau)[i - 1][j - 1]
! 705: .partie_reelle;
! 706: ((complex16 *) matrice_f77)[(i - 1) + (((j - 1) *
! 707: ((nombre_colonnes_a * 2) - j)) / 2)]
! 708: .partie_imaginaire = ((complex16 **)
! 709: (*s_matrice).tableau)[i - 1][j - 1]
! 710: .partie_imaginaire;
! 711: }
! 712: }
! 713: }
! 714:
! 715: zpptrf_(&mode, &nombre_colonnes_a, matrice_f77, &erreur, 1);
! 716:
! 717: if (erreur != 0)
! 718: {
! 719: if (erreur > 0)
! 720: {
! 721: (*s_etat_processus).exception =
! 722: d_ep_matrice_non_definie_positive;
! 723: }
! 724: else
! 725: {
! 726: (*s_etat_processus).erreur_execution =
! 727: d_ex_routines_mathematiques;
! 728: }
! 729:
! 730: free(matrice_f77);
! 731: return;
! 732: }
! 733:
! 734: if (mode == 'U')
! 735: {
! 736: for(j = 1; j <= nombre_colonnes_a; j++)
! 737: {
! 738: for(i = 1; i <= j; i++)
! 739: {
! 740: ((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
! 741: .partie_reelle = ((complex16 *) matrice_f77)
! 742: [i + (((j - 1) * j) / 2) - 1].partie_reelle;
! 743: ((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
! 744: .partie_imaginaire = ((complex16 *) matrice_f77)
! 745: [i + (((j - 1) * j) / 2) - 1].partie_imaginaire;
! 746: }
! 747: }
! 748:
! 749: for(j = 1; j <= nombre_colonnes_a; j++)
! 750: {
! 751: for(i = j + 1; i <= nombre_colonnes_a; i++)
! 752: {
! 753: ((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
! 754: .partie_reelle = 0;
! 755: ((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
! 756: .partie_imaginaire = 0;
! 757: }
! 758: }
! 759: }
! 760: else
! 761: {
! 762: for(j = 1; j <= nombre_colonnes_a; j++)
! 763: {
! 764: for(i = j; i <= nombre_colonnes_a; i++)
! 765: {
! 766: ((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
! 767: .partie_reelle = ((complex16 *)
! 768: matrice_f77)[(i - 1) + (((j - 1) *
! 769: ((nombre_colonnes_a * 2) - j)) / 2)]
! 770: .partie_reelle;
! 771: ((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
! 772: .partie_imaginaire = ((complex16 *)
! 773: matrice_f77)[(i - 1) + (((j - 1) *
! 774: ((nombre_colonnes_a * 2) - j)) / 2)]
! 775: .partie_imaginaire;
! 776: }
! 777: }
! 778:
! 779: for(j = 1; j <= nombre_colonnes_a; j++)
! 780: {
! 781: for(i = 1; i < j; i++)
! 782: {
! 783: ((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
! 784: .partie_reelle = 0;
! 785: ((complex16 **) (*s_matrice).tableau)[i - 1][j - 1]
! 786: .partie_imaginaire = 0;
! 787: }
! 788: }
! 789: }
! 790:
! 791: free(matrice_f77);
! 792: break;
! 793: }
! 794: }
! 795:
! 796: return;
! 797: }
! 798:
! 799: // vim: ts=4
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to algebre_lineaire2.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / src