blas.c

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

/* blas/blas.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001 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 2 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., 675 Mass Ave, Cambridge, MA 02139, 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);
    }
}