A static class for orthogonalization of vector bases. This class acts as a templated namespace, where all member methods are public and static. More...

#include <cu_orthogonalization.h>

Static Public Member Functions
static void	gram_schmidt_process (cublasHandle_t cublas_handle, const DataType V, const LongIndexType vector_size, const IndexType num_vectors, const IndexType last_vector, const FlagType num_ortho, DataType r)
	Modified Gram-Schmidt orthogonalization process to orthogonalize the vector `v` against a subset of the column vectors in the array `V`.

static void	orthogonalize_vectors (cublasHandle_t cublas_handle, DataType *vectors, const LongIndexType vector_size, const IndexType num_vectors)
	Orthogonalizes set of vectors mutually using modified Gram-Schmidt process.

Detailed Description

template<typename DataType>
class cuOrthogonalization< DataType >

A static class for orthogonalization of vector bases. This class acts as a templated namespace, where all member methods are public and static.

See also: RandomVectors

Definition at line 53 of file cu_orthogonalization.h.

Member Function Documentation

◆ gram_schmidt_process()

template<typename DataType >

void cuOrthogonalization< DataType >::gram_schmidt_process	(	cublasHandle_t	cublas_handle,
		const DataType *	V,
		const LongIndexType	vector_size,
		const IndexType	num_vectors,
		const IndexType	last_vector,
		const FlagType	num_ortho,
		DataType *	v
	)

static

Modified Gram-Schmidt orthogonalization process to orthogonalize the vector v against a subset of the column vectors in the array V.

V is 1D array of the length vector_size*num_vectors to represent a 2D array of a set of num_vectors column vectors, each of the length vector_size. The length of v is also vector_size.

v is orthogonalized against the last num_ortho columns of V starting from the column vector of the index last_vector. If the backward indexing from last_vector becomes a negative index, the index wraps around from the last column vector index, i.e., num_vectors-1 .

If num_ortho is zero, or if num_vectors is zero, no orthogonalization is performed.
If num_ortho is negative (usually set to -1), then v is orthogonalized against all column vectors of V.
If num_ortho is larger than num_vectors, then v is orthogonalized against all column vectors of V.
If num_ortho is smaller than num_vectors, then v is orthogonalized against the last num_ortho column vectors of V, starting from the column vector with the index last_vector toward its previous vectors. If the iteration runs into negativen column indices, the column indexing wraps around from the end of the columns from num_vectors-1.

The result of the newer v is written in-place in v.

If vector v is identical to one of the vectors in V, the orthogonalization against the identical vector is skipped.

If one of the column vectors of V is zero (have zero norm), that vector is ignored.

Note: It is assumed that the caller function fills the column vectors of V periodically in a wrapped around order from column index 0,1,... to num_vectors-1, and newer vectors are replaced on the wrapped index starting from index 0,1,... again. Thus, V only stores the last num_vectors column vectors. The index of the last filled vector is indicated by last_vector.

Warning: The vector v can be indeed one of the columns of V itself. However, in this case, vector v must NOT be orthogonalized against itself, rather, it should only be orthogonalized to the other vectors in V. For instance, if num_vectors=10, and v is the 3rd vector of V, and if num_ortho is 6, then we may set last_vector=2. Then v is orthogonalized againts the six columns 2,1,0,9,8,7, where the last three of them are wrapped around the end of the columns.

See also: cu_golub_kahn_bidiagonalizaton, cu_lanczos_bidiagonalization

Parameters

[in]	cublas_handle	The cuBLAS object handle.
[in]	V	1D coalesced array of vectors representing a 2D array. The length of this 1D array is `vector_size*num_vectors`, which indicates a 2D array with the shape (vector_size,num_vectors).
[in]	vector_size	The length of each vector. If we assume `V` indicates a 2D vector, this is the number of rows of `V`.
[in]	num_vectors	The number of column vectors. If we assume `V` indicates a 2D vector, this the number of columns of `V`.
[in]	last_vector	The column vectors of the array `V` are rewritten by the caller function in wrapped-around order. That is, once all the columns (from the zeroth to the `num_vector-1` vector) are filled, the next vector is rewritten in the place of the zeroth vector, and the indices of newer vectors wrap around the columns of `V`. Thus, `V` only retains the last `num_vectors` vectors. The column index of the last written vector is given by `last_vector`. This index is a number between `0` and `num_vectors-1`. The index of the last i-th vector is winding back from the last vector by `last_vector-i+1 mod num_vectors`.
[in]	num_ortho	The number of vectors to be orthogonalized starting from the last vector. `0` indicates no orthogonalization will be performed and the function just returns. A negative value means all vectors will be orthogonalized. A poisitive value will orthogonalize the given number of vectors. This value cannot be larger than the number of vectors.
[in,out]	v	The vector that will be orthogonalized against the columns of `V`. The length of `v` is `vector_size`. This vector is modified in-place.

Definition at line 131 of file cu_orthogonalization.cu.

{
    // Determine how many previous vectors to orthogonalize against
    IndexType num_steps;
    if ((num_ortho == 0) || (num_vectors < 2))
    {
        // No orthogonalization is performed
        return;
    }
    else if ((num_ortho < 0) ||
             (num_ortho > static_cast<FlagType>(num_vectors)))
    {
        // Orthogonalize against all vectors
        num_steps = num_vectors;
    }
    else
    {
        // Orthogonalize against only the last num_ortho vectors
        num_steps = num_ortho;
    }
 
    // Vectors can be orthogonalized at most to the full basis of the vector
    // space. Thus, num_steps cannot be larger than the dimension of vector
    // space, which is vector_size.
    if (num_steps > static_cast<IndexType>(vector_size))
    {
        num_steps = vector_size;
    }
 
    IndexType i;
    DataType inner_prod;
    DataType norm;
    DataType norm_v;
    DataType epsilon = cu_arithmetics::epsilon<DataType>();
    DataType distance2;
 
    // Iterate over vectors
    for (IndexType step=0; step < num_steps; ++step)
    {
        // i is the index of a column vector in V to orthogonalize v against it
        if ((last_vector % num_vectors) >= step)
        {
            i = (last_vector % num_vectors) - step;
        }
        else
        {
            // Wrap around negative indices from the end of column index
            i = (last_vector % num_vectors) - step + num_vectors;
        }
 
        // Norm of j-th vector
        norm = cuVectorOperations<DataType>::euclidean_norm(
                cublas_handle, &V[vector_size*i], vector_size);
 
        // Check norm
        if (norm < \
                cu_arithmetics::mul(
                    epsilon,
                    cu_arithmetics::cast<double, DataType>(
                        std::sqrt(vector_size))
                )
           )
        {
            std::cerr << "WARNING: norm of the given vector is too small. " \
                      << "Cannot orthogonalize against zero vector. " \
                      << "Skipping." << std::endl;
            continue;
        }
 
        // Projection
        inner_prod = cuVectorOperations<DataType>::inner_product(
                cublas_handle, &V[vector_size*i], v, vector_size);
 
        // scale for subtraction
        DataType scale = cu_arithmetics::div(
                inner_prod,
                cu_arithmetics::mul(norm, norm));
        DataType scale_abs = cu_arithmetics::abs(scale);
        DataType one = cu_arithmetics::cast<double, DataType>(1.0);
        DataType two = cu_arithmetics::cast<double, DataType>(2.0);
 
        // If scale is 1, it is possible that vector v and j-th vector are
        // identical (or close).
        if (cu_arithmetics::abs(cu_arithmetics::sub(scale_abs, one)) <= \
            cu_arithmetics::mul(two, epsilon))
        {
            // Norm of the vector v
            norm_v = cuVectorOperations<DataType>::euclidean_norm(
                    cublas_handle, v, vector_size);
 
            // Compute distance between the j-th vector and vector v
            distance2 = \
                cu_arithmetics::add(
                    cu_arithmetics::mul(norm_v, norm_v),
                    cu_arithmetics::mul(-two, inner_prod),
                    cu_arithmetics::mul(norm, norm)
                );
 
            // If distance is zero, do not re-orthogonalize i-th against
            // the j-th vector.
            if (distance2 < \
                    cu_arithmetics::mul(
                        two,
                        epsilon,
                        cu_arithmetics::cast<unsigned long long int, DataType>(
                            static_cast<unsigned long long int>(vector_size))
                    )
               )
            {
                continue;
            }
        }
 
        // Subtraction
        cuVectorOperations<DataType>::subtract_scaled_vector(
                cublas_handle, &V[vector_size*i], vector_size, scale, v);
    }
}

References cu_arithmetics::abs(), cu_arithmetics::div(), cuVectorOperations< DataType >::euclidean_norm(), cuVectorOperations< DataType >::inner_product(), cu_arithmetics::mul(), cu_arithmetics::sub(), and cuVectorOperations< DataType >::subtract_scaled_vector().

Referenced by cu_golub_kahn_bidiagonalization(), and cu_lanczos_tridiagonalization().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ orthogonalize_vectors()

template<typename DataType >

void cuOrthogonalization< DataType >::orthogonalize_vectors	(	cublasHandle_t	cublas_handle,
		DataType *	vectors,
		const LongIndexType	vector_size,
		const IndexType	num_vectors
	)

static

Orthogonalizes set of vectors mutually using modified Gram-Schmidt process.

Note: Let m be the number of vectors (num_vectors), and let n be the size of each vector (vector_size). In general, n is much larger (large matrix size), and m is small, in order of a couple of hundred. But for small matrices (where n could be smaller then m), then each vector can be orthogonalized at most to n other vectors. This is because the dimension of the vector space is n. Thus, if there are extra vectors, each vector is orthogonalized to window of the previous n vector.

If one of the column vectors is identical to one of other column vectors in V, one of the vectors is regenerated by random array and the orthogonalization is repeated.

Note: If two vectors are identical (or the norm of their difference is very small), they cannot be orthogonalized against each other. In this case, one of the vectors is re-generated by new random numbers.

Warning: if num_vectors is larger than vector_size, the orthogonalization fails since not all vectors are independent, and at least one vector becomes zero.

Parameters

