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 calculant le nombre de condition d'une matrice
29: ================================================================================
30: Entrées : struct_matrice
31: --------------------------------------------------------------------------------
32: Sorties : nombre de condition de la matrice
33: --------------------------------------------------------------------------------
34: Effets de bord : néant
35: ================================================================================
36: */
37:
38: static integer4
39: calcul_cond(struct_processus *s_etat_processus, void *matrice_f77,
40: integer4 nombre_lignes_a, integer4 nombre_colonnes_a,
41: integer4 *pivot, unsigned char type, real8 *cond)
42: {
43: integer4 erreur;
44: integer4 *iwork;
45: integer4 longueur;
46: integer4 ordre;
47:
48: real8 anorme;
49: real8 rcond;
50: real8 *rwork;
51:
52: unsigned char norme;
53:
54: void *work;
55:
56: norme = '1';
57: longueur = 1;
58:
59: if (type == 'R')
60: {
61: // work est NULL dans le cas de la norme '1'
62: anorme = dlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a,
63: matrice_f77, &nombre_lignes_a, NULL, longueur);
64:
65: dgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77,
66: &nombre_lignes_a, pivot, &erreur);
67:
68: if (erreur < 0)
69: {
70: (*s_etat_processus).erreur_execution =
71: d_ex_routines_mathematiques;
72:
73: free(matrice_f77);
74: return(-1);
75: }
76:
77: if ((iwork = malloc(((size_t) nombre_colonnes_a) *
78: sizeof(integer4))) == NULL)
79: {
80: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
81: return(-1);
82: }
83:
84: if ((work = malloc(4 * ((size_t) nombre_colonnes_a) *
85: sizeof(real8))) == NULL)
86: {
87: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
88: return(-1);
89: }
90:
91: ordre = (nombre_lignes_a > nombre_colonnes_a)
92: ? nombre_colonnes_a : nombre_lignes_a;
93:
94: dgecon_(&norme, &ordre, matrice_f77,
95: &nombre_lignes_a, &anorme, &rcond, work, iwork, &erreur,
96: longueur);
97:
98: free(work);
99: free(iwork);
100:
101: if (erreur < 0)
102: {
103: (*s_etat_processus).erreur_execution =
104: d_ex_routines_mathematiques;
105:
106: free(matrice_f77);
107: return(-1);
108: }
109: }
110: else
111: {
112: // work est NULL dans le cas de la norme '1'
113: anorme = zlange_(&norme, &nombre_lignes_a, &nombre_colonnes_a,
114: matrice_f77, &nombre_lignes_a, NULL, longueur);
115:
116: zgetrf_(&nombre_lignes_a, &nombre_colonnes_a, matrice_f77,
117: &nombre_lignes_a, pivot, &erreur);
118:
119: if (erreur < 0)
120: {
121: (*s_etat_processus).erreur_execution =
122: d_ex_routines_mathematiques;
123:
124: free(matrice_f77);
125: return(-1);
126: }
127:
128: if ((rwork = malloc(2 * ((size_t) nombre_colonnes_a) * sizeof(real8)))
129: == NULL)
130: {
131: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
132: return(-1);
133: }
134:
135: if ((work = malloc(2 * ((size_t) nombre_colonnes_a) *
136: sizeof(complex16))) == NULL)
137: {
138: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
139: return(-1);
140: }
141:
142: ordre = (nombre_lignes_a > nombre_colonnes_a)
143: ? nombre_colonnes_a : nombre_lignes_a;
144:
145: zgecon_(&norme, &ordre, matrice_f77,
146: &nombre_lignes_a, &anorme, &rcond, work, rwork, &erreur,
147: longueur);
148:
149: free(work);
150: free(rwork);
151:
152: if (erreur < 0)
153: {
154: (*s_etat_processus).erreur_execution =
155: d_ex_routines_mathematiques;
156:
157: free(matrice_f77);
158: return(-1);
159: }
160: }
161:
162: (*cond) = ((real8) 1 / rcond);
163: return(0);
164: }
165:
166:
167: void
168: cond(struct_processus *s_etat_processus,
169: struct_matrice *s_matrice, real8 *condition)
170: {
171: integer4 dimension_vecteur_pivot;
172: integer4 nombre_lignes_a;
173: integer4 nombre_colonnes_a;
174: integer4 *pivot;
175: integer4 rang;
176: integer4 taille_matrice_f77;
177:
178: real8 cond;
179:
180: integer8 i;
181: integer8 j;
182: integer8 k;
183:
184: void *matrice_f77;
185:
186: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
187: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
188: dimension_vecteur_pivot = (nombre_lignes_a < nombre_colonnes_a)
189: ? nombre_lignes_a : nombre_colonnes_a;
190: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;
191:
192: if ((pivot = (integer4 *) malloc(((size_t) dimension_vecteur_pivot) *
193: sizeof(integer4))) == NULL)
194: {
195: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
196: return;
197: }
198:
199: switch((*s_matrice).type)
200: {
201: case 'I' :
202: {
203: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
204: sizeof(real8))) == NULL)
205: {
206: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
207: return;
208: }
209:
210: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
211: {
212: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
213: {
214: ((real8 *) matrice_f77)[k++] = (real8) ((integer8 **)
215: (*s_matrice).tableau)[j][i];
216: }
217: }
218:
219: if ((rang = calcul_cond(s_etat_processus, matrice_f77,
220: nombre_lignes_a, nombre_colonnes_a, pivot,
221: 'R', &cond)) < 0)
222: {
223: free(pivot);
224: free(matrice_f77);
225: return;
226: }
227:
228: free(matrice_f77);
229: (*condition) = cond;
230: break;
231: }
232:
233: case 'R' :
234: {
235: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
236: sizeof(real8))) == NULL)
237: {
238: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
239: return;
240: }
241:
242: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
243: {
244: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
245: {
246: ((real8 *) matrice_f77)[k++] = ((real8 **)
247: (*s_matrice).tableau)[j][i];
248: }
249: }
250:
251: if ((rang = calcul_cond(s_etat_processus, matrice_f77,
252: nombre_lignes_a, nombre_colonnes_a, pivot,
253: 'R', &cond)) < 0)
254: {
255: free(pivot);
256: free(matrice_f77);
257: return;
258: }
259:
260: free(matrice_f77);
261: (*condition) = cond;
262: break;
263: }
264:
265: case 'C' :
266: {
267: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
268: sizeof(complex16))) == NULL)
269: {
270: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
271: return;
272: }
273:
274: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
275: {
276: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
277: {
278: ((complex16 *) matrice_f77)[k++] = ((complex16 **)
279: (*s_matrice).tableau)[j][i];
280: }
281: }
282:
283: if ((rang = calcul_cond(s_etat_processus, matrice_f77,
284: nombre_lignes_a, nombre_colonnes_a, pivot,
285: 'C', &cond)) < 0)
286: {
287: free(pivot);
288: free(matrice_f77);
289: return;
290: }
291:
292: free(matrice_f77);
293: (*condition) = cond;
294: break;
295: }
296: }
297:
298: free(pivot);
299:
300: return;
301: }
302:
303:
304: /*
305: ================================================================================
306: Fonction effectuant une décomposition en valeurs singulières
307: ================================================================================
308: Entrées : struct_matrice
309: --------------------------------------------------------------------------------
310: Sorties : valeurs singulières dans le vecteur S. Si les pointeurs sur U
311: et VH ne sont pas nul, les matrices U et VH sont aussi calculées.
312: --------------------------------------------------------------------------------
313: Effets de bord : néant
314: ================================================================================
315: */
316:
317: void valeurs_singulieres(struct_processus *s_etat_processus,
318: struct_matrice *s_matrice, struct_matrice *matrice_u,
319: struct_vecteur *vecteur_s, struct_matrice *matrice_vh)
320: {
321: integer4 erreur;
322: integer4 longueur;
323: integer4 lwork;
324: integer4 nombre_colonnes_a;
325: integer4 nombre_lignes_a;
326: integer4 nombre_valeurs_singulieres;
327: integer4 taille_matrice_f77;
328:
329: integer8 i;
330: integer8 j;
331: integer8 k;
332:
333: real8 *rwork;
334:
335: unsigned char jobu;
336: unsigned char jobvh;
337:
338: void *matrice_f77;
339: void *matrice_f77_u;
340: void *matrice_f77_vh;
341: void *vecteur_f77_s;
342: void *work;
343:
344: longueur = 1;
345:
346: if (matrice_u != NULL)
347: {
348: jobu = 'A';
349: }
350: else
351: {
352: jobu = 'N';
353: }
354:
355: if (matrice_vh != NULL)
356: {
357: jobvh = 'A';
358: }
359: else
360: {
361: jobvh = 'N';
362: }
363:
364: nombre_lignes_a = (integer4) (*s_matrice).nombre_lignes;
365: nombre_colonnes_a = (integer4) (*s_matrice).nombre_colonnes;
366: taille_matrice_f77 = nombre_lignes_a * nombre_colonnes_a;
367: nombre_valeurs_singulieres = (nombre_lignes_a < nombre_colonnes_a)
368: ? nombre_lignes_a : nombre_colonnes_a;
369:
370: switch((*s_matrice).type)
371: {
372: case 'I' :
373: {
374: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
375: sizeof(real8))) == NULL)
376: {
377: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
378: return;
379: }
380:
381: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
382: {
383: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
384: {
385: ((real8 *) matrice_f77)[k++] = (real8) ((integer8 **)
386: (*s_matrice).tableau)[j][i];
387: }
388: }
389:
390: lwork = -1;
391:
392: if ((work = malloc(sizeof(real8))) == NULL)
393: {
394: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
395: return;
396: }
397:
398: if (matrice_u != NULL)
399: {
400: if ((matrice_f77_u = malloc(((size_t) (nombre_lignes_a *
401: nombre_lignes_a)) * sizeof(real8))) == NULL)
402: {
403: (*s_etat_processus).erreur_systeme =
404: d_es_allocation_memoire;
405: return;
406: }
407: }
408: else
409: {
410: matrice_f77_u = NULL;
411: }
412:
413: if ((vecteur_f77_s = malloc(((size_t) nombre_valeurs_singulieres) *
414: sizeof(real8))) == NULL)
415: {
416: (*s_etat_processus).erreur_systeme =
417: d_es_allocation_memoire;
418: return;
419: }
420:
421: if (matrice_vh != NULL)
422: {
423: if ((matrice_f77_vh = malloc(((size_t) (nombre_colonnes_a
424: * nombre_colonnes_a)) * sizeof(real8))) == NULL)
425: {
426: (*s_etat_processus).erreur_systeme =
427: d_es_allocation_memoire;
428: return;
429: }
430: }
431: else
432: {
433: matrice_f77_vh = NULL;
434: }
435:
436: dgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a,
437: matrice_f77, &nombre_lignes_a, vecteur_f77_s,
438: matrice_f77_u, &nombre_lignes_a,
439: matrice_f77_vh, &nombre_colonnes_a,
440: work, &lwork, &erreur, longueur, longueur);
441:
442: lwork = (integer4) ((real8 *) work)[0];
443: free(work);
444:
445: if ((work = malloc(((size_t) lwork) * sizeof(real8))) == NULL)
446: {
447: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
448: return;
449: }
450:
451: dgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a,
452: matrice_f77, &nombre_lignes_a, vecteur_f77_s,
453: matrice_f77_u, &nombre_lignes_a,
454: matrice_f77_vh, &nombre_colonnes_a,
455: work, &lwork, &erreur, longueur, longueur);
456:
457: free(work);
458: free(matrice_f77);
459:
460: if (erreur != 0)
461: {
462: if (erreur > 0)
463: {
464: (*s_etat_processus).exception = d_ep_decomposition_SVD;
465: }
466: else
467: {
468: (*s_etat_processus).erreur_execution =
469: d_ex_routines_mathematiques;
470: }
471:
472: free(matrice_f77_u);
473: free(matrice_f77_vh);
474: free(vecteur_f77_s);
475: return;
476: }
477:
478: if (matrice_u != NULL)
479: {
480: (*matrice_u).nombre_lignes = nombre_lignes_a;
481: (*matrice_u).nombre_colonnes = nombre_lignes_a;
482:
483: if (((*matrice_u).tableau = malloc(((size_t)
484: (*matrice_u).nombre_lignes) * sizeof(real8 *))) == NULL)
485: {
486: (*s_etat_processus).erreur_systeme =
487: d_es_allocation_memoire;
488: return;
489: }
490:
491: for(i = 0; i < (*matrice_u).nombre_lignes; i++)
492: {
493: if ((((real8 **) (*matrice_u).tableau)[i] =
494: malloc(((size_t) (*matrice_u).nombre_colonnes) *
495: sizeof(real8))) == NULL)
496: {
497: (*s_etat_processus).erreur_systeme =
498: d_es_allocation_memoire;
499: return;
500: }
501: }
502:
503: for(k = 0, i = 0; i < (*matrice_u).nombre_colonnes; i++)
504: {
505: for(j = 0; j < (*matrice_u).nombre_lignes; j++)
506: {
507: ((real8 **) (*matrice_u).tableau)[j][i] =
508: ((real8 *) matrice_f77_u)[k++];
509: }
510: }
511:
512: free(matrice_f77_u);
513: }
514:
515: if (matrice_vh != NULL)
516: {
517: (*matrice_vh).nombre_lignes = nombre_colonnes_a;
518: (*matrice_vh).nombre_colonnes = nombre_colonnes_a;
519:
520: if (((*matrice_vh).tableau = malloc(((size_t) (*matrice_vh)
521: .nombre_lignes) * sizeof(real8 *))) == NULL)
522: {
523: (*s_etat_processus).erreur_systeme =
524: d_es_allocation_memoire;
525: return;
526: }
527:
528: for(i = 0; i < (*matrice_vh).nombre_lignes; i++)
529: {
530: if ((((real8 **) (*matrice_vh).tableau)[i] =
531: malloc(((size_t) (*matrice_vh).nombre_colonnes) *
532: sizeof(real8))) == NULL)
533: {
534: (*s_etat_processus).erreur_systeme =
535: d_es_allocation_memoire;
536: return;
537: }
538: }
539:
540: for(k = 0, i = 0; i < (*matrice_vh).nombre_colonnes; i++)
541: {
542: for(j = 0; j < (*matrice_vh).nombre_lignes; j++)
543: {
544: ((real8 **) (*matrice_vh).tableau)[j][i] =
545: ((real8 *) matrice_f77_vh)[k++];
546: }
547: }
548:
549: free(matrice_f77_vh);
550: }
551:
552: (*vecteur_s).taille = nombre_valeurs_singulieres;
553: (*vecteur_s).type = 'R';
554: (*vecteur_s).tableau = vecteur_f77_s;
555:
556: break;
557: }
558:
559: case 'R' :
560: {
561: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
562: sizeof(real8))) == NULL)
563: {
564: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
565: return;
566: }
567:
568: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
569: {
570: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
571: {
572: ((real8 *) matrice_f77)[k++] = ((real8 **)
573: (*s_matrice).tableau)[j][i];
574: }
575: }
576:
577: lwork = -1;
578:
579: if ((work = malloc(sizeof(real8))) == NULL)
580: {
581: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
582: return;
583: }
584:
585: if (matrice_u != NULL)
586: {
587: if ((matrice_f77_u = malloc(((size_t) (nombre_lignes_a *
588: nombre_lignes_a)) * sizeof(real8))) == NULL)
589: {
590: (*s_etat_processus).erreur_systeme =
591: d_es_allocation_memoire;
592: return;
593: }
594: }
595: else
596: {
597: matrice_f77_u = NULL;
598: }
599:
600: if ((vecteur_f77_s = malloc(((size_t) nombre_valeurs_singulieres) *
601: sizeof(real8))) == NULL)
602: {
603: (*s_etat_processus).erreur_systeme =
604: d_es_allocation_memoire;
605: return;
606: }
607:
608: if (matrice_vh != NULL)
609: {
610: if ((matrice_f77_vh = malloc(((size_t) (nombre_colonnes_a
611: * nombre_colonnes_a)) * sizeof(real8))) == NULL)
612: {
613: (*s_etat_processus).erreur_systeme =
614: d_es_allocation_memoire;
615: return;
616: }
617: }
618: else
619: {
620: matrice_f77_vh = NULL;
621: }
622:
623: dgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a,
624: matrice_f77, &nombre_lignes_a, vecteur_f77_s,
625: matrice_f77_u, &nombre_lignes_a,
626: matrice_f77_vh, &nombre_colonnes_a,
627: work, &lwork, &erreur, longueur, longueur);
628:
629: lwork = (integer4) ((real8 *) work)[0];
630: free(work);
631:
632: if ((work = malloc(((size_t) lwork) * sizeof(real8))) == NULL)
633: {
634: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
635: return;
636: }
637:
638: dgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a,
639: matrice_f77, &nombre_lignes_a, vecteur_f77_s,
640: matrice_f77_u, &nombre_lignes_a,
641: matrice_f77_vh, &nombre_colonnes_a,
642: work, &lwork, &erreur, longueur, longueur);
643:
644: free(work);
645: free(matrice_f77);
646:
647: if (erreur != 0)
648: {
649: if (erreur > 0)
650: {
651: (*s_etat_processus).exception = d_ep_decomposition_SVD;
652: }
653: else
654: {
655: (*s_etat_processus).erreur_execution =
656: d_ex_routines_mathematiques;
657: }
658:
659: free(matrice_f77_u);
660: free(matrice_f77_vh);
661: free(vecteur_f77_s);
662: return;
663: }
664:
665: if (matrice_u != NULL)
666: {
667: (*matrice_u).nombre_lignes = nombre_lignes_a;
668: (*matrice_u).nombre_colonnes = nombre_lignes_a;
669:
670: if (((*matrice_u).tableau = malloc(((size_t)
671: (*matrice_u).nombre_lignes) * sizeof(real8 *))) == NULL)
672: {
673: (*s_etat_processus).erreur_systeme =
674: d_es_allocation_memoire;
675: return;
676: }
677:
678: for(i = 0; i < (*matrice_u).nombre_lignes; i++)
679: {
680: if ((((real8 **) (*matrice_u).tableau)[i] =
681: malloc(((size_t) (*matrice_u).nombre_colonnes) *
682: sizeof(real8))) == NULL)
683: {
684: (*s_etat_processus).erreur_systeme =
685: d_es_allocation_memoire;
686: return;
687: }
688: }
689:
690: for(k = 0, i = 0; i < (*matrice_u).nombre_colonnes; i++)
691: {
692: for(j = 0; j < (*matrice_u).nombre_lignes; j++)
693: {
694: ((real8 **) (*matrice_u).tableau)[j][i] =
695: ((real8 *) matrice_f77_u)[k++];
696: }
697: }
698:
699: free(matrice_f77_u);
700: }
701:
702: if (matrice_vh != NULL)
703: {
704: (*matrice_vh).nombre_lignes = nombre_colonnes_a;
705: (*matrice_vh).nombre_colonnes = nombre_colonnes_a;
706:
707: if (((*matrice_vh).tableau = malloc(((size_t) (*matrice_vh)
708: .nombre_lignes) * sizeof(real8 *))) == NULL)
709: {
710: (*s_etat_processus).erreur_systeme =
711: d_es_allocation_memoire;
712: return;
713: }
714:
715: for(i = 0; i < (*matrice_vh).nombre_lignes; i++)
716: {
717: if ((((real8 **) (*matrice_vh).tableau)[i] =
718: malloc(((size_t) (*matrice_vh).nombre_colonnes) *
719: sizeof(real8))) == NULL)
720: {
721: (*s_etat_processus).erreur_systeme =
722: d_es_allocation_memoire;
723: return;
724: }
725: }
726:
727: for(k = 0, i = 0; i < (*matrice_vh).nombre_colonnes; i++)
728: {
729: for(j = 0; j < (*matrice_vh).nombre_lignes; j++)
730: {
731: ((real8 **) (*matrice_vh).tableau)[j][i] =
732: ((real8 *) matrice_f77_vh)[k++];
733: }
734: }
735:
736: free(matrice_f77_vh);
737: }
738:
739: (*vecteur_s).taille = nombre_valeurs_singulieres;
740: (*vecteur_s).type = 'R';
741: (*vecteur_s).tableau = vecteur_f77_s;
742:
743: break;
744: }
745:
746: case 'C' :
747: {
748: if ((matrice_f77 = malloc(((size_t) taille_matrice_f77) *
749: sizeof(complex16))) == NULL)
750: {
751: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
752: return;
753: }
754:
755: for(k = 0, i = 0; i < (*s_matrice).nombre_colonnes; i++)
756: {
757: for(j = 0; j < (*s_matrice).nombre_lignes; j++)
758: {
759: ((complex16 *) matrice_f77)[k++] = ((complex16 **)
760: (*s_matrice).tableau)[j][i];
761: }
762: }
763:
764: lwork = -1;
765:
766: if ((work = malloc(sizeof(complex16))) == NULL)
767: {
768: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
769: return;
770: }
771:
772: if (matrice_u != NULL)
773: {
774: if ((matrice_f77_u = malloc(((size_t) (nombre_lignes_a *
775: nombre_lignes_a)) * sizeof(complex16))) == NULL)
776: {
777: (*s_etat_processus).erreur_systeme =
778: d_es_allocation_memoire;
779: return;
780: }
781: }
782: else
783: {
784: matrice_f77_u = NULL;
785: }
786:
787: if ((vecteur_f77_s = malloc(((size_t) nombre_valeurs_singulieres) *
788: sizeof(real8))) == NULL)
789: {
790: (*s_etat_processus).erreur_systeme =
791: d_es_allocation_memoire;
792: return;
793: }
794:
795: if (matrice_vh != NULL)
796: {
797: if ((matrice_f77_vh = malloc(((size_t) (nombre_colonnes_a
798: * nombre_colonnes_a)) * sizeof(complex16))) == NULL)
799: {
800: (*s_etat_processus).erreur_systeme =
801: d_es_allocation_memoire;
802: return;
803: }
804: }
805: else
806: {
807: matrice_f77_vh = NULL;
808: }
809:
810: if ((rwork = malloc(5 * ((size_t) nombre_valeurs_singulieres)
811: * sizeof(real8))) == NULL)
812: {
813: (*s_etat_processus).erreur_systeme =
814: d_es_allocation_memoire;
815: return;
816: }
817:
818: zgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a,
819: matrice_f77, &nombre_lignes_a, vecteur_f77_s,
820: matrice_f77_u, &nombre_lignes_a,
821: matrice_f77_vh, &nombre_colonnes_a,
822: work, &lwork, rwork, &erreur, longueur, longueur);
823:
824: lwork = (integer4) ((real8 *) work)[0];
825: free(work);
826:
827: if ((work = malloc(((size_t) lwork) * sizeof(real8))) == NULL)
828: {
829: (*s_etat_processus).erreur_systeme = d_es_allocation_memoire;
830: return;
831: }
832:
833: zgesvd_(&jobu, &jobvh, &nombre_lignes_a, &nombre_colonnes_a,
834: matrice_f77, &nombre_lignes_a, vecteur_f77_s,
835: matrice_f77_u, &nombre_lignes_a,
836: matrice_f77_vh, &nombre_colonnes_a,
837: work, &lwork, rwork, &erreur, longueur, longueur);
838:
839: free(work);
840: free(rwork);
841: free(matrice_f77);
842:
843: if (erreur != 0)
844: {
845: if (erreur > 0)
846: {
847: (*s_etat_processus).exception = d_ep_decomposition_SVD;
848: }
849: else
850: {
851: (*s_etat_processus).erreur_execution =
852: d_ex_routines_mathematiques;
853: }
854:
855: free(matrice_f77_u);
856: free(matrice_f77_vh);
857: free(vecteur_f77_s);
858: return;
859: }
860:
861: if (matrice_u != NULL)
862: {
863: (*matrice_u).nombre_lignes = nombre_lignes_a;
864: (*matrice_u).nombre_colonnes = nombre_lignes_a;
865:
866: if (((*matrice_u).tableau = malloc(((size_t)
867: (*matrice_u).nombre_lignes) * sizeof(complex16 *)))
868: == NULL)
869: {
870: (*s_etat_processus).erreur_systeme =
871: d_es_allocation_memoire;
872: return;
873: }
874:
875: for(i = 0; i < (*matrice_u).nombre_lignes; i++)
876: {
877: if ((((complex16 **) (*matrice_u).tableau)[i] =
878: malloc(((size_t) (*matrice_u).nombre_colonnes) *
879: sizeof(complex16))) == NULL)
880: {
881: (*s_etat_processus).erreur_systeme =
882: d_es_allocation_memoire;
883: return;
884: }
885: }
886:
887: for(k = 0, i = 0; i < (*matrice_u).nombre_colonnes; i++)
888: {
889: for(j = 0; j < (*matrice_u).nombre_lignes; j++)
890: {
891: ((complex16 **) (*matrice_u).tableau)[j][i] =
892: ((complex16 *) matrice_f77_u)[k++];
893: }
894: }
895:
896: free(matrice_f77_u);
897: }
898:
899: if (matrice_vh != NULL)
900: {
901: (*matrice_vh).nombre_lignes = nombre_colonnes_a;
902: (*matrice_vh).nombre_colonnes = nombre_colonnes_a;
903:
904: if (((*matrice_vh).tableau = malloc(((size_t) (*matrice_vh)
905: .nombre_lignes) * sizeof(complex16 *))) == NULL)
906: {
907: (*s_etat_processus).erreur_systeme =
908: d_es_allocation_memoire;
909: return;
910: }
911:
912: for(i = 0; i < (*matrice_vh).nombre_lignes; i++)
913: {
914: if ((((complex16 **) (*matrice_vh).tableau)[i] =
915: malloc(((size_t) (*matrice_vh).nombre_colonnes) *
916: sizeof(complex16))) == NULL)
917: {
918: (*s_etat_processus).erreur_systeme =
919: d_es_allocation_memoire;
920: return;
921: }
922: }
923:
924: for(k = 0, i = 0; i < (*matrice_vh).nombre_colonnes; i++)
925: {
926: for(j = 0; j < (*matrice_vh).nombre_lignes; j++)
927: {
928: ((complex16 **) (*matrice_vh).tableau)[j][i] =
929: ((complex16 *) matrice_f77_vh)[k++];
930: }
931: }
932:
933: free(matrice_f77_vh);
934: }
935:
936: (*vecteur_s).taille = nombre_valeurs_singulieres;
937: (*vecteur_s).type = 'R';
938: (*vecteur_s).tableau = vecteur_f77_s;
939:
940: break;
941: }
942: }
943:
944: return;
945: }
946:
947: // vim: ts=4
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to algebre_lineaire4.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / src