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