blas.c

Go to the documentation of this file.

/*
 * 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, see <http://www.gnu.org/licenses/>.
 */
 
#include "blas.h"
 
static void
gemm_nn(const int M, const int N, const int K, const double ALPHA,
        const double *A, const int lda, const double *B, const int ldb,
        double *C, const int ldc)
{
    for (int i = 0; i < M; ++i) {
        for (int k = 0; k < K; ++k) {
            const double A_PART = ALPHA * A[i * lda + k];
            for (int j = 0; j < N; ++j) {
                C[i * ldc + j] += A_PART * B[k * ldb + j];
            }
        }
    }
}
 
static void
gemm_nt(const int M, const int N, const int K, const double ALPHA,
        const double *A, const int lda, const double *B, const int ldb,
        double *C, const int ldc)
{
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            double sum = 0;
            for (int k = 0; k < K; ++k) {
                sum += ALPHA * A[i * lda + k] * B[j * ldb + k];
            }
            C[i * ldc + j] += sum;
        }
    }
}
 
static void
gemm_tn(const int M, const int N, const int K, const double ALPHA,
        const double *A, const int lda, const double *B, const int ldb,
        double *C, const int ldc)
{
    for (int i = 0; i < M; ++i) {
        for (int k = 0; k < K; ++k) {
            const double A_PART = ALPHA * A[k * lda + i];
            for (int j = 0; j < N; ++j) {
                C[i * ldc + j] += A_PART * B[k * ldb + j];
            }
        }
    }
}
 
static void
gemm_tt(const int M, const int N, const int K, const double ALPHA,
        const double *A, const int lda, const double *B, const int ldb,
        double *C, const int ldc)
{
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            double sum = 0;
            for (int k = 0; k < K; ++k) {
                sum += ALPHA * A[i + k * lda] * B[k + j * ldb];
            }
            C[i * ldc + j] += sum;
        }
    }
}
 
void
blas_gemm(const int TA, const int TB, const int M, const int N, const int K,
          const double ALPHA, const double *A, const int lda, const double *B,
          const int ldb, const double BETA, double *C, const int ldc)
{
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            C[i * ldc + j] *= BETA;
        }
    }
    if (!TA && !TB) {
        gemm_nn(M, N, K, ALPHA, A, lda, B, ldb, C, ldc);
    } else if (TA && !TB) {
        gemm_tn(M, N, K, ALPHA, A, lda, B, ldb, C, ldc);
    } else if (!TA && TB) {
        gemm_nt(M, N, K, ALPHA, A, lda, B, ldb, C, ldc);
    } else {
        gemm_tt(M, N, K, ALPHA, A, lda, B, ldb, C, ldc);
    }
}
 
void
blas_axpy(const int N, const double ALPHA, const double *X, const int INCX,
          double *Y, const int INCY)
{
    if (ALPHA != 1) {
        for (int i = 0; i < N; ++i) {
            Y[i * INCY] += ALPHA * X[i * INCX];
        }
    } else {
        for (int i = 0; i < N; ++i) {
            Y[i * INCY] += X[i * INCX];
        }
    }
}
 
void
blas_scal(const int N, const double ALPHA, double *X, const int INCX)
{
    if (ALPHA != 0) {
        for (int i = 0; i < N; ++i) {
            X[i * INCX] *= ALPHA;
        }
    } else {
        for (int i = 0; i < N; ++i) {
            X[i * INCX] = 0;
        }
    }
}
 
void
blas_fill(const int N, const double ALPHA, double *X, const int INCX)
{
    for (int i = 0; i < N; ++i) {
        X[i * INCX] = ALPHA;
    }
}
 
double
blas_dot(const int N, const double *X, const int INCX, const double *Y,
         const int INCY)
{
    double dot = 0;
    for (int i = 0; i < N; ++i) {
        dot += X[i * INCX] * Y[i * INCY];
    }
    return dot;
}
 
void
blas_mul(const int N, const double *X, const int INCX, double *Y,
         const int INCY)
{
    for (int i = 0; i < N; ++i) {
        Y[i * INCY] *= X[i * INCX];
    }
}
 
double
blas_sum(const double *X, const int N)
{
    double sum = 0;
    for (int i = 0; i < N; ++i) {
        sum += X[i];
    }
    return sum;
}