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