A static class for vector operations, similar to level-1 operations of the BLAS library. This class acts as a templated namespace, where all member methods are public and static. More...

#include <c_vector_operations.h>

Static Public Member Functions
static void	copy_vector (const DataType RESTRICT input_vector, const LongIndexType vector_size, DataType RESTRICT output_vector)
	Copies a vector to a new vector. Result is written in-place.

static void	copy_scaled_vector (const DataType RESTRICT input_vector, const LongIndexType vector_size, const DataType scale, DataType RESTRICT output_vector)
	Scales a vector and stores to a new vector.

static void	subtract_scaled_vector (const DataType RESTRICT input_vector, const LongIndexType vector_size, const DataType scale, DataType RESTRICT output_vector)
	Subtracts the scaled input vector from the output vector.

static DataType	inner_product (const DataType RESTRICT vector1, const DataType RESTRICT vector2, const LongIndexType vector_size)
	Computes Euclidean inner product of two vectors.

static DataType	euclidean_norm (const DataType *RESTRICT vector, const LongIndexType vector_size)
	Computes the Euclidean norm of a 1D array.

static DataType	normalize_vector_in_place (DataType *RESTRICT vector, const LongIndexType vector_size)
	Normalizes a vector based on Euclidean 2-norm. The result is written in-place.

static DataType	normalize_vector_and_copy (const DataType RESTRICT vector, const LongIndexType vector_size, DataType RESTRICT output_vector)
	Normalizes a vector based on Euclidean 2-norm. The result is written into another vector.

Detailed Description

template<typename DataType>
class cVectorOperations< DataType >

A static class for vector operations, similar to level-1 operations of the BLAS library. This class acts as a templated namespace, where all member methods are public and static.

See also: MatrixOperations

Definition at line 46 of file c_vector_operations.h.

Member Function Documentation

◆ copy_scaled_vector()

template<typename DataType >

void cVectorOperations< DataType >::copy_scaled_vector	(	const DataType *RESTRICT	input_vector,
		const LongIndexType	vector_size,
		const DataType	scale,
		DataType *RESTRICT	output_vector
	)

static

Scales a vector and stores to a new vector.

Parameters

[in]	input_vector	A 1D array
[in]	vector_size	Length of vector array
[in]	scale	Scale coefficient to the input vector. If this is equal to one, the function effectively becomes the same as copy_vector.
[out]	output_vector	Output vector (written in place).

Definition at line 93 of file c_vector_operations.cpp.

{
    #if defined(USE_ANY_CBLAS) && (USE_ANY_CBLAS == 1)
 
    // Using BLAS
    int incx = 1;
    int incy = 1;
 
    cblas_api::xcopy(vector_size, input_vector, incx, output_vector, incy);
 
    cblas_api::xscal(vector_size, scale, output_vector, incy); 
 
    #else
 
    // Not using BLAS
    #if defined(USE_OPENMP) && (USE_OPENMP == 1)
    #pragma omp parallel for \
        schedule(static) \
        if ((!omp_in_parallel()) && (vector_size >= LARGE_ARRAY_SIZE)) \
        default(none) \
        shared(input_vector, output_vector, vector_size, scale)
    #endif
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        output_vector[i] = scale * input_vector[i];
    }
 
    #endif
}

Referenced by c_lanczos_tridiagonalization(), and cVectorOperations< DataType >::normalize_vector_and_copy().

Here is the caller graph for this function:

◆ copy_vector()

template<typename DataType >

void cVectorOperations< DataType >::copy_vector	(	const DataType *RESTRICT	input_vector,
		const LongIndexType	vector_size,
		DataType *RESTRICT	output_vector
	)

static

Copies a vector to a new vector. Result is written in-place.

Parameters

[in]	input_vector	A 1D array
[in]	vector_size	Length of vector array
[out]	output_vector	Output vector (written in place).

Definition at line 45 of file c_vector_operations.cpp.

{
    #if defined(USE_ANY_CBLAS) && (USE_ANY_CBLAS == 1)
 
    // Using BLAS
    int incx = 1;
    int incy = 1;
 
    cblas_api::xcopy(vector_size, input_vector, incx, output_vector, incy);
 
    #else
 
    // Not using BLAS
    #if defined(USE_OPENMP) && (USE_OPENMP == 1)
    #pragma omp parallel for \
        schedule(static) \
        if ((!omp_in_parallel()) && (vector_size >= LARGE_ARRAY_SIZE)) \
        default(none) \
        shared(input_vector, output_vector, vector_size)
    #endif
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        output_vector[i] = input_vector[i];
    }
 
    #endif
}

