Annotation of rpl/src/algebre_lineaire2.c, revision 1.21

1.1       bertrand    1: /*
                      2: ================================================================================
1.20      bertrand    3:   RPL/2 (R) version 4.1.0.prerelease.1
1.16      bertrand    4:   Copyright (C) 1989-2011 Dr. BERTRAND Joël
1.1       bertrand    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: 
1.11      bertrand   23: #include "rpl-conv.h"
1.1       bertrand   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>