1: /*
2: ================================================================================
3: RPL/2 (R) version 4.1.32
4: Copyright (C) 1989-2020 Dr. BERTRAND Joël
5:
6: This file is part of RPL/2.
7:
8: RPL/2 is free software; you can redistribute it and/or modify it
9: under the terms of the CeCILL V2 License as published by the french
10: CEA, CNRS and INRIA.
11:
12: RPL/2 is distributed in the hope that it will be useful, but WITHOUT
13: ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14: FITNESS FOR A PARTICULAR PURPOSE. See the CeCILL V2 License
15: for more details.
16:
17: You should have received a copy of the CeCILL License
18: along with RPL/2. If not, write to info@cecill.info.
19: ================================================================================
20: */
21:
22:
23: #include "rpl-conv.h"
24:
25:
26: /*
27: ================================================================================
28: Fonction 'qr'
29: ================================================================================
30: Entrées : pointeur sur une structure struct_processus
31: --------------------------------------------------------------------------------
32: Sorties :
33: --------------------------------------------------------------------------------
34: Effets de bord : néant
35: ================================================================================
36: */
37:
38: void
39: instruction_qr(struct_processus *s_etat_processus)
40: {
41: complex16 registre;
42: complex16 *tau_complexe;
43: complex16 *vecteur_complexe;
44:
45: real8 *tau_reel;
46: real8 *vecteur_reel;
47:
48: struct_liste_chainee *registre_pile_last;
49:
50: struct_objet *s_copie_argument;
51: struct_objet *s_matrice_identite;
52: struct_objet *s_objet;
53: struct_objet *s_objet_argument;
54: struct_objet *s_objet_resultat;
55:
56: integer8 i;
57: integer8 j;
58: integer8 k;
59: integer8 nombre_reflecteurs_elementaires;
60:
61: void *tau;
62:
63: (*s_etat_processus).erreur_execution = d_ex;
64:
65: if ((*s_etat_processus).affichage_arguments == 'Y')
66: {
67: printf("\n QR ");
68:
69: if ((*s_etat_processus).langue == 'F')
70: {
71: printf("(décomposition QR)\n\n");
72: }
73: else
74: {
75: printf("(QR décomposition)\n\n");
76: }
77:
78: printf(" 1: %s, %s\n", d_MIN, d_MRL);
79: printf("-> 2: %s\n", d_MRL);
80: printf(" 1: %s\n\n", d_MRL);
81:
82: printf(" 1: %s\n", d_MCX);
83: printf("-> 2: %s\n", d_MCX);
84: printf(" 1: %s\n", d_MCX);
85:
86: return;
87: }
88: else if ((*s_etat_processus).test_instruction == 'Y')
89: {
90: (*s_etat_processus).nombre_arguments = -1;
91: return;
92: }
93:
94: if (test_cfsf(s_etat_processus, 31) == d_vrai)
95: {
96: if (empilement_pile_last(s_etat_processus, 1) == d_erreur)
97: {
98: return;
99: }
100: }
101:
102: if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
103: &s_objet_argument) == d_erreur)
104: {
105: (*s_etat_processus).erreur_execution = d_ex_manque_argument;
106: return;
107: }
108:
109: if (((*s_objet_argument).type == MIN) ||
110: ((*s_objet_argument).type == MRL))
111: {
112: /*
113: * Matrice entière ou réelle
114: */
115:
116: if ((s_copie_argument = copie_objet(s_etat_processus,
117: s_objet_argument, 'Q')) == NULL)
118: {
119: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
120: return;
121: }
122:
123: factorisation_qr(s_etat_processus, (*s_copie_argument).objet, &tau);
124: (*s_copie_argument).type = MRL;
125:
126: tau_reel = (real8 *) tau;
127:
128: if ((*s_etat_processus).erreur_systeme != d_es)
129: {
130: return;
131: }
132:
133: if (((*s_etat_processus).exception != d_ep) ||
134: ((*s_etat_processus).erreur_execution != d_ex))
135: {
136: free(tau);
137: liberation(s_etat_processus, s_objet_argument);
138: liberation(s_etat_processus, s_copie_argument);
139: return;
140: }
141:
142: if ((s_objet_resultat = copie_objet(s_etat_processus,
143: s_copie_argument, 'O')) == NULL)
144: {
145: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
146: return;
147: }
148:
149: // Matrice Q
150:
151: nombre_reflecteurs_elementaires = ((*((struct_matrice *)
152: (*s_copie_argument).objet)).nombre_colonnes <
153: (*((struct_matrice *) (*s_copie_argument).objet))
154: .nombre_lignes) ? (*((struct_matrice *)
155: (*s_copie_argument).objet)).nombre_colonnes
156: : (*((struct_matrice *) (*s_copie_argument).objet))
157: .nombre_lignes;
158:
159: registre_pile_last = NULL;
160:
161: if (test_cfsf(s_etat_processus, 31) == d_vrai)
162: {
163: registre_pile_last = (*s_etat_processus).l_base_pile_last;
164: (*s_etat_processus).l_base_pile_last = NULL;
165: }
166:
167: if ((s_objet = allocation(s_etat_processus, INT)) == NULL)
168: {
169: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
170: return;
171: }
172:
173: (*((integer8 *) (*s_objet).objet)) = (*((struct_matrice *)
174: (*s_copie_argument).objet)).nombre_lignes;
175:
176: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
177: s_objet) == d_erreur)
178: {
179: return;
180: }
181:
182: instruction_idn(s_etat_processus);
183:
184: if (((*s_etat_processus).erreur_systeme != d_es) ||
185: ((*s_etat_processus).erreur_execution != d_ex) ||
186: ((*s_etat_processus).exception != d_ep))
187: {
188: liberation(s_etat_processus, s_copie_argument);
189: free(tau);
190:
191: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
192: {
193: return;
194: }
195:
196: (*s_etat_processus).l_base_pile_last = registre_pile_last;
197: return;
198: }
199:
200: if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
201: &s_matrice_identite) == d_erreur)
202: {
203: (*s_etat_processus).erreur_execution = d_ex_manque_argument;
204: return;
205: }
206:
207: for(i = 0; i < nombre_reflecteurs_elementaires; i++)
208: {
209: // Calcul de H(i) = I - tau * v * v'
210:
211: if ((s_objet = copie_objet(s_etat_processus, s_matrice_identite,
212: 'P')) == NULL)
213: {
214: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
215: return;
216: }
217:
218: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
219: s_objet) == d_erreur)
220: {
221: return;
222: }
223:
224: if ((s_objet = allocation(s_etat_processus, REL)) == NULL)
225: {
226: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
227: return;
228: }
229:
230: (*((real8 *) (*s_objet).objet)) = tau_reel[i];
231:
232: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
233: s_objet) == d_erreur)
234: {
235: return;
236: }
237:
238: if ((s_objet = allocation(s_etat_processus, MRL)) == NULL)
239: {
240: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
241: return;
242: }
243:
244: (*((struct_matrice *) (*s_objet).objet)).nombre_lignes =
245: (*((struct_matrice *) (*s_copie_argument).objet))
246: .nombre_lignes;
247: (*((struct_matrice *) (*s_objet).objet)).nombre_colonnes =
248: (*((struct_matrice *) (*s_copie_argument).objet))
249: .nombre_lignes;
250:
251: if ((vecteur_reel = malloc(((size_t) (*((struct_matrice *)
252: (*s_objet).objet)).nombre_lignes) * sizeof(real8))) == NULL)
253: {
254: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
255: return;
256: }
257:
258: for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
259: .nombre_lignes; j++)
260: {
261: if (j < i)
262: {
263: vecteur_reel[j] = 0;
264: }
265: else if (j == i)
266: {
267: vecteur_reel[j] = 1;
268: }
269: else
270: {
271: vecteur_reel[j] = ((real8 **) (*((struct_matrice *)
272: (*s_copie_argument).objet)).tableau)[j][i];
273: }
274: }
275:
276: if (((*((struct_matrice *) (*s_objet).objet)).tableau =
277: malloc(((size_t) (*((struct_matrice *) (*s_objet).objet))
278: .nombre_lignes) * sizeof(real8 *))) == NULL)
279: {
280: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
281: return;
282: }
283:
284: for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
285: .nombre_lignes; j++)
286: {
287: if ((((real8 **) (*((struct_matrice *) (*s_objet).objet))
288: .tableau)[j] = malloc(((size_t) (*((struct_matrice *)
289: (*s_objet).objet)).nombre_lignes) * sizeof(real8)))
290: == NULL)
291: {
292: (*s_etat_processus).erreur_systeme =
293: d_es_allocation_memoire;
294: return;
295: }
296:
297: for(k = 0; k < (*((struct_matrice *) (*s_objet).objet))
298: .nombre_colonnes; k++)
299: {
300: ((real8 **) (*((struct_matrice *) (*s_objet).objet))
301: .tableau)[j][k] = vecteur_reel[j] * vecteur_reel[k];
302: }
303: }
304:
305: free(vecteur_reel);
306:
307: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
308: s_objet) == d_erreur)
309: {
310: return;
311: }
312:
313: instruction_multiplication(s_etat_processus);
314:
315: if (((*s_etat_processus).erreur_systeme != d_es) ||
316: ((*s_etat_processus).erreur_execution != d_ex) ||
317: ((*s_etat_processus).exception != d_ep))
318: {
319: liberation(s_etat_processus, s_copie_argument);
320: liberation(s_etat_processus, s_matrice_identite);
321: free(tau);
322:
323: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
324: {
325: return;
326: }
327:
328: (*s_etat_processus).l_base_pile_last = registre_pile_last;
329: return;
330: }
331:
332: instruction_moins(s_etat_processus);
333:
334: if (((*s_etat_processus).erreur_systeme != d_es) ||
335: ((*s_etat_processus).erreur_execution != d_ex) ||
336: ((*s_etat_processus).exception != d_ep))
337: {
338: liberation(s_etat_processus, s_copie_argument);
339: liberation(s_etat_processus, s_matrice_identite);
340: free(tau);
341:
342: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
343: {
344: return;
345: }
346:
347: (*s_etat_processus).l_base_pile_last = registre_pile_last;
348: return;
349: }
350:
351: if (i > 0)
352: {
353: instruction_multiplication(s_etat_processus);
354:
355: if (((*s_etat_processus).erreur_systeme != d_es) ||
356: ((*s_etat_processus).erreur_execution != d_ex) ||
357: ((*s_etat_processus).exception != d_ep))
358: {
359: liberation(s_etat_processus, s_copie_argument);
360: liberation(s_etat_processus, s_matrice_identite);
361: free(tau);
362:
363: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
364: {
365: return;
366: }
367:
368: (*s_etat_processus).l_base_pile_last = registre_pile_last;
369: return;
370: }
371: }
372: }
373:
374: if (test_cfsf(s_etat_processus, 31) == d_vrai)
375: {
376: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
377: {
378: return;
379: }
380:
381: (*s_etat_processus).l_base_pile_last = registre_pile_last;
382: }
383:
384: // Matrice R
385:
386: for(i = 0; i < (*((struct_matrice *) (*s_objet_resultat).objet))
387: .nombre_lignes; i++)
388: {
389: for(j = 0; j < i; j++)
390: {
391: ((real8 **) (*((struct_matrice *) (*s_objet_resultat).objet))
392: .tableau)[i][j] = 0;
393: }
394: }
395:
396: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
397: s_objet_resultat) == d_erreur)
398: {
399: return;
400: }
401:
402: liberation(s_etat_processus, s_matrice_identite);
403: liberation(s_etat_processus, s_copie_argument);
404: free(tau);
405: }
406: else if ((*s_objet_argument).type == MCX)
407: {
408: /*
409: * Matrice complexe
410: */
411:
412: if ((s_copie_argument = copie_objet(s_etat_processus,
413: s_objet_argument, 'Q')) == NULL)
414: {
415: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
416: return;
417: }
418:
419: factorisation_qr(s_etat_processus, (*s_copie_argument).objet, &tau);
420:
421: tau_complexe = (complex16 *) tau;
422:
423: if ((*s_etat_processus).erreur_systeme != d_es)
424: {
425: return;
426: }
427:
428: if (((*s_etat_processus).exception != d_ep) ||
429: ((*s_etat_processus).erreur_execution != d_ex))
430: {
431: free(tau);
432: liberation(s_etat_processus, s_objet_argument);
433: liberation(s_etat_processus, s_copie_argument);
434: return;
435: }
436:
437: if ((s_objet_resultat = copie_objet(s_etat_processus,
438: s_copie_argument, 'O')) == NULL)
439: {
440: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
441: return;
442: }
443:
444: // Matrice Q
445:
446: nombre_reflecteurs_elementaires = ((*((struct_matrice *)
447: (*s_copie_argument).objet)).nombre_colonnes <
448: (*((struct_matrice *) (*s_copie_argument).objet))
449: .nombre_lignes) ? (*((struct_matrice *)
450: (*s_copie_argument).objet)).nombre_colonnes
451: : (*((struct_matrice *) (*s_copie_argument).objet))
452: .nombre_lignes;
453:
454: registre_pile_last = NULL;
455:
456: if (test_cfsf(s_etat_processus, 31) == d_vrai)
457: {
458: registre_pile_last = (*s_etat_processus).l_base_pile_last;
459: (*s_etat_processus).l_base_pile_last = NULL;
460: }
461:
462: if ((s_objet = allocation(s_etat_processus, INT)) == NULL)
463: {
464: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
465: return;
466: }
467:
468: (*((integer8 *) (*s_objet).objet)) = (*((struct_matrice *)
469: (*s_copie_argument).objet)).nombre_lignes;
470:
471: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
472: s_objet) == d_erreur)
473: {
474: return;
475: }
476:
477: instruction_idn(s_etat_processus);
478:
479: if (((*s_etat_processus).erreur_systeme != d_es) ||
480: ((*s_etat_processus).erreur_execution != d_ex) ||
481: ((*s_etat_processus).exception != d_ep))
482: {
483: liberation(s_etat_processus, s_copie_argument);
484: free(tau);
485:
486: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
487: {
488: return;
489: }
490:
491: (*s_etat_processus).l_base_pile_last = registre_pile_last;
492: return;
493: }
494:
495: if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
496: &s_matrice_identite) == d_erreur)
497: {
498: (*s_etat_processus).erreur_execution = d_ex_manque_argument;
499: return;
500: }
501:
502: for(i = 0; i < nombre_reflecteurs_elementaires; i++)
503: {
504: // Calcul de H(i) = I - tau * v * v'
505:
506: if ((s_objet = copie_objet(s_etat_processus, s_matrice_identite,
507: 'P')) == NULL)
508: {
509: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
510: return;
511: }
512:
513: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
514: s_objet) == d_erreur)
515: {
516: return;
517: }
518:
519: if ((s_objet = allocation(s_etat_processus, CPL)) == NULL)
520: {
521: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
522: return;
523: }
524:
525: (*((complex16 *) (*s_objet).objet)) = tau_complexe[i];
526:
527: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
528: s_objet) == d_erreur)
529: {
530: return;
531: }
532:
533: if ((s_objet = allocation(s_etat_processus, MCX)) == NULL)
534: {
535: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
536: return;
537: }
538:
539: (*((struct_matrice *) (*s_objet).objet)).nombre_lignes =
540: (*((struct_matrice *) (*s_copie_argument).objet))
541: .nombre_lignes;
542: (*((struct_matrice *) (*s_objet).objet)).nombre_colonnes =
543: (*((struct_matrice *) (*s_copie_argument).objet))
544: .nombre_lignes;
545:
546: if ((vecteur_complexe = malloc(((size_t) (*((struct_matrice *)
547: (*s_objet).objet)).nombre_lignes) * sizeof(complex16)))
548: == NULL)
549: {
550: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
551: return;
552: }
553:
554: for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
555: .nombre_lignes; j++)
556: {
557: if (j < i)
558: {
559: vecteur_complexe[j].partie_reelle = 0;
560: vecteur_complexe[j].partie_imaginaire = 0;
561: }
562: else if (j == i)
563: {
564: vecteur_complexe[j].partie_reelle = 1;
565: vecteur_complexe[j].partie_imaginaire = 0;
566: }
567: else
568: {
569: vecteur_complexe[j].partie_reelle = ((complex16 **)
570: (*((struct_matrice *) (*s_copie_argument).objet))
571: .tableau)[j][i].partie_reelle;
572: vecteur_complexe[j].partie_imaginaire = ((complex16 **)
573: (*((struct_matrice *) (*s_copie_argument).objet))
574: .tableau)[j][i].partie_imaginaire;
575: }
576: }
577:
578: if (((*((struct_matrice *) (*s_objet).objet)).tableau =
579: malloc(((size_t) (*((struct_matrice *) (*s_objet).objet))
580: .nombre_lignes) * sizeof(complex16 *))) == NULL)
581: {
582: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
583: return;
584: }
585:
586: for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
587: .nombre_lignes; j++)
588: {
589: if ((((complex16 **) (*((struct_matrice *) (*s_objet).objet))
590: .tableau)[j] = malloc(((size_t) (*((struct_matrice *)
591: (*s_objet).objet)).nombre_lignes) * sizeof(complex16)))
592: == NULL)
593: {
594: (*s_etat_processus).erreur_systeme =
595: d_es_allocation_memoire;
596: return;
597: }
598:
599: for(k = 0; k < (*((struct_matrice *) (*s_objet).objet))
600: .nombre_colonnes; k++)
601: {
602: registre = vecteur_complexe[k];
603: registre.partie_imaginaire =
604: -vecteur_complexe[k].partie_imaginaire;
605:
606: f77multiplicationcc_(&(vecteur_complexe[j]),
607: ®istre, &(((complex16 **)
608: (*((struct_matrice *) (*s_objet).objet)).tableau)
609: [j][k]));
610: }
611: }
612:
613: free(vecteur_complexe);
614:
615: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
616: s_objet) == d_erreur)
617: {
618: return;
619: }
620:
621: instruction_multiplication(s_etat_processus);
622:
623: if (((*s_etat_processus).erreur_systeme != d_es) ||
624: ((*s_etat_processus).erreur_execution != d_ex) ||
625: ((*s_etat_processus).exception != d_ep))
626: {
627: liberation(s_etat_processus, s_copie_argument);
628: liberation(s_etat_processus, s_matrice_identite);
629: free(tau);
630:
631: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
632: {
633: return;
634: }
635:
636: (*s_etat_processus).l_base_pile_last = registre_pile_last;
637: return;
638: }
639:
640: instruction_moins(s_etat_processus);
641:
642: if (((*s_etat_processus).erreur_systeme != d_es) ||
643: ((*s_etat_processus).erreur_execution != d_ex) ||
644: ((*s_etat_processus).exception != d_ep))
645: {
646: liberation(s_etat_processus, s_copie_argument);
647: liberation(s_etat_processus, s_matrice_identite);
648: free(tau);
649:
650: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
651: {
652: return;
653: }
654:
655: (*s_etat_processus).l_base_pile_last = registre_pile_last;
656: return;
657: }
658:
659: if (i > 0)
660: {
661: instruction_multiplication(s_etat_processus);
662:
663: if (((*s_etat_processus).erreur_systeme != d_es) ||
664: ((*s_etat_processus).erreur_execution != d_ex) ||
665: ((*s_etat_processus).exception != d_ep))
666: {
667: liberation(s_etat_processus, s_copie_argument);
668: liberation(s_etat_processus, s_matrice_identite);
669: free(tau);
670:
671: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
672: {
673: return;
674: }
675:
676: (*s_etat_processus).l_base_pile_last = registre_pile_last;
677: return;
678: }
679: }
680: }
681:
682: if (test_cfsf(s_etat_processus, 31) == d_vrai)
683: {
684: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
685: {
686: return;
687: }
688:
689: (*s_etat_processus).l_base_pile_last = registre_pile_last;
690: }
691:
692: // Matrice R
693:
694: for(i = 0; i < (*((struct_matrice *) (*s_objet_resultat).objet))
695: .nombre_lignes; i++)
696: {
697: for(j = 0; j < i; j++)
698: {
699: ((complex16 **) (*((struct_matrice *) (*s_objet_resultat)
700: .objet)).tableau)[i][j].partie_reelle = 0;
701: ((complex16 **) (*((struct_matrice *) (*s_objet_resultat)
702: .objet)).tableau)[i][j].partie_imaginaire = 0;
703: }
704: }
705:
706: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
707: s_objet_resultat) == d_erreur)
708: {
709: return;
710: }
711:
712: liberation(s_etat_processus, s_matrice_identite);
713: liberation(s_etat_processus, s_copie_argument);
714: free(tau);
715: }
716:
717: /*
718: * Type d'argument invalide
719: */
720:
721: else
722: {
723: liberation(s_etat_processus, s_objet_argument);
724:
725: (*s_etat_processus).erreur_execution = d_ex_erreur_type_argument;
726: return;
727: }
728:
729: liberation(s_etat_processus, s_objet_argument);
730:
731: return;
732: }
733:
734: // vim: ts=4
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to instructions_q1.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / src