c_orthogonalization.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.
 */
 
 
// =======
// Imports
// =======
 
#include "./c_orthogonalization.h"
#include <cstdlib>  // abort
#include <iostream>  // std::cerr, std::endl
#include <cmath>  // std::sqrt, std::fabs
#include <limits>  // std::numeric_limits
#include "../_c_basic_algebra/c_vector_operations.h"  // cVectorOperations
#include "../_random_generator/random_array_generator.h"  // RandomArrayGene...
#include "../_random_generator/random_number_generator.h"  // RandomNumberGe...
 
 
// ====================
// gram schmidt process
// ====================
 
 
template <typename DataType>
void cOrthogonalization<DataType>::gram_schmidt_process(
        const DataType* V,
        const LongIndexType vector_size,
        const IndexType num_vectors,
        const IndexType last_vector,
        const FlagType num_ortho,
        DataType* v)
{
    // 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 = std::numeric_limits<DataType>::epsilon();
    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 = cVectorOperations<DataType>::euclidean_norm(
                &V[vector_size*i], vector_size);
 
        // Check norm
        if (norm < epsilon * static_cast<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 = cVectorOperations<DataType>::inner_product(
                &V[vector_size*i], v, vector_size);
 
        // scale for subtraction
        DataType scale = inner_prod / (norm * norm);
 
        // If scale is 1, it is possible that vector v and j-th vector are
        // identical (or close).
        if (static_cast<DataType>(std::abs(
                static_cast<DataType>(std::abs(scale)) - \
                static_cast<DataType>(1.0))) <= \
                static_cast<DataType>(2.0) * epsilon)
        {
            // Norm of the vector v
            norm_v = cVectorOperations<DataType>::euclidean_norm(
                    v, vector_size);
 
            // Compute distance between the j-th vector and vector v
            distance2 = norm_v*norm_v - \
                    static_cast<DataType>(2.0)*inner_prod + norm*norm;
 
            // If distance is zero, do not reorthogonalize i-th against
            // the j-th vector.
            if (distance2 < \
                static_cast<DataType>(2.0) * epsilon * \
                static_cast<DataType>(vector_size))
            {
                continue;
            }
        }
 
        // Subtraction
        cVectorOperations<DataType>::subtract_scaled_vector(
                &V[vector_size*i], vector_size, scale, v);
    }
}
 
 
// =====================
// orthogonalize vectors
// =====================
 
 
template <typename DataType>
void cOrthogonalization<DataType>::orthogonalize_vectors(
        DataType* vectors,
        const LongIndexType vector_size,
        const IndexType num_vectors)
{
    // 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 = std::numeric_limits<DataType>::epsilon();
    IndexType success = 1;
    IndexType max_num_trials = 20;
    IndexType num_trials = 0;
    IndexType num_threads = 1;
    RandomNumberGenerator random_number_generator(num_threads);
 
    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;
        }
 
        // Reorthogonalize against previous vectors
        for (j=start_j; j < i; ++j)
        {
            // Norm of the j-th vector
            norm_j = cVectorOperations<DataType>::euclidean_norm(
                    &vectors[j*vector_size], vector_size);
 
            // Check norm
            if (norm_j < \
                    epsilon * static_cast<DataType>(std::sqrt(vector_size)))
            {
                std::cerr << "WARNING: norm of the given vector is too " \
                          << " small. Cannot reorthogonalize against zero" \
                          << "vector. Skipping."
                          << std::endl;
                continue;
            }
 
            // Projecting i-th vector to j-th vector
            inner_prod = cVectorOperations<DataType>::inner_product(
                    &vectors[i*vector_size], &vectors[j*vector_size],
                    vector_size);
 
            // Scale of subtraction
            DataType scale = inner_prod / (norm_j * norm_j);
 
            // Subtraction
            cVectorOperations<DataType>::subtract_scaled_vector(
                    &vectors[vector_size*j], vector_size, scale,
                    &vectors[vector_size*i]);
 
            // Norm of the i-th vector
            norm_i = cVectorOperations<DataType>::euclidean_norm(
                    &vectors[i*vector_size], vector_size);
 
            // If the norm is too small, regenerate the i-th vector randomly
            if (norm_i < \
                    epsilon * static_cast<DataType>(std::sqrt(vector_size)))
            {
                // Regenerate new random vector for i-th vector
                RandomArrayGenerator<DataType>::generate_random_array(
                        random_number_generator, &vectors[i*vector_size],
                        vector_size, num_threads);
 
                // Repeat the reorthogonalization 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;
        }
    }
}
 
 
// ===============================
// Explicit template instantiation
// ===============================
 
template class cOrthogonalization<float>;
template class cOrthogonalization<double>;
template class cOrthogonalization<long double>;