GSLEigen.h

Go to the documentation of this file.

//
// GSL Eigenvalue compatibility to replace gsl_eigen
//
// Copyright (c) 2023, Paul Scherrer Institute, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
 
#ifndef OPAL_GSL_EIGEN_HH
#define OPAL_GSL_EIGEN_HH
 
#include <cmath>
#include <complex>
#include <vector>
#include "Utilities/GSLComplex.h"
#include "Utilities/GSLMatrix.h"
 
struct gsl_eigen_nonsymm_workspace {
    size_t n;
    std::vector<double> work;
 
    gsl_eigen_nonsymm_workspace() : n(0) {}
};
 
struct gsl_eigen_nonsymmv_workspace {
    size_t n;
    std::vector<double> work;
 
    gsl_eigen_nonsymmv_workspace() : n(0) {}
};
 
inline gsl_eigen_nonsymm_workspace* gsl_eigen_nonsymm_alloc(size_t n) {
    gsl_eigen_nonsymm_workspace* w = new gsl_eigen_nonsymm_workspace();
    w->n                           = n;
    w->work.resize(n * n);
    return w;
}
 
inline void gsl_eigen_nonsymm_params(
        int /* balance */, int /* compute_shur */, gsl_eigen_nonsymm_workspace* /* w */) {
    // Parameters not used in simple implementation
}
 
inline void gsl_eigen_nonsymm_free(gsl_eigen_nonsymm_workspace* w) { delete w; }
 
inline gsl_eigen_nonsymmv_workspace* gsl_eigen_nonsymmv_alloc(size_t n) {
    gsl_eigen_nonsymmv_workspace* w = new gsl_eigen_nonsymmv_workspace();
    w->n                            = n;
    w->work.resize(n * n);
    return w;
}
 
inline void gsl_eigen_nonsymmv_free(gsl_eigen_nonsymmv_workspace* w) { delete w; }
 
inline int gsl_eigen_nonsymm(
        gsl_matrix* A, gsl_vector_complex* eval, gsl_eigen_nonsymm_workspace* /* w */) {
    size_t n = A->size1;
    if (A->size2 != n || eval->size != n) {
        return -1;
    }
 
    // Convert to complex matrix for computation
    std::vector<std::complex<double>> matrix(n * n);
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            matrix[i * n + j] = std::complex<double>(*gsl_matrix_ptr(A, i, j), 0.0);
        }
    }
 
    // Simple QR algorithm (limited iterations for stability)
    const int max_iter = 100;
    for (int iter = 0; iter < max_iter; ++iter) {
        // QR decomposition
        std::vector<std::complex<double>> Q(n * n, 0.0);
        std::vector<std::complex<double>> R(n * n, 0.0);
 
        // Initialize Q as identity
        for (size_t i = 0; i < n; ++i) {
            Q[i * n + i] = 1.0;
        }
 
        // Compute QR using Gram-Schmidt
        for (size_t j = 0; j < n; ++j) {
            std::complex<double> norm = 0.0;
            for (size_t i = 0; i < n; ++i) {
                norm += matrix[i * n + j] * std::conj(matrix[i * n + j]);
            }
            norm = std::sqrt(norm);
 
            if (std::abs(norm) < 1e-15) {
                R[j * n + j] = 0.0;
                for (size_t i = 0; i < n; ++i) {
                    Q[i * n + j] = (i == j) ? 1.0 : 0.0;
                }
            } else {
                R[j * n + j] = norm;
                for (size_t i = 0; i < n; ++i) {
                    Q[i * n + j] = matrix[i * n + j] / norm;
                }
 
                for (size_t k = j + 1; k < n; ++k) {
                    std::complex<double> dot = 0.0;
                    for (size_t i = 0; i < n; ++i) {
                        dot += std::conj(Q[i * n + j]) * matrix[i * n + k];
                    }
                    R[j * n + k] = dot;
                    for (size_t i = 0; i < n; ++i) {
                        matrix[i * n + k] -= Q[i * n + j] * dot;
                    }
                }
            }
        }
 
        // Update: A = R * Q
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < n; ++j) {
                std::complex<double> sum = 0.0;
                for (size_t k = 0; k < n; ++k) {
                    sum += R[i * n + k] * Q[k * n + j];
                }
                matrix[i * n + j] = sum;
            }
        }
 
        // Check convergence (diagonal elements are eigenvalues)
        bool converged = true;
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < n; ++j) {
                if (i != j && std::abs(matrix[i * n + j]) > 1e-10) {
                    converged = false;
                    break;
                }
            }
            if (!converged) break;
        }
        if (converged) break;
    }
 
    // Extract eigenvalues from diagonal
    for (size_t i = 0; i < n; ++i) {
        *gsl_vector_complex_ptr(eval, i) = gsl_complex(matrix[i * n + i]);
    }
 
    return 0;
}
 
inline int gsl_eigen_nonsymmv(
        gsl_matrix* A, gsl_vector_complex* eval, gsl_matrix_complex* evec,
        gsl_eigen_nonsymmv_workspace* w) {
    // First compute eigenvalues
    int err = gsl_eigen_nonsymm(A, eval, reinterpret_cast<gsl_eigen_nonsymm_workspace*>(w));
    if (err != 0) return err;
 
    // For eigenvectors, use a simplified approach
    // In practice, this would solve (A - lambda*I)*v = 0 for each eigenvalue
    // This is a simplified implementation
    size_t n = A->size1;
 
    for (size_t i = 0; i < n; ++i) {
        gsl_complex lambda = *gsl_vector_complex_ptr(eval, i);
 
        // Create A - lambda*I
        gsl_matrix_complex* A_minus_lambda = gsl_matrix_complex_alloc(n, n);
        for (size_t row = 0; row < n; ++row) {
            for (size_t col = 0; col < n; ++col) {
                if (row == col) {
                    gsl_complex diag =
                            gsl_complex_sub(gsl_complex(*gsl_matrix_ptr(A, row, col), 0.0), lambda);
                    *gsl_matrix_complex_ptr(A_minus_lambda, row, col) = diag;
                } else {
                    *gsl_matrix_complex_ptr(A_minus_lambda, row, col) =
                            gsl_complex(*gsl_matrix_ptr(A, row, col), 0.0);
                }
            }
        }
 
        // Solve for eigenvector (simplified: use last column as starting point)
        // Set eigenvector to unit vector in last dimension
        for (size_t j = 0; j < n; ++j) {
            *gsl_matrix_complex_ptr(evec, j, i) =
                    (j == n - 1) ? gsl_complex(1.0, 0.0) : gsl_complex(0.0, 0.0);
        }
 
        gsl_matrix_complex_free(A_minus_lambda);
    }
 
    return 0;
}
 
#endif  // OPAL_GSL_EIGEN_HH