![[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|
Return to algebre_lineaire3.c CVS log | Up to [local] / rpl / src |
1.13 bertrand 1: /*
2: ================================================================================
1.60 ! bertrand 3: RPL/2 (R) version 4.1.27
1.59 bertrand 4: Copyright (C) 1989-2017 Dr. BERTRAND Joël
1.13 bertrand 5:
6: This file is part of RPL/2.
7:
8: RPL/2 is free software; you can redistribute it and/or modify it
9: under the terms of the CeCILL V2 License as published by the french
10: CEA, CNRS and INRIA.
11:
12: RPL/2 is distributed in the hope that it will be useful, but WITHOUT
13: ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14: FITNESS FOR A PARTICULAR PURPOSE. See the CeCILL V2 License
15: for more details.
16:
17: You should have received a copy of the CeCILL License
18: along with RPL/2. If not, write to info@cecill.info.
19: ================================================================================
20: */
21:
22:
1.14 bertrand 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.
1.15 bertrand 33: La matrice en entrée est écrasée. La matrice de sortie est
1.14 bertrand 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:
1.43 bertrand 62: integer8 i;
63: integer8 j;
64: integer8 k;
65: integer8 taille_matrice_f77;
1.14 bertrand 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:
1.43 bertrand 85: for(i = 0; i < nombre_lignes_a; i++)
1.14 bertrand 86: {
87: tampon = (*s_matrice).tableau[i];
88:
89: if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *)
1.43 bertrand 90: malloc(((size_t) nombre_colonnes_a) * sizeof(real8)))
91: == NULL)
1.14 bertrand 92: {
93: (*s_etat_processus).erreur_systeme =
94: d_es_allocation_memoire;
95: return;
96: }
97:
1.43 bertrand 98: for(j = 0; j < nombre_colonnes_a; j++)
1.14 bertrand 99: {
1.43 bertrand 100: ((real8 **) (*s_matrice).tableau)[i][j] = (real8)
1.14 bertrand 101: ((integer8 *) tampon)[j];
102: }
103:
104: free(tampon);
105: }
106:
107: (*s_matrice).type = 'R';
108: }
109:
110: case 'R' :
111: {
1.43 bertrand 112: if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 113: sizeof(real8))) == NULL)
114: {
115: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
116: return;
117: }
118:
1.43 bertrand 119: if ((matrice_vs_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 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:
1.43 bertrand 135: if ((wr = (real8 *) malloc(((size_t) nombre_lignes_a) *
136: sizeof(real8))) == NULL)
1.14 bertrand 137: {
138: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
139: return;
140: }
141:
1.43 bertrand 142: if ((wi = (real8 *) malloc(((size_t) nombre_lignes_a) *
143: sizeof(real8))) == NULL)
1.14 bertrand 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:
1.43 bertrand 163: lwork = (integer4) ((real8 *) work)[0];
1.14 bertrand 164: free(work);
165:
1.43 bertrand 166: if ((work = (real8 *) malloc(((size_t) lwork) * sizeof(real8)))
167: == NULL)
1.14 bertrand 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:
1.43 bertrand 213: if (((**s_schur).tableau = malloc(((size_t) (**s_schur)
214: .nombre_lignes) * sizeof(real8 *))) == NULL)
1.14 bertrand 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 *)
1.43 bertrand 223: malloc(((size_t) (**s_schur).nombre_colonnes) *
1.14 bertrand 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: {
1.43 bertrand 249: if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 250: sizeof(complex16))) == NULL)
251: {
252: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
253: return;
254: }
255:
1.43 bertrand 256: if ((matrice_vs_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 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:
1.43 bertrand 276: if ((w = (complex16 *) malloc(((size_t) nombre_lignes_a) *
277: sizeof(complex16))) == NULL)
1.14 bertrand 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:
1.43 bertrand 291: if ((rwork = (real8 *) malloc(((size_t) nombre_lignes_a) *
292: sizeof(real8))) == NULL)
1.14 bertrand 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:
1.43 bertrand 304: lwork = (integer4) ((complex16 *) work)[0].partie_reelle;
1.14 bertrand 305: free(work);
306:
1.43 bertrand 307: if ((work = (complex16 *) malloc(((size_t) lwork) *
308: sizeof(complex16))) == NULL)
1.14 bertrand 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:
1.43 bertrand 358: if (((**s_schur).tableau = malloc(((size_t) (**s_schur)
359: .nombre_lignes) * sizeof(complex16 *))) == NULL)
1.14 bertrand 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 *)
1.43 bertrand 368: malloc(((size_t) (**s_schur).nombre_colonnes) *
1.14 bertrand 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:
1.43 bertrand 421: integer8 i;
422: integer8 j;
423: integer8 k;
424: integer8 taille_matrice_f77;
1.14 bertrand 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:
1.43 bertrand 440: for(i = 0; i < nombre_lignes_a; i++)
1.14 bertrand 441: {
442: tampon = (*s_matrice).tableau[i];
443:
444: if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *)
1.43 bertrand 445: malloc(((size_t) nombre_colonnes_a) * sizeof(real8)))
446: == NULL)
1.14 bertrand 447: {
448: (*s_etat_processus).erreur_systeme =
449: d_es_allocation_memoire;
450: return;
451: }
452:
1.43 bertrand 453: for(j = 0; j < nombre_colonnes_a; j++)
1.14 bertrand 454: {
1.43 bertrand 455: ((real8 **) (*s_matrice).tableau)[i][j] = (real8)
1.14 bertrand 456: ((integer8 *) tampon)[j];
457: }
458:
459: free(tampon);
460: }
461:
462: (*s_matrice).type = 'R';
463: }
464:
465: case 'R' :
466: {
1.43 bertrand 467: if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 468: sizeof(real8))) == NULL)
469: {
470: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
471: return;
472: }
473:
1.43 bertrand 474: if (((*tau) = malloc(((size_t) ((nombre_colonnes_a <
475: nombre_lignes_a) ? nombre_colonnes_a : nombre_lignes_a)) *
476: sizeof(real8))) == NULL)
1.14 bertrand 477: {
478: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
479: return;
480: }
481:
1.43 bertrand 482: if ((work = malloc(((size_t) nombre_lignes_a) * sizeof(real8)))
483: == NULL)
1.14 bertrand 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:
1.43 bertrand 512: for(k = 0, i = 0; i < nombre_colonnes_a; i++)
1.14 bertrand 513: {
1.43 bertrand 514: for(j = 0; j < nombre_lignes_a; j++)
1.14 bertrand 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: {
1.43 bertrand 529: if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 530: sizeof(complex16))) == NULL)
531: {
532: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
533: return;
534: }
535:
1.43 bertrand 536: if (((*tau) = malloc(((size_t) ((nombre_colonnes_a <
537: nombre_lignes_a) ? nombre_colonnes_a : nombre_lignes_a)) *
1.14 bertrand 538: sizeof(complex16))) == NULL)
539: {
540: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
541: return;
542: }
543:
1.43 bertrand 544: if ((work = malloc(((size_t) nombre_lignes_a) * sizeof(complex16)))
545: == NULL)
1.14 bertrand 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:
1.43 bertrand 578: for(k = 0, i = 0; i < nombre_colonnes_a; i++)
1.14 bertrand 579: {
1.43 bertrand 580: for(j = 0; j < nombre_lignes_a; j++)
1.14 bertrand 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:
1.43 bertrand 626: integer8 i;
627: integer8 j;
628: integer8 k;
629: integer8 taille_matrice_f77;
1.14 bertrand 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:
1.43 bertrand 646: for(i = 0; i < nombre_lignes_a; i++)
1.14 bertrand 647: {
648: tampon = (*s_matrice).tableau[i];
649:
650: if ((((real8 **) (*s_matrice).tableau)[i] = (real8 *)
1.43 bertrand 651: malloc(((size_t) nombre_colonnes_a) * sizeof(real8)))
652: == NULL)
1.14 bertrand 653: {
654: (*s_etat_processus).erreur_systeme =
655: d_es_allocation_memoire;
656: return;
657: }
658:
1.43 bertrand 659: for(j = 0; j < nombre_colonnes_a; j++)
1.14 bertrand 660: {
1.43 bertrand 661: ((real8 **) (*s_matrice).tableau)[i][j] = (real8)
1.14 bertrand 662: ((integer8 *) tampon)[j];
663: }
664:
665: free(tampon);
666: }
667:
668: (*s_matrice).type = 'R';
669: }
670:
671: case 'R' :
672: {
1.43 bertrand 673: if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 674: sizeof(real8))) == NULL)
675: {
676: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
677: return;
678: }
679:
1.43 bertrand 680: if (((*tau) = malloc(((size_t) ((nombre_colonnes_a <
681: nombre_lignes_a) ? nombre_colonnes_a : nombre_lignes_a)) *
682: sizeof(real8))) == NULL)
1.14 bertrand 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:
1.43 bertrand 697: if ((pivot = malloc(((size_t) nombre_colonnes_a) *
698: sizeof(integer4))) == NULL)
1.14 bertrand 699: {
700: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
701: return;
702: }
703:
1.43 bertrand 704: for(i = 0; i < nombre_colonnes_a; pivot[i++] = 0);
1.14 bertrand 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:
1.43 bertrand 720: lwork = (integer4) ((real8 *) work)[0];
1.14 bertrand 721:
722: free(work);
723:
1.43 bertrand 724: if ((work = malloc(((size_t) lwork) * sizeof(real8))) == NULL)
1.14 bertrand 725: {
726: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
727: return;
728: }
729:
730: // Calcul de la permutation
731:
1.43 bertrand 732: if ((registre = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 733: sizeof(real8))) == NULL)
734: {
735: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
736: return;
737: }
738:
1.43 bertrand 739: memcpy(registre, matrice_a_f77, ((size_t) taille_matrice_f77) *
740: sizeof(real8));
1.14 bertrand 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:
1.43 bertrand 778: for(i = 0; i < nombre_colonnes_a; i++)
1.14 bertrand 779: {
1.43 bertrand 780: if ((i + 1) != pivot[i])
1.14 bertrand 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:
1.43 bertrand 796: for(k = 0, i = 0; i < nombre_colonnes_a; i++)
1.14 bertrand 797: {
1.43 bertrand 798: for(j = 0; j < nombre_lignes_a; j++)
1.14 bertrand 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: {
1.43 bertrand 813: if ((matrice_a_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 814: sizeof(complex16))) == NULL)
815: {
816: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
817: return;
818: }
819:
1.43 bertrand 820: if (((*tau) = malloc(((size_t) ((nombre_colonnes_a <
821: nombre_lignes_a) ? nombre_colonnes_a : nombre_lignes_a)) *
822: sizeof(complex16))) == NULL)
1.14 bertrand 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:
1.43 bertrand 841: if ((pivot = malloc(((size_t) nombre_colonnes_a) *
842: sizeof(integer4))) == NULL)
1.14 bertrand 843: {
844: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
845: return;
846: }
847:
1.43 bertrand 848: if ((rwork = malloc(2 * ((size_t) nombre_colonnes_a) *
849: sizeof(real8))) == NULL)
1.14 bertrand 850: {
851: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
852: return;
853: }
854:
1.43 bertrand 855: for(i = 0; i < nombre_colonnes_a; pivot[i++] = 0);
1.14 bertrand 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:
1.43 bertrand 886: lwork = (integer4) ((complex16 *) work)[0].partie_reelle;
1.14 bertrand 887:
888: free(work);
889:
1.43 bertrand 890: if ((work = malloc(((size_t) lwork) * sizeof(complex16))) == NULL)
1.14 bertrand 891: {
892: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
893: return;
894: }
895:
896: // Calcul de la permutation
897:
1.43 bertrand 898: if ((registre = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 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,
1.43 bertrand 906: ((size_t) taille_matrice_f77) * sizeof(complex16));
1.14 bertrand 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:
1.43 bertrand 952: for(i = 0; i < nombre_colonnes_a; i++)
1.14 bertrand 953: {
1.43 bertrand 954: if ((i + 1) != pivot[i])
1.14 bertrand 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:
1.43 bertrand 970: for(k = 0, i = 0; i < nombre_colonnes_a; i++)
1.14 bertrand 971: {
1.43 bertrand 972: for(j = 0; j < nombre_lignes_a; j++)
1.14 bertrand 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:
1.43 bertrand 1027: integer8 i;
1.14 bertrand 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:
1.43 bertrand 1037: if ((work = malloc(((size_t) nombre_lignes_a) * sizeof(real8))) == NULL)
1.14 bertrand 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:
1.43 bertrand 1058: if ((matrice_c = malloc(((size_t) (nombre_lignes_a * nombre_colonnes_a))
1059: * sizeof(real8))) == NULL)
1.14 bertrand 1060: {
1061: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1062: return(-1);
1063: }
1064:
1.43 bertrand 1065: memcpy(matrice_c, matrice_f77, ((size_t) (nombre_lignes_a *
1066: nombre_colonnes_a)) * sizeof(real8));
1.14 bertrand 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:
1.43 bertrand 1080: if ((iwork = malloc(((size_t) nombre_colonnes_a) * sizeof(integer4)))
1081: == NULL)
1.14 bertrand 1082: {
1083: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1084: return(-1);
1085: }
1086:
1.43 bertrand 1087: if ((work = malloc(4 * ((size_t) nombre_colonnes_a) * sizeof(real8)))
1088: == NULL)
1.14 bertrand 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:
1.43 bertrand 1113: if ((jpvt = malloc(((size_t) nombre_colonnes_a) * sizeof(integer4)))
1114: == NULL)
1.14 bertrand 1115: {
1116: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1117: return(-1);
1118: }
1119:
1.43 bertrand 1120: for(i = 0; i < nombre_colonnes_a; i++)
1.14 bertrand 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:
1.43 bertrand 1137: if ((matrice_b = malloc(((size_t) nombre_lignes_b) * sizeof(real8)))
1138: == NULL)
1.14 bertrand 1139: {
1140: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1141: return(-1);
1142: }
1143:
1.43 bertrand 1144: for(i = 0; i < nombre_lignes_b; i++)
1.14 bertrand 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:
1.43 bertrand 1154: lwork = (integer4) ((real8 *) work)[0];
1.14 bertrand 1155: free(work);
1156:
1.43 bertrand 1157: if ((work = malloc(((size_t) lwork) * sizeof(real8))) == NULL)
1.14 bertrand 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:
1.43 bertrand 1191: if ((matrice_c = malloc(((size_t) (nombre_lignes_a * nombre_colonnes_a))
1192: * sizeof(complex16))) == NULL)
1.14 bertrand 1193: {
1194: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1195: return(-1);
1196: }
1197:
1.43 bertrand 1198: memcpy(matrice_c, matrice_f77, ((size_t) (nombre_lignes_a *
1199: nombre_colonnes_a)) * sizeof(complex16));
1.14 bertrand 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:
1.43 bertrand 1213: if ((rwork = malloc(2 * ((size_t) nombre_colonnes_a) * sizeof(real8)))
1214: == NULL)
1.14 bertrand 1215: {
1216: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1217: return(-1);
1218: }
1219:
1.43 bertrand 1220: if ((work = malloc(2 * ((size_t) nombre_colonnes_a) *
1221: sizeof(complex16))) == NULL)
1.14 bertrand 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:
1.43 bertrand 1245: if ((jpvt = malloc(((size_t) nombre_colonnes_a) * sizeof(integer4)))
1246: == NULL)
1.14 bertrand 1247: {
1248: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1249: return(-1);
1250: }
1251:
1.43 bertrand 1252: for(i = 0; i < nombre_colonnes_a; i++)
1.14 bertrand 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:
1.43 bertrand 1269: if ((matrice_b = malloc(((size_t) nombre_lignes_b) *
1270: sizeof(complex16))) == NULL)
1.14 bertrand 1271: {
1272: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1273: return(-1);
1274: }
1275:
1.43 bertrand 1276: for(i = 0; i < nombre_lignes_b; i++)
1.14 bertrand 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:
1.43 bertrand 1287: lwork = (integer4) ((complex16 *) work)[0].partie_reelle;
1.14 bertrand 1288: free(work);
1289:
1.43 bertrand 1290: if ((work = malloc(((size_t) lwork) * sizeof(complex16))) == NULL)
1.14 bertrand 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:
1.43 bertrand 1333: integer8 i;
1334: integer8 j;
1335: integer8 k;
1.14 bertrand 1336: integer8 signe;
1.43 bertrand 1337: integer8 taille_matrice_f77;
1.14 bertrand 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: {
1.43 bertrand 1353: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 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: {
1.43 bertrand 1364: ((real8 *) matrice_f77)[k++] = (real8) ((integer8 **)
1.14 bertrand 1365: (*s_matrice).tableau)[j][i];
1366: }
1367: }
1368:
1.43 bertrand 1369: if ((pivot = (integer4 *) malloc(((size_t) dimension_vecteur_pivot)
1370: * sizeof(integer4))) == NULL)
1.14 bertrand 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: {
1.43 bertrand 1389: if ((vecteur_reel = malloc(((size_t) ((*s_matrice)
1390: .nombre_colonnes)) * sizeof(real8))) == NULL)
1.14 bertrand 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: {
1.43 bertrand 1401: if (pivot[i] != (i + 1))
1.14 bertrand 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:
1.43 bertrand 1415: (*((real8 *) valeur)) = vecteur_reel[0] * ((real8) signe);
1.14 bertrand 1416: free(vecteur_reel);
1417: }
1418:
1419: free(matrice_f77);
1420: free(pivot);
1421:
1422: break;
1423: }
1424:
1425: case 'R' :
1426: {
1.43 bertrand 1427: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 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:
1.43 bertrand 1443: if ((pivot = (integer4 *) malloc(((size_t) dimension_vecteur_pivot)
1444: * sizeof(integer4))) == NULL)
1.14 bertrand 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: {
1.43 bertrand 1463: if ((vecteur_reel = malloc(((size_t) (*s_matrice)
1464: .nombre_colonnes) * sizeof(real8))) == NULL)
1.14 bertrand 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: {
1.43 bertrand 1475: if (pivot[i] != (i + 1))
1.14 bertrand 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:
1.43 bertrand 1489: (*((real8 *) valeur)) = vecteur_reel[0] * ((real8) signe);
1.14 bertrand 1490:
1491: free(vecteur_reel);
1492: }
1493:
1494: free(matrice_f77);
1495: free(pivot);
1496:
1497: break;
1498: }
1499:
1500: case 'C' :
1501: {
1.43 bertrand 1502: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 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:
1.43 bertrand 1518: if ((pivot = (integer4 *) malloc(((size_t) dimension_vecteur_pivot)
1519: * sizeof(integer4))) == NULL)
1.14 bertrand 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: {
1.43 bertrand 1539: if ((vecteur_complexe = malloc(((size_t) (*s_matrice)
1540: .nombre_colonnes) * sizeof(complex16))) == NULL)
1.14 bertrand 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: {
1.43 bertrand 1551: if (pivot[i] != (i + 1))
1.14 bertrand 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:
1.43 bertrand 1594: integer8 i;
1595: integer8 j;
1596: integer8 k;
1.14 bertrand 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:
1.43 bertrand 1606: if ((pivot = (integer4 *) malloc(((size_t) dimension_vecteur_pivot) *
1.14 bertrand 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: {
1.43 bertrand 1617: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 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: {
1.43 bertrand 1628: ((real8 *) matrice_f77)[k++] = (real8) ((integer8 **)
1.14 bertrand 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: {
1.43 bertrand 1649: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 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: {
1.43 bertrand 1681: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
1.14 bertrand 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