1 #include <stdio.h>
2 #include "blas_extended.h"
3 #include "blas_extended_test.h"
4
chpmv_copy_row(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a,void * a_vec,int row)5 void chpmv_copy_row(enum blas_order_type order, enum blas_uplo_type uplo,
6 int n, void *a, void *a_vec, int row)
7
8 /*
9 * Copies the given row of packed hermitianmatrix a into the supplied vector.
10 */
11 {
12 int i, ind, inc = 1;
13 const float *a_i = (float *) a;
14 float *a_vec_i = (float *) a_vec;
15 float tmp[2];
16
17 inc *= 2;
18
19 if (((order == blas_rowmajor) && (uplo == blas_upper)) ||
20 ((order == blas_colmajor) && (uplo == blas_lower))) {
21 /* Pretend it is colmajor/lower. We can do this in the
22 symmetric case. */
23 ind = row * inc;
24 for (i = 0; i < row; i++) {
25 tmp[0] = a_i[ind];
26 tmp[1] = a_i[ind + 1];
27 if (uplo == blas_upper) {
28 tmp[1] = -tmp[1];
29 }
30 a_vec_i[i * inc] = tmp[0];
31 a_vec_i[i * inc + 1] = tmp[1];
32 ind += (n - i - 1) * inc;
33 }
34
35 ind = (row + ((2 * n - row - 1) * row) / 2) * inc;
36 for (i = row; i < n; i++) {
37 tmp[0] = a_i[ind];
38 tmp[1] = a_i[ind + 1];
39 if (uplo == blas_lower) {
40 tmp[1] = -tmp[1];
41 }
42 if (i == row) {
43 tmp[1] = 0.0;
44 }
45 a_vec_i[i * inc] = tmp[0];
46 a_vec_i[i * inc + 1] = tmp[1];
47 ind += inc;
48 }
49 } else {
50 /* Pretend it is rowmajor/lower. We can do this in the
51 symmetric case. */
52 ind = row * (row + 1) * inc / 2;
53 for (i = 0; i < row; i++) {
54 tmp[0] = a_i[ind];
55 tmp[1] = a_i[ind + 1];
56 if (uplo == blas_upper) {
57 tmp[1] = -tmp[1];
58 }
59 a_vec_i[i * inc] = tmp[0];
60 a_vec_i[i * inc + 1] = tmp[1];
61 ind += inc;
62 }
63
64 ind = (row + (row * (row + 1)) / 2) * inc;
65 for (i = row; i < n; i++) {
66 tmp[0] = a_i[ind];
67 tmp[1] = a_i[ind + 1];
68 if (uplo == blas_lower) {
69 tmp[1] = -tmp[1];
70 }
71 if (i == row) {
72 tmp[1] = 0.0;;
73 }
74 a_vec_i[i * inc] = tmp[0];
75 a_vec_i[i * inc + 1] = tmp[1];
76 ind += (i + 1) * inc;
77 }
78 }
79 }
80
zhpmv_copy_row(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a,void * a_vec,int row)81 void zhpmv_copy_row(enum blas_order_type order, enum blas_uplo_type uplo,
82 int n, void *a, void *a_vec, int row)
83
84 /*
85 * Copies the given row of packed hermitianmatrix a into the supplied vector.
86 */
87 {
88 int i, ind, inc = 1;
89 const double *a_i = (double *) a;
90 double *a_vec_i = (double *) a_vec;
91 double tmp[2];
92
93 inc *= 2;
94
95 if (((order == blas_rowmajor) && (uplo == blas_upper)) ||
96 ((order == blas_colmajor) && (uplo == blas_lower))) {
97 /* Pretend it is colmajor/lower. We can do this in the
98 symmetric case. */
99 ind = row * inc;
100 for (i = 0; i < row; i++) {
101 tmp[0] = a_i[ind];
102 tmp[1] = a_i[ind + 1];
103 if (uplo == blas_upper) {
104 tmp[1] = -tmp[1];
105 }
106 a_vec_i[i * inc] = tmp[0];
107 a_vec_i[i * inc + 1] = tmp[1];
108 ind += (n - i - 1) * inc;
109 }
110
111 ind = (row + ((2 * n - row - 1) * row) / 2) * inc;
112 for (i = row; i < n; i++) {
113 tmp[0] = a_i[ind];
114 tmp[1] = a_i[ind + 1];
115 if (uplo == blas_lower) {
116 tmp[1] = -tmp[1];
117 }
118 if (i == row) {
119 tmp[1] = 0.0;
120 }
121 a_vec_i[i * inc] = tmp[0];
122 a_vec_i[i * inc + 1] = tmp[1];
123 ind += inc;
124 }
125 } else {
126 /* Pretend it is rowmajor/lower. We can do this in the
127 symmetric case. */
128 ind = row * (row + 1) * inc / 2;
129 for (i = 0; i < row; i++) {
130 tmp[0] = a_i[ind];
131 tmp[1] = a_i[ind + 1];
132 if (uplo == blas_upper) {
133 tmp[1] = -tmp[1];
134 }
135 a_vec_i[i * inc] = tmp[0];
136 a_vec_i[i * inc + 1] = tmp[1];
137 ind += inc;
138 }
139
140 ind = (row + (row * (row + 1)) / 2) * inc;
141 for (i = row; i < n; i++) {
142 tmp[0] = a_i[ind];
143 tmp[1] = a_i[ind + 1];
144 if (uplo == blas_lower) {
145 tmp[1] = -tmp[1];
146 }
147 if (i == row) {
148 tmp[1] = 0.0;;
149 }
150 a_vec_i[i * inc] = tmp[0];
151 a_vec_i[i * inc + 1] = tmp[1];
152 ind += (i + 1) * inc;
153 }
154 }
155 }
156
157
chpmv_commit_row(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a,void * a_vec,int row)158 void chpmv_commit_row(enum blas_order_type order, enum blas_uplo_type uplo,
159 int n, void *a, void *a_vec, int row)
160
161 /*
162 * Commits the supplied vector a_vec to the given row of
163 * packed hermitian matrix a.
164 */
165 {
166 int i, ind, inc = 1;
167 float *a_i = (float *) a;
168 float *a_vec_i = (float *) a_vec;
169 float tmp[2];
170
171 inc *= 2;
172
173 if (((order == blas_rowmajor) && (uplo == blas_upper)) ||
174 ((order == blas_colmajor) && (uplo == blas_lower))) {
175 /* Pretend it is colmajor/lower. We can do this in the
176 symmetric case. */
177 ind = row * inc;
178 for (i = 0; i < row; i++) {
179 tmp[0] = a_vec_i[i * inc];
180 tmp[1] = a_vec_i[i * inc + 1];
181 if (uplo == blas_upper) {
182 tmp[1] = -tmp[1];
183 }
184 a_i[ind] = tmp[0];
185 a_i[ind + 1] = tmp[1];
186 ind += (n - i - 1) * inc;
187 }
188
189 ind = (row + ((2 * n - row - 1) * row) / 2) * inc;
190 for (i = row; i < n; i++) {
191 tmp[0] = a_vec_i[i * inc];
192 tmp[1] = a_vec_i[i * inc + 1];
193 if (uplo == blas_lower) {
194 tmp[1] = -tmp[1];
195 }
196 a_i[ind] = tmp[0];
197 a_i[ind + 1] = tmp[1];
198 ind += inc;
199 }
200 } else {
201 /* Pretend it is rowmajor/lower. We can do this in the
202 symmetric case. */
203 ind = row * (row + 1) * inc / 2;
204 for (i = 0; i < row; i++) {
205 tmp[0] = a_vec_i[i * inc];
206 tmp[1] = a_vec_i[i * inc + 1];
207 if (uplo == blas_upper) {
208 tmp[1] = -tmp[1];
209 }
210 a_i[ind] = tmp[0];
211 a_i[ind + 1] = tmp[1];
212 ind += inc;
213 }
214
215 ind = (row + (row * (row + 1)) / 2) * inc;
216 for (i = row; i < n; i++) {
217 tmp[0] = a_vec_i[i * inc];
218 tmp[1] = a_vec_i[i * inc + 1];
219 if (uplo == blas_lower) {
220 tmp[1] = -tmp[1];
221 }
222 a_i[ind] = tmp[0];
223 a_i[ind + 1] = tmp[1];
224 ind += (i + 1) * inc;
225 }
226 }
227 }
228
zhpmv_commit_row(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a,void * a_vec,int row)229 void zhpmv_commit_row(enum blas_order_type order, enum blas_uplo_type uplo,
230 int n, void *a, void *a_vec, int row)
231
232 /*
233 * Commits the supplied vector a_vec to the given row of
234 * packed hermitian matrix a.
235 */
236 {
237 int i, ind, inc = 1;
238 double *a_i = (double *) a;
239 double *a_vec_i = (double *) a_vec;
240 double tmp[2];
241
242 inc *= 2;
243
244 if (((order == blas_rowmajor) && (uplo == blas_upper)) ||
245 ((order == blas_colmajor) && (uplo == blas_lower))) {
246 /* Pretend it is colmajor/lower. We can do this in the
247 symmetric case. */
248 ind = row * inc;
249 for (i = 0; i < row; i++) {
250 tmp[0] = a_vec_i[i * inc];
251 tmp[1] = a_vec_i[i * inc + 1];
252 if (uplo == blas_upper) {
253 tmp[1] = -tmp[1];
254 }
255 a_i[ind] = tmp[0];
256 a_i[ind + 1] = tmp[1];
257 ind += (n - i - 1) * inc;
258 }
259
260 ind = (row + ((2 * n - row - 1) * row) / 2) * inc;
261 for (i = row; i < n; i++) {
262 tmp[0] = a_vec_i[i * inc];
263 tmp[1] = a_vec_i[i * inc + 1];
264 if (uplo == blas_lower) {
265 tmp[1] = -tmp[1];
266 }
267 a_i[ind] = tmp[0];
268 a_i[ind + 1] = tmp[1];
269 ind += inc;
270 }
271 } else {
272 /* Pretend it is rowmajor/lower. We can do this in the
273 symmetric case. */
274 ind = row * (row + 1) * inc / 2;
275 for (i = 0; i < row; i++) {
276 tmp[0] = a_vec_i[i * inc];
277 tmp[1] = a_vec_i[i * inc + 1];
278 if (uplo == blas_upper) {
279 tmp[1] = -tmp[1];
280 }
281 a_i[ind] = tmp[0];
282 a_i[ind + 1] = tmp[1];
283 ind += inc;
284 }
285
286 ind = (row + (row * (row + 1)) / 2) * inc;
287 for (i = row; i < n; i++) {
288 tmp[0] = a_vec_i[i * inc];
289 tmp[1] = a_vec_i[i * inc + 1];
290 if (uplo == blas_lower) {
291 tmp[1] = -tmp[1];
292 }
293 a_i[ind] = tmp[0];
294 a_i[ind + 1] = tmp[1];
295 ind += (i + 1) * inc;
296 }
297 }
298 }
299
300
chpmv_pack_matrix(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a_packed,void * a_full,int lda)301 void chpmv_pack_matrix(enum blas_order_type order, enum blas_uplo_type uplo,
302 int n, void *a_packed, void *a_full, int lda)
303
304 /*
305 * Packs the he matrix a_full into packed form a.
306 */
307 {
308 int row;
309 float *a_row;;
310
311 a_row = (float *) blas_malloc(n * sizeof(float) * 2);
312 if (n > 0 && a_row == NULL) {
313 BLAS_error("blas_malloc", 0, 0, "malloc failed.\n");
314 };
315 for (row = 0; row < n; row++) {
316 che_copy_row(order, uplo, blas_left_side, n, a_full, lda, a_row, row);
317 chpmv_commit_row(order, uplo, n, a_packed, a_row, row);
318 }
319
320 blas_free(a_row);
321 }
zhpmv_pack_matrix(enum blas_order_type order,enum blas_uplo_type uplo,int n,void * a_packed,void * a_full,int lda)322 void zhpmv_pack_matrix(enum blas_order_type order, enum blas_uplo_type uplo,
323 int n, void *a_packed, void *a_full, int lda)
324
325 /*
326 * Packs the he matrix a_full into packed form a.
327 */
328 {
329 int row;
330 double *a_row;;
331
332 a_row = (double *) blas_malloc(n * sizeof(double) * 2);
333 if (n > 0 && a_row == NULL) {
334 BLAS_error("blas_malloc", 0, 0, "malloc failed.\n");
335 };
336 for (row = 0; row < n; row++) {
337 zhe_copy_row(order, uplo, blas_left_side, n, a_full, lda, a_row, row);
338 zhpmv_commit_row(order, uplo, n, a_packed, a_row, row);
339 }
340
341 blas_free(a_row);
342 }
343
cprint_hpmv_matrix(void * a,int n,enum blas_order_type order,enum blas_uplo_type uplo)344 void cprint_hpmv_matrix(void *a, int n,
345 enum blas_order_type order, enum blas_uplo_type uplo)
346 {
347
348 {
349 int row;
350 float *x;
351 x = (float *) blas_malloc(n * sizeof(float) * 2);
352 if (n > 0 && x == NULL) {
353 BLAS_error("blas_malloc", 0, 0, "malloc failed.\n");
354 }
355
356 for (row = 0; row < n; row++) {
357 chpmv_copy_row(order, uplo, n, a, x, row);
358 cprint_vector(x, n, 1, NULL);
359 }
360 printf("\n");
361 blas_free(x);
362 }
363
364 }
zprint_hpmv_matrix(void * a,int n,enum blas_order_type order,enum blas_uplo_type uplo)365 void zprint_hpmv_matrix(void *a, int n,
366 enum blas_order_type order, enum blas_uplo_type uplo)
367 {
368
369 {
370 int row;
371 double *x;
372 x = (double *) blas_malloc(n * sizeof(double) * 2);
373 if (n > 0 && x == NULL) {
374 BLAS_error("blas_malloc", 0, 0, "malloc failed.\n");
375 }
376
377 for (row = 0; row < n; row++) {
378 zhpmv_copy_row(order, uplo, n, a, x, row);
379 zprint_vector(x, n, 1, NULL);
380 }
381 printf("\n");
382 blas_free(x);
383 }
384
385 }
386