blas_old.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <jwsc/config.h>

#include <stdlib.h>
#include <assert.h>

#include "jwsc/math/complex.h"
#include "jwsc/math/blas.h"


/* ===============  LEVEL 1 EXTERNS ==============*/

extern int sswap_(int*, float*, int*, float*, int*);
extern int dswap_(int*, double*, int*, double*, int*);
extern int cswap_(int*, Complex_f*, int*, Complex_f*, int*);
extern int zswap_(int*, Complex_d*, int*, Complex_d*, int*);
extern int sscal_(int*, float*, float*, int*);
extern int dscal_(int*, double*, double*, int*);
extern int cscal_(int*, Complex_f*, Complex_f*, int*);
extern int zscal_(int*, Complex_d*, Complex_d*, int*);
extern int scopy_(int*, float*, int*, float*, int*);
extern int dcopy_(int*, double*, int*, double*, int*);
extern int ccopy_(int*, Complex_f*, int*, Complex_f*, int*);
extern int zcopy_(int*, Complex_d*, int*, Complex_d*, int*);
extern int saxpy_(int*, float*, float*, int*, float*, int*);
extern int daxpy_(int*, double*, double*, int*, double*, int*);
extern int caxpy_(int*, Complex_f*, Complex_f*, int*, Complex_f*, int*);
extern int zaxpy_(int*, Complex_d*, Complex_d*, int*, Complex_d*, int*);
extern float sdot_(int*, float*, int*, float*, int*);
extern double ddot_(int*, double*, int*, double*, int*);
extern float snrm2(int*, float*, int*);
extern double dnrm2_(int*, double*, int*);
extern float sasum_(int*, float*, int*);
extern double dasum_(int*, double*, int*);
extern int isamax_(int*, float*, int*);
extern int idamax_(int*, double*, int*);
extern int icamax_(int*, Complex_f*, int*);
extern int izamax_(int*, Complex_d*, int*);


/* ===============  LEVEL 2 EXTERNS ==============*/

extern int sgemv_(char*, int*, int*, float*, float*, int*,
                  float*, int*, float*, float*, int*);
extern int dgemv_(char*, int*, int*, double*, double*, int*,
                  double*, int*, double*, double*, int*);
extern int cgemv_(char*, int*, int*, Complex_f*, Complex_f*, int*,
                  Complex_f*, int*, Complex_f*, Complex_f*, int*);
extern int zgemv_(char*, int*, int*, Complex_d*, Complex_d*, int*,
                  Complex_d*, int*, Complex_d*, Complex_d*, int*);


/* ===============  LEVEL 3 EXTERNS ==============*/

extern int sgemm_(char*, char*, int*, int*, int*, float*, float*, int*,
                  float*, int*, float*, float*, int*);
extern int dgemm_(char*, char*, int*, int*, int*, double*, double*, int*,
                  double*, int*, double*, double*, int*);
extern int cgemm_(char*, char*, int*, int*, int*, Complex_f*, Complex_f*, int*,
                  Complex_f*, int*, Complex_f*, Complex_f*, int*);
extern int zgemm_(char*, char*, int*, int*, int*, Complex_d*, Complex_d*, int*,
                  Complex_d*, int*, Complex_d*, Complex_d*, int*);




/* ===============  LEVEL 1 ==============*/


void blas_sswap(int n, float* x, int inc_x, float* y, int inc_y)
{
    sswap_(&n, x, &inc_x, y, &inc_y);
}

void blas_dswap(int n, double* x, int inc_x, double* y, int inc_y)
{
    dswap_(&n, x, &inc_x, y, &inc_y);
}

void blas_cswap(int n, Complex_f* x, int inc_x, Complex_f* y, int inc_y)
{
    cswap_(&n, x, &inc_x, y, &inc_y);
}

void blas_zswap(int n, Complex_d* x, int inc_x, Complex_d* y, int inc_y)
{
    zswap_(&n, x, &inc_x, y, &inc_y);
}

