Annotated Source Code

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
44int
45gsl_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
61int
62gsl_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
78int
79gsl_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
95int
96gsl_blas_ddot (const gsl_vector * X, const gsl_vector * Y, double *result)
5Enter function.    gsl_blas_ddot(p, gradient, &pg)
97{
98  if (X->size == Y->size)
6Take the false branch.    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);
7Call a function.    gsl_error("invalid length", "/home/xingming/workplace/experiment/fp/gsl-2.1/blas/blas.c", 107, GSL_EBADLEN)
108    }
109}
11Exit function.    gsl_blas_ddot(p, gradient, &pg)
110
111
112int
113gsl_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
129int
130gsl_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
146int
147gsl_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
163int
164gsl_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
181float
182gsl_blas_snrm2 (const gsl_vector_float * X)
183{
184  return cblas_snrm2 (INT (X->size), X->data, INT (X->stride));
185}
186
187double
188gsl_blas_dnrm2 (const gsl_vector * X)
189{
190  return cblas_dnrm2 (INT (X->size), X->data, INT (X->stride));
191}
192
193float
194gsl_blas_scnrm2 (const gsl_vector_complex_float * X)
195{
196  return cblas_scnrm2 (INT (X->size), X->data, INT (X->stride));
197}
198
199double
200gsl_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
207float
208gsl_blas_sasum (const gsl_vector_float * X)
209{
210  return cblas_sasum (INT (X->size), X->data, INT (X->stride));
211}
212
213double
214gsl_blas_dasum (const gsl_vector * X)
215{
216  return cblas_dasum (INT (X->size), X->data, INT (X->stride));
217}
218
219float
220gsl_blas_scasum (const gsl_vector_complex_float * X)
221{
222  return cblas_scasum (INT (X->size), X->data, INT (X->stride));
223}
224
225double
226gsl_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
233CBLAS_INDEX_t
234gsl_blas_isamax (const gsl_vector_float * X)
235{
236  return cblas_isamax (INT (X->size), X->data, INT (X->stride));
237}
238
239CBLAS_INDEX_t
240gsl_blas_idamax (const gsl_vector * X)
241{
242  return cblas_idamax (INT (X->size), X->data, INT (X->stride));
243}
244
245CBLAS_INDEX_t
246gsl_blas_icamax (const gsl_vector_complex_float * X)
247{
248  return cblas_icamax (INT (X->size), X->data, INT (X->stride));
249}
250
251CBLAS_INDEX_t
252gsl_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
260int
261gsl_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
275int
276gsl_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
290int
291gsl_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
305int
306gsl_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
323int
324gsl_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
338int
339gsl_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
353int
354gsl_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
369int
370gsl_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
387int
388gsl_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
402int
403gsl_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
417int
418gsl_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
434int
435gsl_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
452int
453gsl_blas_srotg (float a[], float b[], float c[], float s[])
454{
455  cblas_srotg (a, b, c, s);
456  return GSL_SUCCESS;
457}
458
459int
460gsl_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
468int
469gsl_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
483int
484gsl_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
501int
502gsl_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
508int
509gsl_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
518int
519gsl_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
533int
534gsl_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
551void
552gsl_blas_sscal (float alpha, gsl_vector_float * X)
553{
554  cblas_sscal (INT (X->size), alpha, X->data, INT (X->stride));
555}
556
557void
558gsl_blas_dscal (double alpha, gsl_vector * X)
559{
560  cblas_dscal (INT (X->size), alpha, X->data, INT (X->stride));
561}
562
563void
564gsl_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
570void
571gsl_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
577void
578gsl_blas_csscal (float alpha, gsl_vector_complex_float * X)
579{
580  cblas_csscal (INT (X->size), alpha, X->data, INT (X->stride));
581}
582
583void
584gsl_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
596int
597gsl_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
619int
620gsl_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
641int
642gsl_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
667int
668gsl_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
695int
696gsl_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
719int
720gsl_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
745int
746gsl_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
766int
767gsl_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
790int
791gsl_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
813int
814gsl_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
835int
836gsl_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
858int
859gsl_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
883int
884gsl_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
906int
907gsl_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
928int
929gsl_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
951int
952gsl_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
976int
977gsl_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
997int
998gsl_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
1020int
1021gsl_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
1042int
1043gsl_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
1065int
1066gsl_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
1088int
1089gsl_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
1110int
1111gsl_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
1133int
1134gsl_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
1157int
1158gsl_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
1182int
1183gsl_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
1208int
1209gsl_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
1230int
1231gsl_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
1254int
1255gsl_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
1276int
1277gsl_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
1307int
1308gsl_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
1333int
1334gsl_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
1359int
1360gsl_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
1388int
1389gsl_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
1418int
1419gsl_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
1451int
1452gsl_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
1483int
1484gsl_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
1517int
1518gsl_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
1553int
1554gsl_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
1589int
1590gsl_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
1624int
1625gsl_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
1648int
1649gsl_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
1673int
1674gsl_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
1700int
1701gsl_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
1727int
1728gsl_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
1752int
1753gsl_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
1778int
1779gsl_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
1806int
1807gsl_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
1834int
1835gsl_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
1865int
1866gsl_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
1895int
1896gsl_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
1926int
1927gsl_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
1957int
1958gsl_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
1985int
1986gsl_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
2013int
2014gsl_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
2044int
2045gsl_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
2076int
2077gsl_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
2104int
2105gsl_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
2132int
2133gsl_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
2163int
2164gsl_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
Events list
Event 5
Event 6
Event 7
Event 11