1: /*
2: ================================================================================
3: RPL/2 (R) version 4.1.36
4: Copyright (C) 1989-2025 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 l'inversion d'une matrice carrée
29: ================================================================================
30: Entrées : struct_matrice
31: --------------------------------------------------------------------------------
32: Sorties : inverse de la matrice d'entrée et drapeau d'erreur. La matrice
33: en entrée est écrasée. La matrice de sortie est réelle si
34: la matrice d'entrée est entière ou réelle, et complexe
35: si la matrice d'entrée est complexe.
36: --------------------------------------------------------------------------------
37: Effets de bord : néant
38: ================================================================================
39: */
40:
41: void
42: inversion_matrice(struct_processus *s_etat_processus,
43: struct_matrice *s_matrice)
44: {
45: real8 *work;
46:
47: integer4 dim_matrice;
48: integer4 dim_work;
49: integer4 erreur;
50: integer4 *pivot;
51:
52: integer8 i;
53: integer8 j;
54: integer8 k;
55: integer8 rang_estime;
56: integer8 taille_matrice_f77;
57:
58: struct_complexe16 *c_work;
59:
60: void *matrice_f77;
61:
62: rang(s_etat_processus, s_matrice, &rang_estime);
63:
64: if ((*s_etat_processus).erreur_systeme != d_es)
65: {
66: return;
67: }
68:
69: if (((*s_etat_processus).erreur_execution != d_ex) ||
70: ((*s_etat_processus).exception != d_ep))
71: {
72: return;
73: }
74:
75: if (rang_estime < (integer8) (*s_matrice).nombre_lignes)
76: {
77: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
78: return;
79: }
80:
81: taille_matrice_f77 = (*s_matrice).nombre_lignes *
82: (*s_matrice).nombre_colonnes;
83: dim_matrice = (integer4) (*s_matrice).nombre_lignes;
84:
85: switch((*s_matrice).type)
86: {
87: case 'I' :
88: {
89: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
90: sizeof(real8))) == NULL)
91: {
92: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
93: return;
94: }
95:
96: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
97: {
98: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
99: {
100: ((real8 *) matrice_f77)[k++] = (real8) ((integer8 **)
101: (*s_matrice).tableau)[j][i];
102: }
103: }
104:
105: for(i = 0; i < (*s_matrice).nombre_lignes; i++)
106: {
107: free(((integer8 **) (*s_matrice).tableau)[i]);
108: }
109:
110: free((integer8 **) (*s_matrice).tableau);
111:
112: if (((*s_matrice).tableau = malloc(((size_t) (*s_matrice)
113: .nombre_lignes) * sizeof(real8 *))) == NULL)
114: {
115: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
116: return;
117: }
118:
119: for(i = 0; i < (*s_matrice).nombre_lignes; i++)
120: {
121: if ((((*s_matrice).tableau)[i] =
122: (real8 *) malloc(((size_t) (*s_matrice)
123: .nombre_colonnes) * sizeof(real8))) == NULL)
124: {
125: (*s_etat_processus).erreur_systeme =
126: d_es_allocation_memoire;
127: return;
128: }
129: }
130:
131: (*s_matrice).type = 'R';
132:
133: if ((pivot = (integer4 *) malloc(((size_t) (*s_matrice)
134: .nombre_lignes) * sizeof(integer4))) == NULL)
135: {
136: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
137: return;
138: }
139:
140: if ((work = (real8 *) malloc(sizeof(real8))) == NULL)
141: {
142: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
143: return;
144: }
145:
146: dim_work = -1;
147:
148: dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work,
149: &dim_work, &erreur);
150:
151: if (erreur != 0)
152: {
153: if (erreur > 0)
154: {
155: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
156: }
157: else
158: {
159: (*s_etat_processus).erreur_execution =
160: d_ex_routines_mathematiques;
161: }
162:
163: free(pivot);
164: free(work);
165: free(matrice_f77);
166: return;
167: }
168:
169: dim_work = (integer4) work[0];
170:
171: free(work);
172:
173: if ((work = (real8 *) malloc(((size_t) dim_work) * sizeof(real8)))
174: == NULL)
175: {
176: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
177: return;
178: }
179:
180: dgetrf_(&dim_matrice, &dim_matrice, matrice_f77,
181: &dim_matrice, pivot, &erreur);
182:
183: if (erreur != 0)
184: {
185: if (erreur > 0)
186: {
187: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
188: }
189: else
190: {
191: (*s_etat_processus).erreur_execution =
192: d_ex_routines_mathematiques;
193: }
194:
195: free(pivot);
196: free(work);
197: free(matrice_f77);
198: return;
199: }
200:
201: dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work,
202: &dim_work, &erreur);
203:
204: if (erreur != 0)
205: {
206: if (erreur > 0)
207: {
208: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
209: }
210: else
211: {
212: (*s_etat_processus).erreur_execution =
213: d_ex_routines_mathematiques;
214: }
215:
216: free(pivot);
217: free(work);
218: free(matrice_f77);
219: return;
220: }
221:
222: free(work);
223: free(pivot);
224:
225: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
226: {
227: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
228: {
229: ((real8 **) (*s_matrice).tableau)[j][i] =
230: ((real8 *) matrice_f77)[k++];
231: }
232: }
233:
234: free(matrice_f77);
235:
236: break;
237: }
238:
239: case 'R' :
240: {
241: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
242: sizeof(real8))) == NULL)
243: {
244: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
245: return;
246: }
247:
248: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
249: {
250: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
251: {
252: ((real8 *) matrice_f77)[k++] = ((real8 **)
253: (*s_matrice).tableau)[j][i];
254: }
255: }
256:
257: if ((pivot = (integer4 *) malloc(((size_t) (*s_matrice)
258: .nombre_lignes) * sizeof(integer4))) == NULL)
259: {
260: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
261: return;
262: }
263:
264: if ((work = (real8 *) malloc(sizeof(real8))) == NULL)
265: {
266: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
267: return;
268: }
269:
270: dim_work = -1;
271:
272: dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work,
273: &dim_work, &erreur);
274:
275: if (erreur != 0)
276: {
277: if (erreur > 0)
278: {
279: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
280: }
281: else
282: {
283: (*s_etat_processus).erreur_execution =
284: d_ex_routines_mathematiques;
285: }
286:
287: free(pivot);
288: free(work);
289: free(matrice_f77);
290: return;
291: }
292:
293: dim_work = (integer4) work[0];
294:
295: free(work);
296:
297: if ((work = (real8 *) malloc(((size_t) dim_work) * sizeof(real8)))
298: == NULL)
299: {
300: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
301: return;
302: }
303:
304: dgetrf_(&dim_matrice, &dim_matrice, matrice_f77,
305: &dim_matrice, pivot, &erreur);
306:
307: if (erreur != 0)
308: {
309: if (erreur > 0)
310: {
311: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
312: }
313: else
314: {
315: (*s_etat_processus).erreur_execution =
316: d_ex_routines_mathematiques;
317: }
318:
319: free(pivot);
320: free(work);
321: free(matrice_f77);
322: return;
323: }
324:
325: dgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, work,
326: &dim_work, &erreur);
327:
328: if (erreur != 0)
329: {
330: if (erreur > 0)
331: {
332: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
333: }
334: else
335: {
336: (*s_etat_processus).erreur_execution =
337: d_ex_routines_mathematiques;
338: }
339:
340: free(pivot);
341: free(work);
342: free(matrice_f77);
343: return;
344: }
345:
346: free(work);
347: free(pivot);
348:
349: for(k = 0, i = 0; i < (*s_matrice).nombre_lignes; i++)
350: {
351: for(j = 0; j < (*s_matrice).nombre_colonnes; j++)
352: {
353: ((real8 **) (*s_matrice).tableau)[j][i] =
354: ((real8 *) matrice_f77)[k++];
355: }
356: }
357:
358: free(matrice_f77);
359:
360: break;
361: }
362:
363: case 'C' :
364: {
365: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
366: sizeof(struct_complexe16))) == NULL)
367: {
368: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
369: return;
370: }
371:
372: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
373: {
374: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
375: {
376: (((struct_complexe16 *) matrice_f77)[k]).partie_reelle =
377: (((struct_complexe16 **) (*s_matrice).tableau)
378: [j][i]).partie_reelle;
379: (((struct_complexe16 *) matrice_f77)[k]).partie_imaginaire =
380: (((struct_complexe16 **) (*s_matrice).tableau)
381: [j][i]).partie_imaginaire;
382: k++;
383: }
384: }
385:
386: if ((pivot = (integer4 *) malloc(((size_t) (*s_matrice)
387: .nombre_lignes) * sizeof(integer4))) == NULL)
388: {
389: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
390: return;
391: }
392:
393: if ((c_work = (struct_complexe16 *) malloc(
394: sizeof(struct_complexe16))) == NULL)
395: {
396: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
397: return;
398: }
399:
400: dim_work = -1;
401:
402: zgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, c_work,
403: &dim_work, &erreur);
404:
405: if (erreur != 0)
406: {
407: if (erreur > 0)
408: {
409: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
410: }
411: else
412: {
413: (*s_etat_processus).erreur_execution =
414: d_ex_routines_mathematiques;
415: }
416:
417: free(pivot);
418: free(c_work);
419: free(matrice_f77);
420: return;
421: }
422:
423: dim_work = (integer4) c_work[0].partie_reelle;
424:
425: free(c_work);
426:
427: if ((c_work = (struct_complexe16 *) malloc(((size_t) dim_work) *
428: sizeof(struct_complexe16))) == NULL)
429: {
430: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
431: return;
432: }
433:
434: zgetrf_(&dim_matrice, &dim_matrice, matrice_f77,
435: &dim_matrice, pivot, &erreur);
436:
437: if (erreur != 0)
438: {
439: if (erreur > 0)
440: {
441: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
442: }
443: else
444: {
445: (*s_etat_processus).erreur_execution =
446: d_ex_routines_mathematiques;
447: }
448:
449: free(pivot);
450: free(c_work);
451: free(matrice_f77);
452: return;
453: }
454:
455: zgetri_(&dim_matrice, matrice_f77, &dim_matrice, pivot, c_work,
456: &dim_work, &erreur);
457:
458: if (erreur != 0)
459: {
460: if (erreur > 0)
461: {
462: (*s_etat_processus).exception = d_ep_matrice_non_inversible;
463: }
464: else
465: {
466: (*s_etat_processus).erreur_execution =
467: d_ex_routines_mathematiques;
468: }
469:
470: free(pivot);
471: free(c_work);
472: free(matrice_f77);
473: return;
474: }
475:
476: free(c_work);
477: free(pivot);
478:
479: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
480: {
481: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
482: {
483: (((struct_complexe16 **) (*s_matrice).tableau)
484: [j][i]).partie_reelle = (((struct_complexe16 *)
485: matrice_f77)[k]).partie_reelle;
486: (((struct_complexe16 **) (*s_matrice).tableau)
487: [j][i]).partie_imaginaire = (((struct_complexe16 *)
488: matrice_f77)[k]).partie_imaginaire;
489: k++;
490: }
491: }
492:
493: free(matrice_f77);
494:
495: break;
496: }
497: }
498:
499: return;
500: }
501:
502:
503: /*
504: ================================================================================
505: Fonction calculant les vecteurs propres d'une matrice carrée
506: ================================================================================
507: Entrées : struct_matrice, pointeur sur la matrice d'entrée, ainsi que
508: struct_matrice, pointeur sur le vecteur des valeurs propres
509: (vecteur de complexes) et sur deux matrices complexes
510: contenant les différents vecteurs propres gauches et droits.
511: Les pointeurs sur les vecteurs propres peuvent être nuls,
512: et seules les valeurs propres sont calculées.
513: L'allocation des tableaux de sortie est effectuée dans la
514: routine, mais les structures s_matrice et s_vecteurs
515: doivent être passées en paramètre.
516: La matrice présente en entrée n'est pas libérée.
517: --------------------------------------------------------------------------------
518: Sorties : vecteur contenant les valeurs propres, matrice contenant
519: les vecteurs propres gauches et matrice contenant les vecteurs
520: propres droits. Toutes les sorties sont complexes.
521: --------------------------------------------------------------------------------
522: Effets de bord : néant
523: ================================================================================
524: */
525:
526: void
527: valeurs_propres(struct_processus *s_etat_processus,
528: struct_matrice *s_matrice,
529: struct_vecteur *s_valeurs_propres,
530: struct_matrice *s_vecteurs_propres_gauches,
531: struct_matrice *s_vecteurs_propres_droits)
532: {
533: real8 *rwork;
534:
535: integer4 dim_matrice;
536: integer4 lwork;
537: integer4 erreur;
538:
539: integer8 i;
540: integer8 j;
541: integer8 k;
542: integer8 taille_matrice_f77;
543:
544: struct_complexe16 *matrice_f77;
545: struct_complexe16 *vpd_f77;
546: struct_complexe16 *vpg_f77;
547: struct_complexe16 *work;
548:
549: unsigned char calcul_vp_droits;
550: unsigned char calcul_vp_gauches;
551: unsigned char negatif;
552:
553: taille_matrice_f77 = (*s_matrice).nombre_lignes *
554: (*s_matrice).nombre_colonnes;
555: dim_matrice = (integer4) (*s_matrice).nombre_lignes;
556:
557: /*
558: * Allocation de la matrice complexe
559: */
560:
561: if ((matrice_f77 = (struct_complexe16 *) malloc(((size_t)
562: taille_matrice_f77) * sizeof(struct_complexe16))) == NULL)
563: {
564: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
565: return;
566: }
567:
568: /*
569: * Garniture de la matrice de compatibilité Fortran
570: */
571:
572: switch((*s_matrice).type)
573: {
574: /*
575: * La matrice en entrée est une matrice d'entiers.
576: */
577:
578: case 'I' :
579: {
580: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
581: {
582: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
583: {
584: ((struct_complexe16 *) matrice_f77)[k]
585: .partie_reelle = (real8) ((integer8 **)
586: (*s_matrice).tableau)[j][i];
587: ((struct_complexe16 *) matrice_f77)[k++]
588: .partie_imaginaire = (real8) 0;
589: }
590: }
591:
592: break;
593: }
594:
595: /*
596: * La matrice en entrée est une matrice de réels.
597: */
598:
599: case 'R' :
600: {
601: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
602: {
603: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
604: {
605: ((struct_complexe16 *) matrice_f77)[k]
606: .partie_reelle = ((real8 **)
607: (*s_matrice).tableau)[j][i];
608: ((struct_complexe16 *) matrice_f77)[k++]
609: .partie_imaginaire = (real8) 0;
610: }
611: }
612:
613: break;
614: }
615:
616: /*
617: * La matrice en entrée est une matrice de complexes.
618: */
619:
620: case 'C' :
621: {
622: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
623: {
624: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
625: {
626: ((struct_complexe16 *) matrice_f77)[k]
627: .partie_reelle = ((struct_complexe16 **)
628: (*s_matrice).tableau)[j][i].partie_reelle;
629: ((struct_complexe16 *) matrice_f77)[k++]
630: .partie_imaginaire = ((struct_complexe16 **)
631: (*s_matrice).tableau)[j][i].partie_imaginaire;
632: }
633: }
634:
635: break;
636: }
637: }
638:
639: (*s_valeurs_propres).type = 'C';
640: (*s_valeurs_propres).taille = (*s_matrice).nombre_lignes;
641:
642: if (((*s_valeurs_propres).tableau = (struct_complexe16 *)
643: malloc(((size_t) (*s_valeurs_propres).taille) *
644: sizeof(struct_complexe16))) == NULL)
645: {
646: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
647: return;
648: }
649:
650: if (s_vecteurs_propres_gauches != NULL)
651: {
652: (*s_vecteurs_propres_gauches).type = 'C';
653: calcul_vp_gauches = 'V';
654:
655: if ((vpg_f77 = (struct_complexe16 *) malloc(((size_t)
656: taille_matrice_f77) * sizeof(struct_complexe16))) == NULL)
657: {
658: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
659: return;
660: }
661: }
662: else
663: {
664: vpg_f77 = NULL;
665: calcul_vp_gauches = 'N';
666: }
667:
668: if (s_vecteurs_propres_droits != NULL)
669: {
670: (*s_vecteurs_propres_droits).type = 'C';
671: calcul_vp_droits = 'V';
672:
673: if ((vpd_f77 = (struct_complexe16 *) malloc(((size_t)
674: taille_matrice_f77) * sizeof(struct_complexe16))) == NULL)
675: {
676: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
677: return;
678: }
679: }
680: else
681: {
682: vpd_f77 = NULL;
683: calcul_vp_droits = 'N';
684: }
685:
686: negatif = 'N';
687: lwork = -1;
688:
689: if ((rwork = (real8 *) malloc(2 * ((size_t) (*s_matrice).nombre_lignes) *
690: sizeof(real8))) == NULL)
691: {
692: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
693: return;
694: }
695:
696: if ((work = (struct_complexe16 *) malloc(sizeof(struct_complexe16)))
697: == NULL)
698: {
699: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
700: return;
701: }
702:
703: zgeev_(&negatif, &negatif, &dim_matrice, matrice_f77, &dim_matrice,
704: (struct_complexe16 *) (*s_valeurs_propres).tableau,
705: vpg_f77, &dim_matrice, vpd_f77, &dim_matrice,
706: work, &lwork, rwork, &erreur, 1, 1);
707:
708: if (erreur != 0)
709: {
710: if (erreur > 0)
711: {
712: (*s_etat_processus).exception = d_ep_decomposition_QR;
713: }
714: else
715: {
716: (*s_etat_processus).erreur_execution =
717: d_ex_routines_mathematiques;
718: }
719:
720: free(work);
721: free(rwork);
722: free(matrice_f77);
723:
724: if (calcul_vp_gauches == 'V')
725: {
726: free(vpg_f77);
727: }
728:
729: if (calcul_vp_droits == 'V')
730: {
731: free(vpd_f77);
732: }
733:
734: return;
735: }
736:
737: lwork = (integer4) work[0].partie_reelle;
738: free(work);
739:
740: if ((work = (struct_complexe16 *) malloc(((size_t) lwork) *
741: sizeof(struct_complexe16))) == NULL)
742: {
743: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
744: return;
745: }
746:
747: zgeev_(&calcul_vp_gauches, &calcul_vp_droits, &dim_matrice,
748: matrice_f77, &dim_matrice,
749: (struct_complexe16 *) (*s_valeurs_propres).tableau,
750: vpg_f77, &dim_matrice, vpd_f77, &dim_matrice,
751: work, &lwork, rwork, &erreur, 1, 1);
752:
753: if (erreur != 0)
754: {
755: if (erreur > 0)
756: {
757: (*s_etat_processus).exception = d_ep_decomposition_QR;
758: }
759: else
760: {
761: (*s_etat_processus).erreur_execution =
762: d_ex_routines_mathematiques;
763: }
764:
765: free(work);
766: free(rwork);
767: free(matrice_f77);
768:
769: if (calcul_vp_gauches == 'V')
770: {
771: free(vpg_f77);
772: }
773:
774: if (calcul_vp_droits == 'V')
775: {
776: free(vpd_f77);
777: }
778:
779: return;
780: }
781:
782: free(work);
783: free(rwork);
784:
785: if (calcul_vp_gauches == 'V')
786: {
787: (*s_vecteurs_propres_gauches).type = 'C';
788: (*s_vecteurs_propres_gauches).nombre_lignes =
789: (*s_matrice).nombre_lignes;
790: (*s_vecteurs_propres_gauches).nombre_colonnes =
791: (*s_matrice).nombre_colonnes;
792:
793: if (((*s_vecteurs_propres_gauches).tableau = malloc(
794: ((size_t) (*s_vecteurs_propres_gauches).nombre_lignes) *
795: sizeof(struct_complexe16 *))) == NULL)
796: {
797: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
798: return;
799: }
800:
801: for(i = 0; i < (*s_vecteurs_propres_gauches).nombre_lignes; i++)
802: {
803: if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches)
804: .tableau)[i] = (struct_complexe16 *) malloc(((size_t)
805: (*s_vecteurs_propres_gauches).nombre_colonnes) *
806: sizeof(struct_complexe16))) == NULL)
807: {
808: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
809: return;
810: }
811: }
812:
813: for(k = 0, i = 0; i < (*s_vecteurs_propres_gauches).nombre_colonnes;
814: i++)
815: {
816: for(j = 0; j < (*s_vecteurs_propres_gauches).nombre_lignes; j++)
817: {
818: ((struct_complexe16 **) (*s_vecteurs_propres_gauches)
819: .tableau)[j][i].partie_reelle =
820: vpg_f77[k].partie_reelle;
821: ((struct_complexe16 **) (*s_vecteurs_propres_gauches)
822: .tableau)[j][i].partie_imaginaire =
823: vpg_f77[k++].partie_imaginaire;
824: }
825: }
826:
827: free(vpg_f77);
828: }
829:
830: if (calcul_vp_droits == 'V')
831: {
832: (*s_vecteurs_propres_droits).type = 'C';
833: (*s_vecteurs_propres_droits).nombre_lignes =
834: (*s_matrice).nombre_lignes;
835: (*s_vecteurs_propres_droits).nombre_colonnes =
836: (*s_matrice).nombre_colonnes;
837:
838: if (((*s_vecteurs_propres_droits).tableau = malloc(
839: ((size_t) (*s_vecteurs_propres_droits).nombre_lignes) *
840: sizeof(struct_complexe16 *))) == NULL)
841: {
842: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
843: return;
844: }
845:
846: for(i = 0; i < (*s_vecteurs_propres_droits).nombre_lignes; i++)
847: {
848: if ((((struct_complexe16 **) (*s_vecteurs_propres_droits)
849: .tableau)[i] = (struct_complexe16 *) malloc(((size_t)
850: (*s_vecteurs_propres_droits).nombre_colonnes) *
851: sizeof(struct_complexe16))) == NULL)
852: {
853: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
854: return;
855: }
856: }
857:
858: for(k = 0, i = 0; i < (*s_vecteurs_propres_droits).nombre_colonnes; i++)
859: {
860: for(j = 0; j < (*s_vecteurs_propres_droits).nombre_lignes; j++)
861: {
862: ((struct_complexe16 **) (*s_vecteurs_propres_droits)
863: .tableau)[j][i].partie_reelle =
864: vpd_f77[k].partie_reelle;
865: ((struct_complexe16 **) (*s_vecteurs_propres_droits)
866: .tableau)[j][i].partie_imaginaire =
867: vpd_f77[k++].partie_imaginaire;
868: }
869: }
870:
871: free(vpd_f77);
872: }
873:
874: free(matrice_f77);
875:
876: return;
877: }
878:
879:
880: /*
881: ================================================================================
882: Fonction calculant les vecteurs propres généralisés d'une matrice carrée
883: dans une métrique.
884: ================================================================================
885: Entrées : struct_matrice, pointeur sur la matrice d'entrée, ainsi que
886: struct_matrice, pointeur sur la métrique et
887: struct_matrice, pointeur sur le vecteur des valeurs propres
888: (vecteur de complexes) et sur deux matrices complexes
889: contenant les différents vecteurs propres gauches et droits.
890: Les pointeurs sur les vecteurs propres peuvent être nuls,
891: et seules les valeurs propres sont calculées.
892: L'allocation des tableaux de sortie est effectuée dans la
893: routine, mais les structures s_matrice et s_vecteurs
894: doivent être passées en paramètre.
895: La matrice présente en entrée n'est pas libérée.
896: --------------------------------------------------------------------------------
897: Sorties : vecteur contenant les valeurs propres, matrice contenant
898: les vecteurs propres gauches et matrice contenant les vecteurs
899: propres droits. Toutes les sorties sont complexes.
900: --------------------------------------------------------------------------------
901: Effets de bord : néant
902: ================================================================================
903: */
904:
905: void
906: valeurs_propres_generalisees(struct_processus *s_etat_processus,
907: struct_matrice *s_matrice,
908: struct_matrice *s_metrique,
909: struct_vecteur *s_valeurs_propres,
910: struct_matrice *s_vecteurs_propres_gauches,
911: struct_matrice *s_vecteurs_propres_droits)
912: {
913: real8 *rwork;
914:
915: integer4 dim_matrice;
916: integer4 lwork;
917: integer4 erreur;
918:
919: integer8 i;
920: integer8 j;
921: integer8 k;
922: integer8 taille_matrice_f77;
923:
924: struct_complexe16 *alpha;
925: struct_complexe16 *beta;
926: struct_complexe16 *matrice_f77;
927: struct_complexe16 *metrique_f77;
928: struct_complexe16 *vpd_f77;
929: struct_complexe16 *vpg_f77;
930: struct_complexe16 *work;
931:
932: unsigned char calcul_vp_droits;
933: unsigned char calcul_vp_gauches;
934: unsigned char negatif;
935:
936: taille_matrice_f77 = (*s_matrice).nombre_lignes *
937: (*s_matrice).nombre_colonnes;
938: dim_matrice = (integer4) (*s_matrice).nombre_lignes;
939:
940: /*
941: * Allocation de la matrice complexe
942: */
943:
944: if ((matrice_f77 = (struct_complexe16 *) malloc(((size_t)
945: taille_matrice_f77) * sizeof(struct_complexe16))) == NULL)
946: {
947: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
948: return;
949: }
950:
951: /*
952: * Garniture de la matrice de compatibilité Fortran
953: */
954:
955: switch((*s_matrice).type)
956: {
957: /*
958: * La matrice en entrée est une matrice d'entiers.
959: */
960:
961: case 'I' :
962: {
963: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
964: {
965: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
966: {
967: ((struct_complexe16 *) matrice_f77)[k]
968: .partie_reelle = (real8) ((integer8 **)
969: (*s_matrice).tableau)[j][i];
970: ((struct_complexe16 *) matrice_f77)[k++]
971: .partie_imaginaire = (real8) 0;
972: }
973: }
974:
975: break;
976: }
977:
978: /*
979: * La matrice en entrée est une matrice de réels.
980: */
981:
982: case 'R' :
983: {
984: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
985: {
986: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
987: {
988: ((struct_complexe16 *) matrice_f77)[k]
989: .partie_reelle = ((real8 **)
990: (*s_matrice).tableau)[j][i];
991: ((struct_complexe16 *) matrice_f77)[k++]
992: .partie_imaginaire = (real8) 0;
993: }
994: }
995:
996: break;
997: }
998:
999: /*
1000: * La matrice en entrée est une matrice de complexes.
1001: */
1002:
1003: case 'C' :
1004: {
1005: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
1006: {
1007: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
1008: {
1009: ((struct_complexe16 *) matrice_f77)[k]
1010: .partie_reelle = ((struct_complexe16 **)
1011: (*s_matrice).tableau)[j][i].partie_reelle;
1012: ((struct_complexe16 *) matrice_f77)[k++]
1013: .partie_imaginaire = ((struct_complexe16 **)
1014: (*s_matrice).tableau)[j][i].partie_imaginaire;
1015: }
1016: }
1017:
1018: break;
1019: }
1020: }
1021:
1022: /*
1023: * Allocation de la metrique complexe
1024: */
1025:
1026: if ((metrique_f77 = (struct_complexe16 *) malloc(((size_t)
1027: taille_matrice_f77) * sizeof(struct_complexe16))) == NULL)
1028: {
1029: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1030: return;
1031: }
1032:
1033: /*
1034: * Garniture de la matrice de compatibilité Fortran
1035: */
1036:
1037: switch((*s_metrique).type)
1038: {
1039: /*
1040: * La matrice en entrée est une matrice d'entiers.
1041: */
1042:
1043: case 'I' :
1044: {
1045: for(k = 0, i = 0; i < (*s_metrique).nombre_colonnes; i++)
1046: {
1047: for(j = 0; j < (*s_metrique).nombre_lignes; j++)
1048: {
1049: ((struct_complexe16 *) metrique_f77)[k]
1050: .partie_reelle = (real8) ((integer8 **)
1051: (*s_metrique).tableau)[j][i];
1052: ((struct_complexe16 *) metrique_f77)[k++]
1053: .partie_imaginaire = (real8) 0;
1054: }
1055: }
1056:
1057: break;
1058: }
1059:
1060: /*
1061: * La matrice en entrée est une matrice de réels.
1062: */
1063:
1064: case 'R' :
1065: {
1066: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
1067: {
1068: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
1069: {
1070: ((struct_complexe16 *) metrique_f77)[k]
1071: .partie_reelle = ((real8 **)
1072: (*s_metrique).tableau)[j][i];
1073: ((struct_complexe16 *) metrique_f77)[k++]
1074: .partie_imaginaire = (real8) 0;
1075: }
1076: }
1077:
1078: break;
1079: }
1080:
1081: /*
1082: * La matrice en entrée est une matrice de complexes.
1083: */
1084:
1085: case 'C' :
1086: {
1087: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
1088: {
1089: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
1090: {
1091: ((struct_complexe16 *) metrique_f77)[k]
1092: .partie_reelle = ((struct_complexe16 **)
1093: (*s_metrique).tableau)[j][i].partie_reelle;
1094: ((struct_complexe16 *) metrique_f77)[k++]
1095: .partie_imaginaire = ((struct_complexe16 **)
1096: (*s_metrique).tableau)[j][i].partie_imaginaire;
1097: }
1098: }
1099:
1100: break;
1101: }
1102: }
1103:
1104: (*s_valeurs_propres).type = 'C';
1105: (*s_valeurs_propres).taille = (*s_matrice).nombre_lignes;
1106:
1107: if (((*s_valeurs_propres).tableau = (struct_complexe16 *)
1108: malloc(((size_t) (*s_valeurs_propres).taille)
1109: * sizeof(struct_complexe16))) == NULL)
1110: {
1111: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1112: return;
1113: }
1114:
1115: if (s_vecteurs_propres_gauches != NULL)
1116: {
1117: (*s_vecteurs_propres_gauches).type = 'C';
1118: calcul_vp_gauches = 'V';
1119:
1120: if ((vpg_f77 = (struct_complexe16 *) malloc(((size_t)
1121: taille_matrice_f77) * sizeof(struct_complexe16))) == NULL)
1122: {
1123: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1124: return;
1125: }
1126: }
1127: else
1128: {
1129: vpg_f77 = NULL;
1130: calcul_vp_gauches = 'N';
1131: }
1132:
1133: if (s_vecteurs_propres_droits != NULL)
1134: {
1135: (*s_vecteurs_propres_droits).type = 'C';
1136: calcul_vp_droits = 'V';
1137:
1138: if ((vpd_f77 = (struct_complexe16 *) malloc(((size_t)
1139: taille_matrice_f77) * sizeof(struct_complexe16))) == NULL)
1140: {
1141: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1142: return;
1143: }
1144: }
1145: else
1146: {
1147: vpd_f77 = NULL;
1148: calcul_vp_droits = 'N';
1149: }
1150:
1151: negatif = 'N';
1152: lwork = -1;
1153:
1154: if ((rwork = (real8 *) malloc(8 * ((size_t) (*s_matrice).nombre_lignes) *
1155: sizeof(real8))) == NULL)
1156: {
1157: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1158: return;
1159: }
1160:
1161: if ((work = (struct_complexe16 *) malloc(sizeof(struct_complexe16)))
1162: == NULL)
1163: {
1164: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1165: return;
1166: }
1167:
1168: if ((alpha = (struct_complexe16 *) malloc(((size_t) (*s_valeurs_propres)
1169: .taille) * sizeof(struct_complexe16))) == NULL)
1170: {
1171: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1172: return;
1173: }
1174:
1175: if ((beta = (struct_complexe16 *) malloc(((size_t) (*s_valeurs_propres)
1176: .taille) * sizeof(struct_complexe16))) == NULL)
1177: {
1178: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1179: return;
1180: }
1181:
1182: zggev_(&negatif, &negatif, &dim_matrice, matrice_f77, &dim_matrice,
1183: metrique_f77, &dim_matrice, alpha, beta,
1184: vpg_f77, &dim_matrice, vpd_f77, &dim_matrice,
1185: work, &lwork, rwork, &erreur, 1, 1);
1186:
1187: if (erreur != 0)
1188: {
1189: if (erreur > 0)
1190: {
1191: (*s_etat_processus).exception = d_ep_decomposition_QZ;
1192: }
1193: else
1194: {
1195: (*s_etat_processus).erreur_execution =
1196: d_ex_routines_mathematiques;
1197: }
1198:
1199: free(work);
1200: free(rwork);
1201: free(matrice_f77);
1202: free(metrique_f77);
1203:
1204: if (calcul_vp_gauches == 'V')
1205: {
1206: free(vpg_f77);
1207:
1208: (*s_vecteurs_propres_gauches).type = 'C';
1209: (*s_vecteurs_propres_gauches).nombre_lignes = 1;
1210: (*s_vecteurs_propres_gauches).nombre_colonnes = 1;
1211:
1212: if (((*s_vecteurs_propres_gauches).tableau = malloc(((size_t)
1213: (*s_vecteurs_propres_gauches).nombre_lignes) *
1214: sizeof(struct_complexe16 *))) == NULL)
1215: {
1216: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1217: return;
1218: }
1219:
1220: if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches)
1221: .tableau)[0] = (struct_complexe16 *) malloc(((size_t)
1222: (*s_vecteurs_propres_gauches).nombre_colonnes) *
1223: sizeof(struct_complexe16))) == NULL)
1224: {
1225: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1226: return;
1227: }
1228: }
1229:
1230: if (calcul_vp_droits == 'V')
1231: {
1232: free(vpd_f77);
1233:
1234: (*s_vecteurs_propres_droits).type = 'C';
1235: (*s_vecteurs_propres_droits).nombre_lignes = 1;
1236: (*s_vecteurs_propres_droits).nombre_colonnes = 1;
1237:
1238: if (((*s_vecteurs_propres_droits).tableau = malloc(((size_t)
1239: (*s_vecteurs_propres_droits).nombre_lignes) *
1240: sizeof(struct_complexe16 *))) == NULL)
1241: {
1242: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1243: return;
1244: }
1245:
1246: if ((((struct_complexe16 **) (*s_vecteurs_propres_droits)
1247: .tableau)[0] = (struct_complexe16 *) malloc(((size_t)
1248: (*s_vecteurs_propres_droits).nombre_colonnes) *
1249: sizeof(struct_complexe16))) == NULL)
1250: {
1251: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1252: return;
1253: }
1254: }
1255:
1256: return;
1257: }
1258:
1259: lwork = (integer4) work[0].partie_reelle;
1260: free(work);
1261:
1262: if ((work = (struct_complexe16 *) malloc(((size_t) lwork) *
1263: sizeof(struct_complexe16))) == NULL)
1264: {
1265: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1266: return;
1267: }
1268:
1269: zggev_(&calcul_vp_gauches, &calcul_vp_droits, &dim_matrice,
1270: matrice_f77, &dim_matrice,
1271: metrique_f77, &dim_matrice, alpha, beta,
1272: vpg_f77, &dim_matrice, vpd_f77, &dim_matrice,
1273: work, &lwork, rwork, &erreur, 1, 1);
1274:
1275: if (erreur != 0)
1276: {
1277: if (erreur > 0)
1278: {
1279: (*s_etat_processus).exception = d_ep_decomposition_QZ;
1280: }
1281: else
1282: {
1283: (*s_etat_processus).erreur_execution =
1284: d_ex_routines_mathematiques;
1285: }
1286:
1287: free(work);
1288: free(rwork);
1289: free(matrice_f77);
1290: free(metrique_f77);
1291:
1292: if (calcul_vp_gauches == 'V')
1293: {
1294: free(vpg_f77);
1295:
1296: (*s_vecteurs_propres_gauches).type = 'C';
1297: (*s_vecteurs_propres_gauches).nombre_lignes = 1;
1298: (*s_vecteurs_propres_gauches).nombre_colonnes = 1;
1299:
1300: if (((*s_vecteurs_propres_gauches).tableau = malloc(((size_t)
1301: (*s_vecteurs_propres_gauches).nombre_lignes) *
1302: sizeof(struct_complexe16 *))) == NULL)
1303: {
1304: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1305: return;
1306: }
1307:
1308: if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches)
1309: .tableau)[0] = (struct_complexe16 *) malloc(((size_t)
1310: (*s_vecteurs_propres_gauches).nombre_colonnes) *
1311: sizeof(struct_complexe16))) == NULL)
1312: {
1313: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1314: return;
1315: }
1316: }
1317:
1318: if (calcul_vp_droits == 'V')
1319: {
1320: free(vpd_f77);
1321:
1322: (*s_vecteurs_propres_droits).type = 'C';
1323: (*s_vecteurs_propres_droits).nombre_lignes = 1;
1324: (*s_vecteurs_propres_droits).nombre_colonnes = 1;
1325:
1326: if (((*s_vecteurs_propres_droits).tableau = malloc(((size_t)
1327: (*s_vecteurs_propres_droits).nombre_lignes) *
1328: sizeof(struct_complexe16 *))) == NULL)
1329: {
1330: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1331: return;
1332: }
1333:
1334: if ((((struct_complexe16 **) (*s_vecteurs_propres_droits)
1335: .tableau)[0] = (struct_complexe16 *) malloc(((size_t)
1336: (*s_vecteurs_propres_droits).nombre_colonnes) *
1337: sizeof(struct_complexe16))) == NULL)
1338: {
1339: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1340: return;
1341: }
1342: }
1343:
1344: return;
1345: }
1346:
1347: for(i = 0; i < (*s_valeurs_propres).taille; i++)
1348: {
1349: if ((beta[i].partie_reelle != 0) || (beta[i].partie_imaginaire != 0))
1350: {
1351: f77divisioncc_(&(alpha[i]), &(beta[i]), &(((struct_complexe16 *)
1352: (*s_valeurs_propres).tableau)[i]));
1353: }
1354: else
1355: {
1356: ((struct_complexe16 *) (*s_valeurs_propres).tableau)[i]
1357: .partie_reelle = 0;
1358: ((struct_complexe16 *) (*s_valeurs_propres).tableau)[i]
1359: .partie_imaginaire = 0;
1360:
1361: (*s_etat_processus).exception = d_ep_division_par_zero;
1362: }
1363: }
1364:
1365: free(alpha);
1366: free(beta);
1367:
1368: free(work);
1369: free(rwork);
1370:
1371: if (calcul_vp_gauches == 'V')
1372: {
1373: (*s_vecteurs_propres_gauches).type = 'C';
1374: (*s_vecteurs_propres_gauches).nombre_lignes =
1375: (*s_matrice).nombre_lignes;
1376: (*s_vecteurs_propres_gauches).nombre_colonnes =
1377: (*s_matrice).nombre_colonnes;
1378:
1379: if (((*s_vecteurs_propres_gauches).tableau = malloc(((size_t)
1380: (*s_vecteurs_propres_gauches).nombre_lignes) *
1381: sizeof(struct_complexe16 *))) == NULL)
1382: {
1383: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1384: return;
1385: }
1386:
1387: for(i = 0; i < (*s_vecteurs_propres_gauches).nombre_lignes; i++)
1388: {
1389: if ((((struct_complexe16 **) (*s_vecteurs_propres_gauches)
1390: .tableau)[i] = (struct_complexe16 *) malloc(((size_t)
1391: (*s_vecteurs_propres_gauches).nombre_colonnes) *
1392: sizeof(struct_complexe16))) == NULL)
1393: {
1394: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1395: return;
1396: }
1397: }
1398:
1399: for(k = 0, i = 0; i < (*s_vecteurs_propres_gauches).nombre_colonnes;
1400: i++)
1401: {
1402: for(j = 0; j < (*s_vecteurs_propres_gauches).nombre_lignes; j++)
1403: {
1404: ((struct_complexe16 **) (*s_vecteurs_propres_gauches)
1405: .tableau)[j][i].partie_reelle =
1406: vpg_f77[k].partie_reelle;
1407: ((struct_complexe16 **) (*s_vecteurs_propres_gauches)
1408: .tableau)[j][i].partie_imaginaire =
1409: vpg_f77[k++].partie_imaginaire;
1410: }
1411: }
1412:
1413: free(vpg_f77);
1414: }
1415:
1416: if (calcul_vp_droits == 'V')
1417: {
1418: (*s_vecteurs_propres_droits).type = 'C';
1419: (*s_vecteurs_propres_droits).nombre_lignes =
1420: (*s_matrice).nombre_lignes;
1421: (*s_vecteurs_propres_droits).nombre_colonnes =
1422: (*s_matrice).nombre_colonnes;
1423:
1424: if (((*s_vecteurs_propres_droits).tableau = malloc(((size_t)
1425: (*s_vecteurs_propres_droits).nombre_lignes) *
1426: sizeof(struct_complexe16 *))) == NULL)
1427: {
1428: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1429: return;
1430: }
1431:
1432: for(i = 0; i < (*s_vecteurs_propres_droits).nombre_lignes; i++)
1433: {
1434: if ((((struct_complexe16 **) (*s_vecteurs_propres_droits)
1435: .tableau)[i] = (struct_complexe16 *) malloc(((size_t)
1436: (*s_vecteurs_propres_droits).nombre_colonnes) *
1437: sizeof(struct_complexe16))) == NULL)
1438: {
1439: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1440: return;
1441: }
1442: }
1443:
1444: for(k = 0, i = 0; i < (*s_vecteurs_propres_droits).nombre_colonnes; i++)
1445: {
1446: for(j = 0; j < (*s_vecteurs_propres_droits).nombre_lignes; j++)
1447: {
1448: ((struct_complexe16 **) (*s_vecteurs_propres_droits)
1449: .tableau)[j][i].partie_reelle =
1450: vpd_f77[k].partie_reelle;
1451: ((struct_complexe16 **) (*s_vecteurs_propres_droits)
1452: .tableau)[j][i].partie_imaginaire =
1453: vpd_f77[k++].partie_imaginaire;
1454: }
1455: }
1456:
1457: free(vpd_f77);
1458: }
1459:
1460: free(matrice_f77);
1461: free(metrique_f77);
1462:
1463: return;
1464: }
1465:
1466:
1467: /*
1468: ================================================================================
1469: Fonction réalisation le calcul d'un moindres carrés
1470: (minimise [B-AX]^2)
1471: ================================================================================
1472: Entrées : struct_matrice
1473: --------------------------------------------------------------------------------
1474: Sorties : coefficients dans une struct_matrice allouée au vol
1475: --------------------------------------------------------------------------------
1476: Effets de bord : néant
1477: ================================================================================
1478: */
1479:
1480: void
1481: moindres_carres(struct_processus *s_etat_processus,
1482: struct_matrice *s_matrice_a, struct_matrice *s_matrice_b,
1483: struct_matrice *s_matrice_x)
1484: {
1485: real8 *work;
1486: real8 rcond;
1487: real8 *rwork;
1488: real8 *vecteur_s;
1489:
1490: integer4 erreur;
1491: integer4 *iwork;
1492: integer4 lrwork;
1493: integer4 lwork;
1494: integer4 nlvl;
1495: integer4 rank;
1496: integer4 registre_1;
1497: integer4 registre_2;
1498: integer4 registre_3;
1499: integer4 registre_4;
1500: integer4 registre_5;
1501: integer4 smlsiz;
1502:
1503: complex16 *cwork;
1504:
1505: integer8 i;
1506: integer8 j;
1507: integer8 k;
1508: integer8 taille_matrice_a_f77;
1509: integer8 taille_matrice_b_f77;
1510: integer8 taille_matrice_x_f77;
1511:
1512: void *matrice_a_f77;
1513: void *matrice_b_f77;
1514:
1515: taille_matrice_a_f77 = (*s_matrice_a).nombre_lignes *
1516: (*s_matrice_a).nombre_colonnes;
1517: taille_matrice_b_f77 = (*s_matrice_b).nombre_lignes *
1518: (*s_matrice_b).nombre_colonnes;
1519: taille_matrice_x_f77 = (*s_matrice_b).nombre_colonnes *
1520: (*s_matrice_a).nombre_colonnes;
1521:
1522: /*
1523: * Résultat réel
1524: */
1525:
1526: if (((*s_matrice_a).type != 'C') && ((*s_matrice_b).type != 'C'))
1527: {
1528:
1529: /*
1530: * Garniture de la matrice A
1531: */
1532:
1533: if ((matrice_a_f77 = malloc(((size_t) taille_matrice_a_f77) *
1534: sizeof(real8))) == NULL)
1535: {
1536: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1537: return;
1538: }
1539:
1540: if ((*s_matrice_a).type == 'I')
1541: {
1542: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
1543: {
1544: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
1545: {
1546: ((real8 *) matrice_a_f77)[k++] = (real8) ((integer8 **)
1547: (*s_matrice_a).tableau)[j][i];
1548: }
1549: }
1550: }
1551: else
1552: {
1553: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
1554: {
1555: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
1556: {
1557: ((real8 *) matrice_a_f77)[k++] = ((real8 **)
1558: (*s_matrice_a).tableau)[j][i];
1559: }
1560: }
1561: }
1562:
1563: /*
1564: * Garniture de la matrice B
1565: */
1566:
1567: if ((matrice_b_f77 = malloc(((size_t) ((taille_matrice_b_f77
1568: < taille_matrice_x_f77) ? taille_matrice_x_f77
1569: : taille_matrice_b_f77)) * sizeof(real8))) == NULL)
1570: {
1571: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1572: return;
1573: }
1574:
1575: if ((*s_matrice_b).type == 'I')
1576: {
1577: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
1578: {
1579: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
1580: {
1581: ((real8 *) matrice_b_f77)[k++] = (real8) ((integer8 **)
1582: (*s_matrice_b).tableau)[j][i];
1583: }
1584:
1585: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
1586: }
1587: }
1588: else
1589: {
1590: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
1591: {
1592: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
1593: {
1594: ((real8 *) matrice_b_f77)[k++] = ((real8 **)
1595: (*s_matrice_b).tableau)[j][i];
1596: }
1597:
1598: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
1599: }
1600: }
1601:
1602: rcond = -1;
1603:
1604: registre_1 = 9;
1605: registre_2 = 0;
1606:
1607: smlsiz = ilaenv_(®istre_1, "DGELSD", " ", ®istre_2,
1608: ®istre_2, ®istre_2, ®istre_2, 6, 1);
1609:
1610: nlvl = 1 + ((integer4) (log(((real8) (((*s_matrice_a).nombre_lignes <
1611: (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes
1612: : (*s_matrice_a).nombre_colonnes)) / (smlsiz + 1)) /
1613: log((real8) 2)));
1614:
1615: if ((iwork = (integer4 *) malloc(((size_t) ((((*s_matrice_a)
1616: .nombre_lignes < (*s_matrice_a).nombre_colonnes)
1617: ? (*s_matrice_a).nombre_lignes
1618: : (*s_matrice_a).nombre_colonnes) * (11 + (3 * nlvl)))) *
1619: sizeof(integer4))) == NULL)
1620: {
1621: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1622: return;
1623: }
1624:
1625: registre_1 = (integer4) (*s_matrice_a).nombre_lignes;
1626: registre_2 = (integer4) (*s_matrice_a).nombre_colonnes;
1627: registre_3 = (integer4) (*s_matrice_b).nombre_colonnes;
1628: registre_4 = (registre_1 > registre_2) ? registre_1 : registre_2;
1629: registre_5 = (registre_1 > registre_2) ? registre_2 : registre_1;
1630:
1631: if ((work = malloc(sizeof(real8))) == NULL)
1632: {
1633: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1634: return;
1635: }
1636:
1637: if ((vecteur_s = (real8 *) malloc(((size_t) registre_5) *
1638: sizeof(real8))) == NULL)
1639: {
1640: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1641: return;
1642: }
1643:
1644: lwork = -1;
1645:
1646: dgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77,
1647: ®istre_1, matrice_b_f77, ®istre_4, vecteur_s,
1648: &rcond, &rank, work, &lwork, iwork, &erreur);
1649:
1650: if (erreur != 0)
1651: {
1652: if (erreur > 0)
1653: {
1654: (*s_etat_processus).exception = d_ep_decomposition_SVD;
1655: }
1656: else
1657: {
1658: (*s_etat_processus).erreur_execution =
1659: d_ex_routines_mathematiques;
1660: }
1661:
1662: free(work);
1663: free(iwork);
1664: free(matrice_a_f77);
1665: free(matrice_b_f77);
1666: return;
1667: }
1668:
1669: lwork = (integer4) work[0];
1670: free(work);
1671:
1672: if ((work = malloc(((size_t) lwork) * sizeof(real8))) == NULL)
1673: {
1674: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1675: return;
1676: }
1677:
1678: dgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77,
1679: ®istre_1, matrice_b_f77, ®istre_4, vecteur_s,
1680: &rcond, &rank, work, &lwork, iwork, &erreur);
1681:
1682: free(vecteur_s);
1683:
1684: if (erreur != 0)
1685: {
1686: if (erreur > 0)
1687: {
1688: (*s_etat_processus).exception = d_ep_decomposition_SVD;
1689: }
1690: else
1691: {
1692: (*s_etat_processus).erreur_execution =
1693: d_ex_routines_mathematiques;
1694: }
1695:
1696: free(work);
1697: free(iwork);
1698: free(matrice_a_f77);
1699: free(matrice_b_f77);
1700: return;
1701: }
1702:
1703: free(work);
1704: free(iwork);
1705:
1706: (*s_matrice_x).nombre_lignes = (*s_matrice_a).nombre_colonnes;
1707: (*s_matrice_x).nombre_colonnes = (*s_matrice_b).nombre_colonnes;
1708:
1709: if (((*s_matrice_x).tableau = malloc(((size_t) (*s_matrice_x)
1710: .nombre_lignes) * sizeof(real8 *))) == NULL)
1711: {
1712: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1713: return;
1714: }
1715:
1716: for(j = 0; j < (*s_matrice_x).nombre_lignes; j++)
1717: {
1718: if ((((real8 **) (*s_matrice_x).tableau)[j] = (real8 *)
1719: malloc(((size_t) (*s_matrice_x).nombre_colonnes) *
1720: sizeof(real8)))== NULL)
1721: {
1722: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1723: return;
1724: }
1725: }
1726:
1727: for(k = 0, i = 0; i < (*s_matrice_x).nombre_colonnes; i++)
1728: {
1729: for(j = 0; j < (*s_matrice_x).nombre_lignes; j++)
1730: {
1731: (((real8 **) (*s_matrice_x).tableau)[j][i]) = (((real8 *)
1732: matrice_b_f77)[k++]);
1733: }
1734:
1735: for(; j < (*s_matrice_b).nombre_lignes; j++, k++);
1736: }
1737: }
1738:
1739: /*
1740: * Résultat complexe
1741: */
1742:
1743: else
1744: {
1745:
1746: /*
1747: * Garniture de la matrice A
1748: */
1749:
1750: if ((matrice_a_f77 = malloc(((size_t) taille_matrice_a_f77) *
1751: sizeof(struct_complexe16))) == NULL)
1752: {
1753: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1754: return;
1755: }
1756:
1757: if ((*s_matrice_a).type == 'I')
1758: {
1759: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
1760: {
1761: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
1762: {
1763: ((struct_complexe16 *) matrice_a_f77)[k].partie_reelle =
1764: (real8) ((integer8 **) (*s_matrice_a)
1765: .tableau)[j][i];
1766: ((struct_complexe16 *) matrice_a_f77)[k++]
1767: .partie_imaginaire = (real8) 0;
1768: }
1769: }
1770: }
1771: else if ((*s_matrice_a).type == 'R')
1772: {
1773: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
1774: {
1775: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
1776: {
1777: ((struct_complexe16 *) matrice_a_f77)[k].partie_reelle =
1778: ((real8 **) (*s_matrice_a).tableau)[j][i];
1779: ((struct_complexe16 *) matrice_a_f77)[k++]
1780: .partie_imaginaire = (real8) 0;
1781: }
1782: }
1783: }
1784: else
1785: {
1786: for(k = 0, i = 0; i < (*s_matrice_a).nombre_colonnes; i++)
1787: {
1788: for(j = 0; j < (*s_matrice_a).nombre_lignes; j++)
1789: {
1790: ((struct_complexe16 *) matrice_a_f77)[k].partie_reelle =
1791: ((struct_complexe16 **) (*s_matrice_a).tableau)
1792: [j][i].partie_reelle;
1793: ((struct_complexe16 *) matrice_a_f77)[k++]
1794: .partie_imaginaire = ((struct_complexe16 **)
1795: (*s_matrice_a).tableau)[j][i].partie_imaginaire;
1796: }
1797: }
1798: }
1799:
1800: /*
1801: * Garniture de la matrice B
1802: */
1803:
1804: if ((matrice_b_f77 = malloc(((size_t) ((taille_matrice_b_f77
1805: < taille_matrice_x_f77) ? taille_matrice_x_f77
1806: : taille_matrice_b_f77)) * sizeof(struct_complexe16))) == NULL)
1807: {
1808: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1809: return;
1810: }
1811:
1812: if ((*s_matrice_b).type == 'I')
1813: {
1814: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
1815: {
1816: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
1817: {
1818: ((struct_complexe16 *) matrice_b_f77)[k].partie_reelle =
1819: (real8) ((integer8 **) (*s_matrice_b)
1820: .tableau)[j][i];
1821: ((struct_complexe16 *) matrice_b_f77)[k++]
1822: .partie_imaginaire = (real8) 0;
1823: }
1824:
1825: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
1826: }
1827: }
1828: else if ((*s_matrice_b).type == 'R')
1829: {
1830: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
1831: {
1832: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
1833: {
1834: ((struct_complexe16 *) matrice_b_f77)[k].partie_reelle =
1835: ((real8 **) (*s_matrice_b).tableau)[j][i];
1836: ((struct_complexe16 *) matrice_b_f77)[k++]
1837: .partie_imaginaire = (real8) 0;
1838: }
1839:
1840: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
1841: }
1842: }
1843: else
1844: {
1845: for(k = 0, i = 0; i < (*s_matrice_b).nombre_colonnes; i++)
1846: {
1847: for(j = 0; j < (*s_matrice_b).nombre_lignes; j++)
1848: {
1849: ((struct_complexe16 *) matrice_b_f77)[k].partie_reelle =
1850: ((struct_complexe16 **) (*s_matrice_b).tableau)
1851: [j][i].partie_reelle;
1852: ((struct_complexe16 *) matrice_b_f77)[k++]
1853: .partie_imaginaire = ((struct_complexe16 **)
1854: (*s_matrice_b).tableau)[j][i].partie_imaginaire;
1855: }
1856:
1857: for(; j < (*s_matrice_a).nombre_lignes; j++, k++);
1858: }
1859: }
1860:
1861: rcond = -1;
1862:
1863: registre_1 = 9;
1864: registre_2 = 0;
1865:
1866: smlsiz = ilaenv_(®istre_1, "ZGELSD", " ", ®istre_2,
1867: ®istre_2, ®istre_2, ®istre_2, 6, 1);
1868:
1869: nlvl = 1 + ((integer4) (log(((real8) (((*s_matrice_a).nombre_lignes <
1870: (*s_matrice_a).nombre_colonnes) ? (*s_matrice_a).nombre_lignes
1871: : (*s_matrice_a).nombre_colonnes)) / (smlsiz + 1))
1872: / log((real8) 2)));
1873:
1874: if ((*s_matrice_a).nombre_lignes >= (*s_matrice_a).nombre_colonnes)
1875: {
1876: lrwork = (integer4) ((*s_matrice_a).nombre_colonnes * (8 + (2 *
1877: smlsiz) + (8 * nlvl) + (*s_matrice_b).nombre_colonnes));
1878: }
1879: else
1880: {
1881: lrwork = (integer4) ((*s_matrice_a).nombre_lignes * (8 + (2 *
1882: smlsiz) + (8 * nlvl) + (*s_matrice_b).nombre_colonnes));
1883: }
1884:
1885: if ((rwork = (real8 *) malloc(((size_t) lrwork) * sizeof(real8)))
1886: == NULL)
1887: {
1888: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1889: return;
1890: }
1891:
1892: if ((iwork = (integer4 *) malloc(((size_t) ((((*s_matrice_a)
1893: .nombre_lignes < (*s_matrice_a).nombre_colonnes)
1894: ? (*s_matrice_a).nombre_lignes
1895: : (*s_matrice_a).nombre_colonnes) * (11 + (3 * nlvl)))) *
1896: sizeof(integer4))) == NULL)
1897: {
1898: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1899: return;
1900: }
1901:
1902: registre_1 = (integer4) (*s_matrice_a).nombre_lignes;
1903: registre_2 = (integer4) (*s_matrice_a).nombre_colonnes;
1904: registre_3 = (integer4) (*s_matrice_b).nombre_colonnes;
1905: registre_4 = (registre_1 > registre_2) ? registre_1 : registre_2;
1906: registre_5 = (registre_1 > registre_2) ? registre_2 : registre_1;
1907:
1908: if ((cwork = malloc(sizeof(struct_complexe16))) == NULL)
1909: {
1910: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1911: return;
1912: }
1913:
1914: if ((vecteur_s = (real8 *) malloc(((size_t) registre_5) *
1915: sizeof(real8))) == NULL)
1916: {
1917: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1918: return;
1919: }
1920:
1921: lwork = -1;
1922:
1923: zgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77,
1924: ®istre_1, matrice_b_f77, ®istre_4, vecteur_s,
1925: &rcond, &rank, cwork, &lwork, rwork, iwork, &erreur);
1926:
1927: if (erreur != 0)
1928: {
1929: if (erreur > 0)
1930: {
1931: (*s_etat_processus).exception = d_ep_decomposition_SVD;
1932: }
1933: else
1934: {
1935: (*s_etat_processus).erreur_execution =
1936: d_ex_routines_mathematiques;
1937: }
1938:
1939: free(cwork);
1940: free(iwork);
1941: free(rwork);
1942: free(matrice_a_f77);
1943: free(matrice_b_f77);
1944: return;
1945: }
1946:
1947: lwork = (integer4) cwork[0].partie_reelle;
1948: free(cwork);
1949:
1950: if ((cwork = malloc(((size_t) lwork) * sizeof(struct_complexe16)))
1951: == NULL)
1952: {
1953: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1954: return;
1955: }
1956:
1957: zgelsd_(®istre_1, ®istre_2, ®istre_3, matrice_a_f77,
1958: ®istre_1, matrice_b_f77, ®istre_4, vecteur_s,
1959: &rcond, &rank, cwork, &lwork, rwork, iwork, &erreur);
1960:
1961: free(vecteur_s);
1962:
1963: if (erreur != 0)
1964: {
1965: if (erreur > 0)
1966: {
1967: (*s_etat_processus).exception = d_ep_decomposition_SVD;
1968: }
1969: else
1970: {
1971: (*s_etat_processus).erreur_execution =
1972: d_ex_routines_mathematiques;
1973: }
1974:
1975: free(cwork);
1976: free(iwork);
1977: free(rwork);
1978: free(matrice_a_f77);
1979: free(matrice_b_f77);
1980: return;
1981: }
1982:
1983: free(cwork);
1984: free(iwork);
1985: free(rwork);
1986:
1987: (*s_matrice_x).nombre_lignes = (*s_matrice_a).nombre_colonnes;
1988: (*s_matrice_x).nombre_colonnes = (*s_matrice_b).nombre_colonnes;
1989:
1990: if (((*s_matrice_x).tableau = malloc(((size_t) (*s_matrice_x)
1991: .nombre_lignes) * sizeof(struct_complexe16 *))) == NULL)
1992: {
1993: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
1994: return;
1995: }
1996:
1997: for(j = 0; j < (*s_matrice_x).nombre_lignes; j++)
1998: {
1999: if ((((struct_complexe16 **) (*s_matrice_x).tableau)[j] =
2000: (struct_complexe16 *) malloc(((size_t)
2001: (*s_matrice_x).nombre_colonnes)
2002: * sizeof(struct_complexe16)))== NULL)
2003: {
2004: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
2005: return;
2006: }
2007: }
2008:
2009: for(k = 0, i = 0; i < (*s_matrice_x).nombre_colonnes; i++)
2010: {
2011: for(j = 0; j < (*s_matrice_x).nombre_lignes; j++)
2012: {
2013: (((struct_complexe16 **) (*s_matrice_x).tableau)[j][i])
2014: .partie_reelle = (((struct_complexe16 *)
2015: matrice_b_f77)[k]).partie_reelle;
2016: (((struct_complexe16 **) (*s_matrice_x).tableau)[j][i])
2017: .partie_imaginaire = (((struct_complexe16 *)
2018: matrice_b_f77)[k++]).partie_imaginaire;
2019: }
2020:
2021: for(; j < (*s_matrice_b).nombre_lignes; j++, k++);
2022: }
2023: }
2024:
2025: free(matrice_a_f77);
2026: free(matrice_b_f77);
2027:
2028: return;
2029: }
2030:
2031: // vim: ts=4
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to algebre_lineaire1.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / src