00001
00005 #include <bpm/bpm_messages.h>
00006 #include <bpm/bpm_nr.h>
00007
00035 int gsl_matrix_swap_columns(gsl_matrix * m, const size_t i, const size_t j) {
00036
00037 const size_t size1 = m->size1;
00038 const size_t size2 = m->size2;
00039
00040 if (i >= size2)
00041 {
00042 bpm_error("first column index is out of range in gsl_matrix_swap_columns(...)",
00043 __FILE__, __LINE__ );
00044 return BPM_FAILURE;
00045 }
00046
00047 if (j >= size2)
00048 {
00049 bpm_error("second column index is out of range in gsl_matrix_swap_columns(...)",
00050 __FILE__, __LINE__ );
00051 return BPM_FAILURE;
00052 }
00053
00054 if (i != j)
00055 {
00056 double *col1 = m->data + i;
00057 double *col2 = m->data + j;
00058
00059 size_t p;
00060
00061 for (p = 0; p < size1; p++)
00062 {
00063 size_t k;
00064 size_t n = p * m->tda;
00065
00066 for (k = 0; k < 1; k++)
00067 {
00068 double tmp = col1[n+k] ;
00069 col1[n+k] = col2[n+k] ;
00070 col2[n+k] = tmp ;
00071 }
00072 }
00073 }
00074
00075 return BPM_SUCCESS;
00076 }
00077
00078
00079
00080
00090 _gsl_vector_view gsl_matrix_column (gsl_matrix * m, const size_t j) {
00091
00092 _gsl_vector_view view = NULL_VECTOR_VIEW;
00093
00094 if (j >= m->size2)
00095 {
00096 bpm_error("column index is out of range in gsl_matrix_column",
00097 __FILE__, __LINE__ );
00098 return view;
00099 }
00100
00101 {
00102 gsl_vector v = NULL_VECTOR;
00103
00104 v.data = m->data + j;
00105 v.size = m->size1;
00106 v.stride = m->tda;
00107 v.block = m->block;
00108 v.owner = 0;
00109
00110 view.vector = v;
00111 return view;
00112 }
00113 }
00114
00124 double gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j) {
00125
00126 return *(double*) (m->data + (i * m->tda + j));
00127 }
00128
00129
00130
00131
00141 void gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x) {
00142
00143 *(double*) (m->data + (i * m->tda + j)) = x;
00144
00145 }
00146
00147
00148
00152 _gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * m, const size_t i, const size_t j,
00153 const size_t n1, const size_t n2) {
00154
00155 _gsl_matrix_view view = NULL_MATRIX_VIEW;
00156
00157 if (i >= m->size1)
00158 {
00159 bpm_error("row index is out of range in gsl_matrix_submatrix(...)",
00160 __FILE__, __LINE__ );
00161 return view;
00162 }
00163 else if (j >= m->size2)
00164 {
00165 bpm_error("column index is out of range in gsl_matrix_submatrix(...)",
00166 __FILE__, __LINE__ );
00167 return view;
00168 }
00169 else if (n1 == 0)
00170 {
00171 bpm_error("first dimension must be non-zero in gsl_matrix_submatrix(...)",
00172 __FILE__, __LINE__ );
00173 return view;
00174 }
00175 else if (n2 == 0)
00176 {
00177 bpm_error("second dimension must be non-zero in gsl_matrix_submatrix(...)",
00178 __FILE__, __LINE__ );
00179 return view;
00180 }
00181 else if (i + n1 > m->size1)
00182 {
00183 bpm_error("first dimension overflows matrix in gsl_matrix_submatrix(...)",
00184 __FILE__, __LINE__ );
00185 return view;
00186 }
00187 else if (j + n2 > m->size2)
00188 {
00189 bpm_error("second dimension overflows matrix in gsl_matrix_submatrix(...)",
00190 __FILE__, __LINE__ );
00191 return view;
00192 }
00193
00194 {
00195 gsl_matrix s = NULL_MATRIX;
00196
00197 s.data = m->data + (i * m->tda + j);
00198 s.size1 = n1;
00199 s.size2 = n2;
00200 s.tda = m->tda;
00201 s.block = m->block;
00202 s.owner = 0;
00203
00204 view.matrix = s;
00205 return view;
00206 }
00207 }
00208
00209 gsl_matrix * gsl_matrix_alloc (const size_t n1, const size_t n2)
00210 {
00211
00212 gsl_block * block;
00213 gsl_matrix * m;
00214
00215 if (n1 == 0)
00216 {
00217 bpm_error("matrix dimension n1 must be positive integer in gsl_matrix_alloc(...)",
00218 __FILE__, __LINE__);
00219 return NULL;
00220 }
00221 else if (n2 == 0)
00222 {
00223 bpm_error("matrix dimension n2 must be positive integer in gsl_matrix_alloc(...)",
00224 __FILE__, __LINE__);
00225 return NULL;
00226 }
00227
00228 m = (gsl_matrix *) malloc (sizeof(gsl_matrix));
00229
00230 if (m == 0)
00231 {
00232 bpm_error("failed to allocate space for matrix struct in gsl_matrix_alloc(...)",
00233 __FILE__, __LINE__);
00234 return NULL;
00235 }
00236
00237
00238
00239 block = gsl_block_alloc (n1 * n2) ;
00240
00241 if (block == 0)
00242 {
00243 bpm_error("failed to allocate space for block in gsl_matrix_alloc(...)",
00244 __FILE__, __LINE__);
00245 return NULL;
00246 }
00247
00248 m->data = block->data;
00249 m->size1 = n1;
00250 m->size2 = n2;
00251 m->tda = n2;
00252 m->block = block;
00253 m->owner = 1;
00254
00255 return m;
00256 }
00257
00258 gsl_matrix * gsl_matrix_calloc (const size_t n1, const size_t n2)
00259 {
00260 size_t i;
00261
00262 gsl_matrix * m = gsl_matrix_alloc( n1, n2);
00263
00264 if (m == 0)
00265 return 0;
00266
00267
00268
00269 for (i = 0; i < n1 * n2; i++)
00270 {
00271 m->data[i] = 0;
00272 }
00273
00274 return m;
00275 }
00276
00277 _gsl_vector_const_view gsl_matrix_const_row (const gsl_matrix * m, const size_t i)
00278 {
00279
00280 _gsl_vector_const_view view = NULL_VECTOR_VIEW;
00281
00282 if (i >= m->size1)
00283 {
00284 bpm_error("row index is out of range in gsl_matrix_const_row(...)",
00285 __FILE__, __LINE__);
00286 return view;
00287 }
00288
00289 {
00290 gsl_vector v = NULL_VECTOR;
00291
00292 v.data = m->data + i * m->tda;
00293 v.size = m->size2;
00294 v.stride = 1;
00295 v.block = m->block;
00296 v.owner = 0;
00297
00298 view.vector = v;
00299 return view;
00300 }
00301 }
00302
00303 _gsl_vector_view gsl_matrix_row (gsl_matrix * m, const size_t i)
00304 {
00305
00306 _gsl_vector_view view = NULL_VECTOR_VIEW;
00307
00308 if (i >= m->size1)
00309 {
00310 bpm_error("row index is out of range in gsl_matrix_row(...)",
00311 __FILE__, __LINE__);
00312 return view;
00313 }
00314
00315 {
00316 gsl_vector v = NULL_VECTOR;
00317
00318 v.data = m->data + i * m->tda;
00319 v.size = m->size2;
00320 v.stride = 1;
00321 v.block = m->block;
00322 v.owner = 0;
00323
00324 view.vector = v;
00325 return view;
00326 }
00327 }
00328
00329 _gsl_vector_const_view gsl_matrix_const_column (const gsl_matrix * m, const size_t i)
00330 {
00331
00332 _gsl_vector_const_view view = NULL_VECTOR_VIEW;
00333
00334 if (i >= m->size1)
00335 {
00336 bpm_error("row index is out of range", __FILE__, __LINE__ );
00337 return view;
00338 }
00339
00340 {
00341 gsl_vector v = NULL_VECTOR;
00342
00343 v.data = m->data + i * m->tda;
00344 v.size = m->size2;
00345 v.stride = 1;
00346 v.block = m->block;
00347 v.owner = 0;
00348
00349 view.vector = v;
00350 return view;
00351 }
00352 }
00353
00354 void gsl_matrix_set_identity (gsl_matrix * m)
00355 {
00356
00357 size_t i, j;
00358 double * const data = m->data;
00359 const size_t p = m->size1 ;
00360 const size_t q = m->size2 ;
00361 const size_t tda = m->tda ;
00362
00363 const double zero = 0.;
00364 const double one = 1.;
00365
00366 for (i = 0; i < p; i++)
00367 {
00368 for (j = 0; j < q; j++)
00369 {
00370 *(double *) (data + (i * tda + j)) = ((i == j) ? 1. : 0.);
00371 }
00372 }
00373 }
00374