source: Sophya/trunk/ArchTOIPipe/Processors/nrutil.c@ 3263

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

decorrelateur wiener

File size: 16.0 KB
Line 
1#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
2
3#include <stdio.h>
4#include <stddef.h>
5#include <stdlib.h>
6#define NR_END 1
7#define FREE_ARG char*
8
9void nrerror(char error_text[]);
10
11void nrerror(char error_text[])
12/* Numerical Recipes standard error handler */
13{
14 fprintf(stderr,"Numerical Recipes run-time error...\n");
15 fprintf(stderr,"%s\n",error_text);
16 fprintf(stderr,"...now exiting to system...\n");
17 exit(1);
18}
19
20double *dvector(long nl, long nh)
21/* allocate a double vector with subscript range v[nl..nh] */
22{
23 double *v;
24
25 v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
26 if (!v) nrerror("allocation failure in vector()");
27 return v-nl+NR_END;
28}
29
30int *ivector(long nl, long nh)
31/* allocate an int vector with subscript range v[nl..nh] */
32{
33 int *v;
34
35 v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
36 if (!v) nrerror("allocation failure in ivector()");
37 return v-nl+NR_END;
38}
39
40unsigned char *cvector(long nl, long nh)
41/* allocate an unsigned char vector with subscript range v[nl..nh] */
42{
43 unsigned char *v;
44
45 v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
46 if (!v) nrerror("allocation failure in cvector()");
47 return v-nl+NR_END;
48}
49
50unsigned long *lvector(long nl, long nh)
51/* allocate an unsigned long vector with subscript range v[nl..nh] */
52{
53 unsigned long *v;
54
55 v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
56 if (!v) nrerror("allocation failure in lvector()");
57 return v-nl+NR_END;
58}
59
60#if 0
61double *dvector(long nl, long nh)
62/* allocate a double vector with subscript range v[nl..nh] */
63{
64 double *v;
65
66 v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
67 if (!v) nrerror("allocation failure in dvector()");
68 return v-nl+NR_END;
69}
70#endif
71
72double **matrix(long nrl, long nrh, long ncl, long nch)
73/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
74{
75 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
76 double **m;
77
78 /* allocate pointers to rows */
79 m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
80 if (!m) nrerror("allocation failure 1 in matrix()");
81 m += NR_END;
82 m -= nrl;
83
84 /* allocate rows and set pointers to them */
85 m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
86 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
87 m[nrl] += NR_END;
88 m[nrl] -= ncl;
89
90 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
91
92 /* return pointer to array of pointers to rows */
93 return m;
94}
95
96double **dmatrix(long nrl, long nrh, long ncl, long nch)
97/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
98{
99 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
100 double **m;
101
102 /* allocate pointers to rows */
103 m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
104 if (!m) nrerror("allocation failure 1 in matrix()");
105 m += NR_END;
106 m -= nrl;
107
108 /* allocate rows and set pointers to them */
109 m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
110 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
111 m[nrl] += NR_END;
112 m[nrl] -= ncl;
113
114 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
115
116 /* return pointer to array of pointers to rows */
117 return m;
118}
119
120int **imatrix(long nrl, long nrh, long ncl, long nch)
121/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
122{
123 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
124 int **m;
125
126 /* allocate pointers to rows */
127 m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
128 if (!m) nrerror("allocation failure 1 in matrix()");
129 m += NR_END;
130 m -= nrl;
131
132
133 /* allocate rows and set pointers to them */
134 m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
135 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
136 m[nrl] += NR_END;
137 m[nrl] -= ncl;
138
139 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
140
141 /* return pointer to array of pointers to rows */
142 return m;
143}
144
145double **submatrix(double **a, long oldrl, long oldrh, long oldcl, long oldch,
146 long newrl, long newcl)
147/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
148{
149 long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
150 double **m;
151
152 /* allocate array of pointers to rows */
153 m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*)));
154 if (!m) nrerror("allocation failure in submatrix()");
155 m += NR_END;
156 m -= newrl;
157
158 /* set pointers to rows */
159 for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
160
161 /* return pointer to array of pointers to rows */
162 return m;
163}
164
165double **convert_matrix(double *a, long nrl, long nrh, long ncl, long nch)
166/* allocate a double matrix m[nrl..nrh][ncl..nch] that points to the matrix
167declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
168and ncol=nch-ncl+1. The routine should be called with the address
169&a[0][0] as the first argument. */
170{
171 long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
172 double **m;
173
174 /* allocate pointers to rows */
175 m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*)));
176 if (!m) nrerror("allocation failure in convert_matrix()");
177 m += NR_END;
178 m -= nrl;
179
180 /* set pointers to rows */
181 m[nrl]=a-ncl;
182 for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
183 /* return pointer to array of pointers to rows */
184 return m;
185}
186
187double ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
188/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
189{
190 long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
191 double ***t;
192
193 /* allocate pointers to pointers to rows */
194 t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
195 if (!t) nrerror("allocation failure 1 in f3tensor()");
196 t += NR_END;
197 t -= nrl;
198
199 /* allocate pointers to rows and set pointers to them */
200 t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
201 if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
202 t[nrl] += NR_END;
203 t[nrl] -= ncl;
204
205 /* allocate rows and set pointers to them */
206 t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
207 if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
208 t[nrl][ncl] += NR_END;
209 t[nrl][ncl] -= ndl;
210
211 for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
212 for(i=nrl+1;i<=nrh;i++) {
213 t[i]=t[i-1]+ncol;
214 t[i][ncl]=t[i-1][ncl]+ncol*ndep;
215 for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
216 }
217
218 /* return pointer to array of pointers to rows */
219 return t;
220}
221
222void free_vector(double *v, long nl, long nh)
223/* free a double vector allocated with vector() */
224{
225 free((FREE_ARG) (v+nl-NR_END));
226}
227
228void free_ivector(int *v, long nl, long nh)
229/* free an int vector allocated with ivector() */
230{
231 free((FREE_ARG) (v+nl-NR_END));
232}
233
234void free_cvector(unsigned char *v, long nl, long nh)
235/* free an unsigned char vector allocated with cvector() */
236{
237 free((FREE_ARG) (v+nl-NR_END));
238}
239
240void free_lvector(unsigned long *v, long nl, long nh)
241/* free an unsigned long vector allocated with lvector() */
242{
243 free((FREE_ARG) (v+nl-NR_END));
244}
245
246void free_dvector(double *v, long nl, long nh)
247/* free a double vector allocated with dvector() */
248{
249 free((FREE_ARG) (v+nl-NR_END));
250}
251
252void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
253/* free a double matrix allocated by matrix() */
254{
255 free((FREE_ARG) (m[nrl]+ncl-NR_END));
256 free((FREE_ARG) (m+nrl-NR_END));
257}
258
259void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
260/* free a double matrix allocated by dmatrix() */
261{
262 free((FREE_ARG) (m[nrl]+ncl-NR_END));
263 free((FREE_ARG) (m+nrl-NR_END));
264}
265
266void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
267/* free an int matrix allocated by imatrix() */
268{
269 free((FREE_ARG) (m[nrl]+ncl-NR_END));
270 free((FREE_ARG) (m+nrl-NR_END));
271}
272
273void free_submatrix(double **b, long nrl, long nrh, long ncl, long nch)
274/* free a submatrix allocated by submatrix() */
275{
276 free((FREE_ARG) (b+nrl-NR_END));
277}
278
279void free_convert_matrix(double **b, long nrl, long nrh, long ncl, long nch)
280/* free a matrix allocated by convert_matrix() */
281{
282 free((FREE_ARG) (b+nrl-NR_END));
283}
284
285void free_f3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
286 long ndl, long ndh)
287/* free a double f3tensor allocated by f3tensor() */
288{
289 free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
290 free((FREE_ARG) (t[nrl]+ncl-NR_END));
291 free((FREE_ARG) (t+nrl-NR_END));
292}
293
294#else /* ANSI */
295/* traditional - K&R */
296
297#include <stdio.h>
298#define NR_END 1
299#define FREE_ARG char*
300
301void nrerror(error_text)
302char error_text[];
303/* Numerical Recipes standard error handler */
304{
305 void exit();
306
307 fprintf(stderr,"Numerical Recipes run-time error...\n");
308 fprintf(stderr,"%s\n",error_text);
309 fprintf(stderr,"...now exiting to system...\n");
310 exit(1);
311}
312
313double *vector(nl,nh)
314long nh,nl;
315/* allocate a double vector with subscript range v[nl..nh] */
316{
317 double *v;
318
319 v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
320 if (!v) nrerror("allocation failure in vector()");
321 return v-nl+NR_END;
322}
323
324int *ivector(nl,nh)
325long nh,nl;
326/* allocate an int vector with subscript range v[nl..nh] */
327{
328 int *v;
329
330 v=(int *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(int)));
331 if (!v) nrerror("allocation failure in ivector()");
332 return v-nl+NR_END;
333}
334
335unsigned char *cvector(nl,nh)
336long nh,nl;
337/* allocate an unsigned char vector with subscript range v[nl..nh] */
338{
339 unsigned char *v;
340
341 v=(unsigned char *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
342 if (!v) nrerror("allocation failure in cvector()");
343 return v-nl+NR_END;
344}
345
346unsigned long *lvector(nl,nh)
347long nh,nl;
348/* allocate an unsigned long vector with subscript range v[nl..nh] */
349{
350 unsigned long *v;
351
352 v=(unsigned long *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(long)));
353 if (!v) nrerror("allocation failure in lvector()");
354 return v-nl+NR_END;
355}
356
357double *dvector(nl,nh)
358long nh,nl;
359/* allocate a double vector with subscript range v[nl..nh] */
360{
361 double *v;
362
363 v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
364 if (!v) nrerror("allocation failure in dvector()");
365 return v-nl+NR_END;
366}
367
368double **matrix(nrl,nrh,ncl,nch)
369long nch,ncl,nrh,nrl;
370/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
371{
372 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
373 double **m;
374
375 /* allocate pointers to rows */
376 m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
377 if (!m) nrerror("allocation failure 1 in matrix()");
378 m += NR_END;
379 m -= nrl;
380
381 /* allocate rows and set pointers to them */
382 m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
383 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
384 m[nrl] += NR_END;
385 m[nrl] -= ncl;
386
387 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
388
389 /* return pointer to array of pointers to rows */
390 return m;
391}
392
393double **dmatrix(nrl,nrh,ncl,nch)
394long nch,ncl,nrh,nrl;
395/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
396{
397 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
398 double **m;
399
400 /* allocate pointers to rows */
401 m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
402 if (!m) nrerror("allocation failure 1 in matrix()");
403 m += NR_END;
404 m -= nrl;
405
406 /* allocate rows and set pointers to them */
407 m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
408 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
409 m[nrl] += NR_END;
410 m[nrl] -= ncl;
411
412 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
413
414 /* return pointer to array of pointers to rows */
415 return m;
416}
417
418int **imatrix(nrl,nrh,ncl,nch)
419long nch,ncl,nrh,nrl;
420/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
421{
422 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
423 int **m;
424
425 /* allocate pointers to rows */
426 m=(int **) malloc((unsigned int)((nrow+NR_END)*sizeof(int*)));
427 if (!m) nrerror("allocation failure 1 in matrix()");
428 m += NR_END;
429 m -= nrl;
430
431
432 /* allocate rows and set pointers to them */
433 m[nrl]=(int *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(int)));
434 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
435 m[nrl] += NR_END;
436 m[nrl] -= ncl;
437
438 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
439
440 /* return pointer to array of pointers to rows */
441 return m;
442}
443
444double **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
445double **a;
446long newcl,newrl,oldch,oldcl,oldrh,oldrl;
447/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
448{
449 long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
450 double **m;
451
452 /* allocate array of pointers to rows */
453 m=(double **) malloc((unsigned int) ((nrow+NR_END)*sizeof(double*)));
454 if (!m) nrerror("allocation failure in submatrix()");
455 m += NR_END;
456 m -= newrl;
457
458 /* set pointers to rows */
459 for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
460
461 /* return pointer to array of pointers to rows */
462 return m;
463}
464
465double **convert_matrix(a,nrl,nrh,ncl,nch)
466double *a;
467long nch,ncl,nrh,nrl;
468/* allocate a double matrix m[nrl..nrh][ncl..nch] that points to the matrix
469declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
470and ncol=nch-ncl+1. The routine should be called with the address
471&a[0][0] as the first argument. */
472{
473 long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
474 double **m;
475
476 /* allocate pointers to rows */
477 m=(double **) malloc((unsigned int) ((nrow+NR_END)*sizeof(double*)));
478 if (!m) nrerror("allocation failure in convert_matrix()");
479 m += NR_END;
480 m -= nrl;
481
482 /* set pointers to rows */
483 m[nrl]=a-ncl;
484 for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
485 /* return pointer to array of pointers to rows */
486 return m;
487}
488
489double ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
490long nch,ncl,ndh,ndl,nrh,nrl;
491/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
492{
493 long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
494 double ***t;
495
496 /* allocate pointers to pointers to rows */
497 t=(double ***) malloc((unsigned int)((nrow+NR_END)*sizeof(double**)));
498 if (!t) nrerror("allocation failure 1 in f3tensor()");
499 t += NR_END;
500 t -= nrl;
501
502 /* allocate pointers to rows and set pointers to them */
503 t[nrl]=(double **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double*)));
504 if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
505 t[nrl] += NR_END;
506 t[nrl] -= ncl;
507
508 /* allocate rows and set pointers to them */
509 t[nrl][ncl]=(double *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(double)));
510 if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
511 t[nrl][ncl] += NR_END;
512 t[nrl][ncl] -= ndl;
513
514 for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
515 for(i=nrl+1;i<=nrh;i++) {
516 t[i]=t[i-1]+ncol;
517 t[i][ncl]=t[i-1][ncl]+ncol*ndep;
518 for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
519 }
520
521 /* return pointer to array of pointers to rows */
522 return t;
523}
524
525void free_vector(v,nl,nh)
526double *v;
527long nh,nl;
528/* free a double vector allocated with vector() */
529{
530 free((FREE_ARG) (v+nl-NR_END));
531}
532
533void free_ivector(v,nl,nh)
534int *v;
535long nh,nl;
536/* free an int vector allocated with ivector() */
537{
538 free((FREE_ARG) (v+nl-NR_END));
539}
540
541void free_cvector(v,nl,nh)
542long nh,nl;
543unsigned char *v;
544/* free an unsigned char vector allocated with cvector() */
545{
546 free((FREE_ARG) (v+nl-NR_END));
547}
548
549void free_lvector(v,nl,nh)
550long nh,nl;
551unsigned long *v;
552/* free an unsigned long vector allocated with lvector() */
553{
554 free((FREE_ARG) (v+nl-NR_END));
555}
556
557void free_dvector(v,nl,nh)
558double *v;
559long nh,nl;
560/* free a double vector allocated with dvector() */
561{
562 free((FREE_ARG) (v+nl-NR_END));
563}
564
565void free_matrix(m,nrl,nrh,ncl,nch)
566double **m;
567long nch,ncl,nrh,nrl;
568/* free a double matrix allocated by matrix() */
569{
570 free((FREE_ARG) (m[nrl]+ncl-NR_END));
571 free((FREE_ARG) (m+nrl-NR_END));
572}
573
574void free_dmatrix(m,nrl,nrh,ncl,nch)
575double **m;
576long nch,ncl,nrh,nrl;
577/* free a double matrix allocated by dmatrix() */
578{
579 free((FREE_ARG) (m[nrl]+ncl-NR_END));
580 free((FREE_ARG) (m+nrl-NR_END));
581}
582
583void free_imatrix(m,nrl,nrh,ncl,nch)
584int **m;
585long nch,ncl,nrh,nrl;
586/* free an int matrix allocated by imatrix() */
587{
588 free((FREE_ARG) (m[nrl]+ncl-NR_END));
589 free((FREE_ARG) (m+nrl-NR_END));
590}
591
592void free_submatrix(b,nrl,nrh,ncl,nch)
593double **b;
594long nch,ncl,nrh,nrl;
595/* free a submatrix allocated by submatrix() */
596{
597 free((FREE_ARG) (b+nrl-NR_END));
598}
599
600void free_convert_matrix(b,nrl,nrh,ncl,nch)
601double **b;
602long nch,ncl,nrh,nrl;
603/* free a matrix allocated by convert_matrix() */
604{
605 free((FREE_ARG) (b+nrl-NR_END));
606}
607
608void free_f3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
609double ***t;
610long nch,ncl,ndh,ndl,nrh,nrl;
611/* free a double f3tensor allocated by f3tensor() */
612{
613 free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
614 free((FREE_ARG) (t[nrl]+ncl-NR_END));
615 free((FREE_ARG) (t+nrl-NR_END));
616}
617
618#endif /* ANSI */
Note: See TracBrowser for help on using the repository browser.