1 /**
2 @file
3 @ingroup nr
4 */
5 #include <bpm/bpm_messages.h>
6 #include <bpm/bpm_nr.h>
7
8 /**
9 Swap two matrix columns
10
11 *
12 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman, Brian Gough
13 *
14 * This program is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License as published by
16 * the Free Software Foundation; either version 2 of the License, or (at
17 * your option) any later version.
18 *
19 * This program is distributed in the hope that it will be useful, but
20 * WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 * General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with this program; if not, write to the Free Software
26 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
27
28 @param m The matrix
29 @param i index of column one
30 @param j index of column two
31
32 @return BPM_SUCCESS if everything was OK, BPM_FAILURE if not
33 */
34
gsl_matrix_swap_columns(gsl_matrix * m,const size_t i,const size_t j)35 int gsl_matrix_swap_columns(gsl_matrix * m, const size_t i, const size_t j) {
36
37 const size_t size1 = m->size1;
38 const size_t size2 = m->size2;
39
40 if (i >= size2)
41 {
42 bpm_error("first column index is out of range in gsl_matrix_swap_columns(...)",
43 __FILE__, __LINE__ );
44 return BPM_FAILURE;
45 }
46
47 if (j >= size2)
48 {
49 bpm_error("second column index is out of range in gsl_matrix_swap_columns(...)",
50 __FILE__, __LINE__ );
51 return BPM_FAILURE;
52 }
53
54 if (i != j)
55 {
56 double *col1 = m->data + i;
57 double *col2 = m->data + j;
58
59 size_t p;
60
61 for (p = 0; p < size1; p++)
62 {
63 size_t k;
64 size_t n = p * m->tda;
65
66 for (k = 0; k < 1; k++)
67 {
68 double tmp = col1[n+k] ;
69 col1[n+k] = col2[n+k] ;
70 col2[n+k] = tmp ;
71 }
72 }
73 }
74
75 return BPM_SUCCESS;
76 }
77
78
79 // -------------------------------------------------------------------
80
81 /**
82 Retrieve a column of a matrix
83
84 @param m The matrix
85 @param j index of the column
86
87 @return BPM_SUCCESS if everything was OK, BPM_FAILURE if not
88 */
89
gsl_matrix_column(gsl_matrix * m,const size_t j)90 _gsl_vector_view gsl_matrix_column (gsl_matrix * m, const size_t j) {
91
92 _gsl_vector_view view = NULL_VECTOR_VIEW;
93
94 if (j >= m->size2)
95 {
96 bpm_error("column index is out of range in gsl_matrix_column",
97 __FILE__, __LINE__ );
98 return view;
99 }
100
101 {
102 gsl_vector v = NULL_VECTOR;
103
104 v.data = m->data + j;
105 v.size = m->size1;
106 v.stride = m->tda;
107 v.block = m->block;
108 v.owner = 0;
109
110 view.vector = v;
111 return view;
112 }
113 }
114
115 /**
116 Get the matrix value associated with the given row and column
117
118 @param m The matrix
119 @param i The row number
120 @param j The column number
121
122 @return The value of the matrix element
123 */
gsl_matrix_get(const gsl_matrix * m,const size_t i,const size_t j)124 double gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j) {
125
126 return *(double*) (m->data + (i * m->tda + j));
127 }
128
129
130
131
132 /**
133 Set the matrix value associated with the given row and column
134
135 @param m The matrix
136 @param i The row number
137 @param j The column number
138 @param x the value to set
139
140 */
gsl_matrix_set(gsl_matrix * m,const size_t i,const size_t j,const double x)141 void gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x) {
142
143 *(double*) (m->data + (i * m->tda + j)) = x;
144
145 }
146
147
148
149 /**
150 Retrieve a submatrix of the given matrix
151 */
gsl_matrix_submatrix(gsl_matrix * m,const size_t i,const size_t j,const size_t n1,const size_t n2)152 _gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * m, const size_t i, const size_t j,
153 const size_t n1, const size_t n2) {
154
155 _gsl_matrix_view view = NULL_MATRIX_VIEW;
156
157 if (i >= m->size1)
158 {
159 bpm_error("row index is out of range in gsl_matrix_submatrix(...)",
160 __FILE__, __LINE__ );
161 return view;
162 }
163 else if (j >= m->size2)
164 {
165 bpm_error("column index is out of range in gsl_matrix_submatrix(...)",
166 __FILE__, __LINE__ );
167 return view;
168 }
169 else if (n1 == 0)
170 {
171 bpm_error("first dimension must be non-zero in gsl_matrix_submatrix(...)",
172 __FILE__, __LINE__ );
173 return view;
174 }
175 else if (n2 == 0)
176 {
177 bpm_error("second dimension must be non-zero in gsl_matrix_submatrix(...)",
178 __FILE__, __LINE__ );
179 return view;
180 }
181 else if (i + n1 > m->size1)
182 {
183 bpm_error("first dimension overflows matrix in gsl_matrix_submatrix(...)",
184 __FILE__, __LINE__ );
185 return view;
186 }
187 else if (j + n2 > m->size2)
188 {
189 bpm_error("second dimension overflows matrix in gsl_matrix_submatrix(...)",
190 __FILE__, __LINE__ );
191 return view;
192 }
193
194 {
195 gsl_matrix s = NULL_MATRIX;
196
197 s.data = m->data + (i * m->tda + j);
198 s.size1 = n1;
199 s.size2 = n2;
200 s.tda = m->tda;
201 s.block = m->block;
202 s.owner = 0;
203
204 view.matrix = s;
205 return view;
206 }
207 }
208
gsl_matrix_alloc(const size_t n1,const size_t n2)209 gsl_matrix * gsl_matrix_alloc (const size_t n1, const size_t n2)
210 {
211
212 gsl_block * block;
213 gsl_matrix * m;
214
215 if (n1 == 0)
216 {
217 bpm_error("matrix dimension n1 must be positive integer in gsl_matrix_alloc(...)",
218 __FILE__, __LINE__);
219 return NULL;
220 }
221 else if (n2 == 0)
222 {
223 bpm_error("matrix dimension n2 must be positive integer in gsl_matrix_alloc(...)",
224 __FILE__, __LINE__);
225 return NULL;
226 }
227
228 m = (gsl_matrix *) malloc (sizeof(gsl_matrix));
229
230 if (m == 0)
231 {
232 bpm_error("failed to allocate space for matrix struct in gsl_matrix_alloc(...)",
233 __FILE__, __LINE__);
234 return NULL;
235 }
236
237 /* FIXME: n1*n2 could overflow for large dimensions */
238
239 block = gsl_block_alloc (n1 * n2) ;
240
241 if (block == 0)
242 {
243 bpm_error("failed to allocate space for block in gsl_matrix_alloc(...)",
244 __FILE__, __LINE__);
245 return NULL;
246 }
247
248 m->data = block->data;
249 m->size1 = n1;
250 m->size2 = n2;
251 m->tda = n2;
252 m->block = block;
253 m->owner = 1;
254
255 return m;
256 }
257
gsl_matrix_calloc(const size_t n1,const size_t n2)258 gsl_matrix * gsl_matrix_calloc (const size_t n1, const size_t n2)
259 {
260 size_t i;
261
262 gsl_matrix * m = gsl_matrix_alloc( n1, n2);
263
264 if (m == 0)
265 return 0;
266
267 /* initialize matrix to zero */
268
269 for (i = 0; i < n1 * n2; i++)
270 {
271 m->data[i] = 0;
272 }
273
274 return m;
275 }
276
gsl_matrix_const_row(const gsl_matrix * m,const size_t i)277 _gsl_vector_const_view gsl_matrix_const_row (const gsl_matrix * m, const size_t i)
278 {
279
280 _gsl_vector_const_view view = NULL_VECTOR_VIEW;
281
282 if (i >= m->size1)
283 {
284 bpm_error("row index is out of range in gsl_matrix_const_row(...)",
285 __FILE__, __LINE__);
286 return view;
287 }
288
289 {
290 gsl_vector v = NULL_VECTOR;
291
292 v.data = m->data + i * m->tda;
293 v.size = m->size2;
294 v.stride = 1;
295 v.block = m->block;
296 v.owner = 0;
297
298 view.vector = v;
299 return view;
300 }
301 }
302
gsl_matrix_row(gsl_matrix * m,const size_t i)303 _gsl_vector_view gsl_matrix_row (gsl_matrix * m, const size_t i)
304 {
305
306 _gsl_vector_view view = NULL_VECTOR_VIEW;
307
308 if (i >= m->size1)
309 {
310 bpm_error("row index is out of range in gsl_matrix_row(...)",
311 __FILE__, __LINE__);
312 return view;
313 }
314
315 {
316 gsl_vector v = NULL_VECTOR;
317
318 v.data = m->data + i * m->tda;
319 v.size = m->size2;
320 v.stride = 1;
321 v.block = m->block;
322 v.owner = 0;
323
324 view.vector = v;
325 return view;
326 }
327 }
328
gsl_matrix_const_column(const gsl_matrix * m,const size_t i)329 _gsl_vector_const_view gsl_matrix_const_column (const gsl_matrix * m, const size_t i)
330 {
331
332 _gsl_vector_const_view view = NULL_VECTOR_VIEW;
333
334 if (i >= m->size1)
335 {
336 bpm_error("row index is out of range", __FILE__, __LINE__ );
337 return view;
338 }
339
340 {
341 gsl_vector v = NULL_VECTOR;
342
343 v.data = m->data + i * m->tda;
344 v.size = m->size2;
345 v.stride = 1;
346 v.block = m->block;
347 v.owner = 0;
348
349 view.vector = v;
350 return view;
351 }
352 }
353
gsl_matrix_set_identity(gsl_matrix * m)354 void gsl_matrix_set_identity (gsl_matrix * m)
355 {
356
357 size_t i, j;
358 double * const data = m->data;
359 const size_t p = m->size1 ;
360 const size_t q = m->size2 ;
361 const size_t tda = m->tda ;
362
363 const double zero = 0.;
364 const double one = 1.;
365
366 for (i = 0; i < p; i++)
367 {
368 for (j = 0; j < q; j++)
369 {
370 *(double *) (data + (i * tda + j)) = ((i == j) ? 1. : 0.);
371 }
372 }
373 }
374
375