gauss_matrix.hpp

Go to the documentation of this file.

/* 
 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
 * Copyright (c) 2022 Lodz University of Technology
 * 
 * This program 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, version 3.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 */
#ifndef PLASK__SOLVER__OPTICAL_EFFECTIVE_GAUSS_MATRIX_H
#define PLASK__SOLVER__OPTICAL_EFFECTIVE_GAUSS_MATRIX_H
 
#include <cstddef>
#include <plask/plask.hpp>
using plask::dcomplex;
 
// BLAS routine to multiply matrix by vector
#define zgbmv F77_GLOBAL(zgbmv,ZGBMV)
F77SUB zgbmv(const char& trans, const int& m, const int& n, const int& kl, const int& ku, const dcomplex& alpha, dcomplex* a, const int& lda,
             const dcomplex* x, int incx, const dcomplex& beta, dcomplex* y, int incy);
 
 
// LAPACK routines to solve set of linear equations
#define zgbtrf F77_GLOBAL(zgbtrf,ZGBTRF)
F77SUB zgbtrf(const int& m, const int& n, const int& kl, const int& ku, dcomplex* ab, const int& ldab, int* ipiv, int& info);
 
#define zgbtrs F77_GLOBAL(zgbtrs,ZGBTRS)
F77SUB zgbtrs(const char& trans, const int& n, const int& kl, const int& ku, const int& nrhs, dcomplex* ab, const int& ldab, int* ipiv, dcomplex* b, const int& ldb, int& info);
 
namespace plask { namespace optical { namespace effective {
 
constexpr const int LD = 7;
 
struct ZgbMatrix {
 
    const size_t size;              
    dcomplex* data;                 
 
    ZgbMatrix(size_t rank): size(rank), data(aligned_malloc<dcomplex>(LD*rank)) {}
 
    ZgbMatrix(const ZgbMatrix&) = delete; // this object is non-copyable
 
    ~ZgbMatrix() { aligned_free(data); }
 
    size_t index(size_t r, size_t c) {
        assert(r < size && c < size);
        assert(abs(int(c)-int(r)) < 3);
        // AB(kl+ku+1+i-j,j) = A(i,j)
        return (LD-1)*c + r + 4;
    }
 
    dcomplex& operator()(size_t r, size_t c) {
        return data[index(r,c)];
    }
 
    void clear() {
        std::fill_n(data, LD*size, 0.);
    }
 
    void mult(const DataVector<const dcomplex>& vector, DataVector<dcomplex>& result) {
        zgbmv('N', int(size), int(size), 2, 2, 1., data, LD, vector.data(), 1, 0., result.data(), 1);
    }
 
    void addmult(const DataVector<const dcomplex>& vector, DataVector<dcomplex>& result) {
        zgbmv('N', int(size), int(size), 2, 2, 1., data, LD, vector.data(), 1, 1., result.data(), 1);
    }
 
    dcomplex determinant() {
        int info = 0;
        aligned_unique_ptr<int> upiv(aligned_malloc<int>(size));
        int* ipiv = upiv.get();
        zgbtrf(int(size), int(size), 2, 2, data, LD, ipiv, info);
        assert(info >= 0);
 
        dcomplex det = 1.;
        for (std::size_t i = 0; i < size; ++i) {
            det *= data[LD*i + 4];
            if (ipiv[i] != int(i+1)) det = -det;
        }
        return det;
    }
};
 
}}} // # namespace plask::gain::freecarrier
 
#endif // PLASK__SOLVER__OPTICAL_EFFECTIVE_GAUSS_MATRIX_H