[in]	cublas_handle	The cuBLAS object handle.
[in,out]	vectors	2D array of size `vector_sizenum_vectors`. This array will be modified in-place and will be output of this function. Note that this is Fortran ordering, meaning that the first index is contiguous. Hence, to call the j-th element of the i-th vector, use `&vectors`[ivector_size + j].
[in]	num_vectors	Number of columns of vectors array.
[in]	vector_size	Number of rows of vectors array.

Definition at line 302 of file cu_orthogonalization.cu.

{
    // Do nothing if there is only one vector
    if (num_vectors < 2)
    {
        return;
    }
 
    IndexType i = 0;
    IndexType j;
    IndexType start_j;
    DataType inner_prod;
    DataType norm_j;
    DataType norm_i;
    DataType epsilon = cu_arithmetics::epsilon<DataType>();
    IndexType success = 1;
    IndexType max_num_trials = 20;
    IndexType num_trials = 0;
    IndexType num_threads = 1;
    RandomNumberGenerator random_number_generator(num_threads);
    DataType* buffer = NULL;
 
    while (i < num_vectors)
    {
        if ((success == 0) && (num_trials >= max_num_trials))
        {
            std::cerr << "ERROR: Cannot orthogonalize vectors after " \
                      << num_trials << " trials. Aborting." \
                      << std::endl;
            abort();
        }
 
        // Reset on new trial (if it was set to 0 before to start a new trial)
        success = 1;
 
        // j iterates on previous vectors in a window of at most vector_size
        if (static_cast<LongIndexType>(i) > vector_size)
        {
            // When vector_size is smaller than i, it is fine to cast to signed
            start_j = i - static_cast<IndexType>(vector_size);
        }
        else
        {
            start_j = 0;
        }
 
        // Re-orthogonalize against previous vectors
        for (j=start_j; j < i; ++j)
        {
            // Norm of the j-th vector
            norm_j = cuVectorOperations<DataType>::euclidean_norm(
                    cublas_handle, &vectors[j*vector_size], vector_size);
 
            // Check norm
            if (norm_j < \
                    cu_arithmetics::mul(
                        epsilon,
                        cu_arithmetics::cast<double, DataType>(
                            std::sqrt(vector_size))
                    )
               )
            {
                std::cerr << "WARNING: norm of the given vector is too " \
                          << " small. Cannot re-orthogonalize against zero" \
                          << "vector. Skipping."
                          << std::endl;
                continue;
            }
 
            // Projecting i-th vector to j-th vector
            inner_prod = cuVectorOperations<DataType>::inner_product(
                    cublas_handle, &vectors[i*vector_size],
                    &vectors[j*vector_size], vector_size);
 
            // Scale of subtraction
            DataType scale = cu_arithmetics::div(
                    inner_prod,
                    cu_arithmetics::mul(norm_j, norm_j));
 
            // Subtraction
            cuVectorOperations<DataType>::subtract_scaled_vector(
                    cublas_handle, &vectors[vector_size*j], vector_size, scale,
                    &vectors[vector_size*i]);
 
            // Norm of the i-th vector
            norm_i = cuVectorOperations<DataType>::euclidean_norm(
                    cublas_handle, &vectors[i*vector_size], vector_size);
 
            // If the norm is too small, regenerate the i-th vector randomly
            if (norm_i < \
                    cu_arithmetics::mul(
                        epsilon,
                        cu_arithmetics::cast<double, DataType>(
                            std::sqrt(vector_size))
                    )
                )
            {
                // Allocate buffer
                if (buffer == NULL)
                {
                    buffer = new DataType[vector_size];
                }
 
                // Regenerate new random vector on buffer
                RandomArrayGenerator<DataType>::generate_random_array(
                        random_number_generator, buffer,
                        vector_size, num_threads);
 
                // Copy buffer to the i-th vector on device
                CudaAPI<DataType>::copy_to_device(
                        buffer, vector_size, &vectors[i*vector_size]);
 
                // Repeat the re-orthogonalization for i-th vector against
                // all previous vectors again.
                success = 0;
                ++num_trials;
                break;
            }
        }
 
        if (success == 1)
        {
            ++i;
 
            // Reset if num_trials was incremented before.
            num_trials = 0;
        }
    }
 
    // Deallocate buffer
    if (buffer != NULL)
    {
        delete[] buffer;
        buffer = NULL;
    }
}

References cu_arithmetics::abs(), CudaAPI< ArrayType >::copy_to_device(), cu_arithmetics::div(), cuVectorOperations< DataType >::euclidean_norm(), RandomArrayGenerator< DataType >::generate_random_array(), cuVectorOperations< DataType >::inner_product(), cu_arithmetics::mul(), and cuVectorOperations< DataType >::subtract_scaled_vector().

Here is the call graph for this function:

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

/home/runner/work/imate/imate/imate/_cu_trace_estimator/cu_orthogonalization.h
/home/runner/work/imate/imate/imate/_cu_trace_estimator/cu_orthogonalization.cu

Static Public Member Functions

Detailed Description

Member Function Documentation

◆ gram_schmidt_process()

◆ orthogonalize_vectors()