syr2k.hh

// Copyright (c) 2017-2023, University of Tennessee. All rights reserved.
// SPDX-License-Identifier: BSD-3-Clause
// This program is free software: you can redistribute it and/or modify it under
// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
 
#ifndef BLAS_SYR2K_HH
#define BLAS_SYR2K_HH
 
#include "blas/util.hh"
 
#include <limits>
 
namespace blas {
 
// =============================================================================
 
template <typename TA, typename TB, typename TC>
void syr2k(
    blas::Layout layout,
    blas::Uplo uplo,
    blas::Op trans,
    int64_t n, int64_t k,
    scalar_type<TA, TB, TC> alpha,
    TA const *A, int64_t lda,
    TB const *B, int64_t ldb,
    scalar_type<TA, TB, TC> beta,
    TC       *C, int64_t ldc )
{
    typedef blas::scalar_type<TA, TB, TC> scalar_t;
 
    #define A(i_, j_) A[ (i_) + (j_)*lda ]
    #define B(i_, j_) B[ (i_) + (j_)*ldb ]
    #define C(i_, j_) C[ (i_) + (j_)*ldc ]
 
    // constants
    const scalar_t zero = 0;
    const scalar_t one  = 1;
 
    // check arguments
    blas_error_if( layout != Layout::ColMajor &&
                   layout != Layout::RowMajor );
    blas_error_if( uplo != Uplo::Lower &&
                   uplo != Uplo::Upper &&
                   uplo != Uplo::General );
    blas_error_if( n < 0 );
    blas_error_if( k < 0 );
 
    // check and interpret argument trans
    if (trans == Op::ConjTrans) {
        blas_error_if_msg(
                ( blas::is_complex<TA>::value ||
                  blas::is_complex<TB>::value ),
                "trans == Op::ConjTrans && "
                "( blas::is_complex<TA>::value ||"
                "  blas::is_complex<TB>::value )" );
        trans = Op::Trans;
    }
    else {
        blas_error_if( trans != Op::NoTrans &&
                       trans != Op::Trans );
    }
 
    // adapt if row major
    if (layout == Layout::RowMajor) {
        if (uplo == Uplo::Lower)
            uplo = Uplo::Upper;
        else if (uplo == Uplo::Upper)
            uplo = Uplo::Lower;
        trans = (trans == Op::NoTrans)
                ? Op::Trans
                : Op::NoTrans;
    }
 
    // check remaining arguments
    blas_error_if( lda < ((trans == Op::NoTrans) ? n : k) );
    blas_error_if( ldb < ((trans == Op::NoTrans) ? n : k) );
    blas_error_if( ldc < n );
 
    // quick return
    if (n == 0 || k == 0)
        return;
 
    // alpha == zero
    if (alpha == zero) {
        if (beta == zero) {
            if (uplo != Uplo::Upper) {
                for (int64_t j = 0; j < n; ++j) {
                    for (int64_t i = 0; i <= j; ++i)
                        C(i, j) = zero;
                }
            }
            else if (uplo != Uplo::Lower) {
                for (int64_t j = 0; j < n; ++j) {
                    for (int64_t i = j; i < n; ++i)
                        C(i, j) = zero;
                }
            }
            else {
                for (int64_t j = 0; j < n; ++j) {
                    for (int64_t i = 0; i < n; ++i)
                        C(i, j) = zero;
                }
            }
        }
        else if (beta != one) {
            if (uplo != Uplo::Upper) {
                for (int64_t j = 0; j < n; ++j) {
                    for (int64_t i = 0; i <= j; ++i)
                        C(i, j) *= beta;
                }
            }
            else if (uplo != Uplo::Lower) {
                for (int64_t j = 0; j < n; ++j) {
                    for (int64_t i = j; i < n; ++i)
                        C(i, j) *= beta;
                }
            }
            else {
                for (int64_t j = 0; j < n; ++j) {
                    for (int64_t i = 0; i < n; ++i)
                        C(i, j) *= beta;
                }
            }
        }
        return;
    }
 
    // alpha != zero
    if (trans == Op::NoTrans) {
        if (uplo != Uplo::Lower) {
            // uplo == Uplo::Upper or uplo == Uplo::General
            for (int64_t j = 0; j < n; ++j) {
 
                for (int64_t i = 0; i <= j; ++i)
                    C(i, j) *= beta;
 
                for (int64_t l = 0; l < k; ++l) {
                    scalar_t alpha_Bjl = alpha*B(j, l);
                    scalar_t alpha_Ajl = alpha*A(j, l);
                    for (int64_t i = 0; i <= j; ++i)
                        C(i, j) += A(i, l)*alpha_Bjl + B(i, l)*alpha_Ajl;
                }
            }
        }
        else { // uplo == Uplo::Lower
            for (int64_t j = 0; j < n; ++j) {
 
                for (int64_t i = j; i < n; ++i)
                    C(i, j) *= beta;
 
                for (int64_t l = 0; l < k; ++l) {
                    scalar_t alpha_Bjl = alpha*B(j, l);
                    scalar_t alpha_Ajl = alpha*A(j, l);
                    for (int64_t i = j; i < n; ++i)
                        C(i, j) += A(i, l)*alpha_Bjl + B(i, l)*alpha_Ajl;
                }
            }
        }
    }
    else { // trans == Op::Trans
        if (uplo != Uplo::Lower) {
            // uplo == Uplo::Upper or uplo == Uplo::General
            for (int64_t j = 0; j < n; ++j) {
                for (int64_t i = 0; i <= j; ++i) {
                    scalar_t sum1 = zero;
                    scalar_t sum2 = zero;
                    for (int64_t l = 0; l < k; ++l) {
                        sum1 += A(l, i) * B(l, j);
                        sum2 += B(l, i) * A(l, j);
                    }
                    C(i, j) = alpha*sum1 + alpha*sum2 + beta*C(i, j);
                }
            }
        }
        else { // uplo == Uplo::Lower
            for (int64_t j = 0; j < n; ++j) {
                for (int64_t i = j; i < n; ++i) {
                    scalar_t sum1 = zero;
                    scalar_t sum2 = zero;
                    for (int64_t l = 0; l < k; ++l) {
                        sum1 +=  A(l, i) * B(l, j);
                        sum2 +=  B(l, i) * A(l, j);
                    }
                    C(i, j) = alpha*sum1 + alpha*sum2 + beta*C(i, j);
                }
            }
        }
    }
 
    if (uplo == Uplo::General) {
        for (int64_t j = 0; j < n; ++j) {
            for (int64_t i = j+1; i < n; ++i)
                C(i, j) = C(j, i);
        }
    }
 
    #undef A
    #undef B
    #undef C
}
 
}  // namespace blas
 
#endif        //  #ifndef BLAS_SYMM_HH