diff --git a/src/banded_source.c b/src/banded_source.c index fdae32e..7b1809e 100644 --- a/src/banded_source.c +++ b/src/banded_source.c @@ -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); @@ -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; diff --git a/src/banded_source.h b/src/banded_source.h index 6eb5457..2f60ea0 100644 --- a/src/banded_source.h +++ b/src/banded_source.h @@ -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); @@ -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); diff --git a/src/ftinternal.h b/src/ftinternal.h index ff3ae6f..3b8c4f7 100644 --- a/src/ftinternal.h +++ b/src/ftinternal.h @@ -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; diff --git a/src/transforms.c b/src/transforms.c index a75fbb8..3e0836f 100644 --- a/src/transforms.c +++ b/src/transforms.c @@ -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" diff --git a/test/test_transforms.c b/test/test_transforms.c index b0336ef..18a47f9 100644 --- a/test/test_transforms.c +++ b/test/test_transforms.c @@ -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; }