void blas_sscal(int n, float a, float* x, int inc_x)
{
    sscal_(&n, &a, x, &inc_x);
}

void blas_dscal(int n, double a, double* x, int inc_x)
{
    dscal_(&n, &a, x, &inc_x);
}

void blas_cscal(int n, Complex_f a, Complex_f* x, int inc_x)
{
    cscal_(&n, &a, x, &inc_x);
}

void blas_zscal(int n, Complex_d a, Complex_d* x, int inc_x)
{
    zscal_(&n, &a, x, &inc_x);
}

void blas_scopy(int n, const float* x, int inc_x, float* y, int inc_y)
{
    scopy_(&n, (float*)x, &inc_x, y, &inc_y);
}

void blas_dcopy(int n, const double* x, int inc_x, double* y, int inc_y)
{
    dcopy_(&n, (double*)x, &inc_x, y, &inc_y);
}

void blas_ccopy(int n, const Complex_f* x, int inc_x, Complex_f* y, int inc_y)
{
    ccopy_(&n, (Complex_f*)x, &inc_x, y, &inc_y);
}

void blas_zcopy(int n, const Complex_d* x, int inc_x, Complex_d* y, int inc_y)
{
    zcopy_(&n, (Complex_d*)x, &inc_x, y, &inc_y);
}

void blas_saxpy
(
    int          n, 
    float        a, 
    const float* x, 
    int          inc_x, 
    float*       y, 
    int          inc_y
)
{
    saxpy_(&n, &a, (float*)x, &inc_x, y, &inc_y);
}

void blas_daxpy
(
    int           n, 
    double        a, 
    const double* x, 
    int           inc_x, 
    double*       y, 
    int           inc_y
)
{
    daxpy_(&n, &a, (double*)x, &inc_x, y, &inc_y);
}

void blas_caxpy
(
    int              n, 
    Complex_f        a, 
    const Complex_f* x, 
    int              inc_x, 
    Complex_f*       y, 
    int              inc_y
)
{
    caxpy_(&n, &a, (Complex_f*)x, &inc_x, y, &inc_y);
}

void blas_zaxpy
(
    int              n, 
    Complex_d        a, 
    const Complex_d* x, 
    int              inc_x, 
    Complex_d*       y, 
    int              inc_y
)
{
    zaxpy_(&n, &a, (Complex_d*)x, &inc_x, y, &inc_y);
}

float blas_sdot(int n, const float* x, int inc_x, const float* y, int inc_y)
{
    return sdot_(&n, (float*)x, &inc_x, (float*)y, &inc_y);
}

double blas_ddot(int n, const double* x, int inc_x, const double* y, int inc_y)
{
    return ddot_(&n, (double*)x, &inc_x, (double*)y, &inc_y);
}

float blas_snrm2(int n, const float* x, int inc_x)
{
    return snrm2(&n, (float*)x, &inc_x);
}

double blas_dnrm2(int n, const double* x, int inc_x)
{
    return dnrm2_(&n, (double*)x, &inc_x);
}

float blas_sasum(int n, const float* x, int inc_x)
{
    return sasum_(&n, (float*)x, &inc_x);
}

double blas_dasum(int n, const double* x, int inc_x)
{
    return dasum_(&n, (double*)x, &inc_x);
}

int blas_isamax(int n, const float* x, int inc_x)
{
    return isamax_(&n, (float*)x, &inc_x);
}

int blas_idamax(int n, const double* x, int inc_x)
{
    return idamax_(&n, (double*)x, &inc_x);
}

int blas_icamax(int n, const Complex_f* x, int inc_x)
{
    return icamax_(&n, (Complex_f*)x, &inc_x);
}

int blas_izamax(int n, const Complex_d* x, int inc_x)
{
    return izamax_(&n, (Complex_d*)x, &inc_x);
}

/* ===============  LEVEL 2 PROTOTYPES ==============*/


