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