1 /* vector/oper_source.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 int
FUNCTION(gsl_vector,add)21 FUNCTION(gsl_vector, add) (TYPE(gsl_vector) * a, const TYPE(gsl_vector) * b)
22 {
23 const size_t N = a->size;
24
25 if (b->size != N)
26 {
27 GSL_ERROR ("vectors must have same length", GSL_EBADLEN);
28 }
29 else
30 {
31 const size_t stride_a = a->stride;
32 const size_t stride_b = b->stride;
33
34 size_t i;
35
36 for (i = 0; i < N; i++)
37 {
38 a->data[i * stride_a] += b->data[i * stride_b];
39 }
40
41 return GSL_SUCCESS;
42 }
43 }
44
45 int
FUNCTION(gsl_vector,sub)46 FUNCTION(gsl_vector, sub) (TYPE(gsl_vector) * a, const TYPE(gsl_vector) * b)
47 {
48 const size_t N = a->size;
49
50 if (b->size != N)
51 {
52 GSL_ERROR ("vectors must have same length", GSL_EBADLEN);
53 }
54 else
55 {
56 const size_t stride_a = a->stride;
57 const size_t stride_b = b->stride;
58
59 size_t i;
60
61 for (i = 0; i < N; i++)
62 {
63 a->data[i * stride_a] -= b->data[i * stride_b];
64 }
65
66 return GSL_SUCCESS;
67 }
68 }
69
70 int
FUNCTION(gsl_vector,mul)71 FUNCTION(gsl_vector, mul) (TYPE(gsl_vector) * a, const TYPE(gsl_vector) * b)
72 {
73 const size_t N = a->size;
74
75 if (b->size != N)
76 {
77 GSL_ERROR ("vectors must have same length", GSL_EBADLEN);
78 }
79 else
80 {
81 const size_t stride_a = a->stride;
82 const size_t stride_b = b->stride;
83
84 size_t i;
85
86 for (i = 0; i < N; i++)
87 {
88 a->data[i * stride_a] *= b->data[i * stride_b];
89 }
90
91 return GSL_SUCCESS;
92 }
93 }
94
95 int
FUNCTION(gsl_vector,div)96 FUNCTION(gsl_vector, div) (TYPE(gsl_vector) * a, const TYPE(gsl_vector) * b)
97 {
98 const size_t N = a->size;
99
100 if (b->size != N)
101 {
102 GSL_ERROR ("vectors must have same length", GSL_EBADLEN);
103 }
104 else
105 {
106 const size_t stride_a = a->stride;
107 const size_t stride_b = b->stride;
108
109 size_t i;
110
111 for (i = 0; i < N; i++)
112 {
113 a->data[i * stride_a] /= b->data[i * stride_b];
114 }
115
116 return GSL_SUCCESS;
117 }
118 }
119
120 int
FUNCTION(gsl_vector,scale)121 FUNCTION(gsl_vector, scale) (TYPE(gsl_vector) * a, const ATOMIC x)
122 {
123 #if defined(BASE_DOUBLE)
124
125 gsl_blas_dscal(x, a);
126
127 #elif defined(BASE_FLOAT)
128
129 gsl_blas_sscal(x, a);
130
131 #else
132
133 const size_t N = a->size;
134 const size_t stride = a->stride;
135
136 size_t i;
137
138 for (i = 0; i < N; i++)
139 {
140 a->data[i * stride] *= x;
141 }
142
143 #endif
144
145 return GSL_SUCCESS;
146 }
147
148 int
FUNCTION(gsl_vector,add_constant)149 FUNCTION(gsl_vector, add_constant) (TYPE(gsl_vector) * a, const ATOMIC x)
150 {
151 const size_t N = a->size;
152 const size_t stride = a->stride;
153
154 size_t i;
155
156 for (i = 0; i < N; i++)
157 {
158 a->data[i * stride] += x;
159 }
160
161 return GSL_SUCCESS;
162 }
163
164 int
FUNCTION(gsl_vector,axpby)165 FUNCTION (gsl_vector, axpby) (const BASE alpha,
166 const TYPE (gsl_vector) * x,
167 const BASE beta,
168 TYPE (gsl_vector) * y)
169 {
170 const size_t x_size = x->size;
171
172 if (x_size != y->size)
173 {
174 GSL_ERROR ("vector lengths are not equal", GSL_EBADLEN);
175 }
176 else if (beta == (ATOMIC) 0)
177 {
178 const size_t x_stride = x->stride;
179 const size_t y_stride = y->stride;
180 size_t j;
181
182 for (j = 0; j < x_size; j++)
183 {
184 y->data[y_stride * j] = alpha * x->data[x_stride * j];
185 }
186
187 return GSL_SUCCESS;
188 }
189 else
190 {
191 const size_t x_stride = x->stride;
192 const size_t y_stride = y->stride;
193 size_t j;
194
195 for (j = 0; j < x_size; j++)
196 {
197 y->data[y_stride * j] = alpha * x->data[x_stride * j] + beta * y->data[y_stride * j];
198 }
199
200 return GSL_SUCCESS;
201 }
202 }
203
204 BASE
FUNCTION(gsl_vector,sum)205 FUNCTION(gsl_vector, sum) (const TYPE(gsl_vector) * a)
206 {
207 const size_t N = a->size;
208 const size_t stride = a->stride;
209 BASE sum = (BASE) 0;
210 size_t i;
211
212 for (i = 0; i < N; i++)
213 {
214 sum += a->data[i * stride];
215 }
216
217 return sum;
218 }
219