Referenced by c_lanczos_tridiagonalization().

Here is the caller graph for this function:

◆ euclidean_norm()

template<typename DataType >

DataType cVectorOperations< DataType >::euclidean_norm	(	const DataType *RESTRICT	vector,
		const LongIndexType	vector_size
	)

static

Computes the Euclidean norm of a 1D array.

The reduction variable (here, inner_prod ) is of the type long double. This is becase when DataType is float, the summation loses the precision, especially when the vector size is large. It seems that using long double is slightly faster than using double. The advantage of using a type with larger bits for the reduction variable is only sensible if the compiler is optimized with -O2 or -O3 flags.

Using a larger bit type for the reduction variable is very important for this function. If DataType is float, without such consideration, the result of estimation of trace can be completely wrong, just becase of the wrong norm results. For large array sizes, even libraries such as openblas does not compute the dot product accurately.

The chunk computation of the dot product (as seen in the code with chunk=5) improves the preformance with gaining twice speedup. This result is not much dependet on chunk. For example, chunk=10 also yields a similar result.

Parameters

[in]	vector	A pointer to 1D array
[in]	vector_size	Length of the array

Returns: Euclidean norm

Definition at line 332 of file c_vector_operations.cpp.

{
    #if defined(USE_ANY_CBLAS) && (USE_ANY_CBLAS == 1)
 
    // Using BLAS
    int incx = 1;
 
    DataType norm = cblas_api::xnrm2(vector_size, vector, incx);
 
    return norm;
 
    #else
 
    // Compute norm squared
    long double norm2 = 0.0;
    #if defined(USE_LOOP_UNROLLING) && (USE_LOOP_UNROLLING == 1)
    LongIndexType chunk = 5;
    LongIndexType vector_size_chunked = vector_size - (vector_size % chunk);
    #endif
 
    #if defined(USE_LOOP_UNROLLING) && (USE_LOOP_UNROLLING == 1)
        #if defined(USE_OPENMP) && (USE_OPENMP == 1)
        #pragma omp parallel for \
            schedule(static) \
            if ((!omp_in_parallel()) && \
                (vector_size_chunked >= LARGE_ARRAY_SIZE)) \
            default(none) \
            shared(vector, vector_size, vector_size_chunked, chunk) \
            reduction(+: norm2)
        #endif
    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];
    }
    #endif
 
    #if defined(USE_LOOP_UNROLLING) && (USE_LOOP_UNROLLING == 1)
    for (LongIndexType i=vector_size_chunked; i < vector_size; ++i)
    #else
        #if defined(USE_OPENMP) && (USE_OPENMP == 1)
        #pragma omp parallel for \
            schedule(static) \
            if ((!omp_in_parallel()) && (vector_size >= LARGE_ARRAY_SIZE)) \
            default(none) \
            shared(vector, vector_size) \
            reduction(+: norm2)
        #endif
    for (LongIndexType i=0; i < vector_size; ++i)
    #endif
    {
        norm2 += vector[i] * vector[i];
    }
 
    // Norm
    DataType norm = std::sqrt(static_cast<DataType>(norm2));
 
    return norm;
 
    #endif
}

References LARGE_ARRAY_SIZE, and USE_OPENMP.

Referenced by c_lanczos_tridiagonalization(), cOrthogonalization< DataType >::gram_schmidt_process(), cVectorOperations< DataType >::normalize_vector_and_copy(), cVectorOperations< DataType >::normalize_vector_in_place(), and cOrthogonalization< DataType >::orthogonalize_vectors().

Here is the caller graph for this function:

◆ inner_product()

template<typename DataType >

DataType cVectorOperations< DataType >::inner_product	(	const DataType *RESTRICT	vector1,
		const DataType *RESTRICT	vector2,
		const LongIndexType	vector_size
	)

static

Computes Euclidean inner product of two vectors.