void blas_sgemv
(
    char         trans, 
    int          m, 
    int          n, 
    float        alpha, 
    const float* A, 
    int          lda,
    const float* x,
    int          inc_x,
    float        beta,
    float*       y,
    int          inc_y
)
{
    trans = (trans == 'N' || trans == 'n') ? 'T' : 'N';
    sgemv_(&trans, &n, &m, &alpha, (float*)A, &lda, (float*)x, &inc_x,
            &beta, y, &inc_y);
}

void blas_dgemv
(
    char          trans, 
    int           m, 
    int           n, 
    double        alpha, 
    const double* A, 
    int           lda,
    const double* x,
    int           inc_x,
    double        beta,
    double*       y,
    int           inc_y
)
{
    trans = (trans == 'N' || trans == 'n') ? 'T' : 'N';
    dgemv_(&trans, &n, &m, &alpha, (double*)A, &lda, (double*)x, &inc_x,
            &beta, y, &inc_y);
}

void blas_cgemv
(
    char             trans, 
    int              m, 
    int              n, 
    Complex_f        alpha, 
    const Complex_f* A, 
    int              lda,
    const Complex_f* x,
    int              inc_x,
    Complex_f        beta,
    Complex_f*       y,
    int              inc_y
)
{
    trans = (trans == 'N' || trans == 'n') ? 'T' : 'N';
    cgemv_(&trans, &n, &m, &alpha, (Complex_f*)A, &lda, (Complex_f*)x, &inc_x,
            &beta, y, &inc_y);
}

void blas_zgemv
(
    char             trans, 
    int              m, 
    int              n, 
    Complex_d        alpha, 
    const Complex_d* A, 
    int              lda,
    const Complex_d* x,
    int              inc_x,
    Complex_d        beta,
    Complex_d*       y,
    int              inc_y
)
{
    trans = (trans == 'N' || trans == 'n') ? 'T' : 'N';
    zgemv_(&trans, &n, &m, &alpha, (Complex_d*)A, &lda, (Complex_d*)x, &inc_x,
            &beta, y, &inc_y);
}

/* ===============  LEVEL 3 ==============*/


void blas_sgemm
(
    char         trans_a, 
    char         trans_b, 
    int          m, 
    int          n, 
    int          k, 
    float        alpha, 
    const float* A, 
    int          lda,
    const float* B,
    int          ldb,
    float        beta,
    float*       C,
    int          ldc
)
{
    sgemm_(&trans_a, &trans_b, &n, &m, &k, &alpha, (float*) B, &ldb, 
            (float*) A, &lda, &beta, C, &ldc);
}

void blas_dgemm
(
    char          trans_a, 
    char          trans_b, 
    int           m, 
    int           n, 
    int           k, 
    double        alpha, 
    const double* A, 
    int           lda,
    const double* B,
    int           ldb,
    double        beta,
    double*       C,
    int           ldc
)
{
    dgemm_(&trans_a, &trans_b, &n, &m, &k, &alpha, (double*) B, &ldb, 
            (double*) A, &lda, &beta, C, &ldc);
}

void blas_cgemm
(
    char             trans_a, 
    char             trans_b, 
    int              m, 
    int              n, 
    int              k, 
    Complex_f        alpha, 
    const Complex_f* A, 
    int              lda,
    const Complex_f* B,
    int              ldb,
    Complex_f        beta,
    Complex_f*       C,
    int              ldc
)
{
    cgemm_(&trans_a, &trans_b, &n, &m, &k, &alpha, (Complex_f*) B, &ldb, 
            (Complex_f*) A, &lda, &beta, C, &ldc);
}

void blas_zgemm
(
    char             trans_a, 
    char             trans_b, 
    int              m, 
    int              n, 
    int              k, 
    Complex_d        alpha, 
    const Complex_d* A, 
    int              lda,
    const Complex_d* B,
    int              ldb,
    Complex_d        beta,
    Complex_d*       C,
    int              ldc
)
{
    zgemm_(&trans_a, &trans_b, &n, &m, &k, &alpha, (Complex_d*) B, &ldb, 
            (Complex_d*) A, &lda, &beta, C, &ldc);
}