1: /*
2: ================================================================================
3: RPL/2 (R) version 4.1.0
4: Copyright (C) 1989-2011 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: unsigned long i;
57: unsigned long j;
58: unsigned long k;
59: unsigned long 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((*((struct_matrice *) (*s_objet).objet))
252: .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((*((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((*((struct_matrice *) (*s_objet)
289: .objet)).nombre_lignes * sizeof(real8))) == NULL)
290: {
291: (*s_etat_processus).erreur_systeme =
292: d_es_allocation_memoire;
293: return;
294: }
295:
296: for(k = 0; k < (*((struct_matrice *) (*s_objet).objet))
297: .nombre_colonnes; k++)
298: {
299: ((real8 **) (*((struct_matrice *) (*s_objet).objet))
300: .tableau)[j][k] = vecteur_reel[j] * vecteur_reel[k];
301: }
302: }
303:
304: free(vecteur_reel);
305:
306: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
307: s_objet) == d_erreur)
308: {
309: return;
310: }
311:
312: instruction_multiplication(s_etat_processus);
313:
314: if (((*s_etat_processus).erreur_systeme != d_es) ||
315: ((*s_etat_processus).erreur_execution != d_ex) ||
316: ((*s_etat_processus).exception != d_ep))
317: {
318: liberation(s_etat_processus, s_copie_argument);
319: liberation(s_etat_processus, s_matrice_identite);
320: free(tau);
321:
322: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
323: {
324: return;
325: }
326:
327: (*s_etat_processus).l_base_pile_last = registre_pile_last;
328: return;
329: }
330:
331: instruction_moins(s_etat_processus);
332:
333: if (((*s_etat_processus).erreur_systeme != d_es) ||
334: ((*s_etat_processus).erreur_execution != d_ex) ||
335: ((*s_etat_processus).exception != d_ep))
336: {
337: liberation(s_etat_processus, s_copie_argument);
338: liberation(s_etat_processus, s_matrice_identite);
339: free(tau);
340:
341: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
342: {
343: return;
344: }
345:
346: (*s_etat_processus).l_base_pile_last = registre_pile_last;
347: return;
348: }
349:
350: if (i > 0)
351: {
352: instruction_multiplication(s_etat_processus);
353:
354: if (((*s_etat_processus).erreur_systeme != d_es) ||
355: ((*s_etat_processus).erreur_execution != d_ex) ||
356: ((*s_etat_processus).exception != d_ep))
357: {
358: liberation(s_etat_processus, s_copie_argument);
359: liberation(s_etat_processus, s_matrice_identite);
360: free(tau);
361:
362: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
363: {
364: return;
365: }
366:
367: (*s_etat_processus).l_base_pile_last = registre_pile_last;
368: return;
369: }
370: }
371: }
372:
373: if (test_cfsf(s_etat_processus, 31) == d_vrai)
374: {
375: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
376: {
377: return;
378: }
379:
380: (*s_etat_processus).l_base_pile_last = registre_pile_last;
381: }
382:
383: // Matrice R
384:
385: for(i = 0; i < (*((struct_matrice *) (*s_objet_resultat).objet))
386: .nombre_lignes; i++)
387: {
388: for(j = 0; j < i; j++)
389: {
390: ((real8 **) (*((struct_matrice *) (*s_objet_resultat).objet))
391: .tableau)[i][j] = 0;
392: }
393: }
394:
395: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
396: s_objet_resultat) == d_erreur)
397: {
398: return;
399: }
400:
401: liberation(s_etat_processus, s_matrice_identite);
402: liberation(s_etat_processus, s_copie_argument);
403: free(tau);
404: }
405: else if ((*s_objet_argument).type == MCX)
406: {
407: /*
408: * Matrice complexe
409: */
410:
411: if ((s_copie_argument = copie_objet(s_etat_processus,
412: s_objet_argument, 'Q')) == NULL)
413: {
414: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
415: return;
416: }
417:
418: factorisation_qr(s_etat_processus, (*s_copie_argument).objet, &tau);
419:
420: tau_complexe = (complex16 *) tau;
421:
422: if ((*s_etat_processus).erreur_systeme != d_es)
423: {
424: return;
425: }
426:
427: if (((*s_etat_processus).exception != d_ep) ||
428: ((*s_etat_processus).erreur_execution != d_ex))
429: {
430: free(tau);
431: liberation(s_etat_processus, s_objet_argument);
432: liberation(s_etat_processus, s_copie_argument);
433: return;
434: }
435:
436: if ((s_objet_resultat = copie_objet(s_etat_processus,
437: s_copie_argument, 'O')) == NULL)
438: {
439: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
440: return;
441: }
442:
443: // Matrice Q
444:
445: nombre_reflecteurs_elementaires = ((*((struct_matrice *)
446: (*s_copie_argument).objet)).nombre_colonnes <
447: (*((struct_matrice *) (*s_copie_argument).objet))
448: .nombre_lignes) ? (*((struct_matrice *)
449: (*s_copie_argument).objet)).nombre_colonnes
450: : (*((struct_matrice *) (*s_copie_argument).objet))
451: .nombre_lignes;
452:
453: registre_pile_last = NULL;
454:
455: if (test_cfsf(s_etat_processus, 31) == d_vrai)
456: {
457: registre_pile_last = (*s_etat_processus).l_base_pile_last;
458: (*s_etat_processus).l_base_pile_last = NULL;
459: }
460:
461: if ((s_objet = allocation(s_etat_processus, INT)) == NULL)
462: {
463: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
464: return;
465: }
466:
467: (*((integer8 *) (*s_objet).objet)) = (*((struct_matrice *)
468: (*s_copie_argument).objet)).nombre_lignes;
469:
470: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
471: s_objet) == d_erreur)
472: {
473: return;
474: }
475:
476: instruction_idn(s_etat_processus);
477:
478: if (((*s_etat_processus).erreur_systeme != d_es) ||
479: ((*s_etat_processus).erreur_execution != d_ex) ||
480: ((*s_etat_processus).exception != d_ep))
481: {
482: liberation(s_etat_processus, s_copie_argument);
483: free(tau);
484:
485: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
486: {
487: return;
488: }
489:
490: (*s_etat_processus).l_base_pile_last = registre_pile_last;
491: return;
492: }
493:
494: if (depilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
495: &s_matrice_identite) == d_erreur)
496: {
497: (*s_etat_processus).erreur_execution = d_ex_manque_argument;
498: return;
499: }
500:
501: for(i = 0; i < nombre_reflecteurs_elementaires; i++)
502: {
503: // Calcul de H(i) = I - tau * v * v'
504:
505: if ((s_objet = copie_objet(s_etat_processus, s_matrice_identite,
506: 'P')) == NULL)
507: {
508: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
509: return;
510: }
511:
512: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
513: s_objet) == d_erreur)
514: {
515: return;
516: }
517:
518: if ((s_objet = allocation(s_etat_processus, CPL)) == NULL)
519: {
520: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
521: return;
522: }
523:
524: (*((complex16 *) (*s_objet).objet)) = tau_complexe[i];
525:
526: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
527: s_objet) == d_erreur)
528: {
529: return;
530: }
531:
532: if ((s_objet = allocation(s_etat_processus, MCX)) == NULL)
533: {
534: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
535: return;
536: }
537:
538: (*((struct_matrice *) (*s_objet).objet)).nombre_lignes =
539: (*((struct_matrice *) (*s_copie_argument).objet))
540: .nombre_lignes;
541: (*((struct_matrice *) (*s_objet).objet)).nombre_colonnes =
542: (*((struct_matrice *) (*s_copie_argument).objet))
543: .nombre_lignes;
544:
545: if ((vecteur_complexe = malloc((*((struct_matrice *)
546: (*s_objet).objet)).nombre_lignes * sizeof(complex16)))
547: == NULL)
548: {
549: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
550: return;
551: }
552:
553: for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
554: .nombre_lignes; j++)
555: {
556: if (j < i)
557: {
558: vecteur_complexe[j].partie_reelle = 0;
559: vecteur_complexe[j].partie_imaginaire = 0;
560: }
561: else if (j == i)
562: {
563: vecteur_complexe[j].partie_reelle = 1;
564: vecteur_complexe[j].partie_imaginaire = 0;
565: }
566: else
567: {
568: vecteur_complexe[j].partie_reelle = ((complex16 **)
569: (*((struct_matrice *) (*s_copie_argument).objet))
570: .tableau)[j][i].partie_reelle;
571: vecteur_complexe[j].partie_imaginaire = ((complex16 **)
572: (*((struct_matrice *) (*s_copie_argument).objet))
573: .tableau)[j][i].partie_imaginaire;
574: }
575: }
576:
577: if (((*((struct_matrice *) (*s_objet).objet)).tableau =
578: malloc((*((struct_matrice *) (*s_objet).objet))
579: .nombre_lignes * sizeof(complex16 *))) == NULL)
580: {
581: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
582: return;
583: }
584:
585: for(j = 0; j < (*((struct_matrice *) (*s_objet).objet))
586: .nombre_lignes; j++)
587: {
588: if ((((complex16 **) (*((struct_matrice *) (*s_objet).objet))
589: .tableau)[j] = malloc((*((struct_matrice *) (*s_objet)
590: .objet)).nombre_lignes * sizeof(complex16))) == NULL)
591: {
592: (*s_etat_processus).erreur_systeme =
593: d_es_allocation_memoire;
594: return;
595: }
596:
597: for(k = 0; k < (*((struct_matrice *) (*s_objet).objet))
598: .nombre_colonnes; k++)
599: {
600: registre = vecteur_complexe[k];
601: registre.partie_imaginaire =
602: -vecteur_complexe[k].partie_imaginaire;
603:
604: f77multiplicationcc_(&(vecteur_complexe[j]),
605: ®istre, &(((complex16 **)
606: (*((struct_matrice *) (*s_objet).objet)).tableau)
607: [j][k]));
608: }
609: }
610:
611: free(vecteur_complexe);
612:
613: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
614: s_objet) == d_erreur)
615: {
616: return;
617: }
618:
619: instruction_multiplication(s_etat_processus);
620:
621: if (((*s_etat_processus).erreur_systeme != d_es) ||
622: ((*s_etat_processus).erreur_execution != d_ex) ||
623: ((*s_etat_processus).exception != d_ep))
624: {
625: liberation(s_etat_processus, s_copie_argument);
626: liberation(s_etat_processus, s_matrice_identite);
627: free(tau);
628:
629: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
630: {
631: return;
632: }
633:
634: (*s_etat_processus).l_base_pile_last = registre_pile_last;
635: return;
636: }
637:
638: instruction_moins(s_etat_processus);
639:
640: if (((*s_etat_processus).erreur_systeme != d_es) ||
641: ((*s_etat_processus).erreur_execution != d_ex) ||
642: ((*s_etat_processus).exception != d_ep))
643: {
644: liberation(s_etat_processus, s_copie_argument);
645: liberation(s_etat_processus, s_matrice_identite);
646: free(tau);
647:
648: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
649: {
650: return;
651: }
652:
653: (*s_etat_processus).l_base_pile_last = registre_pile_last;
654: return;
655: }
656:
657: if (i > 0)
658: {
659: instruction_multiplication(s_etat_processus);
660:
661: if (((*s_etat_processus).erreur_systeme != d_es) ||
662: ((*s_etat_processus).erreur_execution != d_ex) ||
663: ((*s_etat_processus).exception != d_ep))
664: {
665: liberation(s_etat_processus, s_copie_argument);
666: liberation(s_etat_processus, s_matrice_identite);
667: free(tau);
668:
669: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
670: {
671: return;
672: }
673:
674: (*s_etat_processus).l_base_pile_last = registre_pile_last;
675: return;
676: }
677: }
678: }
679:
680: if (test_cfsf(s_etat_processus, 31) == d_vrai)
681: {
682: if (empilement_pile_last(s_etat_processus, 0) == d_erreur)
683: {
684: return;
685: }
686:
687: (*s_etat_processus).l_base_pile_last = registre_pile_last;
688: }
689:
690: // Matrice R
691:
692: for(i = 0; i < (*((struct_matrice *) (*s_objet_resultat).objet))
693: .nombre_lignes; i++)
694: {
695: for(j = 0; j < i; j++)
696: {
697: ((complex16 **) (*((struct_matrice *) (*s_objet_resultat)
698: .objet)).tableau)[i][j].partie_reelle = 0;
699: ((complex16 **) (*((struct_matrice *) (*s_objet_resultat)
700: .objet)).tableau)[i][j].partie_imaginaire = 0;
701: }
702: }
703:
704: if (empilement(s_etat_processus, &((*s_etat_processus).l_base_pile),
705: s_objet_resultat) == d_erreur)
706: {
707: return;
708: }
709:
710: liberation(s_etat_processus, s_matrice_identite);
711: liberation(s_etat_processus, s_copie_argument);
712: free(tau);
713: }
714:
715: /*
716: * Type d'argument invalide
717: */
718:
719: else
720: {
721: liberation(s_etat_processus, s_objet_argument);
722:
723: (*s_etat_processus).erreur_execution = d_ex_erreur_type_argument;
724: return;
725: }
726:
727: liberation(s_etat_processus, s_objet_argument);
728:
729: return;
730: }
731:
732: // 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