/* blas/blas.c

 * 

 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2009 Gerard Jungman & Brian 

 * Gough

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

/* GSL implementation of BLAS operations for vectors and dense

 * matrices.  Note that GSL native storage is row-major.  */

#include <config.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_cblas.h>

#include <gsl/gsl_cblas.h>

#include <gsl/gsl_blas_types.h>

#include <gsl/gsl_blas.h>

/* ========================================================================

 * Level 1

 * ========================================================================

 */

/* CBLAS defines vector sizes in terms of int. GSL defines sizes in

   terms of size_t, so we need to convert these into integers.  There

   is the possibility of overflow here. FIXME: Maybe this could be

   caught */

#define INT(X) ((int)(X))

int

gsl_blas_sdsdot (float alpha, const gsl_vector_float * X,

                 const gsl_vector_float * Y, float *result)

{

  if (X->size == Y->size)

    {

      *result =

        cblas_sdsdot (INT (X->size), alpha, X->data, INT (X->stride), Y->data,

                      INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_dsdot (const gsl_vector_float * X, const gsl_vector_float * Y,

                double *result)

{

  if (X->size == Y->size)

    {

      *result =

        cblas_dsdot (INT (X->size), X->data, INT (X->stride), Y->data,

                     INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_sdot (const gsl_vector_float * X, const gsl_vector_float * Y,

               float *result)

{

  if (X->size == Y->size)

    {

      *result =

        cblas_sdot (INT (X->size), X->data, INT (X->stride), Y->data,

                    INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_ddot (const gsl_vector * X, const gsl_vector * Y, double *result)

{

  if (X->size == Y->size)

    {

      *result =

        cblas_ddot (INT (X->size), X->data, INT (X->stride), Y->data,

                    INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_cdotu (const gsl_vector_complex_float * X,

                const gsl_vector_complex_float * Y, gsl_complex_float * dotu)

{

  if (X->size == Y->size)

    {

      cblas_cdotu_sub (INT (X->size), X->data, INT (X->stride), Y->data,

                       INT (Y->stride), GSL_COMPLEX_P (dotu));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_cdotc (const gsl_vector_complex_float * X,

                const gsl_vector_complex_float * Y, gsl_complex_float * dotc)

{

  if (X->size == Y->size)

    {

      cblas_cdotc_sub (INT (X->size), X->data, INT (X->stride), Y->data,

                       INT (Y->stride), GSL_COMPLEX_P (dotc));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_zdotu (const gsl_vector_complex * X, const gsl_vector_complex * Y,

                gsl_complex * dotu)

{

  if (X->size == Y->size)

    {

      cblas_zdotu_sub (INT (X->size), X->data, INT (X->stride), Y->data,

                       INT (Y->stride), GSL_COMPLEX_P (dotu));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_zdotc (const gsl_vector_complex * X, const gsl_vector_complex * Y,

                gsl_complex * dotc)

{

  if (X->size == Y->size)

    {

      cblas_zdotc_sub (INT (X->size), X->data, INT (X->stride), Y->data,

                       INT (Y->stride), GSL_COMPLEX_P (dotc));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

/* Norms of vectors */

float

gsl_blas_snrm2 (const gsl_vector_float * X)

{

  return cblas_snrm2 (INT (X->size), X->data, INT (X->stride));

}

double

gsl_blas_dnrm2 (const gsl_vector * X)

{

  return cblas_dnrm2 (INT (X->size), X->data, INT (X->stride));

}

float

gsl_blas_scnrm2 (const gsl_vector_complex_float * X)

{

  return cblas_scnrm2 (INT (X->size), X->data, INT (X->stride));

}

double

gsl_blas_dznrm2 (const gsl_vector_complex * X)

{

  return cblas_dznrm2 (INT (X->size), X->data, INT (X->stride));

}

/* Absolute sums of vectors */

float

gsl_blas_sasum (const gsl_vector_float * X)

{

  return cblas_sasum (INT (X->size), X->data, INT (X->stride));

}

double

gsl_blas_dasum (const gsl_vector * X)

{

  return cblas_dasum (INT (X->size), X->data, INT (X->stride));

}

float

gsl_blas_scasum (const gsl_vector_complex_float * X)

{

  return cblas_scasum (INT (X->size), X->data, INT (X->stride));

}

double

gsl_blas_dzasum (const gsl_vector_complex * X)

{

  return cblas_dzasum (INT (X->size), X->data, INT (X->stride));

}

/* Maximum elements of vectors */

CBLAS_INDEX_t

gsl_blas_isamax (const gsl_vector_float * X)

{

  return cblas_isamax (INT (X->size), X->data, INT (X->stride));

}

CBLAS_INDEX_t

gsl_blas_idamax (const gsl_vector * X)

{

  return cblas_idamax (INT (X->size), X->data, INT (X->stride));

}

CBLAS_INDEX_t

gsl_blas_icamax (const gsl_vector_complex_float * X)

{

  return cblas_icamax (INT (X->size), X->data, INT (X->stride));

}

CBLAS_INDEX_t

gsl_blas_izamax (const gsl_vector_complex * X)

{

  return cblas_izamax (INT (X->size), X->data, INT (X->stride));

}

/* Swap vectors */

int

gsl_blas_sswap (gsl_vector_float * X, gsl_vector_float * Y)

{

  if (X->size == Y->size)

    {

      cblas_sswap (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_dswap (gsl_vector * X, gsl_vector * Y)

{

  if (X->size == Y->size)

    {

      cblas_dswap (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    };

}

int

gsl_blas_cswap (gsl_vector_complex_float * X, gsl_vector_complex_float * Y)

{

  if (X->size == Y->size)

    {

      cblas_cswap (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_zswap (gsl_vector_complex * X, gsl_vector_complex * Y)

{

  if (X->size == Y->size)

    {

      cblas_zswap (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

/* Copy vectors */

int

gsl_blas_scopy (const gsl_vector_float * X, gsl_vector_float * Y)

{

  if (X->size == Y->size)

    {

      cblas_scopy (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_dcopy (const gsl_vector * X, gsl_vector * Y)

{

  if (X->size == Y->size)

    {

      cblas_dcopy (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_ccopy (const gsl_vector_complex_float * X,

                gsl_vector_complex_float * Y)

{

  if (X->size == Y->size)

    {

      cblas_ccopy (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_zcopy (const gsl_vector_complex * X, gsl_vector_complex * Y)

{

  if (X->size == Y->size)

    {

      cblas_zcopy (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

/* Compute Y = alpha X + Y */

int

gsl_blas_saxpy (float alpha, const gsl_vector_float * X, gsl_vector_float * Y)

{

  if (X->size == Y->size)

    {

      cblas_saxpy (INT (X->size), alpha, X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_daxpy (double alpha, const gsl_vector * X, gsl_vector * Y)

{

  if (X->size == Y->size)

    {

      cblas_daxpy (INT (X->size), alpha, X->data, INT (X->stride), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_caxpy (const gsl_complex_float alpha,

                const gsl_vector_complex_float * X,

                gsl_vector_complex_float * Y)

{

  if (X->size == Y->size)

    {

      cblas_caxpy (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,

                   INT (X->stride), Y->data, INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_zaxpy (const gsl_complex alpha, const gsl_vector_complex * X,

                gsl_vector_complex * Y)

{

  if (X->size == Y->size)

    {

      cblas_zaxpy (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,

                   INT (X->stride), Y->data, INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

/* Generate rotation */

int

gsl_blas_srotg (float a[], float b[], float c[], float s[])

{

  cblas_srotg (a, b, c, s);

  return GSL_SUCCESS;

}

int

gsl_blas_drotg (double a[], double b[], double c[], double s[])

{

  cblas_drotg (a, b, c, s);

  return GSL_SUCCESS;

}

/* Apply rotation to vectors */

int

gsl_blas_srot (gsl_vector_float * X, gsl_vector_float * Y, float c, float s)

{

  if (X->size == Y->size)

    {

      cblas_srot (INT (X->size), X->data, INT (X->stride), Y->data,

                  INT (Y->stride), c, s);

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_drot (gsl_vector * X, gsl_vector * Y, const double c, const double s)

{

  if (X->size == Y->size)

    {

      cblas_drot (INT (X->size), X->data, INT (X->stride), Y->data,

                  INT (Y->stride), c, s);

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

/* Generate modified rotation */

int

gsl_blas_srotmg (float d1[], float d2[], float b1[], float b2, float P[])

{

  cblas_srotmg (d1, d2, b1, b2, P);

  return GSL_SUCCESS;

}

int

gsl_blas_drotmg (double d1[], double d2[], double b1[], double b2, double P[])

{

  cblas_drotmg (d1, d2, b1, b2, P);

  return GSL_SUCCESS;

}

/* Apply modified rotation */

int

gsl_blas_srotm (gsl_vector_float * X, gsl_vector_float * Y, const float P[])

{

  if (X->size == Y->size)

    {

      cblas_srotm (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride), P);

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_drotm (gsl_vector * X, gsl_vector * Y, const double P[])

{

  if (X->size == Y->size)

    {

      cblas_drotm (INT (X->size), X->data, INT (X->stride), Y->data,

                   INT (Y->stride), P);

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

/* Scale vector */

void

gsl_blas_sscal (float alpha, gsl_vector_float * X)

{

  cblas_sscal (INT (X->size), alpha, X->data, INT (X->stride));

}

void

gsl_blas_dscal (double alpha, gsl_vector * X)

{

  cblas_dscal (INT (X->size), alpha, X->data, INT (X->stride));

}

void

gsl_blas_cscal (const gsl_complex_float alpha, gsl_vector_complex_float * X)

{

  cblas_cscal (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,

               INT (X->stride));

}

void

gsl_blas_zscal (const gsl_complex alpha, gsl_vector_complex * X)

{

  cblas_zscal (INT (X->size), GSL_COMPLEX_P (&alpha), X->data,

               INT (X->stride));

}

void

gsl_blas_csscal (float alpha, gsl_vector_complex_float * X)

{

  cblas_csscal (INT (X->size), alpha, X->data, INT (X->stride));

}

void

gsl_blas_zdscal (double alpha, gsl_vector_complex * X)

{

  cblas_zdscal (INT (X->size), alpha, X->data, INT (X->stride));

}

/* ===========================================================================

 * Level 2

 * ===========================================================================

 */

/* GEMV */

int

gsl_blas_sgemv (CBLAS_TRANSPOSE_t TransA, float alpha,

                const gsl_matrix_float * A, const gsl_vector_float * X,

                float beta, gsl_vector_float * Y)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if ((TransA == CblasNoTrans && N == X->size && M == Y->size)

      || (TransA == CblasTrans && M == X->size && N == Y->size))

    {

      cblas_sgemv (CblasRowMajor, TransA, INT (M), INT (N), alpha, A->data,

                   INT (A->tda), X->data, INT (X->stride), beta, Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA, double alpha, const gsl_matrix * A,

                const gsl_vector * X, double beta, gsl_vector * Y)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if ((TransA == CblasNoTrans && N == X->size && M == Y->size)

      || (TransA == CblasTrans && M == X->size && N == Y->size))

    {

      cblas_dgemv (CblasRowMajor, TransA, INT (M), INT (N), alpha, A->data,

                   INT (A->tda), X->data, INT (X->stride), beta, Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_cgemv (CBLAS_TRANSPOSE_t TransA, const gsl_complex_float alpha,

                const gsl_matrix_complex_float * A,

                const gsl_vector_complex_float * X,

                const gsl_complex_float beta, gsl_vector_complex_float * Y)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if ((TransA == CblasNoTrans && N == X->size && M == Y->size)

      || (TransA == CblasTrans && M == X->size && N == Y->size)

      || (TransA == CblasConjTrans && M == X->size && N == Y->size))

    {

      cblas_cgemv (CblasRowMajor, TransA, INT (M), INT (N),

                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), X->data,

                   INT (X->stride), GSL_COMPLEX_P (&beta), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

int

gsl_blas_zgemv (CBLAS_TRANSPOSE_t TransA, const gsl_complex alpha,

                const gsl_matrix_complex * A, const gsl_vector_complex * X,

                const gsl_complex beta, gsl_vector_complex * Y)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if ((TransA == CblasNoTrans && N == X->size && M == Y->size)

      || (TransA == CblasTrans && M == X->size && N == Y->size)

      || (TransA == CblasConjTrans && M == X->size && N == Y->size))

    {

      cblas_zgemv (CblasRowMajor, TransA, INT (M), INT (N),

                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), X->data,

                   INT (X->stride), GSL_COMPLEX_P (&beta), Y->data,

                   INT (Y->stride));

      return GSL_SUCCESS;

    }

  else

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

}

/* HEMV */

int

gsl_blas_chemv (CBLAS_UPLO_t Uplo, const gsl_complex_float alpha,

                const gsl_matrix_complex_float * A,

                const gsl_vector_complex_float * X,

                const gsl_complex_float beta, gsl_vector_complex_float * Y)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if (M != N)

    {

      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);

    }

  else if (N != X->size || N != Y->size)

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

  cblas_chemv (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), A->data,

               INT (A->tda), X->data, INT (X->stride), GSL_COMPLEX_P (&beta),

               Y->data, INT (Y->stride));

  return GSL_SUCCESS;

}

int

gsl_blas_zhemv (CBLAS_UPLO_t Uplo, const gsl_complex alpha,

                const gsl_matrix_complex * A, const gsl_vector_complex * X,

                const gsl_complex beta, gsl_vector_complex * Y)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if (M != N)

    {

      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);

    }

  else if (N != X->size || N != Y->size)

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

  cblas_zhemv (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), A->data,

               INT (A->tda), X->data, INT (X->stride), GSL_COMPLEX_P (&beta),

               Y->data, INT (Y->stride));

  return GSL_SUCCESS;

}

/* SYMV */

int

gsl_blas_ssymv (CBLAS_UPLO_t Uplo, float alpha, const gsl_matrix_float * A,

                const gsl_vector_float * X, float beta, gsl_vector_float * Y)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if (M != N)

    {

      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);

    }

  else if (N != X->size || N != Y->size)

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

  cblas_ssymv (CblasRowMajor, Uplo, INT (N), alpha, A->data, INT (A->tda),

               X->data, INT (X->stride), beta, Y->data, INT (Y->stride));

  return GSL_SUCCESS;

}

int

gsl_blas_dsymv (CBLAS_UPLO_t Uplo, double alpha, const gsl_matrix * A,

                const gsl_vector * X, double beta, gsl_vector * Y)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if (M != N)

    {

      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);

    }

  else if (N != X->size || N != Y->size)

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

  cblas_dsymv (CblasRowMajor, Uplo, INT (N), alpha, A->data, INT (A->tda),

               X->data, INT (X->stride), beta, Y->data, INT (Y->stride));

  return GSL_SUCCESS;

}

/* TRMV */

int

gsl_blas_strmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,

                CBLAS_DIAG_t Diag, const gsl_matrix_float * A,

                gsl_vector_float * X)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if (M != N)

    {

      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);

    }

  else if (N != X->size)

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

  cblas_strmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,

               INT (A->tda), X->data, INT (X->stride));

  return GSL_SUCCESS;

}

int

gsl_blas_dtrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,

                CBLAS_DIAG_t Diag, const gsl_matrix * A, gsl_vector * X)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if (M != N)

    {

      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);

    }

  else if (N != X->size)

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

  cblas_dtrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,

               INT (A->tda), X->data, INT (X->stride));

  return GSL_SUCCESS;

}

int

gsl_blas_ctrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,

                CBLAS_DIAG_t Diag, const gsl_matrix_complex_float * A,

                gsl_vector_complex_float * X)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if (M != N)

    {

      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);

    }

  else if (N != X->size)

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

  cblas_ctrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,

               INT (A->tda), X->data, INT (X->stride));

  return GSL_SUCCESS;

}

int

gsl_blas_ztrmv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,

                CBLAS_DIAG_t Diag, const gsl_matrix_complex * A,

                gsl_vector_complex * X)

{

  const size_t M = A->size1;

  const size_t N = A->size2;

  if (M != N)

    {

      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);

    }

  else if (N != X->size)

    {

      GSL_ERROR ("invalid length", GSL_EBADLEN);

    }

  cblas_ztrmv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,

               INT (A->tda), X->data, INT (X->stride));

  return GSL_SUCCESS;

}