1 /* matrix/init_source.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2009 Gerard Jungman, 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
TYPE(gsl_matrix)20 TYPE (gsl_matrix) *
21 FUNCTION (gsl_matrix, alloc) (const size_t n1, const size_t n2)
22 {
23 TYPE (gsl_block) * block;
24 TYPE (gsl_matrix) * m;
25
26 if (n1 == 0)
27 {
28 GSL_ERROR_VAL ("matrix dimension n1 must be positive integer",
29 GSL_EINVAL, 0);
30 }
31 else if (n2 == 0)
32 {
33 GSL_ERROR_VAL ("matrix dimension n2 must be positive integer",
34 GSL_EINVAL, 0);
35 }
36
37 m = (TYPE (gsl_matrix) *) malloc (sizeof (TYPE (gsl_matrix)));
38
39 if (m == 0)
40 {
41 GSL_ERROR_VAL ("failed to allocate space for matrix struct",
42 GSL_ENOMEM, 0);
43 }
44
45 /* FIXME: n1*n2 could overflow for large dimensions */
46
47 block = FUNCTION(gsl_block, alloc) (n1 * n2) ;
48
49 if (block == 0)
50 {
51 GSL_ERROR_VAL ("failed to allocate space for block",
52 GSL_ENOMEM, 0);
53 }
54
55 m->data = block->data;
56 m->size1 = n1;
57 m->size2 = n2;
58 m->tda = n2;
59 m->block = block;
60 m->owner = 1;
61
62 return m;
63 }
64
TYPE(gsl_matrix)65 TYPE (gsl_matrix) *
66 FUNCTION (gsl_matrix, calloc) (const size_t n1, const size_t n2)
67 {
68 size_t i;
69
70 TYPE (gsl_matrix) * m = FUNCTION (gsl_matrix, alloc) (n1, n2);
71
72 if (m == 0)
73 return 0;
74
75 /* initialize matrix to zero */
76
77 for (i = 0; i < MULTIPLICITY * n1 * n2; i++)
78 {
79 m->data[i] = 0;
80 }
81
82 return m;
83 }
84
TYPE(gsl_matrix)85 TYPE (gsl_matrix) *
86 FUNCTION (gsl_matrix, alloc_from_block) (TYPE(gsl_block) * block,
87 const size_t offset,
88 const size_t n1,
89 const size_t n2,
90 const size_t d2)
91 {
92 TYPE (gsl_matrix) * m;
93
94 if (n1 == 0)
95 {
96 GSL_ERROR_VAL ("matrix dimension n1 must be positive integer",
97 GSL_EINVAL, 0);
98 }
99 else if (n2 == 0)
100 {
101 GSL_ERROR_VAL ("matrix dimension n2 must be positive integer",
102 GSL_EINVAL, 0);
103 }
104 else if (d2 < n2)
105 {
106 GSL_ERROR_VAL ("matrix dimension d2 must be greater than n2",
107 GSL_EINVAL, 0);
108 }
109 else if (block->size < offset + n1 * d2)
110 {
111 GSL_ERROR_VAL ("matrix size exceeds available block size",
112 GSL_EINVAL, 0);
113 }
114
115 m = (TYPE (gsl_matrix) *) malloc (sizeof (TYPE (gsl_matrix)));
116
117 if (m == 0)
118 {
119 GSL_ERROR_VAL ("failed to allocate space for matrix struct",
120 GSL_ENOMEM, 0);
121 }
122
123 m->data = block->data + MULTIPLICITY * offset;
124 m->size1 = n1;
125 m->size2 = n2;
126 m->tda = d2;
127 m->block = block;
128 m->owner = 0;
129
130 return m;
131 }
132
133
TYPE(gsl_matrix)134 TYPE (gsl_matrix) *
135 FUNCTION (gsl_matrix, alloc_from_matrix) (TYPE(gsl_matrix) * mm,
136 const size_t k1,
137 const size_t k2,
138 const size_t n1,
139 const size_t n2)
140 {
141 TYPE (gsl_matrix) * m;
142
143 if (n1 == 0)
144 {
145 GSL_ERROR_VAL ("matrix dimension n1 must be positive integer",
146 GSL_EINVAL, 0);
147 }
148 else if (n2 == 0)
149 {
150 GSL_ERROR_VAL ("matrix dimension n2 must be positive integer",
151 GSL_EINVAL, 0);
152 }
153 else if (k1 + n1 > mm->size1)
154 {
155 GSL_ERROR_VAL ("submatrix dimension 1 exceeds size of original",
156 GSL_EINVAL, 0);
157 }
158 else if (k2 + n2 > mm->size2)
159 {
160 GSL_ERROR_VAL ("submatrix dimension 2 exceeds size of original",
161 GSL_EINVAL, 0);
162 }
163
164 m = (TYPE (gsl_matrix) *) malloc (sizeof (TYPE (gsl_matrix)));
165
166 if (m == 0)
167 {
168 GSL_ERROR_VAL ("failed to allocate space for matrix struct",
169 GSL_ENOMEM, 0);
170 }
171
172 m->data = mm->data + k1 * mm->tda + k2 ;
173 m->size1 = n1;
174 m->size2 = n2;
175 m->tda = mm->tda;
176 m->block = mm->block;
177 m->owner = 0;
178
179 return m;
180 }
181
182 void
FUNCTION(gsl_matrix,free)183 FUNCTION (gsl_matrix, free) (TYPE (gsl_matrix) * m)
184 {
185 RETURN_IF_NULL (m);
186
187 if (m->owner)
188 {
189 FUNCTION(gsl_block, free) (m->block);
190 }
191
192 free (m);
193 }
194 void
FUNCTION(gsl_matrix,set_identity)195 FUNCTION (gsl_matrix, set_identity) (TYPE (gsl_matrix) * m)
196 {
197 size_t i, j;
198 ATOMIC * const data = m->data;
199 const size_t p = m->size1 ;
200 const size_t q = m->size2 ;
201 const size_t tda = m->tda ;
202
203 const BASE zero = ZERO;
204 const BASE one = ONE;
205
206 for (i = 0; i < p; i++)
207 {
208 for (j = 0; j < q; j++)
209 {
210 *(BASE *) (data + MULTIPLICITY * (i * tda + j)) = ((i == j) ? one : zero);
211 }
212 }
213 }
214
215 void
FUNCTION(gsl_matrix,set_zero)216 FUNCTION (gsl_matrix, set_zero) (TYPE (gsl_matrix) * m)
217 {
218 size_t i, j;
219 ATOMIC * const data = m->data;
220 const size_t p = m->size1 ;
221 const size_t q = m->size2 ;
222 const size_t tda = m->tda ;
223
224 const BASE zero = ZERO;
225
226 for (i = 0; i < p; i++)
227 {
228 for (j = 0; j < q; j++)
229 {
230 *(BASE *) (data + MULTIPLICITY * (i * tda + j)) = zero;
231 }
232 }
233 }
234
235 void
FUNCTION(gsl_matrix,set_all)236 FUNCTION (gsl_matrix, set_all) (TYPE (gsl_matrix) * m, BASE x)
237 {
238 size_t i, j;
239 ATOMIC * const data = m->data;
240 const size_t p = m->size1 ;
241 const size_t q = m->size2 ;
242 const size_t tda = m->tda ;
243
244 for (i = 0; i < p; i++)
245 {
246 for (j = 0; j < q; j++)
247 {
248 *(BASE *) (data + MULTIPLICITY * (i * tda + j)) = x;
249 }
250 }
251 }
252
253