c_vector_operations.cpp

Go to the documentation of this file.

/*
 *  SPDX-FileCopyrightText: Copyright 2021, Siavash Ameli <sameli@berkeley.edu>
 *  SPDX-License-Identifier: BSD-3-Clause
 *  SPDX-FileType: SOURCE
 *
 *  This program is free software: you can redistribute it and/or modify it
 *  under the terms of the license found in the LICENSE.txt file in the root
 *  directory of this source tree.
 */
 
 
// =======
// Headers
// =======
 
#include "./c_vector_operations.h"
#include <cmath>  // sqrt
#include "../_definitions/definitions.h"  // USE_CBLAS
 
#if (USE_CBLAS == 1)
    #include "./cblas_interface.h"
#endif
 
 
// ===========
// copy vector
// ===========
 
 
template <typename DataType>
void cVectorOperations<DataType>::copy_vector(
        const DataType* input_vector,
        const LongIndexType vector_size,
        DataType* output_vector)
{
    #if (USE_CBLAS == 1)
 
    // Using Openblas
    int incx = 1;
    int incy = 1;
 
    cblas_interface::xcopy(vector_size, input_vector, incx, output_vector,
                           incy);
 
    #else
 
    // Not using OpenBlas
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        output_vector[i] = input_vector[i];
    }
 
    #endif
}
 
// ==================
// copy scaled vector
// ==================
 
 
template <typename DataType>
void cVectorOperations<DataType>::copy_scaled_vector(
        const DataType* input_vector,
        const LongIndexType vector_size,
        const DataType scale,
        DataType* output_vector)
{
    #if (USE_CBLAS == 1)
 
    // Using OpenBlas
    int incx = 1;
    int incy = 1;
 
    cblas_interface::xcopy(vector_size, input_vector, incx, output_vector,
                           incy);
 
    cblas_interface::xscal(vector_size, scale, output_vector, incy); 
 
    #else
 
    // Not using OpenBlas
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        output_vector[i] = scale * input_vector[i];
    }
 
    #endif
}
 
 
// ======================
// subtract scaled vector
// ======================
 
 
template <typename DataType>
void cVectorOperations<DataType>::subtract_scaled_vector(
        const DataType* input_vector,
        const LongIndexType vector_size,
        const DataType scale,
        DataType* output_vector)
{
 
    #if (USE_CBLAS == 1)
 
    // Using OpenBlas
    int incx = 1;
    int incy = 1;
 
    DataType neg_scale = -scale;
    cblas_interface::xaxpy(vector_size, neg_scale, input_vector, incx,
                           output_vector, incy);
 
    #else
 
    // Not using OpenBlas
    if (scale == 0.0)
    {
        return;
    }
 
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        output_vector[i] -= scale * input_vector[i];
    }
 
    #endif
}
 
 
// =============
// inner product
// =============
 
 
template <typename DataType>
DataType cVectorOperations<DataType>::inner_product(
        const DataType* vector1,
        const DataType* vector2,
        const LongIndexType vector_size)
{
    #if (USE_CBLAS == 1)
 
    // Using OpenBlas
    int incx = 1;
    int incy = 1;
 
    DataType inner_prod = cblas_interface::xdot(vector_size, vector1, incx,
                                                vector2, incy);
 
    return inner_prod;
 
    #else
 
    // Not using OpenBlas
    long double inner_prod = 0.0;
    LongIndexType chunk = 5;
    LongIndexType vector_size_chunked = vector_size - (vector_size % chunk);
 
    for (LongIndexType i=0; i < vector_size_chunked; i += chunk)
    {
        inner_prod += vector1[i] * vector2[i] +
                      vector1[i+1] * vector2[i+1] +
                      vector1[i+2] * vector2[i+2] +
                      vector1[i+3] * vector2[i+3] +
                      vector1[i+4] * vector2[i+4];
    }
 
    for (LongIndexType i=vector_size_chunked; i < vector_size; ++i)
    {
        inner_prod += vector1[i] * vector2[i];
    }
 
    return static_cast<DataType>(inner_prod);
 
    #endif
}
 
 
// ==============
// euclidean norm
// ==============
 
 
template <typename DataType>
DataType cVectorOperations<DataType>::euclidean_norm(
        const DataType* vector,
        const LongIndexType vector_size)
{
    #if (USE_CBLAS == 1)
 
    // Using OpenBlas
    int incx = 1;
 
    DataType norm = cblas_interface::xnrm2(vector_size, vector, incx);
 
    return norm;
 
    #else
 
    // Compute norm squared
    long double norm2 = 0.0;
    LongIndexType chunk = 5;
    LongIndexType vector_size_chunked = vector_size - (vector_size % chunk);
 
    for (LongIndexType i=0; i < vector_size_chunked; i += chunk)
    {
        norm2 += vector[i] * vector[i] +
                 vector[i+1] * vector[i+1] +
                 vector[i+2] * vector[i+2] +
                 vector[i+3] * vector[i+3] +
                 vector[i+4] * vector[i+4];
    }
 
    for (LongIndexType i=vector_size_chunked; i < vector_size; ++i)
    {
        norm2 += vector[i] * vector[i];
    }
 
    // Norm
    DataType norm = sqrt(static_cast<DataType>(norm2));
 
    return norm;
 
    #endif
}
 
 
// =========================
// normalize vector in place
// =========================
 
 
template <typename DataType>
DataType cVectorOperations<DataType>::normalize_vector_in_place(
        DataType* vector,
        const LongIndexType vector_size)
{
    #if (USE_CBLAS == 1)
 
    // Norm of vector
    DataType norm = cVectorOperations<DataType>::euclidean_norm(
            vector, vector_size);
 
    // Normalize in place
    DataType scale = 1.0 / norm;
    int incx = 1;
    cblas_interface::xscal(vector_size, scale, vector, incx);
 
    return norm;
 
    #else
 
    // Norm of vector
    DataType norm = cVectorOperations<DataType>::euclidean_norm(vector,
                                                                vector_size);
 
    // Normalize in place
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        vector[i] /= norm;
    }
 
    return norm;
 
    #endif
}
 
 
// =========================
// normalize vector and copy
// =========================
 
 
template <typename DataType>
DataType cVectorOperations<DataType>::normalize_vector_and_copy(
        const DataType* vector,
        const LongIndexType vector_size,
        DataType* output_vector)
{
    #if (USE_CBLAS == 1)
 
    // Norm of vector
    DataType norm = cVectorOperations<DataType>::euclidean_norm(
            vector, vector_size);
 
    // Normalize to output
    DataType scale = 1.0 / norm;
    cVectorOperations<DataType>::copy_scaled_vector(vector, vector_size, scale,
                                                    output_vector);
 
    return norm;
 
    #else
 
    // Norm of vector
    DataType norm = cVectorOperations<DataType>::euclidean_norm(vector,
                                                                vector_size);
 
    // Normalize to output
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        output_vector[i] = vector[i] / norm;
    }
 
    return norm;
 
    #endif
}
 
 
// ===============================
// Explicit template instantiation
// ===============================
 
template class cVectorOperations<float>;
template class cVectorOperations<double>;
template class cVectorOperations<long double>;