blas.c

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

  if (X->size == M && Y->size == N)
    {
      cblas_zgerc (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);
    }
}

/* HER */

int
gsl_blas_cher (CBLAS_UPLO_t Uplo, float alpha,
               const gsl_vector_complex_float * X,
               gsl_matrix_complex_float * A)
{
  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 (X->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_cher (CblasRowMajor, Uplo, INT (M), alpha, X->data, INT (X->stride),
              A->data, INT (A->tda));
  return GSL_SUCCESS;
}


int
gsl_blas_zher (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector_complex * X,
               gsl_matrix_complex * A)
{
  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 (X->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zher (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
              A->data, INT (A->tda));
  return GSL_SUCCESS;
}


/* HER2 */

int
gsl_blas_cher2 (CBLAS_UPLO_t Uplo, const gsl_complex_float alpha,
                const gsl_vector_complex_float * X,
                const gsl_vector_complex_float * Y,
                gsl_matrix_complex_float * A)
{
  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 (X->size != N || Y->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_cher2 (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), X->data,
               INT (X->stride), Y->data, INT (Y->stride), A->data,
               INT (A->tda));
  return GSL_SUCCESS;
}


int
gsl_blas_zher2 (CBLAS_UPLO_t Uplo, const gsl_complex alpha,
                const gsl_vector_complex * X, const gsl_vector_complex * Y,
                gsl_matrix_complex * A)
{
  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 (X->size != N || Y->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zher2 (CblasRowMajor, Uplo, INT (N), GSL_COMPLEX_P (&alpha), X->data,
               INT (X->stride), Y->data, INT (Y->stride), A->data,
               INT (A->tda));
  return GSL_SUCCESS;
}


/* SYR */

int
gsl_blas_ssyr (CBLAS_UPLO_t Uplo, float alpha, const gsl_vector_float * X,
               gsl_matrix_float * A)
{
  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 (X->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ssyr (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
              A->data, INT (A->tda));
  return GSL_SUCCESS;
}


int
gsl_blas_dsyr (CBLAS_UPLO_t Uplo, double alpha, const gsl_vector * X,
               gsl_matrix * A)
{
  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 (X->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dsyr (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
              A->data, INT (A->tda));
  return GSL_SUCCESS;
}


/* SYR2 */

int
gsl_blas_ssyr2 (CBLAS_UPLO_t Uplo, 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 (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N || Y->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ssyr2 (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
               Y->data, INT (Y->stride), A->data, INT (A->tda));
  return GSL_SUCCESS;
}


int
gsl_blas_dsyr2 (CBLAS_UPLO_t Uplo, 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 (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (X->size != N || Y->size != N)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dsyr2 (CblasRowMajor, Uplo, INT (N), alpha, X->data, INT (X->stride),
               Y->data, INT (Y->stride), A->data, INT (A->tda));
  return GSL_SUCCESS;
}


/*
 * ===========================================================================
 * Prototypes for level 3 BLAS
 * ===========================================================================
 */


/* GEMM */

int
gsl_blas_sgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
                float alpha, const gsl_matrix_float * A,
                const gsl_matrix_float * B, float beta, gsl_matrix_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;

  if (M == MA && N == NB && NA == MB)   /* [MxN] = [MAxNA][MBxNB] */
    {
      cblas_sgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
                   alpha, A->data, INT (A->tda), B->data, INT (B->tda), beta,
                   C->data, INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
                double alpha, const gsl_matrix * A, const gsl_matrix * B,
                double beta, gsl_matrix * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;

  if (M == MA && N == NB && NA == MB)   /* [MxN] = [MAxNA][MBxNB] */
    {
      cblas_dgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
                   alpha, A->data, INT (A->tda), B->data, INT (B->tda), beta,
                   C->data, INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_cgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
                const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                const gsl_matrix_complex_float * B,
                const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;

  if (M == MA && N == NB && NA == MB)   /* [MxN] = [MAxNA][MBxNB] */
    {
      cblas_cgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_zgemm (CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB,
                const gsl_complex alpha, const gsl_matrix_complex * A,
                const gsl_matrix_complex * B, const gsl_complex beta,
                gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (TransA == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (TransA == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (TransB == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (TransB == CblasNoTrans) ? B->size2 : B->size1;

  if (M == MA && N == NB && NA == MB)   /* [MxN] = [MAxNA][MBxNB] */
    {
      cblas_zgemm (CblasRowMajor, TransA, TransB, INT (M), INT (N), INT (NA),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


/* SYMM */

int
gsl_blas_ssymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, float alpha,
                const gsl_matrix_float * A, const gsl_matrix_float * B,
                float beta, gsl_matrix_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_ssymm (CblasRowMajor, Side, Uplo, INT (M), INT (N), alpha,
                   A->data, INT (A->tda), B->data, INT (B->tda), beta,
                   C->data, INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

}


int
gsl_blas_dsymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo, double alpha,
                const gsl_matrix * A, const gsl_matrix * B, double beta,
                gsl_matrix * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_dsymm (CblasRowMajor, Side, Uplo, INT (M), INT (N), alpha,
                   A->data, INT (A->tda), B->data, INT (B->tda), beta,
                   C->data, INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_csymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                const gsl_matrix_complex_float * B,
                const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_csymm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

int
gsl_blas_zsymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                const gsl_complex alpha, const gsl_matrix_complex * A,
                const gsl_matrix_complex * B, const gsl_complex beta,
                gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_zsymm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


/* HEMM */

int
gsl_blas_chemm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                const gsl_complex_float alpha,
                const gsl_matrix_complex_float * A,
                const gsl_matrix_complex_float * B,
                const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_chemm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

}


int
gsl_blas_zhemm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                const gsl_complex alpha, const gsl_matrix_complex * A,
                const gsl_matrix_complex * B, const gsl_complex beta,
                gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;
  const size_t MB = B->size1;
  const size_t NB = B->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && (M == MA && N == NB && NA == MB))
      || (Side == CblasRight && (M == MB && N == NA && NB == MA)))
    {
      cblas_zhemm (CblasRowMajor, Side, Uplo, INT (M), INT (N),
                   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
                   INT (B->tda), GSL_COMPLEX_P (&beta), C->data,
                   INT (C->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}

/* SYRK */

int
gsl_blas_ssyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
                const gsl_matrix_float * A, float beta, gsl_matrix_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != K)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ssyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,