1
2 /*
3 * nrutil.c
4 *
5 * Utility functions from Numerical Recipes
6 */
7
8 #include <math.h>
9 #include <stdio.h>
10 #include <stddef.h>
11 #include <stdlib.h>
12 #include <libmisc/misc.h>
13 #define NR_END 1
14 #define FREE_ARG char*
15
16
17 float sqrarg = 0;
18 int iminarg1 = 0;
19 int iminarg2 = 0;
20 int imaxarg1 = 0;
21 int imaxarg2 = 0;
22 long lminarg1 = 0;
23 long lminarg2 = 0;
24 long lmaxarg1 = 0;
25 long lmaxarg2 = 0;
26 float minarg1 = 0;
27 float minarg2 = 0;
28 float maxarg1 = 0;
29 float maxarg2 = 0;
30 double dminarg1 = 0;
31 double dminarg2 = 0;
32 double dmaxarg1 = 0;
33 double dmaxarg2 = 0;
34 double dsqrarg = 0;
35
36 int nrerror_set = 0;
37
38
39 int dummy()
40 {
41 return
42 sqrarg + iminarg1 + iminarg2 + imaxarg1 + imaxarg2 +
43 lminarg1 + lminarg2 + lmaxarg2 + lmaxarg1 + minarg1 + minarg2 +
44 maxarg1 + maxarg2 + dminarg1 + dminarg2 + dmaxarg1 + dmaxarg2 +
45 dsqrarg;
46 }
47
48 /* Numerical Recipes standard error handler */
49 void nrerror(char error_text[])
50 {
51 elog(LOG_CRIT,"Numerical Recipes run-time error: %s",
52 error_text);
53 nrerror_set = 1;
54 }
55
56 /* allocate a float vector with subscript range v[nl..nh] */
57 float *vector(long nl, long nh)
58 {
59 float *v;
60 v = (float *)malloc((size_t)((nh-nl+1+NR_END)*sizeof(float)));
61 if (!v) {
62 elog(LOG_ERR, "Memory allocation failure in vector(): %m");
63 return NULL;
64 }
65 return v-nl+NR_END;
66 }
67
68 /* allocate an int vector with subscript range v[nl..nh] */
69 int *ivector(long nl, long nh)
70 {
71 int *v;
72
73 v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
74 if (!v) nrerror("allocation failure in ivector()");
75 return v-nl+NR_END;
76 }
77
78 /* allocate an unsigned char vector with subscript range v[nl..nh] */
79 unsigned char *cvector(long nl, long nh)
80 {
81 unsigned char *v;
82
83 v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
84 if (!v) nrerror("allocation failure in cvector()");
85 return v-nl+NR_END;
86 }
87
88 /* allocate an unsigned long vector with subscript range v[nl..nh] */
89 unsigned long *lvector(long nl, long nh)
90 {
91 unsigned long *v;
92
93 v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
94 if (!v) nrerror("allocation failure in lvector()");
95 return v-nl+NR_END;
96 }
97
98 /* allocate a double vector with subscript range v[nl..nh] */
99 double *dvector(long nl, long nh)
100 {
101 double *v;
102
103 v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
104 if (!v) nrerror("allocation failure in dvector()");
105 return v-nl+NR_END;
106 }
107
108 /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
109 float **matrix(long nrl, long nrh, long ncl, long nch)
110 {
111 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
112 float **m;
113
114 /* allocate pointers to rows */
115 m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
116 if (!m) nrerror("allocation failure 1 in matrix()");
117 m += NR_END;
118 m -= nrl;
119
120 /* allocate rows and set pointers to them */
121 m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
122 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
123 m[nrl] += NR_END;
124 m[nrl] -= ncl;
125
126 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
127
128 /* return pointer to array of pointers to rows */
129 return m;
130 }
131
132 /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
133 double **dmatrix(long nrl, long nrh, long ncl, long nch)
134 {
135 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
136 double **m;
137
138 /* allocate pointers to rows */
139 m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
140 if (!m) nrerror("allocation failure 1 in matrix()");
141 m += NR_END;
142 m -= nrl;
143
144 /* allocate rows and set pointers to them */
145 m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
146 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
147 m[nrl] += NR_END;
148 m[nrl] -= ncl;
149
150 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
151
152 /* return pointer to array of pointers to rows */
153 return m;
154 }
155
156 /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
157 int **imatrix(long nrl, long nrh, long ncl, long nch)
158 {
159 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
160 int **m;
161
162 /* allocate pointers to rows */
163 m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
164 if (!m) nrerror("allocation failure 1 in matrix()");
165 m += NR_END;
166 m -= nrl;
167
168
169 /* allocate rows and set pointers to them */
170 m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
171 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
172 m[nrl] += NR_END;
173 m[nrl] -= ncl;
174
175 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
176
177 /* return pointer to array of pointers to rows */
178 return m;
179 }
180
181 /* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
182 float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
183 long newrl, long newcl)
184 {
185 long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
186 float **m;
187
188 /* allocate array of pointers to rows */
189 m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
190 if (!m) nrerror("allocation failure in submatrix()");
191 m += NR_END;
192 m -= newrl;
193
194 /* set pointers to rows */
195 for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
196
197 /* return pointer to array of pointers to rows */
198 return m;
199 }
200
201
202 /*
203 allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
204 declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
205 and ncol=nch-ncl+1. The routine should be called with the address
206 &a[0][0] as the first argument.
207 */
208 float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
209 {
210 long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
211 float **m;
212
213 /* allocate pointers to rows */
214 m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
215 if (!m) nrerror("allocation failure in convert_matrix()");
216 m += NR_END;
217 m -= nrl;
218
219 /* set pointers to rows */
220 m[nrl]=a-ncl;
221 for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
222 /* return pointer to array of pointers to rows */
223 return m;
224 }
225
226 /* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
227 float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
228 {
229 long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
230 float ***t;
231
232 /* allocate pointers to pointers to rows */
233 t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
234 if (!t) nrerror("allocation failure 1 in f3tensor()");
235 t += NR_END;
236 t -= nrl;
237
238 /* allocate pointers to rows and set pointers to them */
239 t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
240 if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
241 t[nrl] += NR_END;
242 t[nrl] -= ncl;
243
244 /* allocate rows and set pointers to them */
245 t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
246 if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
247 t[nrl][ncl] += NR_END;
248 t[nrl][ncl] -= ndl;
249
250 for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
251 for(i=nrl+1;i<=nrh;i++) {
252 t[i]=t[i-1]+ncol;
253 t[i][ncl]=t[i-1][ncl]+ncol*ndep;
254 for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
255 }
256
257 /* return pointer to array of pointers to rows */
258 return t;
259 }
260
261 /* free a float vector allocated with vector() */
262 void free_vector(float *v, long nl, long nh)
263 {
264 free((FREE_ARG) (v+nl-NR_END));
265 }
266
267 /* free an int vector allocated with ivector() */
268 void free_ivector(int *v, long nl, long nh)
269 {
270 free((FREE_ARG) (v+nl-NR_END));
271 }
272
273 /* free an unsigned char vector allocated with cvector() */
274 void free_cvector(unsigned char *v, long nl, long nh)
275 {
276 free((FREE_ARG) (v+nl-NR_END));
277 }
278
279 /* free an unsigned long vector allocated with lvector() */
280 void free_lvector(unsigned long *v, long nl, long nh)
281 {
282 free((FREE_ARG) (v+nl-NR_END));
283 }
284
285 /* free a double vector allocated with dvector() */
286 void free_dvector(double *v, long nl, long nh)
287 {
288 free((FREE_ARG) (v+nl-NR_END));
289 }
290
291 /* free a float matrix allocated by matrix() */
292 void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
293 {
294 free((FREE_ARG) (m[nrl]+ncl-NR_END));
295 free((FREE_ARG) (m+nrl-NR_END));
296 }
297
298 /* free a double matrix allocated by dmatrix() */
299 void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
300 {
301 free((FREE_ARG) (m[nrl]+ncl-NR_END));
302 free((FREE_ARG) (m+nrl-NR_END));
303 }
304
305 /* free an int matrix allocated by imatrix() */
306 void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
307 {
308 free((FREE_ARG) (m[nrl]+ncl-NR_END));
309 free((FREE_ARG) (m+nrl-NR_END));
310 }
311
312 /* free a submatrix allocated by submatrix() */
313 void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
314 {
315 free((FREE_ARG) (b+nrl-NR_END));
316 }
317
318 /* free a matrix allocated by convert_matrix() */
319 void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
320 {
321 free((FREE_ARG) (b+nrl-NR_END));
322 }
323
324 /* free a float f3tensor allocated by f3tensor() */
325 void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
326 long ndl, long ndh)
327 {
328 free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
329 free((FREE_ARG) (t[nrl]+ncl-NR_END));
330 free((FREE_ARG) (t+nrl-NR_END));
331 }
332
333
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.