1 | /* nrlin.c **************************************************************** |
---|
2 | |
---|
3 | NR Library Package - Definitions and Linear Algebra Functions |
---|
4 | |
---|
5 | James Trevelyan, University of Western Australia |
---|
6 | Revision 2 January 1996 |
---|
7 | |
---|
8 | **************************************************************************/ |
---|
9 | |
---|
10 | |
---|
11 | #include <stdio.h> |
---|
12 | #include <stddef.h> |
---|
13 | #include <stdlib.h> |
---|
14 | #include <math.h> |
---|
15 | #ifdef __TURBOC__ |
---|
16 | #include <conio.h> |
---|
17 | #endif |
---|
18 | |
---|
19 | #include "nrlin.h" |
---|
20 | |
---|
21 | /* Validity checking - use V_CHECK */ |
---|
22 | #define V_CHECK |
---|
23 | |
---|
24 | /* Keep track of memory used - only for dmatrix and dvector calls at the moment */ |
---|
25 | |
---|
26 | static long mem_t = 0; |
---|
27 | static int mats = 0; |
---|
28 | static int dmats = 0; |
---|
29 | static int vecs = 0; |
---|
30 | static int dvecs = 0; |
---|
31 | |
---|
32 | void using_mem( long space ) |
---|
33 | { |
---|
34 | mem_t += space; |
---|
35 | } |
---|
36 | |
---|
37 | long mem_used( void ) |
---|
38 | { |
---|
39 | return (mem_t); |
---|
40 | } |
---|
41 | |
---|
42 | void reportmemory( FILE *outfile ) |
---|
43 | { |
---|
44 | fprintf( outfile, "Memory: %ld bytes\n", mem_t); |
---|
45 | fprintf( outfile, "%d dmatrices, %d dvectors, %d matrices, %d vectors\n", |
---|
46 | dmats, dvecs, mats, vecs ); |
---|
47 | } |
---|
48 | void choldc(float **a, int n, float p[]) |
---|
49 | { |
---|
50 | nrerror("Obtain \"choldc\" source from NR diskette"); |
---|
51 | } |
---|
52 | |
---|
53 | |
---|
54 | void cholsl(float **a, int n, float p[], float b[], float x[]) |
---|
55 | { |
---|
56 | nrerror("Obtain \"cholsl\" source from NR diskette"); |
---|
57 | } |
---|
58 | |
---|
59 | |
---|
60 | void lubksb(float **a, int n, int *indx, float b[]) |
---|
61 | { |
---|
62 | nrerror("Obtain \"lubksb\" source from NR diskette"); |
---|
63 | } |
---|
64 | |
---|
65 | |
---|
66 | void ludcmp(float **a, int n, int *indx, float *d) |
---|
67 | { |
---|
68 | nrerror("Obtain \"ludcmp\" source from NR diskette"); |
---|
69 | } |
---|
70 | |
---|
71 | void dcholdc(double **a, int n, double p[]) |
---|
72 | { |
---|
73 | nrerror("Obtain \"choldc\" source from NR diskette"); |
---|
74 | } |
---|
75 | |
---|
76 | |
---|
77 | void dcholsl(double **a, int n, double p[], double b[], double x[]) |
---|
78 | { |
---|
79 | nrerror("Obtain \"cholsl\" source from NR diskette"); |
---|
80 | } |
---|
81 | |
---|
82 | |
---|
83 | void dlubksb(double **a, int n, int *indx, double b[]) |
---|
84 | { |
---|
85 | nrerror("Obtain \"lubksb\" source from NR diskette"); |
---|
86 | } |
---|
87 | |
---|
88 | |
---|
89 | void dludcmp(double **a, int n, int *indx, double *d) |
---|
90 | { |
---|
91 | nrerror("Obtain \"ludcmp\" source from NR diskette"); |
---|
92 | } |
---|
93 | |
---|
94 | |
---|
95 | |
---|
96 | |
---|
97 | #define NR_END 1 |
---|
98 | #define NR_TEST 1 |
---|
99 | #define FREE_ARG char* |
---|
100 | static int alt_handler_defined = 0; |
---|
101 | static void (*alt_error_handler)(char error_text[]); |
---|
102 | |
---|
103 | void nrerror(char error_text[]) |
---|
104 | /* Numerical Recipes standard error handler */ |
---|
105 | { |
---|
106 | if (alt_handler_defined) { |
---|
107 | (*alt_error_handler)(error_text); |
---|
108 | } |
---|
109 | else { |
---|
110 | fprintf(stderr,"Numerical Recipes run-time error...\n"); |
---|
111 | fprintf(stderr,"%s\n",error_text); |
---|
112 | fprintf(stderr,"Press <ENTER> to return to system..."); |
---|
113 | getchar(); |
---|
114 | exit(1); |
---|
115 | } |
---|
116 | } |
---|
117 | |
---|
118 | void nrerror_handler( void(*handler)(char error_text[]) ) |
---|
119 | { |
---|
120 | alt_handler_defined = 1; |
---|
121 | alt_error_handler = handler; |
---|
122 | } |
---|
123 | |
---|
124 | float *vector(I_ARG_T nl, I_ARG_T nh) |
---|
125 | /* allocate a float vector with subscript range v[nl..nh] */ |
---|
126 | { |
---|
127 | float *v; |
---|
128 | |
---|
129 | v=(float *)malloc((size_t) ((nh-nl+1+NR_END+NR_TEST)*sizeof(float))); |
---|
130 | using_mem( (nh-nl+1+NR_END+NR_TEST)*sizeof(float) ); |
---|
131 | if (!v) nrerror("allocation failure in dvector()"); |
---|
132 | v -= nl; |
---|
133 | v += NR_END; |
---|
134 | v[nl-1] = -322.0; |
---|
135 | v[nh+1] = -722.0; |
---|
136 | vecs++; |
---|
137 | return (v); |
---|
138 | } |
---|
139 | |
---|
140 | int *ivector(I_ARG_T nl, I_ARG_T nh) |
---|
141 | /* allocate an int vector with subscript range v[nl..nh] */ |
---|
142 | { |
---|
143 | int *v; |
---|
144 | |
---|
145 | v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); |
---|
146 | if (!v) nrerror("allocation failure in ivector()"); |
---|
147 | return v-nl+NR_END; |
---|
148 | } |
---|
149 | |
---|
150 | unsigned char *cvector(I_ARG_T nl, I_ARG_T nh) |
---|
151 | /* allocate an unsigned char vector with subscript range v[nl..nh] */ |
---|
152 | { |
---|
153 | unsigned char *v; |
---|
154 | |
---|
155 | v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char))); |
---|
156 | if (!v) nrerror("allocation failure in cvector()"); |
---|
157 | return v-nl+NR_END; |
---|
158 | } |
---|
159 | |
---|
160 | long *lvector(I_ARG_T nl, I_ARG_T nh) |
---|
161 | /* allocate an unsigned long vector with subscript range v[nl..nh] */ |
---|
162 | { |
---|
163 | long *v; |
---|
164 | |
---|
165 | v=(long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long))); |
---|
166 | if (!v) nrerror("allocation failure in lvector()"); |
---|
167 | return v-nl+NR_END; |
---|
168 | } |
---|
169 | |
---|
170 | void flag_dvector( double *v, I_ARG_T nh) |
---|
171 | { |
---|
172 | v[0] = -322.0; |
---|
173 | v[nh+1] = -722.0; |
---|
174 | } |
---|
175 | |
---|
176 | double *dvector(I_ARG_T nl, I_ARG_T nh) |
---|
177 | /* allocate a double vector with subscript range v[nl..nh] */ |
---|
178 | { |
---|
179 | double *v; |
---|
180 | |
---|
181 | v=(double *)malloc((size_t) ((nh-nl+1+NR_END+NR_TEST)*sizeof(double))); |
---|
182 | using_mem( (nh-nl+1+NR_END+NR_TEST)*sizeof(double) ); |
---|
183 | if (!v) nrerror("allocation failure in dvector()"); |
---|
184 | v -= nl; |
---|
185 | v += NR_END; |
---|
186 | v[nl-1] = -322.0; |
---|
187 | v[nh+1] = -722.0; |
---|
188 | dvecs++; |
---|
189 | return (v); |
---|
190 | } |
---|
191 | |
---|
192 | float **matrix(I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
193 | /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ |
---|
194 | { |
---|
195 | I_ARG_T i, nrow=nrh-nrl+1,ncol=nch-ncl+1; |
---|
196 | float **m; |
---|
197 | |
---|
198 | /* allocate pointers to rows */ |
---|
199 | m=(float **) malloc((size_t)((nrow+NR_END+1)*sizeof(float*))); |
---|
200 | using_mem( (nrow+NR_END+1)*sizeof(float*) ); |
---|
201 | if (!m) nrerror("allocation failure 1 in matrix()"); |
---|
202 | m += NR_END + 1; |
---|
203 | m -= nrl; |
---|
204 | |
---|
205 | /* allocate rows and set pointers to them */ |
---|
206 | m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END+NR_TEST)*sizeof(float))); |
---|
207 | using_mem( (nrow*ncol+NR_END+NR_TEST)*sizeof(float) ); |
---|
208 | if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); |
---|
209 | m[nrl] += NR_END; |
---|
210 | m[nrl] -= ncl; |
---|
211 | m[nrl-2] = (float *)0x555; |
---|
212 | m[nrl-1] = m[nrl]; /* covers m[0][0] mistake */ |
---|
213 | for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; |
---|
214 | m[nrl][0]=-422.0; |
---|
215 | m[nrh][nch+1]=-822.0; |
---|
216 | |
---|
217 | /* return pointer to array of pointers to rows */ |
---|
218 | mats++; |
---|
219 | return m; |
---|
220 | } |
---|
221 | |
---|
222 | /****************************************************************************/ |
---|
223 | /* double **dmatrix(I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
224 | |
---|
225 | Purpose: |
---|
226 | allocate memory for a double matrix |
---|
227 | with subscript range m[nrl..nrh][ncl..nch] |
---|
228 | |
---|
229 | ****************************************************************************/ |
---|
230 | double **dmatrix(I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
231 | { |
---|
232 | I_ARG_T i, nrow=nrh-nrl+1,ncol=nch-ncl+1; |
---|
233 | double **m; |
---|
234 | |
---|
235 | /* allocate pointers to rows */ |
---|
236 | m=(double **) malloc((size_t)((nrow+NR_END+1)*sizeof(double*))); |
---|
237 | using_mem( (nrow+NR_END+1)*sizeof(double*) ); |
---|
238 | if (!m) nrerror("allocation failure 1 in matrix()"); |
---|
239 | m += NR_END + 1; |
---|
240 | m -= nrl; |
---|
241 | |
---|
242 | /* allocate rows and set pointers to them */ |
---|
243 | m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END+NR_TEST)*sizeof(double))); |
---|
244 | using_mem( (nrow*ncol+NR_END+NR_TEST)*sizeof(double) ); |
---|
245 | if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); |
---|
246 | m[nrl] += NR_END; |
---|
247 | m[nrl] -= ncl; |
---|
248 | m[nrl-2] = (double *)0x555; |
---|
249 | m[nrl-1] = m[nrl]; /* covers m[0][0] mistake */ |
---|
250 | for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; |
---|
251 | m[nrl][0]=-422.0; |
---|
252 | m[nrh][nch+1]=-822.0; |
---|
253 | /* return pointer to array of pointers to rows */ |
---|
254 | dmats++; |
---|
255 | return m; |
---|
256 | } |
---|
257 | |
---|
258 | int **imatrix(I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
259 | /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ |
---|
260 | { |
---|
261 | I_ARG_T i, nrow=nrh-nrl+1,ncol=nch-ncl+1; |
---|
262 | int **m; |
---|
263 | |
---|
264 | /* allocate pointers to rows */ |
---|
265 | m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*))); |
---|
266 | if (!m) nrerror("allocation failure 1 in matrix()"); |
---|
267 | m += NR_END; |
---|
268 | m -= nrl; |
---|
269 | |
---|
270 | |
---|
271 | /* allocate rows and set pointers to them */ |
---|
272 | m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END+NR_TEST)*sizeof(int))); |
---|
273 | if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); |
---|
274 | m[nrl] += NR_END; |
---|
275 | m[nrl] -= ncl; |
---|
276 | |
---|
277 | for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; |
---|
278 | |
---|
279 | /* return pointer to array of pointers to rows */ |
---|
280 | return m; |
---|
281 | } |
---|
282 | |
---|
283 | float **submatrix(float **a, I_ARG_T oldrl, I_ARG_T oldrh, I_ARG_T oldcl, I_ARG_T oldch, |
---|
284 | I_ARG_T newrl, I_ARG_T newcl) |
---|
285 | /* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */ |
---|
286 | { |
---|
287 | I_ARG_T i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl; |
---|
288 | float **m; |
---|
289 | |
---|
290 | /* allocate array of pointers to rows */ |
---|
291 | m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); |
---|
292 | if (!m) nrerror("allocation failure in submatrix()"); |
---|
293 | m += NR_END; |
---|
294 | m -= newrl; |
---|
295 | |
---|
296 | /* set pointers to rows */ |
---|
297 | for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol; |
---|
298 | |
---|
299 | /* return pointer to array of pointers to rows */ |
---|
300 | return m; |
---|
301 | } |
---|
302 | |
---|
303 | float **convert_matrix(float *a, I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
304 | /* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix |
---|
305 | declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 |
---|
306 | and ncol=nch-ncl+1. The routine should be called with the address |
---|
307 | &a[0][0] as the first argument. */ |
---|
308 | { |
---|
309 | I_ARG_T i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1; |
---|
310 | float **m; |
---|
311 | |
---|
312 | /* allocate pointers to rows */ |
---|
313 | m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); |
---|
314 | if (!m) nrerror("allocation failure in convert_matrix()"); |
---|
315 | m += NR_END; |
---|
316 | m -= nrl; |
---|
317 | |
---|
318 | /* set pointers to rows */ |
---|
319 | m[nrl]=a-ncl; |
---|
320 | for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol; |
---|
321 | /* return pointer to array of pointers to rows */ |
---|
322 | return m; |
---|
323 | } |
---|
324 | |
---|
325 | float ***f3tensor(I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch, I_ARG_T ndl, I_ARG_T ndh) |
---|
326 | /* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */ |
---|
327 | { |
---|
328 | I_ARG_T i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1; |
---|
329 | float ***t; |
---|
330 | |
---|
331 | /* allocate pointers to pointers to rows */ |
---|
332 | t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**))); |
---|
333 | if (!t) nrerror("allocation failure 1 in f3tensor()"); |
---|
334 | t += NR_END; |
---|
335 | t -= nrl; |
---|
336 | |
---|
337 | /* allocate pointers to rows and set pointers to them */ |
---|
338 | t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*))); |
---|
339 | if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()"); |
---|
340 | t[nrl] += NR_END; |
---|
341 | t[nrl] -= ncl; |
---|
342 | |
---|
343 | /* allocate rows and set pointers to them */ |
---|
344 | t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float))); |
---|
345 | if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()"); |
---|
346 | t[nrl][ncl] += NR_END; |
---|
347 | t[nrl][ncl] -= ndl; |
---|
348 | |
---|
349 | for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep; |
---|
350 | for(i=nrl+1;i<=nrh;i++) { |
---|
351 | t[i]=t[i-1]+ncol; |
---|
352 | t[i][ncl]=t[i-1][ncl]+ncol*ndep; |
---|
353 | for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep; |
---|
354 | } |
---|
355 | |
---|
356 | /* return pointer to array of pointers to rows */ |
---|
357 | return t; |
---|
358 | } |
---|
359 | |
---|
360 | void free_vector(float *v, I_ARG_T nl, I_ARG_T nh) |
---|
361 | /* free a float vector allocated with vector() */ |
---|
362 | { |
---|
363 | if ( valid_vector( v, nl, nh ) ) { |
---|
364 | free((FREE_ARG) (v+nl-NR_END)); |
---|
365 | using_mem( -(nh-nl+1+NR_END+NR_TEST)*sizeof(float) ); |
---|
366 | vecs--; |
---|
367 | } |
---|
368 | else |
---|
369 | nrerror("Invalid vector pointer: free_vector"); |
---|
370 | } |
---|
371 | |
---|
372 | void free_ivector(int *v, I_ARG_T nl, I_ARG_T nh) |
---|
373 | /* free an int vector allocated with ivector() */ |
---|
374 | { |
---|
375 | free((FREE_ARG) (v+nl-NR_END)); |
---|
376 | } |
---|
377 | |
---|
378 | void free_cvector(unsigned char *v, I_ARG_T nl, I_ARG_T nh) |
---|
379 | /* free an unsigned char vector allocated with cvector() */ |
---|
380 | { |
---|
381 | free((FREE_ARG) (v+nl-NR_END)); |
---|
382 | } |
---|
383 | |
---|
384 | void free_lvector(long *v, I_ARG_T nl, I_ARG_T nh) |
---|
385 | /* free an unsigned long vector allocated with lvector() */ |
---|
386 | { |
---|
387 | free((FREE_ARG) (v+nl-NR_END)); |
---|
388 | } |
---|
389 | |
---|
390 | void free_dvector(double *v, I_ARG_T nl, I_ARG_T nh) |
---|
391 | /* free a double vector allocated with dvector() */ |
---|
392 | { |
---|
393 | if ( valid_dvector( v, nl, nh ) ) { |
---|
394 | free((FREE_ARG) (v+nl-NR_END)); |
---|
395 | using_mem( -(long)(nh-nl+1+NR_END+NR_TEST)*sizeof(double) ); |
---|
396 | dvecs--; |
---|
397 | } |
---|
398 | else |
---|
399 | nrerror("Invalid vector pointer: free_vector"); |
---|
400 | } |
---|
401 | |
---|
402 | void free_matrix(float **m, I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
403 | /* free a float matrix allocated by matrix() */ |
---|
404 | { |
---|
405 | I_ARG_T i, nrow=nrh-nrl+1,ncol=nch-ncl+1; |
---|
406 | if ( valid_matrix( m, nrl, nrh, ncl, nch ) ) { |
---|
407 | m[nrl][ncl-1] = 0.0; |
---|
408 | m[nrh][nch+1] = 0.0; |
---|
409 | *(m-1) = (float *)NULL; |
---|
410 | free((FREE_ARG) (m[nrl]+ncl-NR_END)); |
---|
411 | free((FREE_ARG) (m+nrl-NR_END-1)); |
---|
412 | using_mem(-(long)(nrow+NR_END+1)*sizeof(float*) ); |
---|
413 | using_mem(-(long)(nrow*ncol+NR_END+NR_TEST)*sizeof(float) ); |
---|
414 | mats--; |
---|
415 | } |
---|
416 | else { |
---|
417 | nrerror("Invalid pointer to matrix: free_matrix"); |
---|
418 | } |
---|
419 | } |
---|
420 | |
---|
421 | void free_dmatrix(double **m, I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
422 | /* free a double matrix allocated by dmatrix() */ |
---|
423 | { |
---|
424 | I_ARG_T i, nrow=nrh-nrl+1,ncol=nch-ncl+1; |
---|
425 | if ( valid_dmatrix( m, nrl, nrh, ncl, nch ) ) { |
---|
426 | m[nrl][ncl-1] = 0.0; |
---|
427 | m[nrh][nch+1] = 0.0; |
---|
428 | *(m-1) = (double *)NULL; |
---|
429 | free((FREE_ARG) (m[nrl]+ncl-NR_END)); |
---|
430 | free((FREE_ARG) (m+nrl-NR_END-1)); |
---|
431 | using_mem(-(long)(nrow+NR_END+1)*sizeof(double*) ); |
---|
432 | using_mem(-(long)(nrow*ncol+NR_END+NR_TEST)*sizeof(double) ); |
---|
433 | dmats--; |
---|
434 | } |
---|
435 | else { |
---|
436 | nrerror("Invalid pointer to dmatrix: free_dmatrix"); |
---|
437 | } |
---|
438 | } |
---|
439 | |
---|
440 | void free_imatrix(int **m, I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
441 | /* free an int matrix allocated by imatrix() */ |
---|
442 | { |
---|
443 | free((FREE_ARG) (m[nrl]+ncl-NR_END)); |
---|
444 | free((FREE_ARG) (m+nrl-NR_END)); |
---|
445 | } |
---|
446 | |
---|
447 | void free_submatrix(float **b, I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
448 | /* free a submatrix allocated by submatrix() */ |
---|
449 | { |
---|
450 | free((FREE_ARG) (b+nrl-NR_END)); |
---|
451 | } |
---|
452 | |
---|
453 | void free_convert_matrix(float **b, I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch) |
---|
454 | /* free a matrix allocated by convert_matrix() */ |
---|
455 | { |
---|
456 | free((FREE_ARG) (b+nrl-NR_END)); |
---|
457 | } |
---|
458 | |
---|
459 | void free_f3tensor(float ***t, I_ARG_T nrl, I_ARG_T nrh, I_ARG_T ncl, I_ARG_T nch, |
---|
460 | I_ARG_T ndl, I_ARG_T sndh) |
---|
461 | /* free a float f3tensor allocated by f3tensor() */ |
---|
462 | { |
---|
463 | free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END)); |
---|
464 | free((FREE_ARG) (t[nrl]+ncl-NR_END)); |
---|
465 | free((FREE_ARG) (t+nrl-NR_END)); |
---|
466 | } |
---|
467 | |
---|
468 | /* matrix inversion ref - Page 48 */ |
---|
469 | |
---|
470 | /****************************************************************************/ |
---|
471 | /* void dinverse( double **a, int n, double **y ) |
---|
472 | |
---|
473 | Purpose: |
---|
474 | Find inverse of 'a' (decomposed in process!) and return as 'y' |
---|
475 | |
---|
476 | |
---|
477 | ****************************************************************************/ |
---|
478 | void dinverse( double **a, int n, double **y ) |
---|
479 | { |
---|
480 | double d, *col; |
---|
481 | int i, j, *indx /* integer vector */; |
---|
482 | |
---|
483 | #ifdef V_CHECK |
---|
484 | if ( !valid_dmatrix_b( a ) ) |
---|
485 | nrerror("Invalid input matrix: dinverse"); |
---|
486 | if ( !valid_dmatrix_b( y ) ) |
---|
487 | nrerror("Invalid output matrix: dinverse"); |
---|
488 | #endif |
---|
489 | indx = ivector( 1, n ); |
---|
490 | col = dvector( 1, n ); |
---|
491 | dludcmp( a, n, indx, &d); |
---|
492 | for ( j=1; j<=n; j++ ) { |
---|
493 | for ( i=1; i<=n; i++ ) col[i] = 0.0; |
---|
494 | col[j] = 1.0; |
---|
495 | dlubksb( a, n, indx, col ); |
---|
496 | for ( i=1; i<=n; i++ ) y[i][j] = col[i]; |
---|
497 | } |
---|
498 | free_ivector( indx, 1, n ); |
---|
499 | free_dvector( col, 1, n ); |
---|
500 | } |
---|
501 | |
---|
502 | /* if intending to compute inverse(A) * B, use columns of B in 'col' above instead |
---|
503 | of unit vectors as shown - more accurate and faster */ |
---|
504 | |
---|
505 | void dinverse_mult( double **a, int a_rows, double **b, int b_cols, double **y ) |
---|
506 | /* Find inverse of 'a' (decomposed in process!) times 'b' and return result as 'y' */ |
---|
507 | { |
---|
508 | double d, *col; |
---|
509 | int i, j, *indx /* integer vector */; |
---|
510 | |
---|
511 | #ifdef V_CHECK |
---|
512 | if ( !valid_dmatrix_b( a ) ) |
---|
513 | nrerror("Invalid input matrix: dinverse_mult"); |
---|
514 | if ( !valid_dmatrix_b( b ) ) |
---|
515 | nrerror("Invalid output matrix: dinverse_mult"); |
---|
516 | if ( !valid_dmatrix_b( y ) ) |
---|
517 | nrerror("Invalid output matrix: dinverse_mult"); |
---|
518 | #endif |
---|
519 | indx = ivector( 1, a_rows ); |
---|
520 | col = dvector( 1, a_rows ); |
---|
521 | dludcmp( a, a_rows, indx, &d); |
---|
522 | for ( j=1; j<=b_cols; j++ ) { |
---|
523 | for ( i=1; i<=a_rows; i++ ) col[i] = b[i][j]; |
---|
524 | dlubksb( a, a_rows, indx, col ); |
---|
525 | for ( i=1; i<=a_rows; i++ ) y[i][j] = col[i]; |
---|
526 | } |
---|
527 | free_ivector( indx, 1, a_rows ); |
---|
528 | free_dvector( col, 1, a_rows ); |
---|
529 | } |
---|
530 | |
---|
531 | |
---|
532 | |
---|
533 | |
---|
534 | /* positive definite symmetric matrix inversion ref - Page 97, 98 */ |
---|
535 | |
---|
536 | |
---|
537 | void dPDSinverse( double **a, int n, double **y ) |
---|
538 | /* Find inverse of 'a' (decomposed in process!) and return as 'y' */ |
---|
539 | { |
---|
540 | double sum, *p, *col, *yr; |
---|
541 | int i, j; |
---|
542 | |
---|
543 | #ifdef V_CHECK |
---|
544 | if ( !valid_dmatrix_b( a ) ) |
---|
545 | nrerror("Invalid input matrix: dPDSinverse"); |
---|
546 | if ( !valid_dmatrix_b( y ) ) |
---|
547 | nrerror("Invalid output matrix: dPDSinverse"); |
---|
548 | #endif |
---|
549 | p = dvector( 1, n ); |
---|
550 | col = dvector( 1, n ); |
---|
551 | yr = dvector( 1, n ); |
---|
552 | dcholdc( a, n, p); |
---|
553 | for ( j=1; j<=n; j++ ) { |
---|
554 | for ( i=1; i<=n; i++ ) col[i] = 0.0; |
---|
555 | col[j] = 1.0; |
---|
556 | dcholsl( a, n, p, col, yr ); |
---|
557 | for ( i=1; i<=n; i++ ) y[i][j] = yr[i]; |
---|
558 | } |
---|
559 | free_dvector( yr, 1, n ); |
---|
560 | free_dvector( col, 1, n ); |
---|
561 | free_dvector( p, 1, n ); |
---|
562 | } |
---|
563 | |
---|
564 | |
---|
565 | void dPDS_L_inverse( double **a, int n, double **y ) |
---|
566 | /* Find inverse of L (decomposed a) and return as 'y' */ |
---|
567 | { |
---|
568 | double sum, *p; |
---|
569 | int i, j, k; |
---|
570 | #ifdef V_CHECK |
---|
571 | if ( !valid_dmatrix_b( a) ) |
---|
572 | nrerror("Invalid input matrix: dPDS_L_inverse"); |
---|
573 | if ( !valid_dmatrix_b( y ) ) |
---|
574 | nrerror("Invalid output matrix: dPDS_L_inverse"); |
---|
575 | #endif |
---|
576 | p = dvector( 1, n ); |
---|
577 | dcholdc( a, n, p); |
---|
578 | for ( i=1; i<=n; i++ ) { |
---|
579 | a[i][i] = 1.0/p[i]; |
---|
580 | for ( j=i+1; j<=n; j++ ) { |
---|
581 | sum = 0.0; |
---|
582 | for (k=i; k<j; k++) sum -= a[j][k]*a[k][i]; |
---|
583 | a[j][i] = sum/p[j]; |
---|
584 | } |
---|
585 | } |
---|
586 | free_dvector( p, 1, n ); |
---|
587 | } |
---|