source: Sophya/trunk/Poubelle/DPC:FitsIOServer/NTools/matxop.c@ 1084

Last change on this file since 1084 was 658, checked in by ansari, 26 years ago

no message
File size: 19.5 KB
Line 
1#include <stdlib.h>
2#include <unistd.h>
3#include <stdio.h>
4#include <math.h>
5
6
7#include "matxop.h"
8#include "nbmath.h"
9
10/* Fonctions de manipulation de matrices et de vecteurs */
11/* Resolution de systemes lineaires */
12
13/* R. Ansari Juillet 1993 */
14/* LaSilla (Chili) */
15
16
17#define MXXFLOAT 1.e36
18
19
20
21/* Nouvelle-Fonction */
22
23int IVecxVec(int *v1, int *v2, int n)
24
25/* Produit scalaire de deux vecteurs entiers */
26/* Retour (int) = v1(n) . v2(n) */
27
28{
29register int i,rc;
30register int *ip1, *ip2;
31
32ip1 = v1; ip2 = v2;
33rc = 0;
34/* for(i=0; i<n; i++) rc += *(ip1+i) * *(ip2+i); */
35for(i=0; i<n; i++) rc += *ip1++ * *ip2++;
36return(rc);
37}
38
39
40/* Nouvelle-Fonction */
41
42float RVecxVec(float *v1, float *v2, int n)
43
44/* Produit scalaire de deux vecteurs reels */
45/* Retour (float) = v1(n) . v2(n) */
46
47{
48register int i;
49register float rc;
50register float *fp1, *fp2;
51
52fp1 = v1; fp2 = v2;
53rc = 0.0;
54/* for(i=0; i<n; i++) rc += *(fp1+i) * *(fp2+i); */
55for(i=0; i<n; i++) rc += *fp1++ * *fp2++;
56return(rc);
57}
58
59
60/* Nouvelle-Fonction */
61
62double DVecxVec(double *v1, double *v2, int n)
63
64/* Produit scalaire de deux vecteurs double */
65/* Retour (float) = v1(n) . v2(n) */
66
67{
68register int i;
69register double rc;
70register double *fp1, *fp2;
71
72fp1 = v1; fp2 = v2;
73rc = 0.0;
74/* for(i=0; i<n; i++) rc += *(fp1+i) * *(fp2+i); */
75for(i=0; i<n; i++) rc += *fp1++ * *fp2++;
76return(rc);
77}
78
79
80/* Nouvelle-Fonction */
81
82void IMatxVec(int *mx, int *vi, int *vo, int n)
83
84/* Matrice * Vecteur (Entiers) Integer */
85/* Vecteur Vo(n) = Matrice Mx(n,n) . Vi(n) */
86
87{
88register int *vp, *mxp;
89register int s;
90register int j,i;
91
92vp = vi; mxp = mx;
93
94for (i=0; i<n; i++)
95 { s = 0;
96 for(j=0; j<n; j++) s += *(mxp+j) * *(vp+j);
97 *(vo+i) = s; mxp += n; }
98
99return;
100}
101
102
103/* Nouvelle-Fonction */
104
105void RMatxVec(float *mx, float *vi, float *vo, int n)
106
107/* Matrice * Vecteur (Reels) float */
108/* Vecteur Vo(n) = Matrice Mx(n,n) . Vi(n) */
109
110
111/* Matrice * Vecteur (Reels) */
112{
113register float *vp, *mxp;
114register float s;
115register int j,i;
116
117vp = vi; mxp = mx;
118
119for (i=0; i<n; i++)
120 { s = 0;
121 for(j=0; j<n; j++) s += *(mxp+j) * *(vp+j);
122 *(vo+i) = s; mxp += n; }
123
124return;
125}
126
127
128/* Nouvelle-Fonction */
129
130void DMatxVec(double *mx, double *vi, double *vo, int n)
131
132/* Matrice * Vecteur (double) double */
133/* Vecteur Vo(n) = Matrice Mx(n,n) . Vi(n) */
134
135
136/* Matrice * Vecteur (Reels) */
137{
138register double *vp, *mxp;
139register double s;
140register int j,i;
141
142vp = vi; mxp = mx;
143
144for (i=0; i<n; i++)
145 { s = 0;
146 for(j=0; j<n; j++) s += *(mxp+j) * *(vp+j);
147 *(vo+i) = s; mxp += n; }
148
149return;
150}
151
152
153
154
155/* Nouvelle-Fonction */
156
157float SolveRLinSyst(float *mx, float *b, float *x, int n)
158
159/* Resolution de systeme lineaire */
160/* Matrice Mx(n,n) * Vecteur X(n) = Vecteur B(n) */
161/* Inconnu : vecteur X */
162/* Retour = Determinant du systeme (=0 Pb) et X[] */
163/* Note : Au retour le vecteur B(n) et la matrice Mx(n,n) */
164/* sont modifies. Mx(n,n) contient une matrice triangulaire */
165/* superieure */
166
167{
168int i,k;
169register int j;
170register float *fp1, *fp2;
171double det;
172float vt1;
173register float vt2;
174register float six;
175int pivok;
176
177
178#define MINPIVOT 1.e-15
179#define MINDETER 1.e-15
180
181
182
183for (i=0; i<n-1; i++) /* Boucle sur les pivots */
184 {
185
186/* printf("\n Iteration %d : \n",i);
187 for (k=0; k<n; k++)
188 {
189 for (j=0; j<n; j++) printf(" %10g ",*(mx+k*n+j));
190 printf(" b= %10g \n",*(b+k));
191 } */
192
193 fp1 = fp2 = mx+i*n;
194 vt1 = *(fp1+i) ;
195
196 if ((vt1 < MINPIVOT) && (vt1 > -MINPIVOT) ) /* Pivot trop petit */
197 {
198 pivok = 0;
199 for (k=i+1; k<n; k++)
200 {
201 vt1 = *(mx+k*n+i) ;
202 if ((vt1 > MINPIVOT) || (vt1 < -MINPIVOT) )
203 {
204 fp2 = mx+k*n;
205 for (j=i; j<n; j++)
206 {
207 *(fp1+j) += *(fp2+j);
208 }
209 *(b+i) += *(b+k);
210 vt1 = *(fp1+i) ;
211 fp2 = fp1;
212 pivok = 1;
213 break;
214 }
215 }
216 if (!pivok) return(0.0);
217 }
218
219 for (k=i+1; k<n; k++)
220 {
221 fp1 += n;
222 vt2 = *(fp1+i) / vt1 ;
223 for (j=i+1; j<n; j++)
224 *(fp1+j) -= (vt2 * *(fp2+j));
225
226 *(fp1+i) = 0.0;
227 *(b+k) -= (vt2 * *(b+i)) ;
228 }
229
230 }
231
232/* Calcul du determinant */
233det = 1.0;
234for (i=0; i<n; i++) det *= *(mx+i*n+i);
235if ((det < MINDETER) && (det > -MINDETER) ) return(0.0);
236if (det > MXXFLOAT) det = MXXFLOAT;
237if (det < -MXXFLOAT) det = -MXXFLOAT;
238
239
240for(i=n-1; i>=0; i--)
241 {
242 fp1 = mx+i*n;
243 six = *(b+i);
244 for (j=i+1; j<n; j++) six -= (x[j] * *(fp1+j) );
245 x[i] = six / (*(fp1+i));
246
247/* printf(" Solution X[%2d] = %g \n",i,x[i]); */
248 }
249
250return(det);
251}
252
253
254
255
256
257/* Nouvelle-Fonction */
258
259double SolveDLinSyst(double *mx, double *b, double *x, int n)
260
261/* Resolution de systeme lineaire */
262/* Matrice Mx(n,n) * Vecteur X(n) = Vecteur B(n) */
263/* Inconnu : vecteur X */
264/* Retour = Determinant du systeme (=0 Pb) et X[] */
265/* Note : Au retour le vecteur B(n) et la matrice Mx(n,n) */
266/* sont modifies. Mx(n,n) contient une matrice triangulaire */
267/* superieure */
268
269{
270int i,k;
271register int j;
272register double *fp1, *fp2;
273double det;
274double vt1;
275register double vt2;
276register double six;
277int pivok;
278
279
280#define DMINPIVOT 1.e-25
281#define DMINDETER 1.e-30
282
283
284
285for (i=0; i<n-1; i++) /* Boucle sur les pivots */
286 {
287
288/* printf("\n Iteration %d : \n",i);
289 for (k=0; k<n; k++)
290 {
291 for (j=0; j<n; j++) printf(" %10g ",*(mx+k*n+j));
292 printf(" b= %10g \n",*(b+k));
293 } */
294
295 fp1 = fp2 = mx+i*n;
296 vt1 = *(fp1+i) ;
297
298 if ((vt1 < DMINPIVOT) && (vt1 > -DMINPIVOT) ) /* Pivot trop petit */
299 {
300 pivok = 0;
301 for (k=i+1; k<n; k++)
302 {
303 vt1 = *(mx+k*n+i) ;
304 if ((vt1 > DMINPIVOT) || (vt1 < -DMINPIVOT) )
305 {
306 fp2 = mx+k*n;
307 for (j=i; j<n; j++)
308 {
309 *(fp1+j) += *(fp2+j);
310 }
311 *(b+i) += *(b+k);
312 vt1 = *(fp1+i) ;
313 fp2 = fp1;
314 pivok = 1;
315 break;
316 }
317 }
318 if (!pivok) return(0.0);
319 }
320
321 for (k=i+1; k<n; k++)
322 {
323 fp1 += n;
324 vt2 = *(fp1+i) / vt1 ;
325 for (j=i+1; j<n; j++)
326 *(fp1+j) -= (vt2 * *(fp2+j));
327
328 *(fp1+i) = 0.0;
329 *(b+k) -= (vt2 * *(b+i)) ;
330 }
331
332 }
333
334/* Calcul du determinant */
335det = 1.0;
336for (i=0; i<n; i++) det *= *(mx+i*n+i);
337if ((det < DMINDETER) && (det > -DMINDETER) ) return(0.0);
338if (det > MXXFLOAT) det = MXXFLOAT;
339if (det < -MXXFLOAT) det = -MXXFLOAT;
340
341
342for(i=n-1; i>=0; i--)
343 {
344 fp1 = mx+i*n;
345 six = *(b+i);
346 for (j=i+1; j<n; j++) six -= (x[j] * *(fp1+j) );
347 x[i] = six / (*(fp1+i));
348
349/* printf(" Solution X[%2d] = %g \n",i,x[i]); */
350 }
351
352return(det);
353}
354
355
356
357
358/* ============================================================== */
359/* Fonctions de Fit lineaire de Xi2 */
360/* ============================================================== */
361
362/* .......... Fit avec vecteurs depart reels ............... */
363
364static int FFBusy = -1; /* Flag Deja Busy (GetFitVect effectue) */
365static int FNVarMax = 0; /* Nb Maxi de variables a ajuster */
366static int FVLenMax = 0; /* Longueur maxi des vecteurs */
367static float **FVeci = NULL; /* Vecteurs en entree *Vecf = Y */
368 /* *(Vecf+i) = X[i] */
369static float **FVecf = NULL; /* Vecteurs *FVecf = *FVecf+2 = B = M.X */
370 /* *(FVecf+1) = La sortie = Ai */
371static float *FVSpace = NULL; /* Espace pour les *FVeci */
372static float *FVSpacef = NULL; /* " " " *FVecf */
373static double *FFMtx = NULL; /* Matrice a inverser = M */
374static double *FFMtx2 = NULL; /* Copie de M */
375static double *FFSort = NULL; /* RFitLinErr() pour Sortie GausPiv() */
376 /* FNVarMax Lignes * (FNVarMax+1) Col. */
377 /* Col1= Solution, Col2..N+1 : MatInv */
378
379
380int InitRFitLin(int nvarmx, int szv)
381
382/* Initialisation du fit de Xi2 Entier */
383/* Allocation d'espace memoire */
384/* nvarmx = Nb maxi de variable */
385/* szv = Taille maxi des vecteurs de points */
386
387{
388int i;
389
390FFBusy = -1;
391nvarmx ++;
392if (nvarmx < 7) nvarmx = 7;
393if (szv < 25) szv = 25;
394
395if ( (FVeci = malloc(nvarmx*sizeof(float *))) == NULL )
396 { printf("InitRFitLin_Erreur: (Pb malloc(Veci)) \n");
397 return(1);
398 }
399
400if ( (FVSpace = malloc(nvarmx*szv*sizeof(float))) == NULL )
401 { printf("InitRFitLin_Erreur: (Pb malloc(VSpace)) \n");
402 free(FVeci);
403 return(2);
404 }
405
406for(i=0; i<nvarmx; i++) *(FVeci+i) = FVSpace+i*szv;
407
408
409if ( (FVecf = malloc(3*sizeof(float *))) == NULL )
410 { printf("InitRFitLin_Erreur: (Pb malloc(Vecf)) \n");
411 free(FVeci);
412 free(FVSpace);
413 return(3);
414 }
415
416if ( (FVSpacef = malloc(3*szv*sizeof(float))) == NULL )
417 { printf("InitRFitLin_Erreur: (Pb malloc(VSpacef)) \n");
418 free(FVeci);
419 free(FVSpace);
420 free(FVecf);
421 return(4);
422 }
423
424for(i=0; i<3; i++) *(FVecf+i) = FVSpacef+i*szv;
425
426
427if ( (FFMtx = malloc(nvarmx*nvarmx*sizeof(double))) == NULL )
428 { printf("InitRFitLin_Erreur: (Pb malloc(FMtx)) \n");
429 free(FVeci);
430 free(FVSpace);
431 free(FVecf);
432 free(FVSpacef);
433 return(5);
434 }
435
436if ( (FFMtx2 = malloc(nvarmx*nvarmx*sizeof(double))) == NULL )
437 { printf("InitRFitLin_Erreur: (Pb malloc(FMtx2)) \n");
438 free(FVeci);
439 free(FVSpace);
440 free(FVecf);
441 free(FVSpacef);
442 free(FFMtx);
443 return(6);
444 }
445
446if ( (FFSort = malloc((nvarmx+1)*nvarmx*sizeof(double))) == NULL )
447 { printf("InitRFitLin_Erreur: (Pb malloc(FFSort)) \n");
448 free(FVeci);
449 free(FVSpace);
450 free(FVecf);
451 free(FVSpacef);
452 free(FFMtx);
453 free(FFMtx2);
454 return(7);
455 }
456
457FFBusy = 0;
458FNVarMax = nvarmx-1;
459FVLenMax = szv;
460
461return(0);
462}
463
464
465/* Nouvelle-Fonction */
466void EndRFitLin()
467{
468
469if (FVeci != NULL) free(FVeci);
470if (FVSpace != NULL) free(FVSpace);
471if (FVecf != NULL) free(FVecf);
472if (FVSpacef != NULL) free(FVSpacef);
473if (FFMtx != NULL) free(FFMtx);
474if (FFMtx2 != NULL) free(FFMtx2);
475if (FFSort != NULL) free(FFSort);
476
477FVeci = NULL;
478FVSpace = NULL;
479FVecf = NULL;
480FVSpacef = NULL;
481FFMtx = NULL;
482FFMtx2 = NULL;
483
484FFBusy = -1;
485FNVarMax = 0;
486FVLenMax = 0;
487return;
488}
489
490
491/* Nouvelle-Fonction */
492
493float **GetRFitVect(int nvar, int vsz)
494/* reservation de l'espace des vecteurs pour le fit */
495{
496int nvmx, vszmx;
497
498if (FFBusy > 0)
499 { printf("GetRFitVect_Erreur: FFBusy= %d \n",FFBusy);
500 return(NULL); }
501if ((nvar < 1) || (vsz < 1) )
502 { printf("GetRFitVect_Erreur: NVar,VSz= %d %d \n",nvar,vsz);
503 return(NULL); }
504if ((nvar > FNVarMax) || (vsz > FVLenMax) || (FFBusy < 0))
505 {
506 EndRFitLin();
507 if (nvar > FNVarMax) nvmx = nvar;
508 else nvmx = FNVarMax;
509 if (vsz > FVLenMax) vszmx = vsz;
510 else vszmx = FVLenMax;
511 if (InitRFitLin(nvmx, vszmx) != 0) return(NULL);
512 }
513
514FFBusy = 1;
515return(FVeci);
516}
517
518
519/* Nouvelle-Fonction */
520void FreeRFitVect()
521{
522FFBusy = 0;
523return;
524}
525
526/* Nouvelle-Fonction */
527
528float RFitLin(int nd, int vsz, float *XVar)
529
530/* Effectue le fit de Xi2 pour nd inconnu, Taille vecteurs vsz */
531/* Les vecteurs doivent etre remplis auparavant */
532/* Retourne la valeur du Xi2 du fit, Negatif si Pb. */
533
534{
535int i,j;
536register float *ip;
537register float *fp, *fp2;
538register double x;
539register float y;
540register float *ffmtx, *ffmtx2;
541float det, rc;
542
543
544if ((nd > FNVarMax) || (vsz > FVLenMax) ) {rc = -1000.0; goto Fin; }
545
546/* Fabrication matrice du syteme lineaire a resoudre */
547ffmtx = (float *)FFMtx;
548ffmtx2 = (float *)FFMtx2;
549for(i=0; i<nd; i++)
550 {
551 fp = ffmtx+i*nd; fp2 = ffmtx2+i*nd;
552 ip = *(FVeci+i+1);
553
554/* Elements diagonals */
555 *(fp2+i) = *(fp+i) = RVecxVec( ip, ip, vsz) ;
556
557 for(j=i+1; j<nd; j++)
558 { /* Les autres elements */
559 *(fp2+j) = *(fp+j) = RVecxVec( ip, *(FVeci+j+1), vsz) ;
560 /* Matrice symetrique */
561 *(ffmtx2+j*nd+i) = *(ffmtx+j*nd+i) = *(fp+j);
562 }
563
564/* Second membre M.X = B */
565 fp = *FVecf; fp2 = *(FVecf+2);
566 *(fp2+i) = *(fp+i) = RVecxVec( ip, *FVeci, vsz) ;
567
568 }
569
570/* On resoud le systeme lineaire : */
571det = SolveRLinSyst(ffmtx, *FVecf, *(FVecf+1), nd);
572if (det == 0.0) {rc = -500.0; goto Fin; }
573
574/* On calcule le Xi2 */
575x = RVecxVec( *FVeci, *FVeci, vsz);
576fp = *(FVecf+2);
577for(i=0; i<nd; i++)
578 {
579 XVar[i] = y = *(*(FVecf+1)+i);
580 fp2 = ffmtx2+i*nd+i;
581 x += (y * y * *fp2);
582 x -= ( 2.0 * ( y * *(fp+i) ) );
583 for (j=i+1; j<nd; j++)
584 { fp2++;
585 x += 2.0 * y * *(*(FVecf+1)+j) * *fp2 ; }
586 }
587
588rc = (float)x;
589
590Fin:
591FreeRFitVect();
592return(rc);
593}
594
595
596
597/* Nouvelle-Fonction */
598
599float RFitLinErr(int nd, int vsz, float *XVar, float *Err)
600
601/* Effectue le fit de Xi2 pour nd inconnu, Taille vecteurs vsz */
602/* Les vecteurs doivent etre remplis auparavant */
603/* Retourne la valeur du Xi2 du fit, Negatif si Pb. */
604
605{
606int i,j,md;
607register float *ip;
608register double *dp, *dp2;
609register float *fp;
610register double x,y;
611float rc;
612double det;
613
614
615if ((nd > FNVarMax) || (vsz > FVLenMax) ) {rc = -1000.0; goto Fin; }
616
617md = nd+1;
618/* Mise a zero de la matrice de sortie et creation matrice identite */
619for(i=0; i<nd*md; i++) FFSort[i] = 0.;
620for(i=1; i<nd*md; i += md+1) FFSort[i] = 1.;
621
622/* Fabrication matrice du syteme lineaire a resoudre */
623
624for(i=0; i<nd; i++)
625 {
626 dp = FFMtx+i*nd; dp2 = FFMtx2+i*nd;
627 ip = *(FVeci+i+1);
628
629/* Elements diagonals */
630 *(dp+i) = *(dp2+i) = RVecxVec( ip, ip, vsz) ;
631
632 for(j=i+1; j<nd; j++)
633 { /* Les autres elements */
634 *(dp+j) = *(dp2+j) = RVecxVec( ip, *(FVeci+j+1), vsz) ;
635 /* Matrice symetrique */
636 *(FFMtx+j*nd+i) = *(FFMtx2+j*nd+i) = *(dp+j);
637 }
638
639/* Second membre M.X = B */
640 dp = FFSort+i*md; fp = *(FVecf+2);
641 *dp = *(fp+i) = RVecxVec( ip, *FVeci, vsz) ;
642
643 }
644
645
646/* Appel GausPiv pour resoudre le systeme et inverser la matrice */
647det = GausPiv(FFMtx, nd, nd, FFSort, md, md, 0);
648/* printf("Det= %g Resul= %g %g %g \n",det,FFSort[0], FFSort[md], FFSort[2*md]); */
649if (det == 0.0) {rc = -500.0; goto Fin; }
650
651/* On calcule le Xi2 Et remplissage vecteurs sortie */
652
653x = RVecxVec( *FVeci, *FVeci, vsz);
654fp = *(FVecf+2);
655for(i=0; i<nd; i++)
656 {
657 y = FFSort[1+i*(md+1)];
658 Err[i] = (y >= 0.) ? (float)sqrt(y) : -(float)sqrt(-y);
659 XVar[i]= y =FFSort[i*md];
660
661 dp2 = FFMtx2+i*nd+i;
662 x += (y * y * *dp2);
663 x -= ( 2.0 * ( y * (double)*(fp+i) ) );
664 for (j=i+1; j<nd; j++)
665 { dp2++;
666 x += 2.0 * y * FFSort[j*md] * *dp2 ; }
667 }
668rc = (float)x;
669
670Fin:
671FreeRFitVect();
672return(rc);
673}
674
675
676
677
678/* .......... Fit avec vecteurs depart double ............... */
679
680/* Voir description des variables ds IniRFitLin() */
681static int DFBusy = -1;
682static int DNVarMax = 0;
683static int DVLenMax = 0;
684static double **DVeci = NULL;
685static double **DVecf = NULL;
686static double *DVSpace = NULL;
687static double *DVSpacef = NULL;
688static double *DFMtx = NULL;
689static double *DFMtx2 = NULL;
690
691
692int InitDFitLin(int nvarmx, int szv)
693
694/* Initialisation du fit de Xi2 Entier */
695/* Allocation d'espace memoire */
696/* nvarmx = Nb maxi de variable */
697/* szv = Taille maxi des vecteurs de points */
698
699{
700int i;
701
702DFBusy = -1;
703nvarmx ++;
704if (nvarmx < 7) nvarmx = 7;
705if (szv < 25) szv = 25;
706
707if ( (DVeci = malloc(nvarmx*sizeof(double *))) == NULL )
708 { printf("InitDFitLin_Erreur: (Pb malloc(Veci)) \n");
709 return(1);
710 }
711
712if ( (DVSpace = malloc(nvarmx*szv*sizeof(double))) == NULL )
713 { printf("InitDFitLin_Erreur: (Pb malloc(VSpace)) \n");
714 free(DVeci);
715 return(2);
716 }
717
718for(i=0; i<nvarmx; i++) *(DVeci+i) = DVSpace+i*szv;
719
720
721if ( (DVecf = malloc(3*sizeof(double *))) == NULL )
722 { printf("InitDFitLin_Erreur: (Pb malloc(Vecf)) \n");
723 free(DVeci);
724 free(DVSpace);
725 return(3);
726 }
727
728if ( (DVSpacef = malloc(3*szv*sizeof(double))) == NULL )
729 { printf("InitDFitLin_Erreur: (Pb malloc(VSpacef)) \n");
730 free(DVeci);
731 free(DVSpace);
732 free(DVecf);
733 return(4);
734 }
735
736for(i=0; i<3; i++) *(DVecf+i) = DVSpacef+i*szv;
737
738
739if ( (DFMtx = malloc(nvarmx*nvarmx*sizeof(double))) == NULL )
740 { printf("InitDFitLin_Erreur: (Pb malloc(FMtx)) \n");
741 free(DVeci);
742 free(DVSpace);
743 free(DVecf);
744 free(DVSpacef);
745 return(5);
746 }
747
748if ( (DFMtx2 = malloc(nvarmx*nvarmx*sizeof(double))) == NULL )
749 { printf("InitDFitLin_Erreur: (Pb malloc(FMtx2)) \n");
750 free(DVeci);
751 free(DVSpace);
752 free(DVecf);
753 free(DVSpacef);
754 free(DFMtx);
755 return(6);
756 }
757
758DFBusy = 0;
759DNVarMax = nvarmx-1;
760DVLenMax = szv;
761
762return(0);
763}
764
765
766/* Nouvelle-Fonction */
767void EndDFitLin()
768{
769
770if (DVeci != NULL) free(DVeci);
771if (DVSpace != NULL) free(DVSpace);
772if (DVecf != NULL) free(DVecf);
773if (DVSpacef != NULL) free(DVSpacef);
774if (DFMtx != NULL) free(DFMtx);
775if (DFMtx2 != NULL) free(DFMtx2);
776
777DVeci = NULL;
778DVSpace = NULL;
779DVecf = NULL;
780DVSpacef = NULL;
781DFMtx = DFMtx2 = NULL;
782
783DFBusy = -1;
784DNVarMax = 0;
785DVLenMax = 0;
786return;
787}
788
789
790/* Nouvelle-Fonction */
791
792double **GetDFitVect(int nvar, int vsz)
793/* reservation de l'espace des vecteurs pour le fit */
794{
795int nvmx, vszmx;
796
797if (DFBusy > 0) return(NULL);
798if ( (nvar < 1) || (vsz < 1) ) return(NULL);
799if ((nvar > DNVarMax) || (vsz > DVLenMax) || (DFBusy < 0))
800 {
801 EndDFitLin();
802 if (nvar > DNVarMax) nvmx = nvar;
803 else nvmx = DNVarMax;
804 if (vsz > DVLenMax) vszmx = vsz;
805 else vszmx = DVLenMax;
806 InitDFitLin(nvmx, vszmx);
807 }
808
809DFBusy = 1;
810return(DVeci);
811}
812
813
814/* Nouvelle-Fonction */
815void FreeDFitVect()
816{
817DFBusy = 0;
818return;
819}
820
821
822/* Nouvelle-Fonction */
823
824double DFitLin(int nd, int vsz, double *XVar)
825
826/* Effectue le fit de Xi2 pour nd inconnu, Taille vecteurs vsz */
827/* Les vecteurs doivent etre remplis auparavant */
828/* Retourne la valeur du Xi2 du fit, Negatif si Pb. */
829
830{
831int i,j;
832register double *ip;
833register double *fp, *fp2;
834register double x,y;
835double det,rc;
836
837
838if ((nd > DNVarMax) || (vsz > DVLenMax) ) { rc = -1000.; goto Fin; }
839
840/* Impression de debug */
841/*
842for(i=0; i<nd+1; i++)
843 {
844 ip = *(DVeci+i);
845 printf("DFitLin_Debug Veci[%d][0..5]= %lg %lg %lg %lg %lg %lg \n",
846 i,*ip, *(ip+1), *(ip+2), *(ip+3), *(ip+4), *(ip+5));
847 }
848*/
849
850
851/* Fabrication matrice du syteme lineaire a resoudre */
852
853for(i=0; i<nd; i++)
854 {
855 fp = DFMtx+i*nd; fp2 = DFMtx2+i*nd;
856 ip = *(DVeci+i+1);
857
858/* Elements diagonals */
859 *(fp2+i) = *(fp+i) = DVecxVec( ip, ip, vsz) ;
860
861 for(j=i+1; j<nd; j++)
862 { /* Les autres elements */
863 *(fp2+j) = *(fp+j) = DVecxVec( ip, *(DVeci+j+1), vsz) ;
864 /* Matrice symetrique */
865 *(DFMtx2+j*nd+i) = *(DFMtx+j*nd+i) = *(fp+j);
866 }
867
868/* Second membre M.X = B */
869 fp = *DVecf; fp2 = *(DVecf+2);
870 *(fp2+i) = *(fp+i) = DVecxVec( ip, *DVeci, vsz) ;
871
872 }
873
874/*
875printf("DFitLin_Debug B[0..5]= %lg %lg %lg %lg %lg %lg \n",
876 *fp2, *(fp2+1), *(fp2+2), *(fp2+3), *(fp2+4), *(fp2+5));
877*/
878
879/* On resoud le systeme lineaire : */
880fp = *(DVecf+1);
881det = SolveDLinSyst(DFMtx, *DVecf, fp, nd);
882if (det == 0.0) { rc = -500.; goto Fin; }
883
884/* On calcule le Xi2 */
885x = DVecxVec( *DVeci, *DVeci, vsz);
886for(i=0; i<nd; i++) /* Et on remplit le vecteur de sortie */
887 {
888 fp2 = DFMtx2+i*nd+i; y = XVar[i] = *(*(DVecf+1)+i);
889 x += (y * y * *fp2);
890 x -= ( 2.0 * ( y * *(*(DVecf+2)+i) ) );
891 for (j=i+1; j<nd; j++)
892 { fp2++;
893 x += ( 2.0 * ( y * *(*(DVecf+1)+j) * *fp2 ) ); }
894 }
895rc = x;
896
897
898Fin:
899FreeDFitVect();
900return(rc);
901}
902
903
904
Note: See TracBrowser for help on using the repository browser.