Skip to content

Commit

Permalink
half-chebyshev transforms
Browse files Browse the repository at this point in the history
  • Loading branch information
MikaelSlevinsky committed Sep 12, 2024
1 parent abd33bc commit 4c9dc99
Show file tree
Hide file tree
Showing 5 changed files with 275 additions and 0 deletions.
172 changes: 172 additions & 0 deletions src/banded_source.c
Original file line number Diff line number Diff line change
Expand Up @@ -2297,6 +2297,21 @@ void X(create_associated_hermite_to_hermite_diagonal_connection_coefficient)(con
}
}

void X(create_half_chebyshev_to_chebyshev_diagonal_connection_coefficient)(const int n, FLT * D, const int INCD) {
FLT v = ONE(FLT);
for (int i = 0; i < n; i++) {
D[i*INCD] = v;
v /= TWO(FLT);
}
}

void X(create_chebyshev_to_half_chebyshev_diagonal_connection_coefficient)(const int n, FLT * D, const int INCD) {
FLT v = ONE(FLT);
for (int i = 0; i < n; i++) {
D[i*INCD] = v;
v *= TWO(FLT);
}
}

X(triangular_banded) * X(create_A_legendre_to_chebyshev)(const int norm, const int n) {
X(triangular_banded) * A = X(calloc_triangular_banded)(n, 2);
Expand Down Expand Up @@ -2831,6 +2846,163 @@ X(triangular_banded) * X(create_B_associated_hermite_to_hermite)(const int norm,

X(triangular_banded) * X(create_C_associated_hermite_to_hermite)(const int n) {return X(create_I_triangular_banded)(n, 0);}

// [(X2+3I)*(X2-I)*D12 + (X2+I)*R12]*DTU
X(triangular_banded) * X(create_A_half_chebyshev_to_chebyshev)(const int n) {
X(banded) * A = X(calloc_banded)(n, n, 0, 4);

X(banded) * DTU = X(malloc_banded)(n, n, -1, 1);
for (int i = 1; i < n; i++)
X(set_banded_index)(DTU, i, i-1, i);
X(banded) * D12 = X(malloc_banded)(n, n, -1, 1);
for (int i = 1; i < n; i++)
X(set_banded_index)(D12, 2, i-1, i);
X(banded) * R12 = X(calloc_banded)(n, n, 0, 2);
if (n > 0)
X(set_banded_index)(R12, 1, 0, 0);
if (n > 1)
X(set_banded_index)(R12, ONE(FLT)/2, 1, 1);
for (int i = 2; i < n; i++) {
X(set_banded_index)(R12, -ONE(FLT)/(i+1), i-2, i);
X(set_banded_index)(R12, ONE(FLT)/(i+1), i, i);
}
X(banded) * X2 = X(calloc_banded)(n, n, 1, 1);
for (int i = 0; i < n; i++) {
FLT v = (i+((FLT) 3))/(2*i+4);
X(set_banded_index)(X2, v, i-1, i);
v = (i+ONE(FLT))/(2*i+4);
X(set_banded_index)(X2, v, i+1, i);
}

// A2 = (X2+3I)*(X2-I)*D12*DTU
X(banded) * A2 = X(calloc_banded)(n, n, 0, 4);
X(banded) * A2a = X(calloc_banded)(n, n, 0, 2);
for (int i = 0; i < n; i++)
X(set_banded_index)(X2, -1, i, i);
X(gbmm)(1, X2, D12, 0, A2a);
X(banded) * A2b = X(calloc_banded)(n, n, 1, 3);
for (int i = 0; i < n; i++)
X(set_banded_index)(X2, 3, i, i);
X(gbmm)(1, X2, A2a, 0, A2b);
X(gbmm)(1, A2b, DTU, 0, A2);

// A1 = (X2+I)*R12*DTU
X(banded) * A1 = X(calloc_banded)(n, n, 0, 4);
X(banded) * A1a = X(calloc_banded)(n, n, 1, 3);
for (int i = 0; i < n; i++)
X(set_banded_index)(X2, 1, i, i);
X(gbmm)(1, X2, R12, 0, A1a);
X(gbmm)(1, A1a, DTU, 0, A1);

X(banded_add)(1, A1, 1, A2, A);

X(destroy_banded)(DTU);
X(destroy_banded)(D12);
X(destroy_banded)(R12);
X(destroy_banded)(X2);

X(destroy_banded)(A2a);
X(destroy_banded)(A2b);
X(destroy_banded)(A2);
X(destroy_banded)(A1a);
X(destroy_banded)(A1);

return X(convert_banded_to_triangular_banded)(A);
}

// R12*RTU
X(triangular_banded) * X(create_B_half_chebyshev_to_chebyshev)(const int n) {
X(banded) * B = X(calloc_banded)(n, n, 0, 4);

X(banded) * RTU = X(calloc_banded)(n, n, 0, 2);
if (n > 0)
X(set_banded_index)(RTU, 1, 0, 0);
if (n > 1)
X(set_banded_index)(RTU, ONE(FLT)/2, 1, 1);
for (int i = 2; i < n; i++) {
X(set_banded_index)(RTU, -ONE(FLT)/2, i-2, i);
X(set_banded_index)(RTU, ONE(FLT)/2, i, i);
}
X(banded) * R12 = X(calloc_banded)(n, n, 0, 2);
if (n > 0)
X(set_banded_index)(R12, 1, 0, 0);
if (n > 1)
X(set_banded_index)(R12, ONE(FLT)/2, 1, 1);
for (int i = 2; i < n; i++) {
X(set_banded_index)(R12, -ONE(FLT)/(i+1), i-2, i);
X(set_banded_index)(R12, ONE(FLT)/(i+1), i, i);
}

X(gbmm)(1, R12, RTU, 0, B);

X(destroy_banded)(RTU);
X(destroy_banded)(R12);

return X(convert_banded_to_triangular_banded)(B);
}

// [(X2-I)*X2*D12 + (X2-0.5I)*R12]*DTU
X(triangular_banded) * X(create_A_chebyshev_to_half_chebyshev)(const int n) {
X(banded) * A = X(calloc_banded)(n, n, 0, 4);

X(banded) * DTU = X(malloc_banded)(n, n, -1, 1);
for (int i = 1; i < n; i++)
X(set_banded_index)(DTU, i, i-1, i);
X(banded) * D12 = X(malloc_banded)(n, n, -1, 1);
for (int i = 1; i < n; i++)
X(set_banded_index)(D12, 2, i-1, i);
X(banded) * R12 = X(calloc_banded)(n, n, 0, 2);
if (n > 0)
X(set_banded_index)(R12, 1, 0, 0);
if (n > 1)
X(set_banded_index)(R12, ONE(FLT)/2, 1, 1);
for (int i = 2; i < n; i++) {
X(set_banded_index)(R12, -ONE(FLT)/(i+1), i-2, i);
X(set_banded_index)(R12, ONE(FLT)/(i+1), i, i);
}
X(banded) * X2 = X(calloc_banded)(n, n, 1, 1);
for (int i = 0; i < n; i++) {
FLT v = (i+((FLT) 3))/(2*i+4);
X(set_banded_index)(X2, v, i-1, i);
v = (i+ONE(FLT))/(2*i+4);
X(set_banded_index)(X2, v, i+1, i);
}

// A2 = (X2-I)*X2*D12*DTU
X(banded) * A2 = X(calloc_banded)(n, n, 0, 4);
X(banded) * A2a = X(calloc_banded)(n, n, 0, 2);
X(gbmm)(1, X2, D12, 0, A2a);
X(banded) * A2b = X(calloc_banded)(n, n, 1, 3);
for (int i = 0; i < n; i++)
X(set_banded_index)(X2, -1, i, i);
X(gbmm)(1, X2, A2a, 0, A2b);
X(gbmm)(1, A2b, DTU, 0, A2);

// A1 = (X2-0.5I)*R12*DTU
X(banded) * A1 = X(calloc_banded)(n, n, 0, 4);
X(banded) * A1a = X(calloc_banded)(n, n, 1, 3);
for (int i = 0; i < n; i++)
X(set_banded_index)(X2, -ONE(FLT)/2, i, i);
X(gbmm)(1, X2, R12, 0, A1a);
X(gbmm)(1, A1a, DTU, 0, A1);

X(banded_add)(1, A1, 1, A2, A);

X(destroy_banded)(DTU);
X(destroy_banded)(D12);
X(destroy_banded)(R12);
X(destroy_banded)(X2);

X(destroy_banded)(A2a);
X(destroy_banded)(A2b);
X(destroy_banded)(A2);
X(destroy_banded)(A1a);
X(destroy_banded)(A1);

return X(convert_banded_to_triangular_banded)(A);
}

X(triangular_banded) * X(create_B_chebyshev_to_half_chebyshev)(const int n) {return X(create_B_half_chebyshev_to_chebyshev)(n);}

X(banded) * X(operator_normalized_jacobi_clenshaw)(const int n, const int nc, const FLT * c, const int incc, const X(cop_params) params) {
FLT alpha = params.alpha;
FLT beta = params.beta;
Expand Down
8 changes: 8 additions & 0 deletions src/banded_source.h
Original file line number Diff line number Diff line change
Expand Up @@ -218,6 +218,8 @@ void X(create_laguerre_to_laguerre_diagonal_connection_coefficient)(const int no
void X(create_associated_jacobi_to_jacobi_diagonal_connection_coefficient)(const int norm1, const int norm2, const int n, const FLT c, const FLT alpha, const FLT beta, const FLT gamma, const FLT delta, FLT * D, const int INCD);
void X(create_associated_laguerre_to_laguerre_diagonal_connection_coefficient)(const int norm1, const int norm2, const int n, const FLT c, const FLT alpha, const FLT beta, FLT * D, const int INCD);
void X(create_associated_hermite_to_hermite_diagonal_connection_coefficient)(const int norm1, const int norm2, const int n, const FLT c, FLT * D, const int INCD);
void X(create_half_chebyshev_to_chebyshev_diagonal_connection_coefficient)(const int n, FLT * D, const int INCD);
void X(create_chebyshev_to_half_chebyshev_diagonal_connection_coefficient)(const int n, FLT * D, const int INCD);

X(triangular_banded) * X(create_A_legendre_to_chebyshev)(const int norm, const int n);
X(triangular_banded) * X(create_B_legendre_to_chebyshev)(const int norm, const int n);
Expand Down Expand Up @@ -246,6 +248,12 @@ X(triangular_banded) * X(create_A_associated_hermite_to_hermite)(const int norm,
X(triangular_banded) * X(create_B_associated_hermite_to_hermite)(const int norm, const int n);
X(triangular_banded) * X(create_C_associated_hermite_to_hermite)(const int n);

X(triangular_banded) * X(create_A_half_chebyshev_to_chebyshev)(const int n);
X(triangular_banded) * X(create_B_half_chebyshev_to_chebyshev)(const int n);

X(triangular_banded) * X(create_A_chebyshev_to_half_chebyshev)(const int n);
X(triangular_banded) * X(create_B_chebyshev_to_half_chebyshev)(const int n);

X(banded) * X(operator_normalized_jacobi_clenshaw)(const int n, const int nc, const FLT * c, const int incc, const X(cop_params) params);
X(banded) * X(operator_normalized_laguerre_clenshaw)(const int n, const int nc, const FLT * c, const int incc, const X(cop_params) params);
X(banded) * X(operator_normalized_hermite_clenshaw)(const int n, const int nc, const FLT * c, const int incc, const X(cop_params) params);
Expand Down
4 changes: 4 additions & 0 deletions src/ftinternal.h
Original file line number Diff line number Diff line change
Expand Up @@ -362,6 +362,10 @@ double * plan_chebyshev_to_ultraspherical(const int normcheb, const int normultr
double * plan_associated_jacobi_to_jacobi(const int norm1, const int norm2, const int n, const int c, const double alpha, const double beta, const double gamma, const double delta);
double * plan_associated_laguerre_to_laguerre(const int norm1, const int norm2, const int n, const int c, const double alpha, const double beta);
double * plan_associated_hermite_to_hermite(const int norm1, const int norm2, const int n, const int c);
double * plan_half_chebyshev_to_chebyshev(const int n);
float * plan_half_chebyshev_to_chebyshevf(const int n);
double * plan_chebyshev_to_half_chebyshev(const int n);
float * plan_chebyshev_to_half_chebyshevf(const int n);

typedef struct ft_rotation_plan_s {
double * s;
Expand Down
60 changes: 60 additions & 0 deletions src/transforms.c
Original file line number Diff line number Diff line change
Expand Up @@ -264,4 +264,64 @@ double * plan_associated_hermite_to_hermite(const int norm1, const int norm2, co
return V;
}

double * plan_half_chebyshev_to_chebyshev(const int n) {
ft_triangular_bandedl * A = ft_create_A_half_chebyshev_to_chebyshevl(n);
ft_triangular_bandedl * B = ft_create_B_half_chebyshev_to_chebyshevl(n);
long double * Vl = calloc(n*n, sizeof(long double));
ft_create_half_chebyshev_to_chebyshev_diagonal_connection_coefficientl(n, Vl, n+1);
ft_triangular_banded_eigenvectorsl(A, B, Vl);
double * V = malloc(n*n*sizeof(double));
for (int i = 0; i < n*n; i++)
V[i] = Vl[i];
ft_destroy_triangular_bandedl(A);
ft_destroy_triangular_bandedl(B);
free(Vl);
return V;
}

float * plan_half_chebyshev_to_chebyshevf(const int n) {
ft_triangular_banded * A = ft_create_A_half_chebyshev_to_chebyshev(n);
ft_triangular_banded * B = ft_create_B_half_chebyshev_to_chebyshev(n);
double * Vd = calloc(n*n, sizeof(double));
ft_create_half_chebyshev_to_chebyshev_diagonal_connection_coefficient(n, Vd, n+1);
ft_triangular_banded_eigenvectors(A, B, Vd);
float * V = malloc(n*n*sizeof(float));
for (int i = 0; i < n*n; i++)
V[i] = Vd[i];
ft_destroy_triangular_banded(A);
ft_destroy_triangular_banded(B);
free(Vd);
return V;
}

double * plan_chebyshev_to_half_chebyshev(const int n) {
ft_triangular_bandedl * A = ft_create_A_chebyshev_to_half_chebyshevl(n);
ft_triangular_bandedl * B = ft_create_B_chebyshev_to_half_chebyshevl(n);
long double * Vl = calloc(n*n, sizeof(long double));
ft_create_chebyshev_to_half_chebyshev_diagonal_connection_coefficientl(n, Vl, n+1);
ft_triangular_banded_eigenvectorsl(A, B, Vl);
double * V = malloc(n*n*sizeof(double));
for (int i = 0; i < n*n; i++)
V[i] = Vl[i];
ft_destroy_triangular_bandedl(A);
ft_destroy_triangular_bandedl(B);
free(Vl);
return V;
}

float * plan_chebyshev_to_half_chebyshevf(const int n) {
ft_triangular_banded * A = ft_create_A_chebyshev_to_half_chebyshev(n);
ft_triangular_banded * B = ft_create_B_chebyshev_to_half_chebyshev(n);
double * Vd = calloc(n*n, sizeof(double));
ft_create_chebyshev_to_half_chebyshev_diagonal_connection_coefficient(n, Vd, n+1);
ft_triangular_banded_eigenvectors(A, B, Vd);
float * V = malloc(n*n*sizeof(float));
for (int i = 0; i < n*n; i++)
V[i] = Vd[i];
ft_destroy_triangular_banded(A);
ft_destroy_triangular_banded(B);
free(Vd);
return V;
}

#include "transforms_mpfr.c"
31 changes: 31 additions & 0 deletions test/test_transforms.c
Original file line number Diff line number Diff line change
Expand Up @@ -130,5 +130,36 @@ int main(void) {
free(colsum);
}
printf("\n");
printf("\nTesting methods for half-to-full Chebyshev polynomial transforms.\n\n");
printf("\t\t\t Test \t\t\t\t | 2-norm Relative Error\n");
printf("---------------------------------------------------------|----------------------\n");
for (int n = 8; n < 24; n += 2) {
double err = 0;
double * V = plan_half_chebyshev_to_chebyshev(n);
double * iV = plan_chebyshev_to_half_chebyshev(n);
double * Id = calloc(n*n, sizeof(double));
for (int i = 0; i < n; i++)
Id[i+i*n] = 1.0;
ft_trmm('N', n, iV, n, V, n, n);
err = ft_norm_2arg(V, Id, n*n)/ft_norm_1arg(Id, n*n);
printf("V⁻¹V ≈ I \t\t\t\t n = %3i \t |%20.2e ", n, err);
ft_checktest(err, n*n, &checksum);
free(V);
V = plan_half_chebyshev_to_chebyshev(n);
ft_trmm('N', n, V, n, iV, n, n);
err = ft_norm_2arg(iV, Id, n*n)/ft_norm_1arg(Id, n*n);
printf("V V⁻¹ ≈ I \t\t\t\t n = %3i \t |%20.2e ", n, err);
ft_checktest(err, n*n*n*n, &checksum);
free(iV);
iV = plan_chebyshev_to_half_chebyshev(n);
ft_trsm('N', n, V, n, Id, n, n);
err = ft_norm_2arg(Id, iV, n*n)/ft_norm_1arg(iV, n*n);
printf("V\\I ≈ V⁻¹ \t\t\t\t n = %3i \t |%20.2e ", n, err);
ft_checktest(err, n*n, &checksum);
free(V);
free(iV);
free(Id);
}
printf("\n");
return checksum;
}

0 comments on commit 4c9dc99

Please sign in to comment.