The reduction variable (here, inner_prod ) is of the type long double. This is becase when DataType is float, the summation loses the precision, especially when the vector size is large. It seems that using long double is slightly faster than using double. The advantage of using a type with larger bits for the reduction variable is only sensible if the compiler is optimized with -O2 or -O3 flags.

Using a larger bit type for the reduction variable is very important for this function. If DataType is float, without such consideration, the result of estimation of trace can be completely wrong, just becase of the wrong inner product results. For large array sizes, even libraries such as openblas does not compute the dot product accurately.

The chunk computation of the dot product (as seen in the code with chunk=5) improves the preformance with gaining twice speedup. This result is not much dependet on chunk. For example, chunk=10 also yields a similar result.

Parameters

[in]	vector1	1D array
[in]	vector2	1D array
[in]	vector_size	Length of array

Returns: Inner product of two vectors.

Definition at line 230 of file c_vector_operations.cpp.

{
    #if defined(USE_ANY_CBLAS) && (USE_ANY_CBLAS == 1)
 
    // Using BLAS
    int incx = 1;
    int incy = 1;
 
    DataType inner_prod = cblas_api::xdot(vector_size, vector1, incx, vector2,
                                          incy);
 
    return inner_prod;
 
    #else
 
    // Not using BLAS
    long double inner_prod = 0.0;
    #if defined(USE_LOOP_UNROLLING) && (USE_LOOP_UNROLLING == 1)
    LongIndexType chunk = 5;
    LongIndexType vector_size_chunked = vector_size - (vector_size % chunk);
    #endif
 
    #if defined(USE_LOOP_UNROLLING) && (USE_LOOP_UNROLLING == 1)
        #if defined(USE_OPENMP) && (USE_OPENMP == 1)
        #pragma omp parallel for \
            schedule(static) \
            if ((!omp_in_parallel()) && \
                (vector_size_chunked >= LARGE_ARRAY_SIZE)) \
            default(none) \
            shared(vector1, vector2, vector_size, vector_size_chunked, chunk) \
            reduction(+: inner_prod)
        #endif
    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];
    }
    #endif
 
    #if defined(USE_LOOP_UNROLLING) && (USE_LOOP_UNROLLING == 1)
    for (LongIndexType i=vector_size_chunked; i < vector_size; ++i)
    #else
        #if defined(USE_OPENMP) && (USE_OPENMP == 1)
        #pragma omp parallel for \
            schedule(static) \
            if ((!omp_in_parallel()) && (vector_size >= LARGE_ARRAY_SIZE)) \
            default(none) \
            shared(vector1, vector2, vector_size) \
            reduction(+: inner_prod)
        #endif
    for (LongIndexType i=0; i < vector_size; ++i)
    #endif
    {
        inner_prod += vector1[i] * vector2[i];
    }
 
    return static_cast<DataType>(inner_prod);
 
    #endif
}

References LARGE_ARRAY_SIZE, and USE_OPENMP.

Referenced by c_lanczos_tridiagonalization(), cOrthogonalization< DataType >::gram_schmidt_process(), and cOrthogonalization< DataType >::orthogonalize_vectors().

Here is the caller graph for this function:

◆ normalize_vector_and_copy()

template<typename DataType >

DataType cVectorOperations< DataType >::normalize_vector_and_copy	(	const DataType *RESTRICT	vector,
		const LongIndexType	vector_size,
		DataType *RESTRICT	output_vector
	)

static

Normalizes a vector based on Euclidean 2-norm. The result is written into another vector.

Parameters

[in]	vector	Input vector.
[in]	vector_size	Length of the input vector
[out]	output_vector	Output vector, which is the normalization of the input vector.

Returns: 2-norm of the input vector

Definition at line 472 of file c_vector_operations.cpp.

{
    #if defined(USE_ANY_CBLAS) && (USE_ANY_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
    #if defined(USE_OPENMP) && (USE_OPENMP == 1)
    #pragma omp parallel for \
        schedule(static) \
        if ((!omp_in_parallel()) && (vector_size >= LARGE_ARRAY_SIZE)) \
        default(none) \
        shared(vector, output_vector, vector_size, norm)
    #endif
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        output_vector[i] = vector[i] / norm;
    }
 
    return norm;
 
    #endif
}

References cVectorOperations< DataType >::copy_scaled_vector(), and cVectorOperations< DataType >::euclidean_norm().

