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