source: Sophya/trunk/ArchTOIPipe/Processors/dludcmp.c@ 2160

Last change on this file since 2160 was 1944, checked in by aubourg, 24 years ago

decorrelateur wiener
File size: 1.1 KB
RevLine 
[1944]1#include <math.h>
2#define NRANSI
3#define TINY 1.0e-20;
4
5void nrerror(char error_text[]);
6double *dvector(long nl, long nh);
7void free_dvector(double *v, long nl, long nh);
8int dludcmp(double **a, int n, int *indx, double *d);
9
10int dludcmp(double **a, int n, int *indx, double *d)
11{
12 int i,imax,j,k;
13 double big,dum,sum,temp;
14 double *vv;
15
16 vv=dvector(1,n);
17 *d=1.0;
18 for (i=1;i<=n;i++) {
19 big=0.0;
20 for (j=1;j<=n;j++)
21 if ((temp=fabs(a[i][j])) > big) big=temp;
22 if (big == 0.0) return -1;
23 vv[i]=1.0/big;
24 }
25 for (j=1;j<=n;j++) {
26 for (i=1;i<j;i++) {
27 sum=a[i][j];
28 for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
29 a[i][j]=sum;
30 }
31 big=0.0;
32 for (i=j;i<=n;i++) {
33 sum=a[i][j];
34 for (k=1;k<j;k++)
35 sum -= a[i][k]*a[k][j];
36 a[i][j]=sum;
37 if ( (dum=vv[i]*fabs(sum)) >= big) {
38 big=dum;
39 imax=i;
40 }
41 }
42 if (j != imax) {
43 for (k=1;k<=n;k++) {
44 dum=a[imax][k];
45 a[imax][k]=a[j][k];
46 a[j][k]=dum;
47 }
48 *d = -(*d);
49 vv[imax]=vv[j];
50 }
51 indx[j]=imax;
52 if (a[j][j] == 0.0) a[j][j]=TINY;
53 if (j != n) {
54 dum=1.0/(a[j][j]);
55 for (i=j+1;i<=n;i++) a[i][j] *= dum;
56 }
57 }
58 free_dvector(vv,1,n);
59 return 0;
60}
61#undef TINY
62#undef NRANSI
Note: See TracBrowser for help on using the repository browser.