Annotation of rpl/src/algebre_lineaire3.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 de Schur d'une matrice carrée
! 29: ================================================================================
! 30: Entrées : struct_matrice
! 31: --------------------------------------------------------------------------------
! 32: Sorties : décomposition de Schur de la matrice d'entrée et drapeau d'erreur.
! 33: La matrice en entrée est écrasée. La matrice de sortie est
! 34: la forme de Schur.
! 35: La routine renvoie aussi une matrice de complexes correspondant
! 36: aux vecteurs de Schur. Cette matrice est allouée par
! 37: la routine et vaut NULL sinon.
! 38: --------------------------------------------------------------------------------
! 39: Effets de bord : néant
! 40: ================================================================================
! 41: */
! 42:
! 43: void
! 44: factorisation_schur(struct_processus *s_etat_processus,
! 45: struct_matrice *s_matrice, struct_matrice **s_schur)
! 46: {
! 47: complex16 *w;
! 48:
! 49: integer4 info;
! 50: integer4 lwork;
! 51: integer4 nombre_colonnes_a;
! 52: integer4 nombre_lignes_a;
! 53: integer4 sdim;
! 54:
! 55: real8 *rwork;
! 56: real8 *wi;
! 57: real8 *wr;
! 58:
! 59: unsigned char calcul_vecteurs_schur;
! 60: unsigned char tri_vecteurs_schur;
! 61:
! 62: unsigned long i;
! 63: unsigned long j;
! 64: unsigned long k;
! 65: unsigned long taille_matrice_f77;
! 66:
! 67: void *matrice_a_f77;
! 68: void *matrice_vs_f77;
! 69: void *tampon;
! 70: void *work;
! 71:
! 72: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
! 73: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
! 74: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;
! 75:
! 76: calcul_vecteurs_schur = 'V';
! 77: tri_vecteurs_schur = 'N';
! 78:
! 79: switch((*s_matrice).type)
! 80: {
! 81: case 'I' :
! 82: {
! 83: /* Conversion de la matrice en matrice réelle */
! 84:
! 85: for(i = 0; i < (unsigned long) nombre_lignes_a; i++)
! 86: {
! 87: tampon = (*s_matrice).tableau[i];
! 88:
! 89: if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *)
! 90: malloc(nombre_colonnes_a * sizeof(real8))) == NULL)
! 91: {
! 92: (*s_etat_processus).erreur_systeme =
! 93: d_es_allocation_memoire;
! 94: return;
! 95: }
! 96:
! 97: for(j = 0; j < (unsigned long) nombre_colonnes_a; j++)
! 98: {
! 99: ((real8 **) (*s_matrice).tableau)[i][j] =
! 100: ((integer8 *) tampon)[j];
! 101: }
! 102:
! 103: free(tampon);
! 104: }
! 105:
! 106: (*s_matrice).type = 'R';
! 107: }
! 108:
! 109: case 'R' :
! 110: {
! 111: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 *
! 112: sizeof(real8))) == NULL)
! 113: {
! 114: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 115: return;
! 116: }
! 117:
! 118: if ((matrice_vs_f77 = (void *) malloc(taille_matrice_f77 *
! 119: sizeof(real8))) == NULL)
! 120: {
! 121: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 122: return;
! 123: }
! 124:
! 125: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 126: {
! 127: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 128: {
! 129: ((real8 *) matrice_a_f77)[k++] = ((real8 **)
! 130: (*s_matrice).tableau)[j][i];
! 131: }
! 132: }
! 133:
! 134: if ((wr = (real8 *) malloc(nombre_lignes_a * sizeof(real8)))
! 135: == NULL)
! 136: {
! 137: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 138: return;
! 139: }
! 140:
! 141: if ((wi = (real8 *) malloc(nombre_lignes_a * sizeof(real8)))
! 142: == NULL)
! 143: {
! 144: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 145: return;
! 146: }
! 147:
! 148: lwork = -1;
! 149:
! 150: if ((work = (real8 *) malloc(sizeof(real8))) == NULL)
! 151: {
! 152: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 153: return;
! 154: }
! 155:
! 156: dgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur,
! 157: NULL, &nombre_lignes_a, matrice_a_f77,
! 158: &nombre_colonnes_a, &sdim, wr, wi,
! 159: matrice_vs_f77, &nombre_colonnes_a,
! 160: work, &lwork, NULL, &info, 1, 1);
! 161:
! 162: lwork = ((real8 *) work)[0];
! 163: free(work);
! 164:
! 165: if ((work = (real8 *) malloc(lwork * sizeof(real8))) == NULL)
! 166: {
! 167: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 168: return;
! 169: }
! 170:
! 171: dgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur,
! 172: NULL, &nombre_lignes_a, matrice_a_f77,
! 173: &nombre_colonnes_a, &sdim, wr, wi,
! 174: matrice_vs_f77, &nombre_colonnes_a,
! 175: work, &lwork, NULL, &info, 1, 1);
! 176:
! 177: free(work);
! 178: free(wr);
! 179: free(wi);
! 180:
! 181: if (info != 0)
! 182: {
! 183: if (info > 0)
! 184: {
! 185: (*s_etat_processus).exception = d_ep_decomposition_QR;
! 186: }
! 187: else
! 188: {
! 189: (*s_etat_processus).erreur_execution =
! 190: d_ex_routines_mathematiques;
! 191: }
! 192:
! 193: free(matrice_a_f77);
! 194: free(matrice_vs_f77);
! 195: return;
! 196: }
! 197:
! 198: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 199: {
! 200: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 201: {
! 202: ((real8 **) (*s_matrice).tableau)[j][i] =
! 203: ((real8 *) matrice_a_f77)[k++];
! 204: }
! 205: }
! 206:
! 207: (**s_schur).nombre_colonnes = (*s_matrice).nombre_colonnes;
! 208: (**s_schur).nombre_lignes = (*s_matrice).nombre_lignes;
! 209: (**s_schur).type = 'R';
! 210:
! 211: if (((**s_schur).tableau = malloc((**s_schur)
! 212: .nombre_lignes * sizeof(real8 *))) == NULL)
! 213: {
! 214: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 215: return;
! 216: }
! 217:
! 218: for(i = 0; i < (**s_schur).nombre_lignes; i++)
! 219: {
! 220: if ((((real8 **) (**s_schur).tableau)[i] = (real8 *)
! 221: malloc((**s_schur).nombre_colonnes *
! 222: sizeof(real8))) == NULL)
! 223: {
! 224: (*s_etat_processus).erreur_systeme =
! 225: d_es_allocation_memoire;
! 226: return;
! 227: }
! 228: }
! 229:
! 230: for(k = 0, i = 0; i < (**s_schur).nombre_colonnes; i++)
! 231: {
! 232: for(j = 0; j < (**s_schur).nombre_lignes; j++)
! 233: {
! 234: ((real8 **) (**s_schur).tableau)[j][i] = ((real8 *)
! 235: matrice_vs_f77)[k++];
! 236: }
! 237: }
! 238:
! 239: free(matrice_a_f77);
! 240: free(matrice_vs_f77);
! 241:
! 242: break;
! 243: }
! 244:
! 245: case 'C' :
! 246: {
! 247: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 *
! 248: sizeof(complex16))) == NULL)
! 249: {
! 250: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 251: return;
! 252: }
! 253:
! 254: if ((matrice_vs_f77 = (void *) malloc(taille_matrice_f77 *
! 255: sizeof(complex16))) == NULL)
! 256: {
! 257: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 258: return;
! 259: }
! 260:
! 261: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 262: {
! 263: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 264: {
! 265: ((complex16 *) matrice_a_f77)[k].partie_reelle =
! 266: ((complex16 **) (*s_matrice).tableau)[j][i]
! 267: .partie_reelle;
! 268: ((complex16 *) matrice_a_f77)[k++].partie_imaginaire =
! 269: ((complex16 **) (*s_matrice).tableau)[j][i]
! 270: .partie_imaginaire;
! 271: }
! 272: }
! 273:
! 274: if ((w = (complex16 *) malloc(nombre_lignes_a * sizeof(complex16)))
! 275: == NULL)
! 276: {
! 277: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 278: return;
! 279: }
! 280:
! 281: lwork = -1;
! 282:
! 283: if ((work = (complex16 *) malloc(sizeof(complex16))) == NULL)
! 284: {
! 285: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 286: return;
! 287: }
! 288:
! 289: if ((rwork = (real8 *) malloc(nombre_lignes_a * sizeof(real8)))
! 290: == NULL)
! 291: {
! 292: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 293: return;
! 294: }
! 295:
! 296: zgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur,
! 297: NULL, &nombre_lignes_a, matrice_a_f77,
! 298: &nombre_colonnes_a, &sdim, w,
! 299: matrice_vs_f77, &nombre_colonnes_a,
! 300: work, &lwork, rwork, NULL, &info, 1, 1);
! 301:
! 302: lwork = ((complex16 *) work)[0].partie_reelle;
! 303: free(work);
! 304:
! 305: if ((work = (complex16 *) malloc(lwork * sizeof(complex16)))
! 306: == NULL)
! 307: {
! 308: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 309: return;
! 310: }
! 311:
! 312: zgees_(&calcul_vecteurs_schur, &tri_vecteurs_schur,
! 313: NULL, &nombre_lignes_a, matrice_a_f77,
! 314: &nombre_colonnes_a, &sdim, w,
! 315: matrice_vs_f77, &nombre_colonnes_a,
! 316: work, &lwork, rwork, NULL, &info, 1, 1);
! 317:
! 318: free(work);
! 319: free(rwork);
! 320: free(w);
! 321:
! 322: if (info != 0)
! 323: {
! 324: if (info > 0)
! 325: {
! 326: (*s_etat_processus).exception = d_ep_decomposition_QR;
! 327: }
! 328: else
! 329: {
! 330: (*s_etat_processus).erreur_execution =
! 331: d_ex_routines_mathematiques;
! 332: }
! 333:
! 334: free(matrice_a_f77);
! 335: free(matrice_vs_f77);
! 336: return;
! 337: }
! 338:
! 339: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 340: {
! 341: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 342: {
! 343: ((complex16 **) (*s_matrice).tableau)[j][i]
! 344: .partie_reelle = ((complex16 *) matrice_a_f77)[k]
! 345: .partie_reelle;
! 346: ((complex16 **) (*s_matrice).tableau)[j][i]
! 347: .partie_imaginaire = ((complex16 *) matrice_a_f77)
! 348: [k++].partie_imaginaire;
! 349: }
! 350: }
! 351:
! 352: (**s_schur).nombre_colonnes = (*s_matrice).nombre_colonnes;
! 353: (**s_schur).nombre_lignes = (*s_matrice).nombre_lignes;
! 354: (**s_schur).type = 'C';
! 355:
! 356: if (((**s_schur).tableau = malloc((**s_schur)
! 357: .nombre_lignes * sizeof(complex16 *))) == NULL)
! 358: {
! 359: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 360: return;
! 361: }
! 362:
! 363: for(i = 0; i < (**s_schur).nombre_lignes; i++)
! 364: {
! 365: if ((((complex16 **) (**s_schur).tableau)[i] = (complex16 *)
! 366: malloc((**s_schur).nombre_colonnes *
! 367: sizeof(complex16))) == NULL)
! 368: {
! 369: (*s_etat_processus).erreur_systeme =
! 370: d_es_allocation_memoire;
! 371: return;
! 372: }
! 373: }
! 374:
! 375: for(k = 0, i = 0; i < (**s_schur).nombre_colonnes; i++)
! 376: {
! 377: for(j = 0; j < (**s_schur).nombre_lignes; j++)
! 378: {
! 379: ((complex16 **) (**s_schur).tableau)[j][i].partie_reelle =
! 380: ((complex16 *) matrice_vs_f77)[k].partie_reelle;
! 381: ((complex16 **) (**s_schur).tableau)[j][i]
! 382: .partie_imaginaire = ((complex16 *) matrice_vs_f77)
! 383: [k++].partie_imaginaire;
! 384: }
! 385: }
! 386:
! 387: free(matrice_a_f77);
! 388: free(matrice_vs_f77);
! 389:
! 390: break;
! 391: }
! 392: }
! 393:
! 394: return;
! 395: }
! 396:
! 397:
! 398: /*
! 399: ================================================================================
! 400: Fonction réalisation la factorisation LQ d'une matrice quelconque
! 401: ================================================================================
! 402: Entrées : struct_matrice
! 403: --------------------------------------------------------------------------------
! 404: Sorties : décomposition de LQ de la matrice d'entrée. Le tableau tau
! 405: est initialisé par la fonction
! 406: --------------------------------------------------------------------------------
! 407: Effets de bord : néant
! 408: ================================================================================
! 409: */
! 410:
! 411: void
! 412: factorisation_lq(struct_processus *s_etat_processus, struct_matrice *s_matrice,
! 413: void **tau)
! 414: {
! 415: integer4 nombre_colonnes_a;
! 416: integer4 nombre_lignes_a;
! 417: integer4 erreur;
! 418:
! 419: unsigned long i;
! 420: unsigned long j;
! 421: unsigned long k;
! 422: unsigned long taille_matrice_f77;
! 423:
! 424: void *matrice_a_f77;
! 425: void *tampon;
! 426: void *work;
! 427:
! 428: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
! 429: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
! 430: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;
! 431:
! 432: switch((*s_matrice).type)
! 433: {
! 434: case 'I' :
! 435: {
! 436: /* Conversion de la matrice en matrice réelle */
! 437:
! 438: for(i = 0; i < (unsigned long) nombre_lignes_a; i++)
! 439: {
! 440: tampon = (*s_matrice).tableau[i];
! 441:
! 442: if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *)
! 443: malloc(nombre_colonnes_a * sizeof(real8))) == NULL)
! 444: {
! 445: (*s_etat_processus).erreur_systeme =
! 446: d_es_allocation_memoire;
! 447: return;
! 448: }
! 449:
! 450: for(j = 0; j < (unsigned long) nombre_colonnes_a; j++)
! 451: {
! 452: ((real8 **) (*s_matrice).tableau)[i][j] =
! 453: ((integer8 *) tampon)[j];
! 454: }
! 455:
! 456: free(tampon);
! 457: }
! 458:
! 459: (*s_matrice).type = 'R';
! 460: }
! 461:
! 462: case 'R' :
! 463: {
! 464: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 *
! 465: sizeof(real8))) == NULL)
! 466: {
! 467: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 468: return;
! 469: }
! 470:
! 471: if (((*tau) = malloc(((nombre_colonnes_a < nombre_lignes_a)
! 472: ? nombre_colonnes_a : nombre_lignes_a) * sizeof(real8)))
! 473: == NULL)
! 474: {
! 475: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 476: return;
! 477: }
! 478:
! 479: if ((work = malloc(nombre_lignes_a * sizeof(real8))) == NULL)
! 480: {
! 481: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 482: return;
! 483: }
! 484:
! 485: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 486: {
! 487: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 488: {
! 489: ((real8 *) matrice_a_f77)[k++] = ((real8 **)
! 490: (*s_matrice).tableau)[j][i];
! 491: }
! 492: }
! 493:
! 494: dgelq2_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
! 495: &nombre_lignes_a, (*((real8 **) tau)), work, &erreur);
! 496:
! 497: if (erreur != 0)
! 498: {
! 499: // L'erreur ne peut être que négative.
! 500:
! 501: (*s_etat_processus).erreur_execution =
! 502: d_ex_routines_mathematiques;
! 503: free(work);
! 504: free(matrice_a_f77);
! 505: return;
! 506: }
! 507:
! 508: for(k = 0, i = 0; i < (unsigned long) nombre_colonnes_a; i++)
! 509: {
! 510: for(j = 0; j < (unsigned long) nombre_lignes_a; j++)
! 511: {
! 512: ((real8 **) (*s_matrice).tableau)[j][i] = ((real8 *)
! 513: matrice_a_f77)[k++];
! 514: }
! 515: }
! 516:
! 517: free(work);
! 518: free(matrice_a_f77);
! 519:
! 520: break;
! 521: }
! 522:
! 523: case 'C' :
! 524: {
! 525: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 *
! 526: sizeof(complex16))) == NULL)
! 527: {
! 528: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 529: return;
! 530: }
! 531:
! 532: if (((*tau) = malloc(((nombre_colonnes_a < nombre_lignes_a)
! 533: ? nombre_colonnes_a : nombre_lignes_a) *
! 534: sizeof(complex16))) == NULL)
! 535: {
! 536: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 537: return;
! 538: }
! 539:
! 540: if ((work = malloc(nombre_lignes_a * sizeof(complex16))) == NULL)
! 541: {
! 542: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 543: return;
! 544: }
! 545:
! 546: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 547: {
! 548: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 549: {
! 550: ((complex16 *) matrice_a_f77)[k].partie_reelle =
! 551: ((complex16 **) (*s_matrice).tableau)[j][i]
! 552: .partie_reelle;
! 553: ((complex16 *) matrice_a_f77)[k++].partie_imaginaire =
! 554: ((complex16 **) (*s_matrice).tableau)[j][i]
! 555: .partie_imaginaire;
! 556: }
! 557: }
! 558:
! 559: zgelq2_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
! 560: &nombre_lignes_a, (*((complex16 **) tau)), work, &erreur);
! 561:
! 562: if (erreur != 0)
! 563: {
! 564: // L'erreur ne peut être que négative.
! 565:
! 566: (*s_etat_processus).erreur_execution =
! 567: d_ex_routines_mathematiques;
! 568: free(work);
! 569: free(matrice_a_f77);
! 570: return;
! 571: }
! 572:
! 573: for(k = 0, i = 0; i < (unsigned long) nombre_colonnes_a; i++)
! 574: {
! 575: for(j = 0; j < (unsigned long) nombre_lignes_a; j++)
! 576: {
! 577: ((complex16 **) (*s_matrice).tableau)[j][i].partie_reelle =
! 578: ((complex16 *) matrice_a_f77)[k].partie_reelle;
! 579: ((complex16 **) (*s_matrice).tableau)[j][i]
! 580: .partie_imaginaire = ((complex16 *) matrice_a_f77)
! 581: [k++].partie_imaginaire;
! 582: }
! 583: }
! 584:
! 585: free(work);
! 586: free(matrice_a_f77);
! 587:
! 588: break;
! 589: }
! 590: }
! 591:
! 592: return;
! 593: }
! 594:
! 595:
! 596: /*
! 597: ================================================================================
! 598: Fonction réalisation la factorisation QR d'une matrice quelconque
! 599: ================================================================================
! 600: Entrées : struct_matrice
! 601: --------------------------------------------------------------------------------
! 602: Sorties : décomposition de QR de la matrice d'entrée. Le tableau tau
! 603: est initialisé par la fonction
! 604: --------------------------------------------------------------------------------
! 605: Effets de bord : néant
! 606: ================================================================================
! 607: */
! 608:
! 609: void
! 610: factorisation_qr(struct_processus *s_etat_processus, struct_matrice *s_matrice,
! 611: void **tau)
! 612: {
! 613: integer4 lwork;
! 614: integer4 nombre_colonnes_a;
! 615: integer4 nombre_lignes_a;
! 616: integer4 erreur;
! 617: integer4 *pivot;
! 618:
! 619: real8 *rwork;
! 620:
! 621: unsigned long i;
! 622: unsigned long j;
! 623: unsigned long k;
! 624: unsigned long taille_matrice_f77;
! 625:
! 626: void *matrice_a_f77;
! 627: void *registre;
! 628: void *tampon;
! 629: void *work;
! 630:
! 631: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
! 632: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
! 633: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;
! 634:
! 635: switch((*s_matrice).type)
! 636: {
! 637: case 'I' :
! 638: {
! 639: /* Conversion de la matrice en matrice réelle */
! 640:
! 641: for(i = 0; i < (unsigned long) nombre_lignes_a; i++)
! 642: {
! 643: tampon = (*s_matrice).tableau[i];
! 644:
! 645: if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *)
! 646: malloc(nombre_colonnes_a * sizeof(real8))) == NULL)
! 647: {
! 648: (*s_etat_processus).erreur_systeme =
! 649: d_es_allocation_memoire;
! 650: return;
! 651: }
! 652:
! 653: for(j = 0; j < (unsigned long) nombre_colonnes_a; j++)
! 654: {
! 655: ((real8 **) (*s_matrice).tableau)[i][j] =
! 656: ((integer8 *) tampon)[j];
! 657: }
! 658:
! 659: free(tampon);
! 660: }
! 661:
! 662: (*s_matrice).type = 'R';
! 663: }
! 664:
! 665: case 'R' :
! 666: {
! 667: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 *
! 668: sizeof(real8))) == NULL)
! 669: {
! 670: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 671: return;
! 672: }
! 673:
! 674: if (((*tau) = malloc(((nombre_colonnes_a < nombre_lignes_a)
! 675: ? nombre_colonnes_a : nombre_lignes_a) * sizeof(real8)))
! 676: == NULL)
! 677: {
! 678: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 679: return;
! 680: }
! 681:
! 682: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 683: {
! 684: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 685: {
! 686: ((real8 *) matrice_a_f77)[k++] = ((real8 **)
! 687: (*s_matrice).tableau)[j][i];
! 688: }
! 689: }
! 690:
! 691: if ((pivot = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL)
! 692: {
! 693: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 694: return;
! 695: }
! 696:
! 697: for(i = 0; i < (unsigned long) nombre_colonnes_a; pivot[i++] = 0);
! 698:
! 699: lwork = -1;
! 700:
! 701: if ((work = malloc(sizeof(real8))) == NULL)
! 702: {
! 703: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 704: return;
! 705: }
! 706:
! 707: // Calcul de la taille de l'espace
! 708:
! 709: dgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
! 710: &nombre_lignes_a, pivot, (*((real8 **) tau)),
! 711: work, &lwork, &erreur);
! 712:
! 713: lwork = ((real8 *) work)[0];
! 714:
! 715: free(work);
! 716:
! 717: if ((work = malloc(lwork * sizeof(real8))) == NULL)
! 718: {
! 719: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 720: return;
! 721: }
! 722:
! 723: // Calcul de la permutation
! 724:
! 725: if ((registre = (void *) malloc(taille_matrice_f77 *
! 726: sizeof(real8))) == NULL)
! 727: {
! 728: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 729: return;
! 730: }
! 731:
! 732: memcpy(registre, matrice_a_f77, taille_matrice_f77 * sizeof(real8));
! 733:
! 734: dgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, registre,
! 735: &nombre_lignes_a, pivot, (*((real8 **) tau)),
! 736: work, &lwork, &erreur);
! 737:
! 738: free(registre);
! 739:
! 740: if (erreur != 0)
! 741: {
! 742: // L'erreur ne peut être que négative.
! 743:
! 744: (*s_etat_processus).erreur_execution =
! 745: d_ex_routines_mathematiques;
! 746: free(work);
! 747: free(matrice_a_f77);
! 748: free(tau);
! 749: return;
! 750: }
! 751:
! 752: // La permutation doit maintenant être unitaire
! 753:
! 754: dgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
! 755: &nombre_lignes_a, pivot, (*((real8 **) tau)),
! 756: work, &lwork, &erreur);
! 757:
! 758: if (erreur != 0)
! 759: {
! 760: // L'erreur ne peut être que négative.
! 761:
! 762: (*s_etat_processus).erreur_execution =
! 763: d_ex_routines_mathematiques;
! 764: free(work);
! 765: free(matrice_a_f77);
! 766: free(tau);
! 767: return;
! 768: }
! 769:
! 770: for(i = 0; i < (unsigned long) nombre_colonnes_a; i++)
! 771: {
! 772: if ((i + 1) != (unsigned long) pivot[i])
! 773: {
! 774: (*s_etat_processus).erreur_execution =
! 775: d_ex_routines_mathematiques;
! 776:
! 777: free(pivot);
! 778: free(work);
! 779: free(matrice_a_f77);
! 780: free(tau);
! 781:
! 782: return;
! 783: }
! 784: }
! 785:
! 786: free(pivot);
! 787:
! 788: for(k = 0, i = 0; i < (unsigned long) nombre_colonnes_a; i++)
! 789: {
! 790: for(j = 0; j < (unsigned long) nombre_lignes_a; j++)
! 791: {
! 792: ((real8 **) (*s_matrice).tableau)[j][i] = ((real8 *)
! 793: matrice_a_f77)[k++];
! 794: }
! 795: }
! 796:
! 797: free(work);
! 798: free(matrice_a_f77);
! 799:
! 800: break;
! 801: }
! 802:
! 803: case 'C' :
! 804: {
! 805: if ((matrice_a_f77 = (void *) malloc(taille_matrice_f77 *
! 806: sizeof(complex16))) == NULL)
! 807: {
! 808: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 809: return;
! 810: }
! 811:
! 812: if (((*tau) = malloc(((nombre_colonnes_a < nombre_lignes_a)
! 813: ? nombre_colonnes_a : nombre_lignes_a) * sizeof(complex16)))
! 814: == NULL)
! 815: {
! 816: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 817: return;
! 818: }
! 819:
! 820: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 821: {
! 822: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 823: {
! 824: ((complex16 *) matrice_a_f77)[k].partie_reelle =
! 825: ((complex16 **) (*s_matrice).tableau)[j][i]
! 826: .partie_reelle;
! 827: ((complex16 *) matrice_a_f77)[k++].partie_imaginaire =
! 828: ((complex16 **) (*s_matrice).tableau)[j][i]
! 829: .partie_imaginaire;
! 830: }
! 831: }
! 832:
! 833: if ((pivot = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL)
! 834: {
! 835: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 836: return;
! 837: }
! 838:
! 839: if ((rwork = malloc(2 * nombre_colonnes_a * sizeof(real8))) == NULL)
! 840: {
! 841: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 842: return;
! 843: }
! 844:
! 845: for(i = 0; i < (unsigned long) nombre_colonnes_a; pivot[i++] = 0);
! 846:
! 847: lwork = -1;
! 848:
! 849: if ((work = malloc(sizeof(complex16))) == NULL)
! 850: {
! 851: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 852: return;
! 853: }
! 854:
! 855: // Calcul de la taille de l'espace
! 856:
! 857: zgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
! 858: &nombre_lignes_a, pivot, (*((complex16 **) tau)),
! 859: work, &lwork, rwork, &erreur);
! 860:
! 861: if (erreur != 0)
! 862: {
! 863: // L'erreur ne peut être que négative.
! 864:
! 865: (*s_etat_processus).erreur_execution =
! 866: d_ex_routines_mathematiques;
! 867:
! 868: free(work);
! 869: free(rwork);
! 870: free(pivot);
! 871: free(matrice_a_f77);
! 872: free(tau);
! 873: return;
! 874: }
! 875:
! 876: lwork = ((complex16 *) work)[0].partie_reelle;
! 877:
! 878: free(work);
! 879:
! 880: if ((work = malloc(lwork * sizeof(complex16))) == NULL)
! 881: {
! 882: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 883: return;
! 884: }
! 885:
! 886: // Calcul de la permutation
! 887:
! 888: if ((registre = (void *) malloc(taille_matrice_f77 *
! 889: sizeof(complex16))) == NULL)
! 890: {
! 891: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 892: return;
! 893: }
! 894:
! 895: memcpy(registre, matrice_a_f77,
! 896: taille_matrice_f77 * sizeof(complex16));
! 897:
! 898: zgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, registre,
! 899: &nombre_lignes_a, pivot, (*((complex16 **) tau)),
! 900: work, &lwork, rwork, &erreur);
! 901:
! 902: free(registre);
! 903:
! 904: if (erreur != 0)
! 905: {
! 906: // L'erreur ne peut être que négative.
! 907:
! 908: (*s_etat_processus).erreur_execution =
! 909: d_ex_routines_mathematiques;
! 910:
! 911: free(work);
! 912: free(rwork);
! 913: free(pivot);
! 914: free(matrice_a_f77);
! 915: free(tau);
! 916: return;
! 917: }
! 918:
! 919: // La permutation doit maintenant être unitaire
! 920:
! 921: zgeqp3_(&nombre_lignes_a, &nombre_colonnes_a, matrice_a_f77,
! 922: &nombre_lignes_a, pivot, (*((complex16 **) tau)),
! 923: work, &lwork, rwork, &erreur);
! 924:
! 925: if (erreur != 0)
! 926: {
! 927: // L'erreur ne peut être que négative.
! 928:
! 929: (*s_etat_processus).erreur_execution =
! 930: d_ex_routines_mathematiques;
! 931:
! 932: free(work);
! 933: free(rwork);
! 934: free(pivot);
! 935: free(matrice_a_f77);
! 936: free(tau);
! 937: return;
! 938: }
! 939:
! 940: free(rwork);
! 941:
! 942: for(i = 0; i < (unsigned long) nombre_colonnes_a; i++)
! 943: {
! 944: if ((i + 1) != (unsigned long) pivot[i])
! 945: {
! 946: (*s_etat_processus).erreur_execution =
! 947: d_ex_routines_mathematiques;
! 948:
! 949: free(pivot);
! 950: free(work);
! 951: free(matrice_a_f77);
! 952: free(tau);
! 953:
! 954: return;
! 955: }
! 956: }
! 957:
! 958: free(pivot);
! 959:
! 960: for(k = 0, i = 0; i < (unsigned long) nombre_colonnes_a; i++)
! 961: {
! 962: for(j = 0; j < (unsigned long) nombre_lignes_a; j++)
! 963: {
! 964: ((complex16 **) (*s_matrice).tableau)[j][i].partie_reelle =
! 965: ((complex16 *) matrice_a_f77)[k].partie_reelle;
! 966: ((complex16 **) (*s_matrice).tableau)[j][i]
! 967: .partie_imaginaire = ((complex16 *)
! 968: matrice_a_f77)[k++].partie_imaginaire;
! 969: }
! 970: }
! 971:
! 972: free(work);
! 973: free(matrice_a_f77);
! 974:
! 975: break;
! 976: }
! 977: }
! 978:
! 979: return;
! 980: }
! 981:
! 982:
! 983: /*
! 984: ================================================================================
! 985: Fonctions calculant le déterminant ou le rang d'une matrice quelconque
! 986: ================================================================================
! 987: Entrées : struct_matrice
! 988: --------------------------------------------------------------------------------
! 989: Sorties : déterminant
! 990: --------------------------------------------------------------------------------
! 991: Effets de bord : néant
! 992: ================================================================================
! 993: */
! 994:
! 995:
! 996: static integer4
! 997: calcul_rang(struct_processus *s_etat_processus, void *matrice_f77,
! 998: integer4 nombre_lignes_a, integer4 nombre_colonnes_a,
! 999: integer4 *pivot, integer4 dimension_vecteur_pivot, unsigned char type)
! 1000: {
! 1001: integer4 erreur;
! 1002: integer4 *iwork;
! 1003: integer4 *jpvt;
! 1004: integer4 lwork;
! 1005: integer4 longueur;
! 1006: integer4 nombre_colonnes_b;
! 1007: integer4 nombre_lignes_b;
! 1008: integer4 ordre;
! 1009: integer4 rang;
! 1010:
! 1011: real8 anorme;
! 1012: real8 rcond;
! 1013: real8 *rwork;
! 1014:
! 1015: unsigned char norme;
! 1016:
! 1017: unsigned long i;
! 1018:
! 1019: void *matrice_b;
! 1020: void *matrice_c;
! 1021: void *work;
! 1022:
! 1023: #undef NORME_I
! 1024: #ifdef NORME_I
! 1025: norme = 'I';
! 1026:
! 1027: if ((work = malloc(nombre_lignes_a * sizeof(real8))) == NULL)
! 1028: {
! 1029: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1030: return(-1);
! 1031: }
! 1032: #else
! 1033: norme = '1';
! 1034: work = NULL;
! 1035: #endif
! 1036:
! 1037: longueur = 1;
! 1038:
! 1039: if (type == 'R')
! 1040: {
! 1041: anorme = dlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a,
! 1042: matrice_f77, &nombre_lignes_a, work, longueur);
! 1043:
! 1044: #ifndef NORME_I
! 1045: free(work);
! 1046: #endif
! 1047:
! 1048: if ((matrice_c = malloc(nombre_lignes_a * nombre_colonnes_a *
! 1049: sizeof(real8))) == NULL)
! 1050: {
! 1051: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1052: return(-1);
! 1053: }
! 1054:
! 1055: memcpy(matrice_c, matrice_f77, nombre_lignes_a * nombre_colonnes_a *
! 1056: sizeof(real8));
! 1057:
! 1058: dgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77,
! 1059: &nombre_lignes_a, pivot, &erreur);
! 1060:
! 1061: if (erreur < 0)
! 1062: {
! 1063: (*s_etat_processus).erreur_execution =
! 1064: d_ex_routines_mathematiques;
! 1065:
! 1066: free(matrice_f77);
! 1067: return(-1);
! 1068: }
! 1069:
! 1070: if ((iwork = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL)
! 1071: {
! 1072: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1073: return(-1);
! 1074: }
! 1075:
! 1076: if ((work = malloc(4 * nombre_colonnes_a * sizeof(real8))) == NULL)
! 1077: {
! 1078: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1079: return(-1);
! 1080: }
! 1081:
! 1082: ordre = (nombre_lignes_a > nombre_colonnes_a)
! 1083: ? nombre_colonnes_a : nombre_lignes_a;
! 1084:
! 1085: dgecon_(&norme, &ordre, matrice_f77,
! 1086: &nombre_lignes_a, &anorme, &rcond, work, iwork, &erreur,
! 1087: longueur);
! 1088:
! 1089: free(work);
! 1090: free(iwork);
! 1091:
! 1092: if (erreur < 0)
! 1093: {
! 1094: (*s_etat_processus).erreur_execution =
! 1095: d_ex_routines_mathematiques;
! 1096:
! 1097: free(matrice_f77);
! 1098: return(-1);
! 1099: }
! 1100:
! 1101: if ((jpvt = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL)
! 1102: {
! 1103: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1104: return(-1);
! 1105: }
! 1106:
! 1107: for(i = 0; i < (unsigned long) nombre_colonnes_a; i++)
! 1108: {
! 1109: ((integer4 *) jpvt)[i] = 0;
! 1110: }
! 1111:
! 1112: lwork = -1;
! 1113:
! 1114: if ((work = malloc(sizeof(real8))) == NULL)
! 1115: {
! 1116: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1117: return(-1);
! 1118: }
! 1119:
! 1120: nombre_colonnes_b = 1;
! 1121: nombre_lignes_b = (nombre_lignes_a > nombre_colonnes_a)
! 1122: ? nombre_lignes_a : nombre_colonnes_a;
! 1123:
! 1124: if ((matrice_b = malloc(nombre_lignes_b * sizeof(real8))) == NULL)
! 1125: {
! 1126: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1127: return(-1);
! 1128: }
! 1129:
! 1130: for(i = 0; i < (unsigned long) nombre_lignes_b; i++)
! 1131: {
! 1132: ((real8 *) matrice_b)[i] = 0;
! 1133: }
! 1134:
! 1135: dgelsy_(&nombre_lignes_a, &nombre_colonnes_a,
! 1136: &nombre_colonnes_b, matrice_c, &nombre_lignes_a,
! 1137: matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang,
! 1138: work, &lwork, &erreur);
! 1139:
! 1140: lwork = ((real8 *) work)[0];
! 1141: free(work);
! 1142:
! 1143: if ((work = malloc(lwork * sizeof(real8))) == NULL)
! 1144: {
! 1145: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1146: return(-1);
! 1147: }
! 1148:
! 1149: dgelsy_(&nombre_lignes_a, &nombre_colonnes_a,
! 1150: &nombre_colonnes_b, matrice_c, &nombre_lignes_a,
! 1151: matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang,
! 1152: work, &lwork, &erreur);
! 1153:
! 1154: free(matrice_b);
! 1155: free(matrice_c);
! 1156: free(work);
! 1157: free(jpvt);
! 1158:
! 1159: if (erreur < 0)
! 1160: {
! 1161: (*s_etat_processus).erreur_execution =
! 1162: d_ex_routines_mathematiques;
! 1163:
! 1164: free(matrice_f77);
! 1165: return(-1);
! 1166: }
! 1167: }
! 1168: else
! 1169: {
! 1170: anorme = zlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a,
! 1171: matrice_f77, &nombre_lignes_a, work, longueur);
! 1172:
! 1173: #ifndef NORME_I
! 1174: free(work);
! 1175: #endif
! 1176:
! 1177: if ((matrice_c = malloc(nombre_lignes_a * nombre_colonnes_a *
! 1178: sizeof(complex16))) == NULL)
! 1179: {
! 1180: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1181: return(-1);
! 1182: }
! 1183:
! 1184: memcpy(matrice_c, matrice_f77, nombre_lignes_a * nombre_colonnes_a *
! 1185: sizeof(complex16));
! 1186:
! 1187: zgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77,
! 1188: &nombre_lignes_a, pivot, &erreur);
! 1189:
! 1190: if (erreur < 0)
! 1191: {
! 1192: (*s_etat_processus).erreur_execution =
! 1193: d_ex_routines_mathematiques;
! 1194:
! 1195: free(matrice_f77);
! 1196: return(-1);
! 1197: }
! 1198:
! 1199: if ((rwork = malloc(2 * nombre_colonnes_a * sizeof(real8))) == NULL)
! 1200: {
! 1201: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1202: return(-1);
! 1203: }
! 1204:
! 1205: if ((work = malloc(2 * nombre_colonnes_a * sizeof(complex16))) == NULL)
! 1206: {
! 1207: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1208: return(-1);
! 1209: }
! 1210:
! 1211: ordre = (nombre_lignes_a > nombre_colonnes_a)
! 1212: ? nombre_colonnes_a : nombre_lignes_a;
! 1213:
! 1214: zgecon_(&norme, &ordre, matrice_f77,
! 1215: &nombre_lignes_a, &anorme, &rcond, work, rwork, &erreur,
! 1216: longueur);
! 1217:
! 1218: free(work);
! 1219:
! 1220: if (erreur < 0)
! 1221: {
! 1222: (*s_etat_processus).erreur_execution =
! 1223: d_ex_routines_mathematiques;
! 1224:
! 1225: free(matrice_f77);
! 1226: return(-1);
! 1227: }
! 1228:
! 1229: if ((jpvt = malloc(nombre_colonnes_a * sizeof(integer4))) == NULL)
! 1230: {
! 1231: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1232: return(-1);
! 1233: }
! 1234:
! 1235: for(i = 0; i < (unsigned long) nombre_colonnes_a; i++)
! 1236: {
! 1237: ((integer4 *) jpvt)[i] = 0;
! 1238: }
! 1239:
! 1240: lwork = -1;
! 1241:
! 1242: if ((work = malloc(sizeof(complex16))) == NULL)
! 1243: {
! 1244: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1245: return(-1);
! 1246: }
! 1247:
! 1248: nombre_colonnes_b = 1;
! 1249: nombre_lignes_b = (nombre_lignes_a > nombre_colonnes_a)
! 1250: ? nombre_lignes_a : nombre_colonnes_a;
! 1251:
! 1252: if ((matrice_b = malloc(nombre_lignes_b * sizeof(complex16))) == NULL)
! 1253: {
! 1254: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1255: return(-1);
! 1256: }
! 1257:
! 1258: for(i = 0; i < (unsigned long) nombre_lignes_b; i++)
! 1259: {
! 1260: ((complex16 *) matrice_b)[i].partie_reelle = 0;
! 1261: ((complex16 *) matrice_b)[i].partie_imaginaire = 0;
! 1262: }
! 1263:
! 1264: zgelsy_(&nombre_lignes_a, &nombre_colonnes_a,
! 1265: &nombre_colonnes_b, matrice_c, &nombre_lignes_a,
! 1266: matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang,
! 1267: work, &lwork, rwork, &erreur);
! 1268:
! 1269: lwork = ((complex16 *) work)[0].partie_reelle;
! 1270: free(work);
! 1271:
! 1272: if ((work = malloc(lwork * sizeof(complex16))) == NULL)
! 1273: {
! 1274: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1275: return(-1);
! 1276: }
! 1277:
! 1278: zgelsy_(&nombre_lignes_a, &nombre_colonnes_a,
! 1279: &nombre_colonnes_b, matrice_c, &nombre_lignes_a,
! 1280: matrice_b, &nombre_lignes_b, jpvt, &rcond, &rang,
! 1281: work, &lwork, rwork, &erreur);
! 1282:
! 1283: free(rwork);
! 1284: free(matrice_b);
! 1285: free(matrice_c);
! 1286: free(work);
! 1287: free(jpvt);
! 1288:
! 1289: if (erreur < 0)
! 1290: {
! 1291: (*s_etat_processus).erreur_execution =
! 1292: d_ex_routines_mathematiques;
! 1293:
! 1294: free(matrice_f77);
! 1295: return(-1);
! 1296: }
! 1297: }
! 1298:
! 1299: return(rang);
! 1300: }
! 1301:
! 1302:
! 1303: void
! 1304: determinant(struct_processus *s_etat_processus, struct_matrice *s_matrice,
! 1305: void *valeur)
! 1306: {
! 1307: complex16 *vecteur_complexe;
! 1308:
! 1309: integer4 dimension_vecteur_pivot;
! 1310: integer4 nombre_colonnes_a;
! 1311: integer4 nombre_lignes_a;
! 1312: integer4 *pivot;
! 1313: integer4 rang;
! 1314:
! 1315: integer8 signe;
! 1316:
! 1317: real8 *vecteur_reel;
! 1318:
! 1319: unsigned long i;
! 1320: unsigned long j;
! 1321: unsigned long k;
! 1322: unsigned long taille_matrice_f77;
! 1323:
! 1324: void *matrice_f77;
! 1325:
! 1326: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
! 1327: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
! 1328: dimension_vecteur_pivot = (nombre_lignes_a < nombre_colonnes_a)
! 1329: ? nombre_lignes_a : nombre_colonnes_a;
! 1330: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;
! 1331:
! 1332: switch((*s_matrice).type)
! 1333: {
! 1334: case 'I' :
! 1335: {
! 1336: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 1337: sizeof(real8))) == NULL)
! 1338: {
! 1339: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1340: return;
! 1341: }
! 1342:
! 1343: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1344: {
! 1345: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 1346: {
! 1347: ((real8 *) matrice_f77)[k++] = ((integer8 **)
! 1348: (*s_matrice).tableau)[j][i];
! 1349: }
! 1350: }
! 1351:
! 1352: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot *
! 1353: sizeof(integer4))) == NULL)
! 1354: {
! 1355: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1356: return;
! 1357: }
! 1358:
! 1359: if ((rang = calcul_rang(s_etat_processus, matrice_f77,
! 1360: nombre_lignes_a, nombre_colonnes_a, pivot,
! 1361: dimension_vecteur_pivot, 'R')) < 0)
! 1362: {
! 1363: return;
! 1364: }
! 1365:
! 1366: if (rang < nombre_lignes_a)
! 1367: {
! 1368: (*((real8 *) valeur)) = 0;
! 1369: }
! 1370: else
! 1371: {
! 1372: if ((vecteur_reel = malloc((*s_matrice).nombre_colonnes *
! 1373: sizeof(real8))) == NULL)
! 1374: {
! 1375: (*s_etat_processus).erreur_systeme =
! 1376: d_es_allocation_memoire;
! 1377: return;
! 1378: }
! 1379:
! 1380: signe = 1;
! 1381:
! 1382: for(i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1383: {
! 1384: if ((unsigned long) pivot[i] != (i + 1))
! 1385: {
! 1386: signe = (signe == 1) ? -1 : 1;
! 1387: }
! 1388:
! 1389: vecteur_reel[i] = ((real8 *) matrice_f77)
! 1390: [(i * nombre_colonnes_a) + i];
! 1391: }
! 1392:
! 1393: for(i = 1; i < (*s_matrice).nombre_colonnes; i++)
! 1394: {
! 1395: vecteur_reel[0] *= vecteur_reel[i];
! 1396: }
! 1397:
! 1398: (*((real8 *) valeur)) = vecteur_reel[0] * signe;
! 1399:
! 1400: free(vecteur_reel);
! 1401: }
! 1402:
! 1403: free(matrice_f77);
! 1404: free(pivot);
! 1405:
! 1406: break;
! 1407: }
! 1408:
! 1409: case 'R' :
! 1410: {
! 1411: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 1412: sizeof(real8))) == NULL)
! 1413: {
! 1414: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1415: return;
! 1416: }
! 1417:
! 1418: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1419: {
! 1420: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 1421: {
! 1422: ((real8 *) matrice_f77)[k++] = ((real8 **)
! 1423: (*s_matrice).tableau)[j][i];
! 1424: }
! 1425: }
! 1426:
! 1427: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot *
! 1428: sizeof(integer4))) == NULL)
! 1429: {
! 1430: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1431: return;
! 1432: }
! 1433:
! 1434: if ((rang = calcul_rang(s_etat_processus, matrice_f77,
! 1435: nombre_lignes_a, nombre_colonnes_a, pivot,
! 1436: dimension_vecteur_pivot, 'R')) < 0)
! 1437: {
! 1438: return;
! 1439: }
! 1440:
! 1441: if (rang < nombre_lignes_a)
! 1442: {
! 1443: (*((real8 *) valeur)) = 0;
! 1444: }
! 1445: else
! 1446: {
! 1447: if ((vecteur_reel = malloc((*s_matrice).nombre_colonnes *
! 1448: sizeof(real8))) == NULL)
! 1449: {
! 1450: (*s_etat_processus).erreur_systeme =
! 1451: d_es_allocation_memoire;
! 1452: return;
! 1453: }
! 1454:
! 1455: signe = 1;
! 1456:
! 1457: for(i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1458: {
! 1459: if ((unsigned long) pivot[i] != (i + 1))
! 1460: {
! 1461: signe = (signe == 1) ? -1 : 1;
! 1462: }
! 1463:
! 1464: vecteur_reel[i] = ((real8 *) matrice_f77)
! 1465: [(i * nombre_colonnes_a) + i];
! 1466: }
! 1467:
! 1468: for(i = 1; i < (*s_matrice).nombre_colonnes; i++)
! 1469: {
! 1470: vecteur_reel[0] *= vecteur_reel[i];
! 1471: }
! 1472:
! 1473: (*((real8 *) valeur)) = vecteur_reel[0] * signe;
! 1474:
! 1475: free(vecteur_reel);
! 1476: }
! 1477:
! 1478: free(matrice_f77);
! 1479: free(pivot);
! 1480:
! 1481: break;
! 1482: }
! 1483:
! 1484: case 'C' :
! 1485: {
! 1486: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 1487: sizeof(complex16))) == NULL)
! 1488: {
! 1489: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1490: return;
! 1491: }
! 1492:
! 1493: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1494: {
! 1495: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 1496: {
! 1497: ((complex16 *) matrice_f77)[k++] = ((complex16 **)
! 1498: (*s_matrice).tableau)[j][i];
! 1499: }
! 1500: }
! 1501:
! 1502: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot *
! 1503: sizeof(integer4))) == NULL)
! 1504: {
! 1505: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1506: return;
! 1507: }
! 1508:
! 1509: if ((rang = calcul_rang(s_etat_processus, matrice_f77,
! 1510: nombre_lignes_a, nombre_colonnes_a, pivot,
! 1511: dimension_vecteur_pivot, 'C')) < 0)
! 1512: {
! 1513: return;
! 1514: }
! 1515:
! 1516: if (rang < nombre_lignes_a)
! 1517: {
! 1518: (*((complex16 *) valeur)).partie_reelle = 0;
! 1519: (*((complex16 *) valeur)).partie_imaginaire = 0;
! 1520: }
! 1521: else
! 1522: {
! 1523: if ((vecteur_complexe = malloc((*s_matrice).nombre_colonnes *
! 1524: sizeof(complex16))) == NULL)
! 1525: {
! 1526: (*s_etat_processus).erreur_systeme =
! 1527: d_es_allocation_memoire;
! 1528: return;
! 1529: }
! 1530:
! 1531: signe = 1;
! 1532:
! 1533: for(i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1534: {
! 1535: if ((unsigned long) pivot[i] != (i + 1))
! 1536: {
! 1537: signe = (signe == 1) ? -1 : 1;
! 1538: }
! 1539:
! 1540: vecteur_complexe[i] = ((complex16 *) matrice_f77)
! 1541: [(i * nombre_colonnes_a) + i];
! 1542: }
! 1543:
! 1544: for(i = 1; i < (*s_matrice).nombre_colonnes; i++)
! 1545: {
! 1546: f77multiplicationcc_(&(vecteur_complexe[0]),
! 1547: &(vecteur_complexe[i]), &(vecteur_complexe[0]));
! 1548: }
! 1549:
! 1550: f77multiplicationci_(&(vecteur_complexe[0]), &signe,
! 1551: ((complex16 *) valeur));
! 1552:
! 1553: free(vecteur_complexe);
! 1554: }
! 1555:
! 1556: free(matrice_f77);
! 1557: free(pivot);
! 1558:
! 1559: break;
! 1560: }
! 1561: }
! 1562:
! 1563: return;
! 1564: }
! 1565:
! 1566:
! 1567: void
! 1568: rang(struct_processus *s_etat_processus, struct_matrice *s_matrice,
! 1569: integer8 *valeur)
! 1570: {
! 1571: integer4 dimension_vecteur_pivot;
! 1572: integer4 nombre_lignes_a;
! 1573: integer4 nombre_colonnes_a;
! 1574: integer4 *pivot;
! 1575: integer4 rang;
! 1576: integer4 taille_matrice_f77;
! 1577:
! 1578: unsigned long i;
! 1579: unsigned long j;
! 1580: unsigned long k;
! 1581:
! 1582: void *matrice_f77;
! 1583:
! 1584: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
! 1585: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
! 1586: dimension_vecteur_pivot = (nombre_lignes_a < nombre_colonnes_a)
! 1587: ? nombre_lignes_a : nombre_colonnes_a;
! 1588: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;
! 1589:
! 1590: if ((pivot = (integer4 *) malloc(dimension_vecteur_pivot *
! 1591: sizeof(integer4))) == NULL)
! 1592: {
! 1593: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1594: return;
! 1595: }
! 1596:
! 1597: switch((*s_matrice).type)
! 1598: {
! 1599: case 'I' :
! 1600: {
! 1601: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 1602: sizeof(real8))) == NULL)
! 1603: {
! 1604: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1605: return;
! 1606: }
! 1607:
! 1608: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1609: {
! 1610: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 1611: {
! 1612: ((real8 *) matrice_f77)[k++] = ((integer8 **)
! 1613: (*s_matrice).tableau)[j][i];
! 1614: }
! 1615: }
! 1616:
! 1617: if ((rang = calcul_rang(s_etat_processus, matrice_f77,
! 1618: nombre_lignes_a, nombre_colonnes_a, pivot,
! 1619: dimension_vecteur_pivot, 'R')) < 0)
! 1620: {
! 1621: free(pivot);
! 1622: free(matrice_f77);
! 1623: return;
! 1624: }
! 1625:
! 1626: free(matrice_f77);
! 1627: (*valeur) = rang;
! 1628: break;
! 1629: }
! 1630:
! 1631: case 'R' :
! 1632: {
! 1633: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 1634: sizeof(real8))) == NULL)
! 1635: {
! 1636: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1637: return;
! 1638: }
! 1639:
! 1640: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1641: {
! 1642: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 1643: {
! 1644: ((real8 *) matrice_f77)[k++] = ((real8 **)
! 1645: (*s_matrice).tableau)[j][i];
! 1646: }
! 1647: }
! 1648:
! 1649: if ((rang = calcul_rang(s_etat_processus, matrice_f77,
! 1650: nombre_lignes_a, nombre_colonnes_a, pivot,
! 1651: dimension_vecteur_pivot, 'R')) < 0)
! 1652: {
! 1653: free(pivot);
! 1654: free(matrice_f77);
! 1655: return;
! 1656: }
! 1657:
! 1658: free(matrice_f77);
! 1659: (*valeur) = rang;
! 1660: break;
! 1661: }
! 1662:
! 1663: case 'C' :
! 1664: {
! 1665: if ((matrice_f77 = (void *) malloc(taille_matrice_f77 *
! 1666: sizeof(complex16))) == NULL)
! 1667: {
! 1668: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
! 1669: return;
! 1670: }
! 1671:
! 1672: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
! 1673: {
! 1674: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
! 1675: {
! 1676: ((complex16 *) matrice_f77)[k++] = ((complex16 **)
! 1677: (*s_matrice).tableau)[j][i];
! 1678: }
! 1679: }
! 1680:
! 1681: if ((rang = calcul_rang(s_etat_processus, matrice_f77,
! 1682: nombre_lignes_a, nombre_colonnes_a, pivot,
! 1683: dimension_vecteur_pivot, 'C')) < 0)
! 1684: {
! 1685: free(pivot);
! 1686: free(matrice_f77);
! 1687: return;
! 1688: }
! 1689:
! 1690: free(matrice_f77);
! 1691: (*valeur) = rang;
! 1692: break;
! 1693: }
! 1694: }
! 1695:
! 1696: free(pivot);
! 1697:
! 1698: return;
! 1699: }
! 1700:
! 1701: // vim: ts=4
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to algebre_lineaire3.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / src