221 lines
6.5 KiB
C
221 lines
6.5 KiB
C
#include <stdlib.h>
|
|
#include <cblas.h>
|
|
#include "scatsystem.h"
|
|
#include "indexing.h"
|
|
#include <gsl/gsl_spline.h>
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
|
|
qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t));
|
|
if (!t) abort();
|
|
else {
|
|
t->spec = bspec;
|
|
size_t n = bspec->n;
|
|
t->m = calloc(n*n, sizeof(complex double));
|
|
if (!t->m) abort();
|
|
t->owns_m = true;
|
|
}
|
|
return t;
|
|
}
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_copy(const qpms_tmatrix_t *T) {
|
|
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
|
size_t n = T->spec->n;
|
|
for(size_t i = 0; i < n*n; ++i)
|
|
t->m = T->m;
|
|
return t;
|
|
}
|
|
|
|
void qpms_tmatrix_free(qpms_tmatrix_t *t){
|
|
if(t && t->owns_m) free(t->m);
|
|
free(t);
|
|
}
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_apply_symop(
|
|
const qpms_tmatrix_t *T,
|
|
const complex double *M
|
|
)
|
|
{
|
|
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
|
const size_t n = T->spec->n;
|
|
complex double tmp[n][n];
|
|
// tmp = M T
|
|
const complex double one = 1, zero = 0;
|
|
cblas_zgemm(CblasRowMajor,
|
|
CblasNoTrans,
|
|
CblasNoTrans,
|
|
n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
|
|
// t->m = tmp M* = M T M*
|
|
cblas_zgemm(CblasRowMajor,
|
|
CblasNoTrans,
|
|
CblasConjTrans,
|
|
n, n, n, &one, tmp, n, M, n, &zero, t->m, n);
|
|
return t;
|
|
}
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution_inplace(
|
|
qpms_tmatrix_t *T,
|
|
const complex double *M
|
|
)
|
|
{
|
|
qpms_tmatrix_t *t = qpms_tmatrix_apply_symop(T, M);
|
|
const size_t n = T->spec->n;
|
|
for(size_t i = 0; i < n*n; ++i)
|
|
T->m[i] = 0.5 * (t->m[i] + T->m[i]);
|
|
qpms_tmatrix_free(t);
|
|
return T;
|
|
}
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_symmetrise_involution(
|
|
const qpms_tmatrix_t *T,
|
|
const complex double *M
|
|
)
|
|
{
|
|
qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
|
|
const size_t n = T->spec->n;
|
|
complex double tmp[n][n];
|
|
// tmp = M T
|
|
const complex double one = 1, zero = 0;
|
|
cblas_zgemm(CblasRowMajor,
|
|
CblasNoTrans,
|
|
CblasNoTrans,
|
|
n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
|
|
// t->m = tmp M* = M T M*
|
|
cblas_zgemm(CblasRowMajor,
|
|
CblasNoTrans,
|
|
CblasConjTrans,
|
|
n, n, n, &one, tmp, n, M, n, &zero, t->m, n);
|
|
for(size_t i = 0; i < n*n; ++i)
|
|
t->m[i] = 0.5 * (t->m[i] + T->m[i]);
|
|
return t;
|
|
}
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf(const qpms_tmatrix_t *T) {
|
|
qpms_tmatrix_t *t = qpms_tmatrix_copy(T);
|
|
return qpms_tmatrix_symmetrise_C_inf_inplace(t);
|
|
}
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_inf_inplace(qpms_tmatrix_t *T) {
|
|
const size_t n = T->spec->n;
|
|
for (size_t row = 0; row < n; row++) {
|
|
qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
|
|
for (size_t col = 0; col < n; col++) {
|
|
qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
|
|
if (rm == cm)
|
|
;// No-op // t->m[n*row + col] = T->m[n*row + col];
|
|
else
|
|
T->m[n*row + col] = 0;
|
|
}
|
|
}
|
|
return T;
|
|
}
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N(const qpms_tmatrix_t *T, int N) {
|
|
qpms_tmatrix_t *t = qpms_tmatrix_copy(T);
|
|
return qpms_tmatrix_symmetrise_C_N_inplace(t, N);
|
|
}
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) {
|
|
const size_t n = T->spec->n;
|
|
for (size_t row = 0; row < n; row++) {
|
|
qpms_m_t rm = qpms_uvswfi2m(T->spec->ilist[row]);
|
|
for (size_t col = 0; col < n; col++) {
|
|
qpms_m_t cm = qpms_uvswfi2m(T->spec->ilist[col]);
|
|
if (((rm - cm) % N) == 0)
|
|
; // T->m[n*row + col] = T->m[n*row + col];
|
|
else
|
|
T->m[n*row + col] = 0;
|
|
}
|
|
}
|
|
return T;
|
|
}
|
|
|
|
bool qpms_tmatrix_isnear(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B,
|
|
const double rtol, const double atol)
|
|
{
|
|
if (!qpms_vswf_set_spec_isidentical(A->spec, B->spec))
|
|
return false;
|
|
if (A->m == B->m)
|
|
return true;
|
|
const size_t n = A->spec->n;
|
|
for (size_t i = 0; i < n*n; ++i) {
|
|
const double tol = atol + rtol * (cabs(B->m[i]));
|
|
if ( cabs(B->m[i] - A->m[i]) > tol )
|
|
return false;
|
|
}
|
|
return true;
|
|
}
|
|
|
|
|
|
qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(const size_t incount,
|
|
const double *freqs, const qpms_tmatrix_t *ta, const gsl_interp_type *iptype//, const bool copy_bspec
|
|
) {
|
|
if (incount <= 0) return NULL;
|
|
qpms_tmatrix_interpolator_t *ip = malloc(sizeof(qpms_tmatrix_interpolator_t));
|
|
/*
|
|
if (copy_bspec) {
|
|
ip->bspec = qpms_vswf_set_spec_copy(ta[0].spec);
|
|
ip->owns_bspec = true;
|
|
}
|
|
else {
|
|
*/
|
|
ip->bspec = ta[0].spec;
|
|
// ip->owns_bspec = false;
|
|
//}
|
|
const size_t n = ip->bspec->n;
|
|
|
|
// check if all matrices have the same bspec
|
|
for (size_t i = 0; i < incount; ++i)
|
|
if (!qpms_vswf_set_spec_isidentical(ip->bspec, ta[i].spec))
|
|
abort();
|
|
|
|
if (!(ip->splines_real = calloc(n*n,sizeof(gsl_spline *)))) abort();
|
|
if (!(ip->splines_imag = calloc(n*n,sizeof(gsl_spline *)))) abort();
|
|
for (size_t row = 0; row < n; ++row)
|
|
for (size_t col = 0; col < n; ++col) {
|
|
double y_real[incount], y_imag[incount];
|
|
bool n0_real = false, n0_imag = false;
|
|
for (size_t i = 0; i < incount; ++i) {
|
|
complex double telem = ta[i].m[n * row + col];
|
|
if (y_real[i] = creal(telem)) n0_real = true;
|
|
if (y_imag[i] = cimag(telem)) n0_imag = true;
|
|
}
|
|
if (n0_real) {
|
|
gsl_spline *s =
|
|
ip->splines_real[n * row + col] = gsl_spline_alloc(iptype, incount);
|
|
if (gsl_spline_init(s, freqs, y_real, incount) != 0 /*GSL_SUCCESS*/) abort();
|
|
}
|
|
else ip->splines_real[n * row + col] = NULL;
|
|
if (n0_imag) {
|
|
gsl_spline *s =
|
|
ip->splines_imag[n * row + col] = gsl_spline_alloc(iptype, incount);
|
|
if (gsl_spline_init(s, freqs, y_imag, incount) != 0 /*GSL_SUCCESS*/) abort();
|
|
}
|
|
else ip->splines_imag[n * row + col] = NULL;
|
|
}
|
|
return ip;
|
|
}
|
|
|
|
void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *ip) {
|
|
if (ip) {
|
|
const size_t n = ip->bspec->n;
|
|
for (size_t i = 0; i < n*n; ++i) {
|
|
if (ip->splines_real[i]) gsl_spline_free(ip->splines_real[i]);
|
|
if (ip->splines_imag[i]) gsl_spline_free(ip->splines_imag[i]);
|
|
}
|
|
//if (ip->owns_bspec)
|
|
// qpms_vswf_set_spec_free(ip->bspec);
|
|
free(ip);
|
|
}
|
|
}
|
|
|
|
qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *ip, double freq) {
|
|
qpms_tmatrix_t *t = qpms_tmatrix_init(ip->bspec);
|
|
const size_t n = ip->bspec->n;
|
|
for (size_t i = 0; i < n*n; ++i){
|
|
if (ip->splines_real[i]) t->m[i] = gsl_spline_eval(ip->splines_real[i], freq, NULL /*does this work?*/);
|
|
if (ip->splines_imag[i]) t->m[i] += I* gsl_spline_eval(ip->splines_imag[i], freq, NULL /*does this work?*/);
|
|
}
|
|
return t;
|
|
}
|
|
|