make mpfr dependency optional

to remove it, just comment out the definition of FT_USE_MPFR=1 in Make.inc

- add ft_mpfr_trsv and ft_mpfr_trsm

- simplify tests

Make.inc

@@ -49,11 +49,22 @@ ifeq ($(UNAME), Windows)
     LDFLAGS += -L.
 endif
 
-LDLIBS += -lm -lquadmath -lmpfr -fopenmp -lblas -lfftw3
+LDLIBS += -lm -lquadmath -fopenmp -lblas -lfftw3
 ifeq ($(UNAME), Linux)
-    LDLIBS += -lfftw3_threads -lgmp
+    LDLIBS += -lfftw3_threads
 else ifeq ($(UNAME), Darwin)
-    LDLIBS += -lfftw3_threads -lgmp
-else ifeq ($(OS), Windows_NT)
-    LDLIBS += -lmpir
+    LDLIBS += -lfftw3_threads
+endif
+
+FT_USE_MPFR=1
+
+ifdef FT_USE_MPFR
+    LDLIBS += -DFT_USE_MPFR -lmpfr
+    ifeq ($(UNAME), Linux)
+        LDLIBS += -lgmp
+    else ifeq ($(UNAME), Darwin)
+        LDLIBS += -lgmp
+    else ifeq ($(OS), Windows_NT)
+        LDLIBS += -lmpir
+    endif
 endif

Makefile

@@ -52,6 +52,7 @@ clean:
 	rm -f additiontheorem
 	rm -f calculus
 	rm -f holomorphic
+	rm -f nonlocaldiffusion
 	rm -f test_*
 
 .PHONY: all lib assembly examples tests clean

src/fasttransforms.h

