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