00001
00005 #include <bpm/bpm_messages.h>
00006 #include <bpm/bpm_nr.h>
00007
00008 _gsl_vector_view gsl_vector_subvector (gsl_vector *v, size_t offset, size_t n) {
00009
00031 _gsl_vector_view view = NULL_VECTOR_VIEW;
00032
00033 if (n == 0)
00034 {
00035 bpm_error("vector length n must be positive integer in gsl_vector_subvector(...)",
00036 __FILE__, __LINE__);
00037 return view;
00038 }
00039
00040 if (offset + (n - 1) >= v->size)
00041 {
00042 bpm_error("view would extend past end of vector in gsl_vector_subvector(...)",
00043 __FILE__, __LINE__);
00044 return view;
00045 }
00046
00047 {
00048 gsl_vector s = NULL_VECTOR;
00049
00050 s.data = v->data + v->stride * offset ;
00051 s.size = n;
00052 s.stride = v->stride;
00053 s.block = v->block;
00054 s.owner = 0;
00055
00056 view.vector = s;
00057 return view;
00058 }
00059 }
00060
00061 double gsl_vector_get (const gsl_vector * v, const size_t i)
00062 {
00067 return *(double *) (v->data + i * v->stride);
00068 }
00069
00070 void gsl_vector_set (gsl_vector * v, const size_t i, double x)
00071 {
00076 *(double *) (v->data + i * v->stride) = x;
00077 }
00078
00079 int gsl_vector_swap_elements (gsl_vector * v, const size_t i, const size_t j) {
00080
00081
00082 double * data = v->data ;
00083 const size_t size = v->size ;
00084 const size_t stride = v->stride ;
00085
00086 if (i >= size)
00087 {
00088 bpm_error("first index is out of range in gsl_vector_swap_elements(...)",
00089 __FILE__, __LINE__ );
00090 return BPM_FAILURE;
00091 }
00092
00093 if (j >= size)
00094 {
00095 bpm_error("second index is out of range in gsl_vector_swap_elements(...)",
00096 __FILE__, __LINE__ );
00097 return BPM_FAILURE;
00098 }
00099
00100 if (i != j)
00101 {
00102 const size_t s = stride ;
00103 size_t k ;
00104
00105 for (k = 0; k < 1; k++)
00106 {
00107 double tmp = data[j*s + k];
00108 data[j*s+k] = data[i*s + k];
00109 data[i*s+k] = tmp;
00110 }
00111 }
00112
00113 return BPM_SUCCESS;
00114 }
00115
00116 gsl_vector *gsl_vector_alloc (const size_t n)
00117 {
00118
00119 gsl_block * block;
00120 gsl_vector * v;
00121
00122 if (n == 0)
00123 {
00124 bpm_error("vector length n must be positive integer in gsl_vector_alloc(...)",
00125 __FILE__, __LINE__ );
00126 return NULL;
00127 }
00128
00129 v = (gsl_vector*) malloc (sizeof (gsl_vector));
00130
00131 if (v == 0)
00132 {
00133 bpm_error("failed to allocate space for vector struct in gsl_vector_alloc(...)",
00134 __FILE__, __LINE__ );
00135 return NULL;
00136 }
00137
00138 block = gsl_block_alloc(n);
00139
00140 if (block == 0)
00141 {
00142 free (v) ;
00143 bpm_error("failed to allocate space for block in gsl_vector_alloc(...)",
00144 __FILE__, __LINE__ );
00145 return NULL;
00146 }
00147
00148 v->data = block->data ;
00149 v->size = n;
00150 v->stride = 1;
00151 v->block = block;
00152 v->owner = 1;
00153
00154 return v;
00155 }
00156
00157 gsl_vector *gsl_vector_calloc (const size_t n)
00158 {
00159
00160 size_t i;
00161
00162 gsl_vector* v = gsl_vector_alloc (n);
00163
00164 if (v == 0)
00165 return 0;
00166
00167
00168
00169 for (i = 0; i < n; i++)
00170 {
00171 v->data[i] = 0;
00172 }
00173
00174 return v;
00175 }
00176
00177
00178 _gsl_vector_const_view gsl_vector_const_subvector (const gsl_vector *v, size_t offset, size_t n)
00179 {
00180
00181 _gsl_vector_const_view view = NULL_VECTOR_VIEW;
00182
00183 if (n == 0)
00184 {
00185 bpm_error("vector length n must be positive integer in gsl_vector_const_subvector(...)",
00186 __FILE__, __LINE__);
00187 return view;
00188 }
00189
00190 if (offset + (n - 1) >= v->size)
00191 {
00192 bpm_error("view would extend past end of vector in gsl_vector_const_subvector(...)",
00193 __FILE__, __LINE__);
00194 return view;
00195 }
00196
00197 {
00198 gsl_vector s = NULL_VECTOR;
00199
00200 s.data = v->data + v->stride * offset ;
00201 s.size = n;
00202 s.stride = v->stride;
00203 s.block = v->block;
00204 s.owner = 0;
00205
00206 view.vector = s;
00207 return view;
00208 }
00209 }
00210
00211 void gsl_vector_free (gsl_vector * v)
00212 {
00213
00214 if (v->owner)
00215 {
00216 gsl_block_free(v->block) ;
00217 }
00218 free (v);
00219 }