blas.c

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

/* 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;
}


/* TRSV */

int
gsl_blas_strsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix_float * A,
                gsl_vector_float * X)
{
  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_strsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}


int
gsl_blas_dtrsv (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_dtrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}


int
gsl_blas_ctrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix_complex_float * A,
                gsl_vector_complex_float * X)
{
  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_ctrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}


int
gsl_blas_ztrsv (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                CBLAS_DIAG_t Diag, const gsl_matrix_complex * A,
                gsl_vector_complex * X)
{
  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_ztrsv (CblasRowMajor, Uplo, TransA, Diag, INT (N), A->data,
               INT (A->tda), X->data, INT (X->stride));
  return GSL_SUCCESS;
}


/* GER */

int
gsl_blas_sger (float alpha, const gsl_vector_float * X,
               const gsl_vector_float * Y, gsl_matrix_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_sger (CblasRowMajor, INT (M), INT (N), alpha, X->data,
                  INT (X->stride), Y->data, INT (Y->stride), A->data,
                  INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_dger (double alpha, const gsl_vector * X, const gsl_vector * Y,
               gsl_matrix * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_dger (CblasRowMajor, INT (M), INT (N), alpha, X->data,
                  INT (X->stride), Y->data, INT (Y->stride), A->data,
                  INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


/* GERU */

int
gsl_blas_cgeru (const gsl_complex_float alpha,
                const gsl_vector_complex_float * X,
                const gsl_vector_complex_float * Y,
                gsl_matrix_complex_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_cgeru (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
                   X->data, INT (X->stride), Y->data, INT (Y->stride),
                   A->data, INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zgeru (const gsl_complex alpha, const gsl_vector_complex * X,
                const gsl_vector_complex * Y, gsl_matrix_complex * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_zgeru (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
                   X->data, INT (X->stride), Y->data, INT (Y->stride),
                   A->data, INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


/* GERC */

int
gsl_blas_cgerc (const gsl_complex_float alpha,
                const gsl_vector_complex_float * X,
                const gsl_vector_complex_float * Y,
                gsl_matrix_complex_float * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (X->size == M && Y->size == N)
    {
      cblas_cgerc (CblasRowMajor, INT (M), INT (N), GSL_COMPLEX_P (&alpha),
                   X->data, INT (X->stride), Y->data, INT (Y->stride),
                   A->data, INT (A->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_zgerc (const gsl_complex alpha, const gsl_vector_complex * X,
                const gsl_vector_complex * Y, gsl_matrix_complex * A)
{
  const size_t M = A->size1;
  const size_t N = A->size2;