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