1 /* matrix/copy_source.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
4 * Copyright (C) 2019 Patrick Alken
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20
21 int
FUNCTION(gsl_matrix,memcpy)22 FUNCTION (gsl_matrix, memcpy) (TYPE (gsl_matrix) * dest,
23 const TYPE (gsl_matrix) * src)
24 {
25 const size_t src_size1 = src->size1;
26 const size_t src_size2 = src->size2;
27 const size_t dest_size1 = dest->size1;
28 const size_t dest_size2 = dest->size2;
29 size_t i;
30
31 if (src_size1 != dest_size1 || src_size2 != dest_size2)
32 {
33 GSL_ERROR ("matrix sizes are different", GSL_EBADLEN);
34 }
35
36 #if defined(BASE_DOUBLE) || defined(BASE_FLOAT) || defined(BASE_GSL_COMPLEX) || defined(BASE_GSL_COMPLEX_FLOAT)
37
38 for (i = 0; i < src_size1; ++i)
39 {
40 VIEW (gsl_vector, const_view) sv = FUNCTION (gsl_matrix, const_row) (src, i);
41 VIEW (gsl_vector, view) dv = FUNCTION (gsl_matrix, row) (dest, i);
42
43 #if defined(BASE_DOUBLE)
44 gsl_blas_dcopy(&sv.vector, &dv.vector);
45 #elif defined(BASE_FLOAT)
46 gsl_blas_scopy(&sv.vector, &dv.vector);
47 #elif defined(BASE_GSL_COMPLEX)
48 gsl_blas_zcopy(&sv.vector, &dv.vector);
49 #elif defined(BASE_GSL_COMPLEX_FLOAT)
50 gsl_blas_ccopy(&sv.vector, &dv.vector);
51 #endif
52 }
53
54 #else
55
56 {
57 const size_t src_tda = src->tda ;
58 const size_t dest_tda = dest->tda ;
59 size_t j;
60
61 for (i = 0; i < src_size1 ; i++)
62 {
63 for (j = 0; j < MULTIPLICITY * src_size2; j++)
64 {
65 dest->data[MULTIPLICITY * dest_tda * i + j]
66 = src->data[MULTIPLICITY * src_tda * i + j];
67 }
68 }
69 }
70
71 #endif
72
73 return GSL_SUCCESS;
74 }
75
76
77 int
FUNCTION(gsl_matrix,swap)78 FUNCTION (gsl_matrix, swap) (TYPE (gsl_matrix) * dest, TYPE (gsl_matrix) * src)
79 {
80 const size_t src_size1 = src->size1;
81 const size_t src_size2 = src->size2;
82 const size_t dest_size1 = dest->size1;
83 const size_t dest_size2 = dest->size2;
84 size_t i;
85
86 if (src_size1 != dest_size1 || src_size2 != dest_size2)
87 {
88 GSL_ERROR ("matrix sizes are different", GSL_EBADLEN);
89 }
90
91 #if defined(BASE_DOUBLE) || defined(BASE_FLOAT) || defined(BASE_GSL_COMPLEX) || defined(BASE_GSL_COMPLEX_FLOAT)
92
93 for (i = 0; i < src_size1; ++i)
94 {
95 VIEW (gsl_vector, view) sv = FUNCTION (gsl_matrix, row) (src, i);
96 VIEW (gsl_vector, view) dv = FUNCTION (gsl_matrix, row) (dest, i);
97
98 #if defined(BASE_DOUBLE)
99 gsl_blas_dswap(&sv.vector, &dv.vector);
100 #elif defined(BASE_FLOAT)
101 gsl_blas_sswap(&sv.vector, &dv.vector);
102 #elif defined(BASE_GSL_COMPLEX)
103 gsl_blas_zswap(&sv.vector, &dv.vector);
104 #elif defined(BASE_GSL_COMPLEX_FLOAT)
105 gsl_blas_cswap(&sv.vector, &dv.vector);
106 #endif
107 }
108
109 #else
110
111 {
112 const size_t src_tda = src->tda ;
113 const size_t dest_tda = dest->tda ;
114 size_t j;
115
116 for (i = 0; i < src_size1 ; i++)
117 {
118 for (j = 0; j < MULTIPLICITY * src_size2; j++)
119 {
120 ATOMIC tmp = src->data[MULTIPLICITY * src_tda * i + j];
121 src->data[MULTIPLICITY * src_tda * i + j]
122 = dest->data[MULTIPLICITY * dest_tda * i + j];
123 dest->data[MULTIPLICITY * dest_tda * i + j] = tmp ;
124 }
125 }
126 }
127
128 #endif
129
130 return GSL_SUCCESS;
131 }
132
133 int
FUNCTION(gsl_matrix,tricpy)134 FUNCTION (gsl_matrix, tricpy) (CBLAS_UPLO_t Uplo,
135 CBLAS_DIAG_t Diag,
136 TYPE (gsl_matrix) * dest,
137 const TYPE (gsl_matrix) * src)
138 {
139 const size_t M = src->size1;
140 const size_t N = src->size2;
141 size_t i;
142
143 if (M != dest->size1 || N != dest->size2)
144 {
145 GSL_ERROR ("matrix sizes are different", GSL_EBADLEN);
146 }
147
148 #if defined(BASE_DOUBLE) || defined(BASE_FLOAT) || defined(BASE_GSL_COMPLEX) || defined(BASE_GSL_COMPLEX_FLOAT)
149
150 if (Uplo == CblasLower)
151 {
152 for (i = 1; i < M; i++)
153 {
154 size_t k = GSL_MIN (i, N);
155 VIEW (gsl_vector, const_view) a = FUNCTION (gsl_matrix, const_subrow) (src, i, 0, k);
156 VIEW (gsl_vector, view) b = FUNCTION (gsl_matrix, subrow) (dest, i, 0, k);
157
158 #if defined(BASE_DOUBLE)
159 gsl_blas_dcopy(&a.vector, &b.vector);
160 #elif defined(BASE_FLOAT)
161 gsl_blas_scopy(&a.vector, &b.vector);
162 #elif defined(BASE_GSL_COMPLEX)
163 gsl_blas_zcopy(&a.vector, &b.vector);
164 #elif defined(BASE_GSL_COMPLEX_FLOAT)
165 gsl_blas_ccopy(&a.vector, &b.vector);
166 #endif
167 }
168 }
169 else if (Uplo == CblasUpper)
170 {
171 for (i = 0; i < GSL_MIN(M, N - 1); i++)
172 {
173 VIEW (gsl_vector, const_view) a = FUNCTION (gsl_matrix, const_subrow) (src, i, i + 1, N - i - 1);
174 VIEW (gsl_vector, view) b = FUNCTION (gsl_matrix, subrow) (dest, i, i + 1, N - i - 1);
175
176 #if defined(BASE_DOUBLE)
177 gsl_blas_dcopy(&a.vector, &b.vector);
178 #elif defined(BASE_FLOAT)
179 gsl_blas_scopy(&a.vector, &b.vector);
180 #elif defined(BASE_GSL_COMPLEX)
181 gsl_blas_zcopy(&a.vector, &b.vector);
182 #elif defined(BASE_GSL_COMPLEX_FLOAT)
183 gsl_blas_ccopy(&a.vector, &b.vector);
184 #endif
185 }
186 }
187
188 if (Diag == CblasNonUnit)
189 {
190 VIEW (gsl_vector, const_view) a = FUNCTION (gsl_matrix, const_diagonal) (src);
191 VIEW (gsl_vector, view) b = FUNCTION (gsl_matrix, diagonal) (dest);
192
193 #if defined(BASE_DOUBLE)
194 gsl_blas_dcopy(&a.vector, &b.vector);
195 #elif defined(BASE_FLOAT)
196 gsl_blas_scopy(&a.vector, &b.vector);
197 #elif defined(BASE_GSL_COMPLEX)
198 gsl_blas_zcopy(&a.vector, &b.vector);
199 #elif defined(BASE_GSL_COMPLEX_FLOAT)
200 gsl_blas_ccopy(&a.vector, &b.vector);
201 #endif
202 }
203
204 #else
205
206 {
207 const size_t src_tda = src->tda ;
208 const size_t dest_tda = dest->tda ;
209 size_t j, k;
210
211 if (Uplo == CblasLower)
212 {
213 for (i = 1; i < M ; i++)
214 {
215 for (j = 0; j < GSL_MIN(i, N); j++)
216 {
217 for (k = 0; k < MULTIPLICITY; k++)
218 {
219 size_t e1 = (i * dest_tda + j) * MULTIPLICITY + k ;
220 size_t e2 = (i * src_tda + j) * MULTIPLICITY + k ;
221 dest->data[e1] = src->data[e2];
222 }
223 }
224 }
225 }
226 else if (Uplo == CblasUpper)
227 {
228 for (i = 0; i < M ; i++)
229 {
230 for (j = i + 1; j < N; j++)
231 {
232 for (k = 0; k < MULTIPLICITY; k++)
233 {
234 size_t e1 = (i * dest_tda + j) * MULTIPLICITY + k ;
235 size_t e2 = (i * src_tda + j) * MULTIPLICITY + k ;
236 dest->data[e1] = src->data[e2];
237 }
238 }
239 }
240 }
241 else
242 {
243 GSL_ERROR ("invalid Uplo parameter", GSL_EINVAL);
244 }
245
246 if (Diag == CblasNonUnit)
247 {
248 for (i = 0; i < GSL_MIN(M, N); i++)
249 {
250 for (k = 0; k < MULTIPLICITY; k++)
251 {
252 size_t e1 = (i * dest_tda + i) * MULTIPLICITY + k ;
253 size_t e2 = (i * src_tda + i) * MULTIPLICITY + k ;
254 dest->data[e1] = src->data[e2];
255 }
256 }
257 }
258 }
259
260 #endif
261
262 return GSL_SUCCESS;
263 }
264
265
266