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