c_dense_affine_matrix_function.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_dense_affine_matrix_function.h"
#include <cassert>  // assert
#include <cstddef>  // NULL
 
 
// =============
// constructor 1
// =============
 
 
template <typename DataType>
cDenseAffineMatrixFunction<DataType>::cDenseAffineMatrixFunction(
        const DataType* A_,
        const LongIndexType num_rows_,
        const LongIndexType num_columns_,
        const FlagType A_is_row_major_,
        const FlagType A_is_symmetric_):
 
    // Base class constructor
    cLinearOperatorBase(num_rows_, num_columns_),
 
    // Initializer list
    A(A_, num_rows_, num_columns_, A_is_row_major_, A_is_symmetric_)
{
    // This constructor is called assuming B is identity
    this->B_is_identity = true;
 
    // When B is identity, the eigenvalues of A+tB are known for any t
    this->eigenvalue_relation_known = 1;
}
 
 
// =============
// constructor 2
// =============
 
 
template <typename DataType>
cDenseAffineMatrixFunction<DataType>::cDenseAffineMatrixFunction(
        const DataType* A_,
        const LongIndexType num_rows_,
        const LongIndexType num_columns_,
        const FlagType A_is_row_major_,
        const FlagType A_is_symmetric_,
        const DataType* B_,
        const FlagType B_is_row_major_,
        const FlagType B_is_symmetric_):
 
    // Base class constructor
    cLinearOperatorBase(num_rows_, num_columns_),
 
    // Initializer list
    A(A_, num_rows_, num_columns_, A_is_row_major_, A_is_symmetric_),
    B(B_, num_rows_, num_columns_, B_is_row_major_, B_is_symmetric_)
{
    // Matrix B is assumed to be non-zero. Check if it is identity or generic
    if (this->B.is_identity_matrix())
    {
        this->B_is_identity = true;
        this->eigenvalue_relation_known = 1;
    }
}
 
 
// ==========
// destructor
// ==========
 
 
template <typename DataType>
cDenseAffineMatrixFunction<DataType>::~cDenseAffineMatrixFunction()
{
}
 
 
// ============
// set symmetry
// ============
 
 
template <typename DataType>
void cDenseAffineMatrixFunction<DataType>::set_symmetry(
        const FlagType symmetric)
{
    if (symmetric == 1)
    {
        this->A.set_symmetry(1);
        this->B.set_symmetry(1);
    }
    else
    {
        this->A.set_symmetry(0);
        this->B.set_symmetry(0);
    }
}
 
 
// ===
// dot
// ===
 
 
template <typename DataType>
void cDenseAffineMatrixFunction<DataType>::dot(
        const DataType* vector,
        DataType* product)
{
    // Matrix A times vector
    this->A.dot(vector, product);
    LongIndexType min_vector_size;
 
    // Matrix B times vector to be added to the product
    if (this->B_is_identity)
    {
        // Check parameter is set
        assert((this->parameters != NULL) && "Parameter is not set.");
 
        // Find minimum of the number of rows and columns
        min_vector_size = \
            (this->num_rows < this->num_columns) ? \
            this->num_rows : this->num_columns;
 
        // Adding input vector to product
        this->_add_scaled_vector(vector, min_vector_size,
                                 this->parameters[0], product);
    }
    else
    {
        // Check parameter is set
        assert((this->parameters != NULL) && "Parameter is not set.");
 
        // Adding parameter times B times input vector to the product
        this->B.dot_plus(vector, this->parameters[0], product);
    }
}
 
 
// =============
// transpose dot
// =============
 
 
template <typename DataType>
void cDenseAffineMatrixFunction<DataType>::transpose_dot(
        const DataType* vector,
        DataType* product)
{
    // Matrix A times vector
    this->A.transpose_dot(vector, product);
    LongIndexType min_vector_size;
 
    // Matrix B times vector to be added to the product
    if (this->B_is_identity)
    {
        // Check parameter is set
        assert((this->parameters != NULL) && "Parameter is not set.");
 
        // Find minimum of the number of rows and columns
        min_vector_size = \
            (this->num_rows < this->num_columns) ? \
            this->num_rows : this->num_columns;
 
        // Adding input vector to product
        this->_add_scaled_vector(vector, min_vector_size,
                                 this->parameters[0], product);
    }
    else
    {
        // Check parameter is set
        assert((this->parameters != NULL) && "Parameter is not set.");
 
        // Adding "parameter * B * input vector" to the product
        this->B.transpose_dot_plus(vector, this->parameters[0], product);
    }
}
 
 
// ===============================
// Explicit template instantiation
// ===============================
 
template class cDenseAffineMatrixFunction<float>;
template class cDenseAffineMatrixFunction<double>;
template class cDenseAffineMatrixFunction<long double>;