diagonalization.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 "./diagonalization.h"
#include <cassert>  // assert
#include <cstdlib>  // NULL
#include "./lapack_api.h"  // lapack_xstev, lapack_xbdsdc,
 
 
// ================
// eigh tridiagonal
// ================
 
 
template <typename DataType>
int Diagonalization<DataType>::eigh_tridiagonal(
        DataType* diagonals,
        DataType* subdiagonals,
        DataType* eigenvectors,
        IndexType matrix_size)
{
    // 'V' computes both eigenvalues and eigenvectors
    char jobz = 'V';
 
    // Workspace array
    DataType* work = new DataType[2*matrix_size - 2];
 
    // Leading dimension of 2D eigenvectors array
    int ldz = matrix_size;
 
    // matrix size
    int n = static_cast<int>(matrix_size);
 
    // Error code output
    int info;
 
    // Calling Fortran subroutine
    lapack_xstev(&jobz, &n, diagonals, subdiagonals, eigenvectors, &ldz, work,
              &info);
 
    delete[] work;
    assert((info == 0, "?stev subroutine returned non-zero status."));
 
    return info;
}
 
 
// ==============
// svd bidiagonal
// ==============
 
 
template <typename DataType>
int Diagonalization<DataType>::svd_bidiagonal(
        DataType* diagonals,
        DataType* supdiagonals,
        DataType* U,
        DataType* Vt,
        IndexType matrix_size)
{
    // Code 'U' indicates the matrix is upper bi-diagonal
    char UPLO = 'U';
 
    // Code 'I' computes both singular values and singular vectors
    char COMPQ = 'I';
 
    // matrix size
    int n = static_cast<int>(matrix_size);
 
    // Leading dimensions of arrays U and Vt
    int LDU = matrix_size;
    int LDVT = matrix_size;
 
    // There arrays are not referenced when COMPQ is set to 'I'
    DataType* Q = NULL;
    int* IQ = NULL;
 
    // Work arrays
    DataType* work = new DataType[3*matrix_size*matrix_size + 4*matrix_size];
    int* iwork = new int[8 * matrix_size];
 
    // Error code output
    int info;
 
    // Calling Fortran subroutine
    lapack_xbdsdc(&UPLO, &COMPQ, &n, diagonals, supdiagonals, U, &LDU, Vt,
            &LDVT, Q, IQ, work, iwork, &info);
 
    delete[] work;
    delete[] iwork;
 
    assert((info == 0, "?stev subroutine returned non-zero status."));
 
    return info;
}
 
 
// ===============================
// Explicit template instantiation
// ===============================
 
template class Diagonalization<float>;
template class Diagonalization<double>;
template class Diagonalization<long double>;