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