htswanalysis/MACS/lib/gsl/gsl-1.11/blas/blas.c

   1 /* blas/blas.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001 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
  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
  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
  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
  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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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 }