Referenced by c_golub_kahn_bidiagonalization().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ normalize_vector_in_place()

template<typename DataType >

DataType cVectorOperations< DataType >::normalize_vector_in_place	(	DataType *RESTRICT	vector,
		const LongIndexType	vector_size
	)

static

Normalizes a vector based on Euclidean 2-norm. The result is written in-place.

Parameters

[in,out]	vector	Input vector to be normalized in-place.
[in]	vector_size	Length of the input vector

Returns: 2-Norm of the input vector (before normalization)

Definition at line 414 of file c_vector_operations.cpp.

{
    #if defined(USE_ANY_CBLAS) && (USE_ANY_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_api::xscal(vector_size, scale, vector, incx);
 
    return norm;
 
    #else
 
    // Norm of vector
    DataType norm = cVectorOperations<DataType>::euclidean_norm(vector,
                                                                vector_size);
 
    // Normalize in place
    #if defined(USE_OPENMP) && (USE_OPENMP == 1)
    #pragma omp parallel for \
        schedule(static) \
        if ((!omp_in_parallel()) && (vector_size >= LARGE_ARRAY_SIZE)) \
        default(none) \
        shared(vector, vector_size, norm)
    #endif
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        vector[i] /= norm;
    }
 
    return norm;
 
    #endif
}

References cVectorOperations< DataType >::euclidean_norm().

Referenced by c_golub_kahn_bidiagonalization().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ subtract_scaled_vector()

template<typename DataType >

void cVectorOperations< DataType >::subtract_scaled_vector	(	const DataType *RESTRICT	input_vector,
		const LongIndexType	vector_size,
		const DataType	scale,
		DataType *RESTRICT	output_vector
	)

static

Subtracts the scaled input vector from the output vector.

Performs the following operation:

\[ \boldsymbol{b} = \boldsymbol{b} - c \boldsymbol{a}, \]

where

\( \boldsymbol{a} \) is the input vector,
\( c \) is a scalar scale to the input vector, and
\( \boldsymbol{b} \) is the output vector that is written in-place.

Parameters

[in]	input_vector	A 1D array
[in]	vector_size	Length of vector array
[in]	scale	Scale coefficient to the input vector.
[in,out]	output_vector	Output vector (written in place).

Definition at line 153 of file c_vector_operations.cpp.

{
 
    #if defined(USE_ANY_CBLAS) && (USE_ANY_CBLAS == 1)
 
    // Using BLAS
    int incx = 1;
    int incy = 1;
 
    DataType neg_scale = -scale;
    cblas_api::xaxpy(vector_size, neg_scale, input_vector, incx, output_vector,
                     incy);
 
    #else
 
    // Not using BLAS
    DataType zero = 0.0;
    if (c_arithmetics::is_equal(scale, zero))
    {
        return;
    }
 
    #if defined(USE_OPENMP) && (USE_OPENMP == 1)
    #pragma omp parallel for \
        schedule(static) \
        if ((!omp_in_parallel()) && (vector_size >= LARGE_ARRAY_SIZE)) \
        default(none) \
        shared(input_vector, output_vector, vector_size, scale)
    #endif
    for (LongIndexType i=0; i < vector_size; ++i)
    {
        output_vector[i] -= scale * input_vector[i];
    }
 
    #endif
}

References c_arithmetics::is_equal().

Referenced by cAffineMatrixFunction< DataType >::_add_scaled_vector(), c_golub_kahn_bidiagonalization(), c_lanczos_tridiagonalization(), cOrthogonalization< DataType >::gram_schmidt_process(), and cOrthogonalization< DataType >::orthogonalize_vectors().

Here is the call graph for this function:

Here is the caller graph for this function:

The documentation for this class was generated from the following files:

/home/runner/work/imate/imate/imate/_c_basic_algebra/c_vector_operations.h
/home/runner/work/imate/imate/imate/_c_basic_algebra/c_vector_operations.cpp

Static Public Member Functions

Detailed Description

Member Function Documentation

◆ copy_scaled_vector()

◆ copy_vector()

◆ euclidean_norm()

◆ inner_product()

◆ normalize_vector_and_copy()

◆ normalize_vector_in_place()

◆ subtract_scaled_vector()