1 /* spmatrix/swap_source.c
2 *
3 * Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017, 2018 Patrick Alken
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_spmatrix,transpose)21 FUNCTION (gsl_spmatrix, transpose) (TYPE (gsl_spmatrix) * m)
22 {
23 /* swap dimensions - this must be done before the tree_rebuild step */
24 if (m->size1 != m->size2)
25 {
26 size_t tmp = m->size1;
27 m->size1 = m->size2;
28 m->size2 = tmp;
29 }
30
31 if (GSL_SPMATRIX_ISCOO(m))
32 {
33 size_t n;
34
35 /* swap row/column indices */
36 for (n = 0; n < m->nz; ++n)
37 {
38 int tmp = m->p[n];
39 m->p[n] = m->i[n];
40 m->i[n] = tmp;
41 }
42
43 /* need to rebuild binary tree, or element searches won't
44 * work correctly with transposed indices */
45 FUNCTION (gsl_spmatrix, tree_rebuild) (m);
46 }
47 else if (GSL_SPMATRIX_ISCSC(m))
48 {
49 m->sptype = GSL_SPMATRIX_CSR;
50 }
51 else if (GSL_SPMATRIX_ISCSR(m))
52 {
53 m->sptype = GSL_SPMATRIX_CSC;
54 }
55 else
56 {
57 GSL_ERROR("unknown sparse matrix type", GSL_EINVAL);
58 }
59
60 return GSL_SUCCESS;
61 }
62
63 /* FIXME XXX: deprecated function */
64 int
FUNCTION(gsl_spmatrix,transpose2)65 FUNCTION (gsl_spmatrix, transpose2) (TYPE (gsl_spmatrix) * m)
66 {
67 return FUNCTION (gsl_spmatrix, transpose) (m);
68 }
69
70 int
FUNCTION(gsl_spmatrix,transpose_memcpy)71 FUNCTION (gsl_spmatrix, transpose_memcpy) (TYPE (gsl_spmatrix) * dest, const TYPE (gsl_spmatrix) * src)
72 {
73 const size_t M = src->size1;
74 const size_t N = src->size2;
75
76 if (M != dest->size2 || N != dest->size1)
77 {
78 GSL_ERROR("dimensions of dest must be transpose of src matrix",
79 GSL_EBADLEN);
80 }
81 else if (dest->sptype != src->sptype)
82 {
83 GSL_ERROR("cannot copy matrices of different storage formats",
84 GSL_EINVAL);
85 }
86 else
87 {
88 int status = GSL_SUCCESS;
89 const size_t nz = src->nz;
90
91 if (dest->nzmax < src->nz)
92 {
93 status = FUNCTION (gsl_spmatrix, realloc) (src->nz, dest);
94 if (status)
95 return status;
96 }
97
98 if (GSL_SPMATRIX_ISCOO(src))
99 {
100 size_t n, r;
101 void *ptr;
102
103 for (n = 0; n < nz; ++n)
104 {
105 dest->i[n] = src->p[n];
106 dest->p[n] = src->i[n];
107
108 for (r = 0; r < MULTIPLICITY; ++r)
109 dest->data[MULTIPLICITY * n + r] = src->data[MULTIPLICITY * n + r];
110
111 /* copy binary tree data */
112 ptr = gsl_bst_insert(&dest->data[MULTIPLICITY * n], dest->tree);
113 if (ptr != NULL)
114 {
115 GSL_ERROR("detected duplicate entry", GSL_EINVAL);
116 }
117 }
118 }
119 else if (GSL_SPMATRIX_ISCSC(src))
120 {
121 int * Ai = src->i;
122 int * Ap = src->p;
123 ATOMIC * Ad = src->data;
124 int * ATi = dest->i;
125 int * ATp = dest->p;
126 ATOMIC * ATd = dest->data;
127 int * w = dest->work.work_int;
128 int p;
129 size_t j, r;
130
131 /* initialize to 0 */
132 for (j = 0; j < M + 1; ++j)
133 ATp[j] = 0;
134
135 /* compute row counts of A (= column counts for A^T) */
136 for (j = 0; j < nz; ++j)
137 ATp[Ai[j]]++;
138
139 /* compute row pointers for A (= column pointers for A^T) */
140 gsl_spmatrix_cumsum(M, ATp);
141
142 /* make copy of row pointers */
143 for (j = 0; j < M; ++j)
144 w[j] = ATp[j];
145
146 for (j = 0; j < N; ++j)
147 {
148 for (p = Ap[j]; p < Ap[j + 1]; ++p)
149 {
150 int k = w[Ai[p]]++;
151 ATi[k] = j;
152
153 for (r = 0; r < MULTIPLICITY; ++r)
154 ATd[MULTIPLICITY * k + r] = Ad[MULTIPLICITY * p + r];
155 }
156 }
157 }
158 else if (GSL_SPMATRIX_ISCSR(src))
159 {
160 int * Aj = src->i;
161 int * Ap = src->p;
162 ATOMIC * Ad = src->data;
163 int * ATj = dest->i;
164 int * ATp = dest->p;
165 ATOMIC * ATd = dest->data;
166 int * w = dest->work.work_int;
167 int p;
168 size_t i, r;
169
170 /* initialize to 0 */
171 for (i = 0; i < N + 1; ++i)
172 ATp[i] = 0;
173
174 /* compute column counts of A (= row counts for A^T) */
175 for (i = 0; i < nz; ++i)
176 ATp[Aj[i]]++;
177
178 /* compute column pointers for A (= row pointers for A^T) */
179 gsl_spmatrix_cumsum(N, ATp);
180
181 /* make copy of column pointers */
182 for (i = 0; i < N; ++i)
183 w[i] = ATp[i];
184
185 for (i = 0; i < M; ++i)
186 {
187 for (p = Ap[i]; p < Ap[i + 1]; ++p)
188 {
189 size_t k = w[Aj[p]]++;
190 ATj[k] = i;
191
192 for (r = 0; r < MULTIPLICITY; ++r)
193 ATd[MULTIPLICITY * k + r] = Ad[MULTIPLICITY * p + r];
194 }
195 }
196 }
197 else
198 {
199 GSL_ERROR("unknown sparse matrix type", GSL_EINVAL);
200 }
201
202 dest->nz = nz;
203
204 return status;
205 }
206 }
207