Return to algebre_lineaire2.c CVS log

Up to [local] / rpl / src

Correction d'une réinitialisation sauvage de la pile des variables par niveau
dans la copie de la structure de description du processus. Cela corrige
la fonction SPAWN qui échouait sur un segmentation fault car la pile des
variables par niveau était vide alors même que l'arbre des variables contenait
bien les variables. Passage à la prerelease 2.

    1: /*
    2: ================================================================================
    3:   RPL/2 (R) version 4.1.0.prerelease.2
    4:   Copyright (C) 1989-2011 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