@@ -171,39 +171,41 @@ ft_tb_eigen_FMMl * ft_plan_ultraspherical_to_chebyshevl(const int normultra, con
 /// A long double precision version of \ref ft_plan_chebyshev_to_ultraspherical.
 ft_tb_eigen_FMMl * ft_plan_chebyshev_to_ultrasphericall(const int normcheb, const int normultra, const int n, const long double lambda);
 
-#include <mpfr.h>
-
-typedef struct {
-    mpfr_t * data;
-    int n;
-    int b;
-} ft_mpfr_triangular_banded;
-
-void ft_mpfr_destroy_plan(mpfr_t * A, int n);
-void ft_mpfr_trmv(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * x, mpfr_rnd_t rnd);
-void ft_mpfr_trmm(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * B, int LDB, int N, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_legendre_to_chebyshev that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_legendre_to_chebyshev(const int normleg, const int normcheb, const int n, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_chebyshev_to_legendre that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_chebyshev_to_legendre(const int normcheb, const int normleg, const int n, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_ultraspherical_to_ultraspherical that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_ultraspherical_to_ultraspherical(const int norm1, const int norm2, const int n, mpfr_srcptr lambda, mpfr_srcptr mu, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_jacobi_to_jacobi that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_jacobi_to_jacobi(const int norm1, const int norm2, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_srcptr gamma, mpfr_srcptr delta, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_laguerre_to_laguerre that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_laguerre_to_laguerre(const int norm1, const int norm2, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_jacobi_to_ultraspherical that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_jacobi_to_ultraspherical(const int normjac, const int normultra, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_srcptr lambda, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_ultraspherical_to_jacobi that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_ultraspherical_to_jacobi(const int normultra, const int normjac, const int n, mpfr_srcptr lambda, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_jacobi_to_chebyshev that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_jacobi_to_chebyshev(const int normjac, const int normcheb, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_chebyshev_to_jacobi that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_chebyshev_to_jacobi(const int normcheb, const int normjac, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_ultraspherical_to_chebyshev that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_ultraspherical_to_chebyshev(const int normultra, const int normcheb, const int n, mpfr_srcptr lambda, mpfr_prec_t prec, mpfr_rnd_t rnd);
-/// A multi-precision version of \ref ft_plan_chebyshev_to_ultraspherical that returns a dense array of connection coefficients.
-mpfr_t * ft_mpfr_plan_chebyshev_to_ultraspherical(const int normcheb, const int normultra, const int n, mpfr_srcptr lambda, mpfr_prec_t prec, mpfr_rnd_t rnd);
+#ifdef FT_USE_MPFR
+    #include <mpfr.h>
+    typedef struct {
+        mpfr_t * data;
+        int n;
+        int b;
+    } ft_mpfr_triangular_banded;
+    void ft_mpfr_destroy_plan(mpfr_t * A, int n);
+    void ft_mpfr_trmv(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * x, mpfr_rnd_t rnd);
+    void ft_mpfr_trsv(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * x, mpfr_rnd_t rnd);
+    void ft_mpfr_trmm(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * B, int LDB, int N, mpfr_rnd_t rnd);
+    void ft_mpfr_trsm(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * B, int LDB, int N, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_legendre_to_chebyshev that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_legendre_to_chebyshev(const int normleg, const int normcheb, const int n, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_chebyshev_to_legendre that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_chebyshev_to_legendre(const int normcheb, const int normleg, const int n, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_ultraspherical_to_ultraspherical that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_ultraspherical_to_ultraspherical(const int norm1, const int norm2, const int n, mpfr_srcptr lambda, mpfr_srcptr mu, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_jacobi_to_jacobi that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_jacobi_to_jacobi(const int norm1, const int norm2, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_srcptr gamma, mpfr_srcptr delta, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_laguerre_to_laguerre that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_laguerre_to_laguerre(const int norm1, const int norm2, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_jacobi_to_ultraspherical that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_jacobi_to_ultraspherical(const int normjac, const int normultra, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_srcptr lambda, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_ultraspherical_to_jacobi that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_ultraspherical_to_jacobi(const int normultra, const int normjac, const int n, mpfr_srcptr lambda, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_jacobi_to_chebyshev that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_jacobi_to_chebyshev(const int normjac, const int normcheb, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_chebyshev_to_jacobi that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_chebyshev_to_jacobi(const int normcheb, const int normjac, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_ultraspherical_to_chebyshev that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_ultraspherical_to_chebyshev(const int normultra, const int normcheb, const int n, mpfr_srcptr lambda, mpfr_prec_t prec, mpfr_rnd_t rnd);
+    /// A multi-precision version of \ref ft_plan_chebyshev_to_ultraspherical that returns a dense array of connection coefficients.
+    mpfr_t * ft_mpfr_plan_chebyshev_to_ultraspherical(const int normcheb, const int normultra, const int n, mpfr_srcptr lambda, mpfr_prec_t prec, mpfr_rnd_t rnd);
+#endif
 
 /// Set the number of OpenMP threads.
 void ft_set_num_threads(const int n);

src/transforms.c

@@ -285,4 +285,6 @@ double * plan_chebyshev_to_ultraspherical(const int normcheb, const int normultr
     return V;
 }
 
-#include "transforms_mpfr.c"
+#ifdef FT_USE_MPFR
+    #include "transforms_mpfr.c"
+#endif

src/transforms_mpfr.c

@@ -31,13 +31,42 @@ void ft_mpfr_trmv(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * x, mpfr_rnd_t
     }
 }
 
+// x ← A⁻¹*x, x ← A⁻ᵀ*x
+void ft_mpfr_trsv(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * x, mpfr_rnd_t rnd) {
+    if (TRANS == 'N') {
+        for (int j = n-1; j >= 0; j--) {
+            mpfr_div(x[j], x[j], A[j+j*LDA], rnd);
+            for (int i = 0; i < j; i++) {
+                mpfr_fms(x[i], A[i+j*LDA], x[j], x[i], rnd);
+                mpfr_neg(x[i], x[i], rnd);
+            }
+        }
+    }
+    else if (TRANS == 'T') {
+        for (int i = 0; i < n; i++) {
+            for (int j = 0; j < i; j++) {
+                mpfr_fms(x[i], A[j+i*LDA], x[j], x[i], rnd);
+                mpfr_neg(x[i], x[i], rnd);
+            }
+            mpfr_div(x[i], x[i], A[i+i*LDA], rnd);
+        }
+    }
+}
+
 // B ← A*B, B ← Aᵀ*B
 void ft_mpfr_trmm(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * B, int LDB, int N, mpfr_rnd_t rnd) {
     #pragma omp parallel for
     for (int j = 0; j < N; j++)
         ft_mpfr_trmv(TRANS, n, A, LDA, B+j*LDB, rnd);
 }
 
+// B ← A⁻¹*B, B ← A⁻ᵀ*B
+void ft_mpfr_trsm(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * B, int LDB, int N, mpfr_rnd_t rnd) {
+    #pragma omp parallel for
+    for (int j = 0; j < N; j++)
+        ft_mpfr_trsv(TRANS, n, A, LDA, B+j*LDB, rnd);
+}
+
 ft_mpfr_triangular_banded * ft_mpfr_calloc_triangular_banded(const int n, const int b, mpfr_prec_t prec) {
     mpfr_t * data = malloc(n*(b+1)*sizeof(mpfr_t));
     for (int j = 0; j < n; j++)

test/test_banded.c

@@ -1,26 +1,6 @@
 #include "fasttransforms.h"
 #include "ftutilities.h"
 
-void test_bandedf(int * checksum);
-void test_banded (int * checksum);
-void test_bandedl(int * checksum);
-void test_bandedq(int * checksum);
-
-int main(void) {
-    int checksum = 0;
-    printf("\nTesting methods for banded matrices.\n");
-    printf("\n\tSingle precision.\n\n");
-    test_bandedf(&checksum);
-    printf("\n\tDouble precision.\n\n");
-    test_banded(&checksum);
-    printf("\n\tLong double precision.\n\n");
-    test_bandedl(&checksum);
-    printf("\n\tQuadruple precision.\n\n");
-    test_bandedq(&checksum);
-    printf("\n");
-    return checksum;
-}
-
 #define FLT float
 #define X(name) FT_CONCAT(ft_, name, f)
 #define Y(name) FT_CONCAT(, name, f)
@@ -52,3 +32,18 @@ int main(void) {
 #undef FLT
 #undef X
 #undef Y
+
+int main(void) {
+    int checksum = 0;
+    printf("\nTesting methods for banded matrices.\n");
+    printf("\n\tSingle precision.\n\n");
+    test_bandedf(&checksum);
+    printf("\n\tDouble precision.\n\n");
+    test_banded(&checksum);
+    printf("\n\tLong double precision.\n\n");
+    test_bandedl(&checksum);
+    printf("\n\tQuadruple precision.\n\n");
+    test_bandedq(&checksum);
+    printf("\n");
+    return checksum;
+}

test/test_dprk.c

@@ -1,26 +1,6 @@
 #include "fasttransforms.h"
 #include "ftutilities.h"
 
-void test_dprkf(int * checksum);
-void test_dprk (int * checksum);
-void test_dprkl(int * checksum);
-void test_dprkq(int * checksum);
-
-int main(void) {
-    int checksum = 0;
-    printf("\nTesting methods for symmetric diagonal-plus-rank-k matrices.\n");
-    printf("\n\tSingle precision.\n\n");
-    test_dprkf(&checksum);
-    printf("\n\tDouble precision.\n\n");
-    test_dprk(&checksum);
-    printf("\n\tLong double precision.\n\n");
-    test_dprkl(&checksum);
-    printf("\n\tQuadruple precision.\n\n");
-    test_dprkq(&checksum);
-    printf("\n");
-    return checksum;
-}
-
 #define FLT float
 #define X(name) FT_CONCAT(ft_, name, f)
 #define Y(name) FT_CONCAT(, name, f)
@@ -52,3 +32,18 @@ int main(void) {
 #undef FLT
 #undef X
 #undef Y
+
+int main(void) {
+    int checksum = 0;
+    printf("\nTesting methods for symmetric diagonal-plus-rank-k matrices.\n");
+    printf("\n\tSingle precision.\n\n");
+    test_dprkf(&checksum);
+    printf("\n\tDouble precision.\n\n");
+    test_dprk(&checksum);
+    printf("\n\tLong double precision.\n\n");
+    test_dprkl(&checksum);
+    printf("\n\tQuadruple precision.\n\n");
+    test_dprkq(&checksum);
+    printf("\n");
+    return checksum;
+}

test/test_hierarchical.c

@@ -1,26 +1,6 @@
 #include "fasttransforms.h"
 #include "ftutilities.h"
 
-void test_hierarchicalf(int * checksum);
-void test_hierarchical (int * checksum);
-void test_hierarchicall(int * checksum);
-void test_hierarchicalq(int * checksum);
-
-int main(void) {
-    int checksum = 0;
-    printf("\nTesting methods for hierarchical matrices.\n");
-    printf("\n\tSingle precision.\n\n");
-    test_hierarchicalf(&checksum);
-    printf("\n\tDouble precision.\n\n");
-    test_hierarchical(&checksum);
-    printf("\n\tLong double precision.\n\n");
-    test_hierarchicall(&checksum);
-    printf("\n\tQuadruple precision.\n\n");
-    test_hierarchicalq(&checksum);
-    printf("\n");
-    return checksum;
-}
-
 #define FLT float
 #define X(name) FT_CONCAT(ft_, name, f)
 #define Y(name) FT_CONCAT(, name, f)
@@ -52,3 +32,18 @@ int main(void) {
 #undef FLT
 #undef X
 #undef Y
+
+int main(void) {
+    int checksum = 0;
+    printf("\nTesting methods for hierarchical matrices.\n");
+    printf("\n\tSingle precision.\n\n");
+    test_hierarchicalf(&checksum);
+    printf("\n\tDouble precision.\n\n");
+    test_hierarchical(&checksum);
+    printf("\n\tLong double precision.\n\n");
+    test_hierarchicall(&checksum);
+    printf("\n\tQuadruple precision.\n\n");
+    test_hierarchicalq(&checksum);
+    printf("\n");
+    return checksum;
+}

test/test_tdc.c

@@ -1,35 +1,6 @@
 #include "fasttransforms.h"
 #include "ftutilities.h"
 
-void test_tdcf(int * checksum);
-void test_tdc (int * checksum);
-void test_tdcl(int * checksum);
-void test_tdcq(int * checksum);
-
-void test_tdc_drop_precisionf(int * checksum);
-void test_tdc_drop_precision (int * checksum);
-
-int main(void) {
-    int checksum = 0;
-    printf("\nTesting methods for symmetric-definite tridiagonal divide and conquer.\n");
-    printf("\n\tSingle precision.\n\n");
-    test_tdcf(&checksum);
-    printf("\n\tDouble precision.\n\n");
-    test_tdc(&checksum);
-    printf("\n\tLong double precision.\n\n");
-    test_tdcl(&checksum);
-    printf("\n\tQuadruple precision.\n\n");
-    test_tdcq(&checksum);
-    printf("\n");
-    printf("\nTesting methods for dropping the precision.\n");
-    printf("\n\tDouble ↘  single precision.\n\n");
-    test_tdc_drop_precisionf(&checksum);
-    printf("\n\tLong double ↘  double precision.\n\n");
-    test_tdc_drop_precision(&checksum);
-    printf("\n");
-    return 0;
-}
-
 #define FLT quadruple
 #define X(name) FT_CONCAT(ft_, name, q)
 #define Y(name) FT_CONCAT(, name, q)
@@ -67,3 +38,24 @@ int main(void) {
 #undef X
 #undef X2
 #undef Y
+
+int main(void) {
+    int checksum = 0;
+    printf("\nTesting methods for symmetric-definite tridiagonal divide and conquer.\n");
+    printf("\n\tSingle precision.\n\n");
+    test_tdcf(&checksum);
+    printf("\n\tDouble precision.\n\n");
+    test_tdc(&checksum);
+    printf("\n\tLong double precision.\n\n");
+    test_tdcl(&checksum);
+    printf("\n\tQuadruple precision.\n\n");
+    test_tdcq(&checksum);
+    printf("\n");
+    printf("\nTesting methods for dropping the precision.\n");
+    printf("\n\tDouble ↘  single precision.\n\n");
+    test_tdc_drop_precisionf(&checksum);
+    printf("\n\tLong double ↘  double precision.\n\n");
+    test_tdc_drop_precision(&checksum);
+    printf("\n");
+    return 0;
+}