1 /* blas/blas.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2009 Gerard Jungman & Brian
4 * Gough
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 /* GSL implementation of BLAS operations for vectors and dense
22 * matrices. Note that GSL native storage is row-major. */
23
24 #include <config.h>
25 #include <gsl/gsl_math.h>
26 #include <gsl/gsl_errno.h>
27 #include <gsl/gsl_cblas.h>
28 #include <gsl/gsl_cblas.h>
29 #include <gsl/gsl_blas_types.h>
30 #include <gsl/gsl_blas.h>
31
32 /* ========================================================================
33 * Level 1
34 * ========================================================================
35 */
36
37 /* CBLAS defines vector sizes in terms of int. GSL defines sizes in
38 terms of size_t, so we need to convert these into integers. There
39 is the possibility of overflow here. FIXME: Maybe this could be
40 caught */
41
42 #define INT(X) ((int)(X))
43
44 int
gsl_blas_sdsdot(float alpha,const gsl_vector_float * X,const gsl_vector_float * Y,float * result)45 gsl_blas_sdsdot (float alpha, const gsl_vector_float * X,
46 const gsl_vector_float * Y, float *result)
47 {
48 if (X->size == Y->size)
49 {
50 *result =
51 cblas_sdsdot (INT (X->size), alpha, X->data, INT (X->stride), Y->data,
52 INT (Y->stride));
53 return GSL_SUCCESS;
54 }
55 else
56 {
57 GSL_ERROR ("invalid length", GSL_EBADLEN);
58 }
59 }
60
61 int
gsl_blas_dsdot(const gsl_vector_float * X,const gsl_vector_float * Y,double * result)62 gsl_blas_dsdot (const gsl_vector_float * X, const gsl_vector_float * Y,
63 double *result)
64 {
65 if (X->size == Y->size)
66 {
67 *result =
68 cblas_dsdot (INT (X->size), X->data, INT (X->stride), Y->data,
69 INT (Y->stride));
70 return GSL_SUCCESS;
71 }
72 else
73 {
74 GSL_ERROR ("invalid length", GSL_EBADLEN);
75 }
76 }
77
78 int
gsl_blas_sdot(const gsl_vector_float * X,const gsl_vector_float * Y,float * result)79 gsl_blas_sdot (const gsl_vector_float * X, const gsl_vector_float * Y,
80 float *result)
81 {
82 if (X->size == Y->size)
83 {
84 *result =
85 cblas_sdot (INT (X->size), X->data, INT (X->stride), Y->data,
86 INT (Y->stride));
87 return GSL_SUCCESS;
88 }
89 else
90 {
91 GSL_ERROR ("invalid length", GSL_EBADLEN);
92 }
93 }
94
95 int
gsl_blas_ddot(const gsl_vector * X,const gsl_vector * Y,double * result)96 gsl_blas_ddot (const gsl_vector * X, const gsl_vector * Y, double *result)
97 {
98 if (X->size == Y->size)
99 {
100 *result =
101 cblas_ddot (INT (X->size), X->data, INT (X->stride), Y->data,
102 INT (Y->stride));
103 return GSL_SUCCESS;
104 }
105 else
106 {
107 GSL_ERROR ("invalid length", GSL_EBADLEN);
108 }
109 }
110
111
112 int
gsl_blas_cdotu(const gsl_vector_complex_float * X,const gsl_vector_complex_float * Y,gsl_complex_float * dotu)113 gsl_blas_cdotu (const gsl_vector_complex_float * X,
114 const gsl_vector_complex_float * Y, gsl_complex_float * dotu)
115 {
116 if (X->size == Y->size)
117 {
118 cblas_cdotu_sub (INT (X->size), X->data, INT (X->stride), Y->data,
119 INT (Y->stride), GSL_COMPLEX_P (dotu));
120 return GSL_SUCCESS;
121 }
122 else
123 {
124 GSL_ERROR ("invalid length", GSL_EBADLEN);
125 }
126 }
127
128
129 int
gsl_blas_cdotc(const gsl_vector_complex_float * X,const gsl_vector_complex_float * Y,gsl_complex_float * dotc)130 gsl_blas_cdotc (const gsl_vector_complex_float * X,
131 const gsl_vector_complex_float * Y, gsl_complex_float * dotc)
132 {
133 if (X->size == Y->size)
134 {
135 cblas_cdotc_sub (INT (X->size), X->data, INT (X->stride), Y->data,
136 INT (Y->stride), GSL_COMPLEX_P (dotc));
137 return GSL_SUCCESS;
138 }
139 else
140 {
141 GSL_ERROR ("invalid length", GSL_EBADLEN);
142 }
143 }
144
145
146 int
gsl_blas_zdotu(const gsl_vector_complex * X,const gsl_vector_complex * Y,gsl_complex * dotu)147 gsl_blas_zdotu (const gsl_vector_complex * X, const gsl_vector_complex * Y,
148 gsl_complex * dotu)
149 {
150 if (X->size == Y->size)
151 {
152 cblas_zdotu_sub (INT (X->size), X->data, INT (X->stride), Y->data,
153 INT (Y->stride), GSL_COMPLEX_P (dotu));
154 return GSL_SUCCESS;
155 }
156 else
157 {
158 GSL_ERROR ("invalid length", GSL_EBADLEN);
159 }
160 }
161
162
163 int
gsl_blas_zdotc(const gsl_vector_complex * X,const gsl_vector_complex * Y,gsl_complex * dotc)164 gsl_blas_zdotc (const gsl_vector_complex * X, const gsl_vector_complex * Y,
165 gsl_complex * dotc)
166 {
167 if (X->size == Y->size)
168 {
169 cblas_zdotc_sub (INT (X->size), X->data, INT (X->stride), Y->data,
170 INT (Y->stride), GSL_COMPLEX_P (dotc));
171 return GSL_SUCCESS;
172 }
173 else
174 {
175 GSL_ERROR ("invalid length", GSL_EBADLEN);
176 }
177 }
178
179 /* Norms of vectors */
180
181 float
gsl_blas_snrm2(const gsl_vector_float * X)182 gsl_blas_snrm2 (const gsl_vector_float * X)
183 {
184 return cblas_snrm2 (INT (X->size), X->data, INT (X->stride));
185 }
186
187 double
gsl_blas_dnrm2(const gsl_vector * X)188 gsl_blas_dnrm2 (const gsl_vector * X)
189 {
190 return cblas_dnrm2 (INT (X->size), X->data, INT (X->stride));
191 }
192
193 float
gsl_blas_scnrm2(const gsl_vector_complex_float * X)194 gsl_blas_scnrm2 (const gsl_vector_complex_float * X)
195 {
196 return cblas_scnrm2 (INT (X->size), X->data, INT (X->stride));
197 }
198
199 double
gsl_blas_dznrm2(const gsl_vector_complex * X)200 gsl_blas_dznrm2 (const gsl_vector_complex * X)
201 {
202 return cblas_dznrm2 (INT (X->size), X->data, INT (X->stride));
203 }
204
205 /* Absolute sums of vectors */
206
207 float
gsl_blas_sasum(const gsl_vector_float * X)208 gsl_blas_sasum (const gsl_vector_float * X)
209 {
210 return cblas_sasum (INT (X->size), X->data, INT (X->stride));
211 }
212
213 double
gsl_blas_dasum(const gsl_vector * X)214 gsl_blas_dasum (const gsl_vector * X)
215 {
216 return cblas_dasum (INT (X->size), X->data, INT (X->stride));
217 }
218
219 float
gsl_blas_scasum(const gsl_vector_complex_float * X)220 gsl_blas_scasum (const gsl_vector_complex_float * X)
221 {
222 return cblas_scasum (INT (X->size), X->data, INT (X->stride));
223 }
224
225 double
gsl_blas_dzasum(const gsl_vector_complex * X)226 gsl_blas_dzasum (const gsl_vector_complex * X)
227 {
228 return cblas_dzasum (INT (X->size), X->data, INT (X->stride));
229 }
230
231 /* Maximum elements of vectors */
232
233 CBLAS_INDEX_t
gsl_blas_isamax(const gsl_vector_float * X)234 gsl_blas_isamax (const gsl_vector_float * X)
235 {
236 return cblas_isamax (INT (X->size), X->data, INT (X->stride));
237 }
238
239 CBLAS_INDEX_t
gsl_blas_idamax(const gsl_vector * X)240 gsl_blas_idamax (const gsl_vector * X)
241 {
242 return cblas_idamax (INT (X->size), X->data, INT (X->stride));
243 }
244
245 CBLAS_INDEX_t
gsl_blas_icamax(const gsl_vector_complex_float * X)246 gsl_blas_icamax (const gsl_vector_complex_float * X)
247 {
248 return cblas_icamax (INT (X->size), X->data, INT (X->stride));
249 }
250
251 CBLAS_INDEX_t
gsl_blas_izamax(const gsl_vector_complex * X)252 gsl_blas_izamax (const gsl_vector_complex * X)
253 {
254 return cblas_izamax (INT (X->size), X->data, INT (X->stride));
255 }
256
257
258 /* Swap vectors */
259
260 int
gsl_blas_sswap(gsl_vector_float * X,gsl_vector_float * Y)261 gsl_blas_sswap (gsl_vector_float * X, gsl_vector_float * Y)
262 {
263 if (X->size == Y->size)
264 {
265 cblas_sswap (INT (X->size), X->data, INT (X->stride), Y->data,
266 INT (Y->stride));
267 return GSL_SUCCESS;
268 }
269 else
270 {
271 GSL_ERROR ("invalid length", GSL_EBADLEN);
272 }
273 }
274
275 int
gsl_blas_dswap(gsl_vector * X,gsl_vector * Y)276 gsl_blas_dswap (gsl_vector * X, gsl_vector * Y)
277 {
278 if (X->size == Y->size)
279 {
280 cblas_dswap (INT (X->size), X->data, INT (X->stride), Y->data,
281 INT (Y->stride));
282 return GSL_SUCCESS;
283 }
284 else
285 {
286 GSL_ERROR ("invalid length", GSL_EBADLEN);
287 };
288 }
289
290 int
gsl_blas_cswap(gsl_vector_complex_float * X,gsl_vector_complex_float * Y)291 gsl_blas_cswap (gsl_vector_complex_float * X, gsl_vector_complex_float * Y)
292 {
293 if (X->size == Y->size)
294 {
295 cblas_cswap (INT (X->size), X->data, INT (X->stride), Y->data,
296 INT (Y->stride));
297 return GSL_SUCCESS;
298 }
299 else
300 {
301 GSL_ERROR ("invalid length", GSL_EBADLEN);
302 }
303 }
304
305 int
gsl_blas_zswap(gsl_vector_complex * X,gsl_vector_complex * Y)306 gsl_blas_zswap (gsl_vector_complex * X, gsl_vector_complex * Y)
307 {
308 if (X->size == Y->size)
309 {
310 cblas_zswap (INT (X->size), X->data, INT (X->stride), Y->data,
311 INT (Y->stride));
312 return GSL_SUCCESS;
313 }
314 else
315 {
316 GSL_ERROR ("invalid length", GSL_EBADLEN);
317 }
318 }
319
320
321 /* Copy vectors */
322
323 int
gsl_blas_scopy(const gsl_vector_float * X,gsl_vector_float * Y)324 gsl_blas_scopy (const gsl_vector_float * X, gsl_vector_float * Y)
325 {
326 if (X->size == Y->size)
327 {
328 cblas_scopy (INT (X->size), X->data, INT (X->stride), Y->data,
329 INT (Y->stride));
330 return GSL_SUCCESS;
331 }
332 else
333 {
334 GSL_ERROR ("invalid length", GSL_EBADLEN);
335 }
336 }
337
338 int
gsl_blas_dcopy(const gsl_vector * X,gsl_vector * Y)339 gsl_blas_dcopy (const gsl_vector * X, gsl_vector * Y)
340 {
341 if (X->size == Y->size)
342 {
343 cblas_dcopy (INT (X->size), X->data, INT (X->stride), Y->data,
344 INT (Y->stride));
345 return GSL_SUCCESS;
346 }
347 else
348 {
349 GSL_ERROR ("invalid length", GSL_EBADLEN);
350 }
351 }
352
353 int
gsl_blas_ccopy(const gsl_vector_complex_float * X,gsl_vector_complex_float * Y)354 gsl_blas_ccopy (const gsl_vector_complex_float * X,
355 gsl_vector_complex_float * Y)
356 {
357 if (X->size == Y->size)
358 {
359 cblas_ccopy (INT (X->size), X->data, INT (X->stride), Y->data,
360 INT (Y->stride));
361 return GSL_SUCCESS;
362 }
363 else
364 {
365 GSL_ERROR ("invalid length", GSL_EBADLEN);
366 }
367 }
368
369 int
gsl_blas_zcopy(const gsl_vector_complex * X,gsl_vector_complex * Y)370 gsl_blas_zcopy (const gsl_vector_complex * X, gsl_vector_complex * Y)
371 {
372 if (X->size == Y->size)
373 {
374 cblas_zcopy (INT (X->size), X->data, INT (X->stride), Y->data,
375 INT (Y->stride));
376 return GSL_SUCCESS;
377 }
378 else
379 {
380 GSL_ERROR ("invalid length", GSL_EBADLEN);
381 }
382 }
383
384
385 /* Compute Y = alpha X + Y */
386
387 int
gsl_blas_saxpy(float alpha,const gsl_vector_float * X,gsl_vector_float * Y)388 gsl_blas_saxpy (float alpha, const gsl_vector_float * X, gsl_vector_float * Y)
389 {
390 if (X->size == Y->size)
391 {
392 cblas_saxpy (INT (X->size), alpha, X->data, INT (X->stride), Y->data,
393 INT (Y->stride));
394 return GSL_SUCCESS;
395 }
396 else
397 {
398 GSL_ERROR ("invalid length", GSL_EBADLEN);
399 }
400 }
401
402 int
gsl_blas_daxpy(double alpha,const gsl_vector * X,gsl_vector * Y)403 gsl_blas_daxpy (double alpha, const gsl_vector * X, gsl_vector * Y)
404 {
405 if (X->size == Y->size)
406 {
407 cblas_daxpy (INT (X->size), alpha, X->data, INT (X->stride), Y->data,
408 INT (Y->stride));
409 return GSL_SUCCESS;
410 }
411 else
412 {
413 GSL_ERROR ("invalid length", GSL_EBADLEN);
414 }
415 }
416
417 int
gsl_blas_caxpy(const gsl_complex_float alpha,const gsl_vector_complex_float * X,gsl_vector_complex_float * Y)418 gsl_blas_caxpy (const gsl_complex_float alpha,
419 const gsl_vector_complex_float * X,
420 gsl_vector_complex_float * Y)
421 {
422 if (X->size == Y->size)
423 {
424 cblas_caxpy (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
425 INT (X->stride), Y->data, INT (Y->stride));
426 return GSL_SUCCESS;
427 }
428 else
429 {
430 GSL_ERROR ("invalid length", GSL_EBADLEN);
431 }
432 }
433
434 int
gsl_blas_zaxpy(const gsl_complex alpha,const gsl_vector_complex * X,gsl_vector_complex * Y)435 gsl_blas_zaxpy (const gsl_complex alpha, const gsl_vector_complex * X,
436 gsl_vector_complex * Y)
437 {
438 if (X->size == Y->size)
439 {
440 cblas_zaxpy (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
441 INT (X->stride), Y->data, INT (Y->stride));
442 return GSL_SUCCESS;
443 }
444 else
445 {
446 GSL_ERROR ("invalid length", GSL_EBADLEN);
447 }
448 }
449
450 /* Generate rotation */
451
452 int
gsl_blas_srotg(float a[],float b[],float c[],float s[])453 gsl_blas_srotg (float a[], float b[], float c[], float s[])
454 {
455 cblas_srotg (a, b, c, s);
456 return GSL_SUCCESS;
457 }
458
459 int
gsl_blas_drotg(double a[],double b[],double c[],double s[])460 gsl_blas_drotg (double a[], double b[], double c[], double s[])
461 {
462 cblas_drotg (a, b, c, s);
463 return GSL_SUCCESS;
464 }
465
466 /* Apply rotation to vectors */
467
468 int
gsl_blas_srot(gsl_vector_float * X,gsl_vector_float * Y,float c,float s)469 gsl_blas_srot (gsl_vector_float * X, gsl_vector_float * Y, float c, float s)
470 {
471 if (X->size == Y->size)
472 {
473 cblas_srot (INT (X->size), X->data, INT (X->stride), Y->data,
474 INT (Y->stride), c, s);
475 return GSL_SUCCESS;
476 }
477 else
478 {
479 GSL_ERROR ("invalid length", GSL_EBADLEN);
480 }
481 }
482
483 int
gsl_blas_drot(gsl_vector * X,gsl_vector * Y,const double c,const double s)484 gsl_blas_drot (gsl_vector * X, gsl_vector * Y, const double c, const double s)
485 {
486 if (X->size == Y->size)
487 {
488 cblas_drot (INT (X->size), X->data, INT (X->stride), Y->data,
489 INT (Y->stride), c, s);
490 return GSL_SUCCESS;
491 }
492 else
493 {
494 GSL_ERROR ("invalid length", GSL_EBADLEN);
495 }
496 }
497
498
499 /* Generate modified rotation */
500
501 int
gsl_blas_srotmg(float d1[],float d2[],float b1[],float b2,float P[])502 gsl_blas_srotmg (float d1[], float d2[], float b1[], float b2, float P[])
503 {
504 cblas_srotmg (d1, d2, b1, b2, P);
505 return GSL_SUCCESS;
506 }
507
508 int
gsl_blas_drotmg(double d1[],double d2[],double b1[],double b2,double P[])509 gsl_blas_drotmg (double d1[], double d2[], double b1[], double b2, double P[])
510 {
511 cblas_drotmg (d1, d2, b1, b2, P);
512 return GSL_SUCCESS;
513 }
514
515
516 /* Apply modified rotation */
517
518 int
gsl_blas_srotm(gsl_vector_float * X,gsl_vector_float * Y,const float P[])519 gsl_blas_srotm (gsl_vector_float * X, gsl_vector_float * Y, const float P[])
520 {
521 if (X->size == Y->size)
522 {
523 cblas_srotm (INT (X->size), X->data, INT (X->stride), Y->data,
524 INT (Y->stride), P);
525 return GSL_SUCCESS;
526 }
527 else
528 {
529 GSL_ERROR ("invalid length", GSL_EBADLEN);
530 }
531 }
532
533 int
gsl_blas_drotm(gsl_vector * X,gsl_vector * Y,const double P[])534 gsl_blas_drotm (gsl_vector * X, gsl_vector * Y, const double P[])
535 {
536 if (X->size == Y->size)
537 {
538 cblas_drotm (INT (X->size), X->data, INT (X->stride), Y->data,
539 INT (Y->stride), P);
540 return GSL_SUCCESS;
541 }
542 else
543 {
544 GSL_ERROR ("invalid length", GSL_EBADLEN);
545 }
546 }
547
548
549 /* Scale vector */
550
551 void
gsl_blas_sscal(float alpha,gsl_vector_float * X)552 gsl_blas_sscal (float alpha, gsl_vector_float * X)
553 {
554 cblas_sscal (INT (X->size), alpha, X->data, INT (X->stride));
555 }
556
557 void
gsl_blas_dscal(double alpha,gsl_vector * X)558 gsl_blas_dscal (double alpha, gsl_vector * X)
559 {
560 cblas_dscal (INT (X->size), alpha, X->data, INT (X->stride));
561 }
562
563 void
gsl_blas_cscal(const gsl_complex_float alpha,gsl_vector_complex_float * X)564 gsl_blas_cscal (const gsl_complex_float alpha, gsl_vector_complex_float * X)
565 {
566 cblas_cscal (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
567 INT (X->stride));
568 }
569
570 void
gsl_blas_zscal(const gsl_complex alpha,gsl_vector_complex * X)571 gsl_blas_zscal (const gsl_complex alpha, gsl_vector_complex * X)
572 {
573 cblas_zscal (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,
574 INT (X->stride));
575 }
576
577 void
gsl_blas_csscal(float alpha,gsl_vector_complex_float * X)578 gsl_blas_csscal (float alpha, gsl_vector_complex_float * X)
579 {
580 cblas_csscal (INT (X->size), alpha, X->data, INT (X->stride));
581 }
582
583 void
gsl_blas_zdscal(double alpha,gsl_vector_complex * X)584 gsl_blas_zdscal (double alpha, gsl_vector_complex * X)
585 {
586 cblas_zdscal (INT (X->size), alpha, X->data, INT (X->stride));
587 }
588
589 /* ===========================================================================
590 * Level 2
591 * ===========================================================================
592 */
593
594 /* GEMV */
595
596 int
gsl_blas_sgemv(CBLAS_TRANSPOSE_t TransA,float alpha,const gsl_matrix_float * A,const gsl_vector_float * X,float beta,gsl_vector_float * Y)597 gsl_blas_sgemv (CBLAS_TRANSPOSE_t TransA, float alpha,
598 const gsl_matrix_float * A, const gsl_vector_float * X,
599 float beta, gsl_vector_float * Y)
600 {
601 const size_t M = A->size1;
602 const size_t N = A->size2;
603
604 if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
605 || (TransA == CblasTrans && M == X->size && N == Y->size))
606 {
607 cblas_sgemv (CblasRowMajor, TransA, INT (M), INT (N), alpha, A->data,
608 INT (A->tda), X->data, INT (X->stride), beta, Y->data,
609 INT (Y->stride));
610 return GSL_SUCCESS;
611 }
612 else
613 {
614 GSL_ERROR ("invalid length", GSL_EBADLEN);
615 }
616 }
617
618
619 int
gsl_blas_dgemv(CBLAS_TRANSPOSE_t TransA,double alpha,const gsl_matrix * A,const gsl_vector * X,double beta,gsl_vector * Y)620 gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA, double alpha, const gsl_matrix * A,
621 const gsl_vector * X, double beta, gsl_vector * Y)
622 {
623 const size_t M = A->size1;
624 const size_t N = A->size2;
625
626 if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
627 || (TransA == CblasTrans && M == X->size && N == Y->size))
628 {
629 cblas_dgemv (CblasRowMajor, TransA, INT (M), INT (N), alpha, A->data,
630 INT (A->tda), X->data, INT (X->stride), beta, Y->data,
631 INT (Y->stride));
632 return GSL_SUCCESS;
633 }
634 else
635 {
636 GSL_ERROR ("invalid length", GSL_EBADLEN);
637 }
638 }
639
640
641 int
gsl_blas_cgemv(CBLAS_TRANSPOSE_t TransA,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,const gsl_vector_complex_float * X,const gsl_complex_float beta,gsl_vector_complex_float * Y)642 gsl_blas_cgemv (CBLAS_TRANSPOSE_t TransA, const gsl_complex_float alpha,
643 const gsl_matrix_complex_float * A,
644 const gsl_vector_complex_float * X,
645 const gsl_complex_float beta, gsl_vector_complex_float * Y)
646 {
647 const size_t M = A->size1;
648 const size_t N = A->size2;
649
650 if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
651 || (TransA == CblasTrans && M == X->size && N == Y->size)
652 || (TransA == CblasConjTrans && M == X->size && N == Y->size))
653 {
654 cblas_cgemv (CblasRowMajor, TransA, INT (M), INT (N),
655 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), X->data,
656 INT (X->stride), GSL_COMPLEX_P (&beta), Y->data,
657 INT (Y->stride));
658 return GSL_SUCCESS;
659 }
660 else
661 {
662 GSL_ERROR ("invalid length", GSL_EBADLEN);
663 }
664 }
665
666
667 int
gsl_blas_zgemv(CBLAS_TRANSPOSE_t TransA,const gsl_complex alpha,const gsl_matrix_complex * A,const gsl_vector_complex * X,const gsl_complex beta,gsl_vector_complex * Y)668 gsl_blas_zgemv (CBLAS_TRANSPOSE_t TransA, const gsl_complex alpha,
669 const gsl_matrix_complex * A, const gsl_vector_complex * X,
670 const gsl_complex beta, gsl_vector_complex * Y)
671 {
672 const size_t M = A->size1;
673 const size_t N = A->size2;
674
675 if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
676 || (TransA == CblasTrans && M == X->size && N == Y->size)
677 || (TransA == CblasConjTrans && M == X->size && N == Y->size))
678 {
679 cblas_zgemv (CblasRowMajor, TransA, INT (M), INT (N),
680 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), X->data,
681 INT (X->stride), GSL_COMPLEX_P (&beta), Y->data,
682 INT (Y->stride));
683 return GSL_SUCCESS;
684 }
685 else
686 {
687 GSL_ERROR ("invalid length", GSL_EBADLEN);
688 }
689 }
690
691
692
693 /* HEMV */
694
695 int
gsl_blas_chemv(CBLAS_UPLO_t Uplo,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,const gsl_vector_complex_float * X,const gsl_complex_float beta,gsl_vector_complex_float * Y)696 gsl_blas_chemv (CBLAS_UPLO_t Uplo, const gsl_complex_float alpha,
697 const gsl_matrix_complex_float * A,
698 const gsl_vector_complex_float * X,
699 const gsl_complex_float beta, gsl_vector_complex_float * Y)
700 {
701 const size_t M = A->size1;
702 const size_t N = A->size2;
703
704 if (M != N)
705 {
706 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
707 }
708 else if (N != X->size || N != Y->size)
709 {
710 GSL_ERROR ("invalid length", GSL_EBADLEN);
711 }
712
713 cblas_chemv (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), A->data,
714 INT (A->tda), X->data, INT (X->stride), GSL_COMPLEX_P (&beta),
715 Y->data, INT (Y->stride));
716 return GSL_SUCCESS;
717 }
718
719 int
gsl_blas_zhemv(CBLAS_UPLO_t Uplo,const gsl_complex alpha,const gsl_matrix_complex * A,const gsl_vector_complex * X,const gsl_complex beta,gsl_vector_complex * Y)720 gsl_blas_zhemv (CBLAS_UPLO_t Uplo, const gsl_complex alpha,
721 const gsl_matrix_complex * A, const gsl_vector_complex * X,
722 const gsl_complex beta, gsl_vector_complex * Y)
723 {
724 const size_t M = A->size1;
725 const size_t N = A->size2;
726
727 if (M != N)
728 {
729 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
730 }
731 else if (N != X->size || N != Y->size)
732 {
733 GSL_ERROR ("invalid length", GSL_EBADLEN);
734 }
735
736 cblas_zhemv (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), A->data,
737 INT (A->tda), X->data, INT (X->stride), GSL_COMPLEX_P (&beta),
738 Y->data, INT (Y->stride));
739 return GSL_SUCCESS;
740 }
741
742
743 /* SYMV */
744
745 int
gsl_blas_ssymv(CBLAS_UPLO_t Uplo,float alpha,const gsl_matrix_float * A,const gsl_vector_float * X,float beta,gsl_vector_float * Y)746 gsl_blas_ssymv (CBLAS_UPLO_t Uplo, float alpha, const gsl_matrix_float * A,
747 const gsl_vector_float * X, float beta, gsl_vector_float * Y)
748 {
749 const size_t M = A->size1;
750 const size_t N = A->size2;
751
752 if (M != N)
753 {
754 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
755 }
756 else if (N != X->size || N != Y->size)
757 {
758 GSL_ERROR ("invalid length", GSL_EBADLEN);
759 }
760
761 cblas_ssymv (CblasRowMajor, Uplo, INT (N), alpha, A->data, INT (A->tda),
762 X->data, INT (X->stride), beta, Y->data, INT (Y->stride));
763 return GSL_SUCCESS;
764 }
765
766 int
gsl_blas_dsymv(CBLAS_UPLO_t Uplo,double alpha,const gsl_matrix * A,const gsl_vector * X,double beta,gsl_vector * Y)767 gsl_blas_dsymv (CBLAS_UPLO_t Uplo, double alpha, const gsl_matrix * A,
768 const gsl_vector * X, double beta, gsl_vector * Y)
769 {
770 const size_t M = A->size1;
771 const size_t N = A->size2;
772
773 if (M != N)
774 {
775 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
776 }
777 else if (N != X->size || N != Y->size)
778 {
779 GSL_ERROR ("invalid length", GSL_EBADLEN);
780 }
781
782 cblas_dsymv (CblasRowMajor, Uplo, INT (N), alpha, A->data, INT (A->tda),
783 X->data, INT (X->stride), beta, Y->data, INT (Y->stride));
784 return GSL_SUCCESS;
785 }
786
787
788 /* TRMV */
789
790 int
gsl_blas_strmv(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_matrix_float * A,gsl_vector_float * X)791 gsl_blas_strmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
792 CBLAS_DIAG_t Diag, const gsl_matrix_float * A,
793 gsl_vector_float * X)
794 {
795 const size_t M = A->size1;
796 const size_t N = A->size2;
797
798 if (M != N)
799 {
800 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
801 }
802 else if (N != X->size)
803 {
804 GSL_ERROR ("invalid length", GSL_EBADLEN);
805 }
806
807 cblas_strmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
808 INT (A->tda), X->data, INT (X->stride));
809 return GSL_SUCCESS;
810 }
811
812
813 int
gsl_blas_dtrmv(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_matrix * A,gsl_vector * X)814 gsl_blas_dtrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
815 CBLAS_DIAG_t Diag, const gsl_matrix * A, gsl_vector * X)
816 {
817 const size_t M = A->size1;
818 const size_t N = A->size2;
819
820 if (M != N)
821 {
822 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
823 }
824 else if (N != X->size)
825 {
826 GSL_ERROR ("invalid length", GSL_EBADLEN);
827 }
828
829 cblas_dtrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
830 INT (A->tda), X->data, INT (X->stride));
831 return GSL_SUCCESS;
832 }
833
834
835 int
gsl_blas_ctrmv(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_matrix_complex_float * A,gsl_vector_complex_float * X)836 gsl_blas_ctrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
837 CBLAS_DIAG_t Diag, const gsl_matrix_complex_float * A,
838 gsl_vector_complex_float * X)
839 {
840 const size_t M = A->size1;
841 const size_t N = A->size2;
842
843 if (M != N)
844 {
845 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
846 }
847 else if (N != X->size)
848 {
849 GSL_ERROR ("invalid length", GSL_EBADLEN);
850 }
851
852 cblas_ctrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
853 INT (A->tda), X->data, INT (X->stride));
854 return GSL_SUCCESS;
855 }
856
857
858 int
gsl_blas_ztrmv(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_matrix_complex * A,gsl_vector_complex * X)859 gsl_blas_ztrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
860 CBLAS_DIAG_t Diag, const gsl_matrix_complex * A,
861 gsl_vector_complex * X)
862 {
863 const size_t M = A->size1;
864 const size_t N = A->size2;
865
866 if (M != N)
867 {
868 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
869 }
870 else if (N != X->size)
871 {
872 GSL_ERROR ("invalid length", GSL_EBADLEN);
873 }
874
875 cblas_ztrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
876 INT (A->tda), X->data, INT (X->stride));
877 return GSL_SUCCESS;
878 }
879
880
881 /* TRSV */
882
883 int
gsl_blas_strsv(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_matrix_float * A,gsl_vector_float * X)884 gsl_blas_strsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
885 CBLAS_DIAG_t Diag, const gsl_matrix_float * A,
886 gsl_vector_float * X)
887 {
888 const size_t M = A->size1;
889 const size_t N = A->size2;
890
891 if (M != N)
892 {
893 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
894 }
895 else if (N != X->size)
896 {
897 GSL_ERROR ("invalid length", GSL_EBADLEN);
898 }
899
900 cblas_strsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
901 INT (A->tda), X->data, INT (X->stride));
902 return GSL_SUCCESS;
903 }
904
905
906 int
gsl_blas_dtrsv(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_matrix * A,gsl_vector * X)907 gsl_blas_dtrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
908 CBLAS_DIAG_t Diag, const gsl_matrix * A, gsl_vector * X)
909 {
910 const size_t M = A->size1;
911 const size_t N = A->size2;
912
913 if (M != N)
914 {
915 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
916 }
917 else if (N != X->size)
918 {
919 GSL_ERROR ("invalid length", GSL_EBADLEN);
920 }
921
922 cblas_dtrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
923 INT (A->tda), X->data, INT (X->stride));
924 return GSL_SUCCESS;
925 }
926
927
928 int
gsl_blas_ctrsv(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_matrix_complex_float * A,gsl_vector_complex_float * X)929 gsl_blas_ctrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
930 CBLAS_DIAG_t Diag, const gsl_matrix_complex_float * A,
931 gsl_vector_complex_float * X)
932 {
933 const size_t M = A->size1;
934 const size_t N = A->size2;
935
936 if (M != N)
937 {
938 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
939 }
940 else if (N != X->size)
941 {
942 GSL_ERROR ("invalid length", GSL_EBADLEN);
943 }
944
945 cblas_ctrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
946 INT (A->tda), X->data, INT (X->stride));
947 return GSL_SUCCESS;
948 }
949
950
951 int
gsl_blas_ztrsv(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_matrix_complex * A,gsl_vector_complex * X)952 gsl_blas_ztrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
953 CBLAS_DIAG_t Diag, const gsl_matrix_complex * A,
954 gsl_vector_complex * X)
955 {
956 const size_t M = A->size1;
957 const size_t N = A->size2;
958
959 if (M != N)
960 {
961 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
962 }
963 else if (N != X->size)
964 {
965 GSL_ERROR ("invalid length", GSL_EBADLEN);
966 }
967
968 cblas_ztrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
969 INT (A->tda), X->data, INT (X->stride));
970 return GSL_SUCCESS;
971 }
972
973
974 /* GER */
975
976 int
gsl_blas_sger(float alpha,const gsl_vector_float * X,const gsl_vector_float * Y,gsl_matrix_float * A)977 gsl_blas_sger (float alpha, const gsl_vector_float * X,
978 const gsl_vector_float * Y, gsl_matrix_float * A)
979 {
980 const size_t M = A->size1;
981 const size_t N = A->size2;
982
983 if (X->size == M && Y->size == N)
984 {
985 cblas_sger (CblasRowMajor, INT (M), INT (N), alpha, X->data,
986 INT (X->stride), Y->data, INT (Y->stride), A->data,
987 INT (A->tda));
988 return GSL_SUCCESS;
989 }
990 else
991 {
992 GSL_ERROR ("invalid length", GSL_EBADLEN);
993 }
994 }
995
996
997 int
gsl_blas_dger(double alpha,const gsl_vector * X,const gsl_vector * Y,gsl_matrix * A)998 gsl_blas_dger (double alpha, const gsl_vector * X, const gsl_vector * Y,
999 gsl_matrix * A)
1000 {
1001 const size_t M = A->size1;
1002 const size_t N = A->size2;
1003
1004 if (X->size == M && Y->size == N)
1005 {
1006 cblas_dger (CblasRowMajor, INT (M), INT (N), alpha, X->data,
1007 INT (X->stride), Y->data, INT (Y->stride), A->data,
1008 INT (A->tda));
1009 return GSL_SUCCESS;
1010 }
1011 else
1012 {
1013 GSL_ERROR ("invalid length", GSL_EBADLEN);
1014 }
1015 }
1016
1017
1018 /* GERU */
1019
1020 int
gsl_blas_cgeru(const gsl_complex_float alpha,const gsl_vector_complex_float * X,const gsl_vector_complex_float * Y,gsl_matrix_complex_float * A)1021 gsl_blas_cgeru (const gsl_complex_float alpha,
1022 const gsl_vector_complex_float * X,
1023 const gsl_vector_complex_float * Y,
1024 gsl_matrix_complex_float * A)
1025 {
1026 const size_t M = A->size1;
1027 const size_t N = A->size2;
1028
1029 if (X->size == M && Y->size == N)
1030 {
1031 cblas_cgeru (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
1032 X->data, INT (X->stride), Y->data, INT (Y->stride),
1033 A->data, INT (A->tda));
1034 return GSL_SUCCESS;
1035 }
1036 else
1037 {
1038 GSL_ERROR ("invalid length", GSL_EBADLEN);
1039 }
1040 }
1041
1042 int
gsl_blas_zgeru(const gsl_complex alpha,const gsl_vector_complex * X,const gsl_vector_complex * Y,gsl_matrix_complex * A)1043 gsl_blas_zgeru (const gsl_complex alpha, const gsl_vector_complex * X,
1044 const gsl_vector_complex * Y, gsl_matrix_complex * A)
1045 {
1046 const size_t M = A->size1;
1047 const size_t N = A->size2;
1048
1049 if (X->size == M && Y->size == N)
1050 {
1051 cblas_zgeru (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
1052 X->data, INT (X->stride), Y->data, INT (Y->stride),
1053 A->data, INT (A->tda));
1054 return GSL_SUCCESS;
1055 }
1056 else
1057 {
1058 GSL_ERROR ("invalid length", GSL_EBADLEN);
1059 }
1060 }
1061
1062
1063 /* GERC */
1064
1065 int
gsl_blas_cgerc(const gsl_complex_float alpha,const gsl_vector_complex_float * X,const gsl_vector_complex_float * Y,gsl_matrix_complex_float * A)1066 gsl_blas_cgerc (const gsl_complex_float alpha,
1067 const gsl_vector_complex_float * X,
1068 const gsl_vector_complex_float * Y,
1069 gsl_matrix_complex_float * A)
1070 {
1071 const size_t M = A->size1;
1072 const size_t N = A->size2;
1073
1074 if (X->size == M && Y->size == N)
1075 {
1076 cblas_cgerc (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
1077 X->data, INT (X->stride), Y->data, INT (Y->stride),
1078 A->data, INT (A->tda));
1079 return GSL_SUCCESS;
1080 }
1081 else
1082 {
1083 GSL_ERROR ("invalid length", GSL_EBADLEN);
1084 }
1085 }
1086
1087
1088 int
gsl_blas_zgerc(const gsl_complex alpha,const gsl_vector_complex * X,const gsl_vector_complex * Y,gsl_matrix_complex * A)1089 gsl_blas_zgerc (const gsl_complex alpha, const gsl_vector_complex * X,
1090 const gsl_vector_complex * Y, gsl_matrix_complex * A)
1091 {
1092 const size_t M = A->size1;
1093 const size_t N = A->size2;
1094
1095 if (X->size == M && Y->size == N)
1096 {
1097 cblas_zgerc (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
1098 X->data, INT (X->stride), Y->data, INT (Y->stride),
1099 A->data, INT (A->tda));
1100 return GSL_SUCCESS;
1101 }
1102 else
1103 {
1104 GSL_ERROR ("invalid length", GSL_EBADLEN);
1105 }
1106 }
1107
1108 /* HER */
1109
1110 int
gsl_blas_cher(CBLAS_UPLO_t Uplo,float alpha,const gsl_vector_complex_float * X,gsl_matrix_complex_float * A)1111 gsl_blas_cher (CBLAS_UPLO_t Uplo, float alpha,
1112 const gsl_vector_complex_float * X,
1113 gsl_matrix_complex_float * A)
1114 {
1115 const size_t M = A->size1;
1116 const size_t N = A->size2;
1117
1118 if (M != N)
1119 {
1120 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
1121 }
1122 else if (X->size != N)
1123 {
1124 GSL_ERROR ("invalid length", GSL_EBADLEN);
1125 }
1126
1127 cblas_cher (CblasRowMajor, Uplo, INT (M), alpha, X->data, INT (X->stride),
1128 A->data, INT (A->tda));
1129 return GSL_SUCCESS;
1130 }
1131
1132
1133 int
gsl_blas_zher(CBLAS_UPLO_t Uplo,double alpha,const gsl_vector_complex * X,gsl_matrix_complex * A)1134 gsl_blas_zher (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector_complex * X,
1135 gsl_matrix_complex * A)
1136 {
1137 const size_t M = A->size1;
1138 const size_t N = A->size2;
1139
1140 if (M != N)
1141 {
1142 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
1143 }
1144 else if (X->size != N)
1145 {
1146 GSL_ERROR ("invalid length", GSL_EBADLEN);
1147 }
1148
1149 cblas_zher (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
1150 A->data, INT (A->tda));
1151 return GSL_SUCCESS;
1152 }
1153
1154
1155 /* HER2 */
1156
1157 int
gsl_blas_cher2(CBLAS_UPLO_t Uplo,const gsl_complex_float alpha,const gsl_vector_complex_float * X,const gsl_vector_complex_float * Y,gsl_matrix_complex_float * A)1158 gsl_blas_cher2 (CBLAS_UPLO_t Uplo, const gsl_complex_float alpha,
1159 const gsl_vector_complex_float * X,
1160 const gsl_vector_complex_float * Y,
1161 gsl_matrix_complex_float * A)
1162 {
1163 const size_t M = A->size1;
1164 const size_t N = A->size2;
1165
1166 if (M != N)
1167 {
1168 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
1169 }
1170 else if (X->size != N || Y->size != N)
1171 {
1172 GSL_ERROR ("invalid length", GSL_EBADLEN);
1173 }
1174
1175 cblas_cher2 (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), X->data,
1176 INT (X->stride), Y->data, INT (Y->stride), A->data,
1177 INT (A->tda));
1178 return GSL_SUCCESS;
1179 }
1180
1181
1182 int
gsl_blas_zher2(CBLAS_UPLO_t Uplo,const gsl_complex alpha,const gsl_vector_complex * X,const gsl_vector_complex * Y,gsl_matrix_complex * A)1183 gsl_blas_zher2 (CBLAS_UPLO_t Uplo, const gsl_complex alpha,
1184 const gsl_vector_complex * X, const gsl_vector_complex * Y,
1185 gsl_matrix_complex * A)
1186 {
1187 const size_t M = A->size1;
1188 const size_t N = A->size2;
1189
1190 if (M != N)
1191 {
1192 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
1193 }
1194 else if (X->size != N || Y->size != N)
1195 {
1196 GSL_ERROR ("invalid length", GSL_EBADLEN);
1197 }
1198
1199 cblas_zher2 (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), X->data,
1200 INT (X->stride), Y->data, INT (Y->stride), A->data,
1201 INT (A->tda));
1202 return GSL_SUCCESS;
1203 }
1204
1205
1206 /* SYR */
1207
1208 int
gsl_blas_ssyr(CBLAS_UPLO_t Uplo,float alpha,const gsl_vector_float * X,gsl_matrix_float * A)1209 gsl_blas_ssyr (CBLAS_UPLO_t Uplo, float alpha, const gsl_vector_float * X,
1210 gsl_matrix_float * A)
1211 {
1212 const size_t M = A->size1;
1213 const size_t N = A->size2;
1214
1215 if (M != N)
1216 {
1217 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
1218 }
1219 else if (X->size != N)
1220 {
1221 GSL_ERROR ("invalid length", GSL_EBADLEN);
1222 }
1223
1224 cblas_ssyr (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
1225 A->data, INT (A->tda));
1226 return GSL_SUCCESS;
1227 }
1228
1229
1230 int
gsl_blas_dsyr(CBLAS_UPLO_t Uplo,double alpha,const gsl_vector * X,gsl_matrix * A)1231 gsl_blas_dsyr (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector * X,
1232 gsl_matrix * A)
1233 {
1234 const size_t M = A->size1;
1235 const size_t N = A->size2;
1236
1237 if (M != N)
1238 {
1239 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
1240 }
1241 else if (X->size != N)
1242 {
1243 GSL_ERROR ("invalid length", GSL_EBADLEN);
1244 }
1245
1246 cblas_dsyr (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
1247 A->data, INT (A->tda));
1248 return GSL_SUCCESS;
1249 }
1250
1251
1252 /* SYR2 */
1253
1254 int
gsl_blas_ssyr2(CBLAS_UPLO_t Uplo,float alpha,const gsl_vector_float * X,const gsl_vector_float * Y,gsl_matrix_float * A)1255 gsl_blas_ssyr2 (CBLAS_UPLO_t Uplo, float alpha, const gsl_vector_float * X,
1256 const gsl_vector_float * Y, gsl_matrix_float * A)
1257 {
1258 const size_t M = A->size1;
1259 const size_t N = A->size2;
1260
1261 if (M != N)
1262 {
1263 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
1264 }
1265 else if (X->size != N || Y->size != N)
1266 {
1267 GSL_ERROR ("invalid length", GSL_EBADLEN);
1268 }
1269
1270 cblas_ssyr2 (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
1271 Y->data, INT (Y->stride), A->data, INT (A->tda));
1272 return GSL_SUCCESS;
1273 }
1274
1275
1276 int
gsl_blas_dsyr2(CBLAS_UPLO_t Uplo,double alpha,const gsl_vector * X,const gsl_vector * Y,gsl_matrix * A)1277 gsl_blas_dsyr2 (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector * X,
1278 const gsl_vector * Y, gsl_matrix * A)
1279 {
1280 const size_t M = A->size1;
1281 const size_t N = A->size2;
1282
1283 if (M != N)
1284 {
1285 GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
1286 }
1287 else if (X->size != N || Y->size != N)
1288 {
1289 GSL_ERROR ("invalid length", GSL_EBADLEN);
1290 }
1291
1292 cblas_dsyr2 (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
1293 Y->data, INT (Y->stride), A->data, INT (A->tda));
1294 return GSL_SUCCESS;
1295 }
1296
1297
1298 /*
1299 * ===========================================================================
1300 * Prototypes for level 3 BLAS
1301 * ===========================================================================
1302 */
1303
1304
1305 /* GEMM */
1306
1307 int
gsl_blas_sgemm(CBLAS_TRANSPOSE_t TransA,CBLAS_TRANSPOSE_t TransB,float alpha,const gsl_matrix_float * A,const gsl_matrix_float * B,float beta,gsl_matrix_float * C)1308 gsl_blas_sgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
1309 float alpha, const gsl_matrix_float * A,
1310 const gsl_matrix_float * B, float beta, gsl_matrix_float * C)
1311 {
1312 const size_t M = C->size1;
1313 const size_t N = C->size2;
1314 const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
1315 const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
1316 const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
1317 const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;
1318
1319 if (M == MA && N == NB && NA == MB) /* [MxN] = [MAxNA][MBxNB] */
1320 {
1321 cblas_sgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
1322 alpha, A->data, INT (A->tda), B->data, INT (B->tda), beta,
1323 C->data, INT (C->tda));
1324 return GSL_SUCCESS;
1325 }
1326 else
1327 {
1328 GSL_ERROR ("invalid length", GSL_EBADLEN);
1329 }
1330 }
1331
1332
1333 int
gsl_blas_dgemm(CBLAS_TRANSPOSE_t TransA,CBLAS_TRANSPOSE_t TransB,double alpha,const gsl_matrix * A,const gsl_matrix * B,double beta,gsl_matrix * C)1334 gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
1335 double alpha, const gsl_matrix * A, const gsl_matrix * B,
1336 double beta, gsl_matrix * C)
1337 {
1338 const size_t M = C->size1;
1339 const size_t N = C->size2;
1340 const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
1341 const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
1342 const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
1343 const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;
1344
1345 if (M == MA && N == NB && NA == MB) /* [MxN] = [MAxNA][MBxNB] */
1346 {
1347 cblas_dgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
1348 alpha, A->data, INT (A->tda), B->data, INT (B->tda), beta,
1349 C->data, INT (C->tda));
1350 return GSL_SUCCESS;
1351 }
1352 else
1353 {
1354 GSL_ERROR ("invalid length", GSL_EBADLEN);
1355 }
1356 }
1357
1358
1359 int
gsl_blas_cgemm(CBLAS_TRANSPOSE_t TransA,CBLAS_TRANSPOSE_t TransB,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,const gsl_matrix_complex_float * B,const gsl_complex_float beta,gsl_matrix_complex_float * C)1360 gsl_blas_cgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
1361 const gsl_complex_float alpha,
1362 const gsl_matrix_complex_float * A,
1363 const gsl_matrix_complex_float * B,
1364 const gsl_complex_float beta, gsl_matrix_complex_float * C)
1365 {
1366 const size_t M = C->size1;
1367 const size_t N = C->size2;
1368 const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
1369 const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
1370 const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
1371 const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;
1372
1373 if (M == MA && N == NB && NA == MB) /* [MxN] = [MAxNA][MBxNB] */
1374 {
1375 cblas_cgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
1376 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1377 INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
1378 INT (C->tda));
1379 return GSL_SUCCESS;
1380 }
1381 else
1382 {
1383 GSL_ERROR ("invalid length", GSL_EBADLEN);
1384 }
1385 }
1386
1387
1388 int
gsl_blas_zgemm(CBLAS_TRANSPOSE_t TransA,CBLAS_TRANSPOSE_t TransB,const gsl_complex alpha,const gsl_matrix_complex * A,const gsl_matrix_complex * B,const gsl_complex beta,gsl_matrix_complex * C)1389 gsl_blas_zgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
1390 const gsl_complex alpha, const gsl_matrix_complex * A,
1391 const gsl_matrix_complex * B, const gsl_complex beta,
1392 gsl_matrix_complex * C)
1393 {
1394 const size_t M = C->size1;
1395 const size_t N = C->size2;
1396 const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
1397 const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
1398 const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
1399 const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;
1400
1401 if (M == MA && N == NB && NA == MB) /* [MxN] = [MAxNA][MBxNB] */
1402 {
1403 cblas_zgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
1404 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1405 INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
1406 INT (C->tda));
1407 return GSL_SUCCESS;
1408 }
1409 else
1410 {
1411 GSL_ERROR ("invalid length", GSL_EBADLEN);
1412 }
1413 }
1414
1415
1416 /* SYMM */
1417
1418 int
gsl_blas_ssymm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,float alpha,const gsl_matrix_float * A,const gsl_matrix_float * B,float beta,gsl_matrix_float * C)1419 gsl_blas_ssymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, float alpha,
1420 const gsl_matrix_float * A, const gsl_matrix_float * B,
1421 float beta, gsl_matrix_float * C)
1422 {
1423 const size_t M = C->size1;
1424 const size_t N = C->size2;
1425 const size_t MA = A->size1;
1426 const size_t NA = A->size2;
1427 const size_t MB = B->size1;
1428 const size_t NB = B->size2;
1429
1430 if (MA != NA)
1431 {
1432 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
1433 }
1434
1435 if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
1436 || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
1437 {
1438 cblas_ssymm (CblasRowMajor, Side, Uplo, INT (M), INT (N), alpha,
1439 A->data, INT (A->tda), B->data, INT (B->tda), beta,
1440 C->data, INT (C->tda));
1441 return GSL_SUCCESS;
1442 }
1443 else
1444 {
1445 GSL_ERROR ("invalid length", GSL_EBADLEN);
1446 }
1447
1448 }
1449
1450
1451 int
gsl_blas_dsymm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,double alpha,const gsl_matrix * A,const gsl_matrix * B,double beta,gsl_matrix * C)1452 gsl_blas_dsymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, double alpha,
1453 const gsl_matrix * A, const gsl_matrix * B, double beta,
1454 gsl_matrix * C)
1455 {
1456 const size_t M = C->size1;
1457 const size_t N = C->size2;
1458 const size_t MA = A->size1;
1459 const size_t NA = A->size2;
1460 const size_t MB = B->size1;
1461 const size_t NB = B->size2;
1462
1463 if (MA != NA)
1464 {
1465 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
1466 }
1467
1468 if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
1469 || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
1470 {
1471 cblas_dsymm (CblasRowMajor, Side, Uplo, INT (M), INT (N), alpha,
1472 A->data, INT (A->tda), B->data, INT (B->tda), beta,
1473 C->data, INT (C->tda));
1474 return GSL_SUCCESS;
1475 }
1476 else
1477 {
1478 GSL_ERROR ("invalid length", GSL_EBADLEN);
1479 }
1480 }
1481
1482
1483 int
gsl_blas_csymm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,const gsl_matrix_complex_float * B,const gsl_complex_float beta,gsl_matrix_complex_float * C)1484 gsl_blas_csymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
1485 const gsl_complex_float alpha,
1486 const gsl_matrix_complex_float * A,
1487 const gsl_matrix_complex_float * B,
1488 const gsl_complex_float beta, gsl_matrix_complex_float * C)
1489 {
1490 const size_t M = C->size1;
1491 const size_t N = C->size2;
1492 const size_t MA = A->size1;
1493 const size_t NA = A->size2;
1494 const size_t MB = B->size1;
1495 const size_t NB = B->size2;
1496
1497 if (MA != NA)
1498 {
1499 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
1500 }
1501
1502 if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
1503 || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
1504 {
1505 cblas_csymm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
1506 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1507 INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
1508 INT (C->tda));
1509 return GSL_SUCCESS;
1510 }
1511 else
1512 {
1513 GSL_ERROR ("invalid length", GSL_EBADLEN);
1514 }
1515 }
1516
1517 int
gsl_blas_zsymm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,const gsl_complex alpha,const gsl_matrix_complex * A,const gsl_matrix_complex * B,const gsl_complex beta,gsl_matrix_complex * C)1518 gsl_blas_zsymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
1519 const gsl_complex alpha, const gsl_matrix_complex * A,
1520 const gsl_matrix_complex * B, const gsl_complex beta,
1521 gsl_matrix_complex * C)
1522 {
1523 const size_t M = C->size1;
1524 const size_t N = C->size2;
1525 const size_t MA = A->size1;
1526 const size_t NA = A->size2;
1527 const size_t MB = B->size1;
1528 const size_t NB = B->size2;
1529
1530 if (MA != NA)
1531 {
1532 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
1533 }
1534
1535 if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
1536 || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
1537 {
1538 cblas_zsymm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
1539 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1540 INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
1541 INT (C->tda));
1542 return GSL_SUCCESS;
1543 }
1544 else
1545 {
1546 GSL_ERROR ("invalid length", GSL_EBADLEN);
1547 }
1548 }
1549
1550
1551 /* HEMM */
1552
1553 int
gsl_blas_chemm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,const gsl_matrix_complex_float * B,const gsl_complex_float beta,gsl_matrix_complex_float * C)1554 gsl_blas_chemm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
1555 const gsl_complex_float alpha,
1556 const gsl_matrix_complex_float * A,
1557 const gsl_matrix_complex_float * B,
1558 const gsl_complex_float beta, gsl_matrix_complex_float * C)
1559 {
1560 const size_t M = C->size1;
1561 const size_t N = C->size2;
1562 const size_t MA = A->size1;
1563 const size_t NA = A->size2;
1564 const size_t MB = B->size1;
1565 const size_t NB = B->size2;
1566
1567 if (MA != NA)
1568 {
1569 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
1570 }
1571
1572 if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
1573 || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
1574 {
1575 cblas_chemm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
1576 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1577 INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
1578 INT (C->tda));
1579 return GSL_SUCCESS;
1580 }
1581 else
1582 {
1583 GSL_ERROR ("invalid length", GSL_EBADLEN);
1584 }
1585
1586 }
1587
1588
1589 int
gsl_blas_zhemm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,const gsl_complex alpha,const gsl_matrix_complex * A,const gsl_matrix_complex * B,const gsl_complex beta,gsl_matrix_complex * C)1590 gsl_blas_zhemm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
1591 const gsl_complex alpha, const gsl_matrix_complex * A,
1592 const gsl_matrix_complex * B, const gsl_complex beta,
1593 gsl_matrix_complex * C)
1594 {
1595 const size_t M = C->size1;
1596 const size_t N = C->size2;
1597 const size_t MA = A->size1;
1598 const size_t NA = A->size2;
1599 const size_t MB = B->size1;
1600 const size_t NB = B->size2;
1601
1602 if (MA != NA)
1603 {
1604 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
1605 }
1606
1607 if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
1608 || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
1609 {
1610 cblas_zhemm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
1611 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1612 INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
1613 INT (C->tda));
1614 return GSL_SUCCESS;
1615 }
1616 else
1617 {
1618 GSL_ERROR ("invalid length", GSL_EBADLEN);
1619 }
1620 }
1621
1622 /* SYRK */
1623
1624 int
gsl_blas_ssyrk(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,float alpha,const gsl_matrix_float * A,float beta,gsl_matrix_float * C)1625 gsl_blas_ssyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
1626 const gsl_matrix_float * A, float beta, gsl_matrix_float * C)
1627 {
1628 const size_t M = C->size1;
1629 const size_t N = C->size2;
1630 const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1631 const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1632
1633 if (M != N)
1634 {
1635 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1636 }
1637 else if (N != J)
1638 {
1639 GSL_ERROR ("invalid length", GSL_EBADLEN);
1640 }
1641
1642 cblas_ssyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
1643 INT (A->tda), beta, C->data, INT (C->tda));
1644 return GSL_SUCCESS;
1645 }
1646
1647
1648 int
gsl_blas_dsyrk(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,double alpha,const gsl_matrix * A,double beta,gsl_matrix * C)1649 gsl_blas_dsyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
1650 const gsl_matrix * A, double beta, gsl_matrix * C)
1651 {
1652 const size_t M = C->size1;
1653 const size_t N = C->size2;
1654 const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1655 const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1656
1657 if (M != N)
1658 {
1659 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1660 }
1661 else if (N != J)
1662 {
1663 GSL_ERROR ("invalid length", GSL_EBADLEN);
1664 }
1665
1666 cblas_dsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
1667 INT (A->tda), beta, C->data, INT (C->tda));
1668 return GSL_SUCCESS;
1669
1670 }
1671
1672
1673 int
gsl_blas_csyrk(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,const gsl_complex_float beta,gsl_matrix_complex_float * C)1674 gsl_blas_csyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
1675 const gsl_complex_float alpha,
1676 const gsl_matrix_complex_float * A,
1677 const gsl_complex_float beta, gsl_matrix_complex_float * C)
1678 {
1679 const size_t M = C->size1;
1680 const size_t N = C->size2;
1681 const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1682 const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1683
1684 if (M != N)
1685 {
1686 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1687 }
1688 else if (N != J)
1689 {
1690 GSL_ERROR ("invalid length", GSL_EBADLEN);
1691 }
1692
1693 cblas_csyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K),
1694 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda),
1695 GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
1696 return GSL_SUCCESS;
1697 }
1698
1699
1700 int
gsl_blas_zsyrk(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,const gsl_complex alpha,const gsl_matrix_complex * A,const gsl_complex beta,gsl_matrix_complex * C)1701 gsl_blas_zsyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
1702 const gsl_complex alpha, const gsl_matrix_complex * A,
1703 const gsl_complex beta, gsl_matrix_complex * C)
1704 {
1705 const size_t M = C->size1;
1706 const size_t N = C->size2;
1707 const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1708 const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1709
1710 if (M != N)
1711 {
1712 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1713 }
1714 else if (N != J)
1715 {
1716 GSL_ERROR ("invalid length", GSL_EBADLEN);
1717 }
1718
1719 cblas_zsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K),
1720 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda),
1721 GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
1722 return GSL_SUCCESS;
1723 }
1724
1725 /* HERK */
1726
1727 int
gsl_blas_cherk(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,float alpha,const gsl_matrix_complex_float * A,float beta,gsl_matrix_complex_float * C)1728 gsl_blas_cherk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
1729 const gsl_matrix_complex_float * A, float beta,
1730 gsl_matrix_complex_float * C)
1731 {
1732 const size_t M = C->size1;
1733 const size_t N = C->size2;
1734 const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1735 const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1736
1737 if (M != N)
1738 {
1739 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1740 }
1741 else if (N != J)
1742 {
1743 GSL_ERROR ("invalid length", GSL_EBADLEN);
1744 }
1745
1746 cblas_cherk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
1747 INT (A->tda), beta, C->data, INT (C->tda));
1748 return GSL_SUCCESS;
1749 }
1750
1751
1752 int
gsl_blas_zherk(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,double alpha,const gsl_matrix_complex * A,double beta,gsl_matrix_complex * C)1753 gsl_blas_zherk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
1754 const gsl_matrix_complex * A, double beta,
1755 gsl_matrix_complex * C)
1756 {
1757 const size_t M = C->size1;
1758 const size_t N = C->size2;
1759 const size_t J = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1760 const size_t K = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1761
1762 if (M != N)
1763 {
1764 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1765 }
1766 else if (N != J)
1767 {
1768 GSL_ERROR ("invalid length", GSL_EBADLEN);
1769 }
1770
1771 cblas_zherk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
1772 INT (A->tda), beta, C->data, INT (C->tda));
1773 return GSL_SUCCESS;
1774 }
1775
1776 /* SYR2K */
1777
1778 int
gsl_blas_ssyr2k(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,float alpha,const gsl_matrix_float * A,const gsl_matrix_float * B,float beta,gsl_matrix_float * C)1779 gsl_blas_ssyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
1780 const gsl_matrix_float * A, const gsl_matrix_float * B,
1781 float beta, gsl_matrix_float * C)
1782 {
1783 const size_t M = C->size1;
1784 const size_t N = C->size2;
1785 const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1786 const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1787 const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
1788 const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
1789
1790 if (M != N)
1791 {
1792 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1793 }
1794 else if (N != MA || N != MB || NA != NB)
1795 {
1796 GSL_ERROR ("invalid length", GSL_EBADLEN);
1797 }
1798
1799 cblas_ssyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA), alpha, A->data,
1800 INT (A->tda), B->data, INT (B->tda), beta, C->data,
1801 INT (C->tda));
1802 return GSL_SUCCESS;
1803 }
1804
1805
1806 int
gsl_blas_dsyr2k(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,double alpha,const gsl_matrix * A,const gsl_matrix * B,double beta,gsl_matrix * C)1807 gsl_blas_dsyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
1808 const gsl_matrix * A, const gsl_matrix * B, double beta,
1809 gsl_matrix * C)
1810 {
1811 const size_t M = C->size1;
1812 const size_t N = C->size2;
1813 const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1814 const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1815 const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
1816 const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
1817
1818 if (M != N)
1819 {
1820 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1821 }
1822 else if (N != MA || N != MB || NA != NB)
1823 {
1824 GSL_ERROR ("invalid length", GSL_EBADLEN);
1825 }
1826
1827 cblas_dsyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA), alpha, A->data,
1828 INT (A->tda), B->data, INT (B->tda), beta, C->data,
1829 INT (C->tda));
1830 return GSL_SUCCESS;
1831 }
1832
1833
1834 int
gsl_blas_csyr2k(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,const gsl_matrix_complex_float * B,const gsl_complex_float beta,gsl_matrix_complex_float * C)1835 gsl_blas_csyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
1836 const gsl_complex_float alpha,
1837 const gsl_matrix_complex_float * A,
1838 const gsl_matrix_complex_float * B,
1839 const gsl_complex_float beta, gsl_matrix_complex_float * C)
1840 {
1841 const size_t M = C->size1;
1842 const size_t N = C->size2;
1843 const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1844 const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1845 const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
1846 const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
1847
1848 if (M != N)
1849 {
1850 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1851 }
1852 else if (N != MA || N != MB || NA != NB)
1853 {
1854 GSL_ERROR ("invalid length", GSL_EBADLEN);
1855 }
1856
1857 cblas_csyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
1858 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1859 INT (B->tda), GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
1860 return GSL_SUCCESS;
1861 }
1862
1863
1864
1865 int
gsl_blas_zsyr2k(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,const gsl_complex alpha,const gsl_matrix_complex * A,const gsl_matrix_complex * B,const gsl_complex beta,gsl_matrix_complex * C)1866 gsl_blas_zsyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
1867 const gsl_complex alpha, const gsl_matrix_complex * A,
1868 const gsl_matrix_complex * B, const gsl_complex beta,
1869 gsl_matrix_complex * C)
1870 {
1871 const size_t M = C->size1;
1872 const size_t N = C->size2;
1873 const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1874 const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1875 const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
1876 const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
1877
1878 if (M != N)
1879 {
1880 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1881 }
1882 else if (N != MA || N != MB || NA != NB)
1883 {
1884 GSL_ERROR ("invalid length", GSL_EBADLEN);
1885 }
1886
1887 cblas_zsyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
1888 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1889 INT (B->tda), GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
1890 return GSL_SUCCESS;
1891 }
1892
1893 /* HER2K */
1894
1895 int
gsl_blas_cher2k(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,const gsl_matrix_complex_float * B,float beta,gsl_matrix_complex_float * C)1896 gsl_blas_cher2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
1897 const gsl_complex_float alpha,
1898 const gsl_matrix_complex_float * A,
1899 const gsl_matrix_complex_float * B, float beta,
1900 gsl_matrix_complex_float * C)
1901 {
1902 const size_t M = C->size1;
1903 const size_t N = C->size2;
1904 const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1905 const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1906 const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
1907 const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
1908
1909 if (M != N)
1910 {
1911 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1912 }
1913 else if (N != MA || N != MB || NA != NB)
1914 {
1915 GSL_ERROR ("invalid length", GSL_EBADLEN);
1916 }
1917
1918 cblas_cher2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
1919 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1920 INT (B->tda), beta, C->data, INT (C->tda));
1921 return GSL_SUCCESS;
1922
1923 }
1924
1925
1926 int
gsl_blas_zher2k(CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t Trans,const gsl_complex alpha,const gsl_matrix_complex * A,const gsl_matrix_complex * B,double beta,gsl_matrix_complex * C)1927 gsl_blas_zher2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
1928 const gsl_complex alpha, const gsl_matrix_complex * A,
1929 const gsl_matrix_complex * B, double beta,
1930 gsl_matrix_complex * C)
1931 {
1932 const size_t M = C->size1;
1933 const size_t N = C->size2;
1934 const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
1935 const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
1936 const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
1937 const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;
1938
1939 if (M != N)
1940 {
1941 GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
1942 }
1943 else if (N != MA || N != MB || NA != NB)
1944 {
1945 GSL_ERROR ("invalid length", GSL_EBADLEN);
1946 }
1947
1948 cblas_zher2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
1949 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
1950 INT (B->tda), beta, C->data, INT (C->tda));
1951 return GSL_SUCCESS;
1952
1953 }
1954
1955 /* TRMM */
1956
1957 int
gsl_blas_strmm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,float alpha,const gsl_matrix_float * A,gsl_matrix_float * B)1958 gsl_blas_strmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
1959 CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, float alpha,
1960 const gsl_matrix_float * A, gsl_matrix_float * B)
1961 {
1962 const size_t M = B->size1;
1963 const size_t N = B->size2;
1964 const size_t MA = A->size1;
1965 const size_t NA = A->size2;
1966
1967 if (MA != NA)
1968 {
1969 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
1970 }
1971
1972 if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
1973 {
1974 cblas_strmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
1975 alpha, A->data, INT (A->tda), B->data, INT (B->tda));
1976 return GSL_SUCCESS;
1977 }
1978 else
1979 {
1980 GSL_ERROR ("invalid length", GSL_EBADLEN);
1981 }
1982 }
1983
1984
1985 int
gsl_blas_dtrmm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,double alpha,const gsl_matrix * A,gsl_matrix * B)1986 gsl_blas_dtrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
1987 CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, double alpha,
1988 const gsl_matrix * A, gsl_matrix * B)
1989 {
1990 const size_t M = B->size1;
1991 const size_t N = B->size2;
1992 const size_t MA = A->size1;
1993 const size_t NA = A->size2;
1994
1995 if (MA != NA)
1996 {
1997 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
1998 }
1999
2000 if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
2001 {
2002 cblas_dtrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
2003 alpha, A->data, INT (A->tda), B->data, INT (B->tda));
2004 return GSL_SUCCESS;
2005 }
2006 else
2007 {
2008 GSL_ERROR ("invalid length", GSL_EBADLEN);
2009 }
2010 }
2011
2012
2013 int
gsl_blas_ctrmm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,gsl_matrix_complex_float * B)2014 gsl_blas_ctrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
2015 CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
2016 const gsl_complex_float alpha,
2017 const gsl_matrix_complex_float * A,
2018 gsl_matrix_complex_float * B)
2019 {
2020 const size_t M = B->size1;
2021 const size_t N = B->size2;
2022 const size_t MA = A->size1;
2023 const size_t NA = A->size2;
2024
2025 if (MA != NA)
2026 {
2027 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
2028 }
2029
2030 if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
2031 {
2032 cblas_ctrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
2033 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
2034 INT (B->tda));
2035 return GSL_SUCCESS;
2036 }
2037 else
2038 {
2039 GSL_ERROR ("invalid length", GSL_EBADLEN);
2040 }
2041 }
2042
2043
2044 int
gsl_blas_ztrmm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_complex alpha,const gsl_matrix_complex * A,gsl_matrix_complex * B)2045 gsl_blas_ztrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
2046 CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
2047 const gsl_complex alpha, const gsl_matrix_complex * A,
2048 gsl_matrix_complex * B)
2049 {
2050 const size_t M = B->size1;
2051 const size_t N = B->size2;
2052 const size_t MA = A->size1;
2053 const size_t NA = A->size2;
2054
2055 if (MA != NA)
2056 {
2057 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
2058 }
2059
2060 if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
2061 {
2062 cblas_ztrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
2063 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
2064 INT (B->tda));
2065 return GSL_SUCCESS;
2066 }
2067 else
2068 {
2069 GSL_ERROR ("invalid length", GSL_EBADLEN);
2070 }
2071 }
2072
2073
2074 /* TRSM */
2075
2076 int
gsl_blas_strsm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,float alpha,const gsl_matrix_float * A,gsl_matrix_float * B)2077 gsl_blas_strsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
2078 CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, float alpha,
2079 const gsl_matrix_float * A, gsl_matrix_float * B)
2080 {
2081 const size_t M = B->size1;
2082 const size_t N = B->size2;
2083 const size_t MA = A->size1;
2084 const size_t NA = A->size2;
2085
2086 if (MA != NA)
2087 {
2088 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
2089 }
2090
2091 if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
2092 {
2093 cblas_strsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
2094 alpha, A->data, INT (A->tda), B->data, INT (B->tda));
2095 return GSL_SUCCESS;
2096 }
2097 else
2098 {
2099 GSL_ERROR ("invalid length", GSL_EBADLEN);
2100 }
2101 }
2102
2103
2104 int
gsl_blas_dtrsm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,double alpha,const gsl_matrix * A,gsl_matrix * B)2105 gsl_blas_dtrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
2106 CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, double alpha,
2107 const gsl_matrix * A, gsl_matrix * B)
2108 {
2109 const size_t M = B->size1;
2110 const size_t N = B->size2;
2111 const size_t MA = A->size1;
2112 const size_t NA = A->size2;
2113
2114 if (MA != NA)
2115 {
2116 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
2117 }
2118
2119 if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
2120 {
2121 cblas_dtrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
2122 alpha, A->data, INT (A->tda), B->data, INT (B->tda));
2123 return GSL_SUCCESS;
2124 }
2125 else
2126 {
2127 GSL_ERROR ("invalid length", GSL_EBADLEN);
2128 }
2129 }
2130
2131
2132 int
gsl_blas_ctrsm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_complex_float alpha,const gsl_matrix_complex_float * A,gsl_matrix_complex_float * B)2133 gsl_blas_ctrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
2134 CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
2135 const gsl_complex_float alpha,
2136 const gsl_matrix_complex_float * A,
2137 gsl_matrix_complex_float * B)
2138 {
2139 const size_t M = B->size1;
2140 const size_t N = B->size2;
2141 const size_t MA = A->size1;
2142 const size_t NA = A->size2;
2143
2144 if (MA != NA)
2145 {
2146 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
2147 }
2148
2149 if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
2150 {
2151 cblas_ctrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
2152 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
2153 INT (B->tda));
2154 return GSL_SUCCESS;
2155 }
2156 else
2157 {
2158 GSL_ERROR ("invalid length", GSL_EBADLEN);
2159 }
2160 }
2161
2162
2163 int
gsl_blas_ztrsm(CBLAS_SIDE_t Side,CBLAS_UPLO_t Uplo,CBLAS_TRANSPOSE_t TransA,CBLAS_DIAG_t Diag,const gsl_complex alpha,const gsl_matrix_complex * A,gsl_matrix_complex * B)2164 gsl_blas_ztrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
2165 CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
2166 const gsl_complex alpha, const gsl_matrix_complex * A,
2167 gsl_matrix_complex * B)
2168 {
2169 const size_t M = B->size1;
2170 const size_t N = B->size2;
2171 const size_t MA = A->size1;
2172 const size_t NA = A->size2;
2173
2174 if (MA != NA)
2175 {
2176 GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
2177 }
2178
2179 if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
2180 {
2181 cblas_ztrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
2182 GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
2183 INT (B->tda));
2184 return GSL_SUCCESS;
2185 }
2186 else
2187 {
2188 GSL_ERROR ("invalid length", GSL_EBADLEN);
2189 }
2190 }
2191