fix associated transforms with eigenvalues with algebraic and geometric multiplicity greater than 1

Use the true eigenvectors to correct the singularities in the low-rank block V_{12} = \hat{V}_{12} + S_{12}, where S_{12} is the sparse correction

If the geometric multiplicity != algebraic multiplicity, there is a problem!

Author: MikaelSlevinsky

examples/holomorphic.c

@@ -98,7 +98,7 @@ int main(void) {
 
     printf("The Zernike coefficients are useful for integration. The integral\n");
     printf("of "MAGENTA("f(x,y)")" over the disk should be "MAGENTA("π/2")" by harmonicity.\n");
-    printf("The coefficient of "MAGENTA("Z_0^0")" multiplied by "MAGENTA("√π")" is: ");
+    printf("The coefficient of "MAGENTA("Z_{0,0}")" multiplied by "MAGENTA("√π")" is: ");
     printf(FMT, F[0]*sqrt(M_PI));
     printf(".\n\n");
     printf("Using an orthonormal basis, the integral of "MAGENTA("[f(x,y)]^2")" over the disk is\n");

src/banded_source.c

@@ -1,3 +1,10 @@
+void X(destroy_sparse)(X(sparse) * A) {
+    free(A->p);
+    free(A->q);
+    free(A->v);
+    free(A);
+}
+
 void X(destroy_banded)(X(banded) * A) {
     free(A->data);
     free(A);
@@ -23,6 +30,7 @@ void X(destroy_tb_eigen_FMM)(X(tb_eigen_FMM) * F) {
         X(destroy_hierarchicalmatrix)(F->F0);
         X(destroy_tb_eigen_FMM)(F->F1);
         X(destroy_tb_eigen_FMM)(F->F2);
+        X(destroy_sparse)(F->S);
         free(F->X);
         free(F->Y);
         free(F->t1);
@@ -75,6 +83,28 @@ size_t X(summary_size_tb_eigen_ADI)(X(tb_eigen_ADI) * F) {
     return S;
 }
 
+X(sparse) * X(malloc_sparse)(const int m, const int n, const int nnz) {
+    X(sparse) * A = malloc(sizeof(X(sparse)));
+    A->p = malloc(nnz*sizeof(int));
+    A->q = malloc(nnz*sizeof(int));
+    A->v = malloc(nnz*sizeof(FLT));
+    A->m = m;
+    A->n = n;
+    A->nnz = nnz;
+    return A;
+}
+
+X(sparse) * X(calloc_sparse)(const int m, const int n, const int nnz) {
+    X(sparse) * A = malloc(sizeof(X(sparse)));
+    A->p = calloc(nnz, sizeof(int));
+    A->q = calloc(nnz, sizeof(int));
+    A->v = calloc(nnz, sizeof(FLT));
+    A->m = m;
+    A->n = n;
+    A->nnz = nnz;
+    return A;
+}
+
 X(banded) * X(malloc_banded)(const int m, const int n, const int l, const int u) {
     FLT * data = malloc(n*(l+u+1)*sizeof(FLT));
     X(banded) * A = malloc(sizeof(X(banded)));
@@ -488,7 +518,7 @@ void X(triangular_banded_eigenvectors)(X(triangular_banded) * A, X(triangular_ba
             }
             d = lam*X(get_triangular_banded_index)(B, i, i) - X(get_triangular_banded_index)(A, i, i);
             kd = Y(fabs)(lam*X(get_triangular_banded_index)(B, i, i)) + Y(fabs)(X(get_triangular_banded_index)(A, i, i));
-            if (Y(fabs)(d) < 4*kd*Y(eps)() && Y(fabs)(t) < 4*kt*Y(eps)())
+            if (Y(fabs)(d) < 4*kd*Y(eps)() || Y(fabs)(t) < 4*kt*Y(eps)())
                 V[i+j*n] = 0;
             else
                 V[i+j*n] = t/d;
@@ -533,6 +563,83 @@ void X(triangular_banded_quadratic_eigenvectors)(X(triangular_banded) * A, X(tri
     }
 }
 
+// Assumptions: x, y are non-decreasing.
+static inline int X(count_intersections)(const int m, const FLT * x, const int n, const FLT * y, const FLT epsilon) {
+    int istart = 0, idx = 0;
+    for (int j = 0; j < n; j++) {
+        int i = istart;
+        int thefirst = 1;
+        while (i < m) {
+            if (Y(fabs)(x[i] - y[j]) < epsilon*MAX(Y(fabs)(x[i]), Y(fabs)(y[j]))) {
+                idx++;
+                if (thefirst) {
+                    istart = i;
+                    thefirst--;
+                }
+            }
+            else if (x[i] > y[j])
+                break;
+            i++;
+        }
+    }
+    return idx;
+}
+
+// Assumptions: p and q have been malloc'ed with `idx` integers.
+static inline void X(produce_intersection_indices)(const int m, const FLT * x, const int n, const FLT * y, const FLT epsilon, int * p, int * q) {
+    int istart = 0, idx = 0;
+    for (int j = 0; j < n; j++) {
+        int i = istart;
+        int thefirst = 1;
+        while (i < m) {
+            if (Y(fabs)(x[i] - y[j]) < epsilon*MAX(Y(fabs)(x[i]), Y(fabs)(y[j]))) {
+                p[idx] = i;
+                q[idx] = j;
+                idx++;
+                if (thefirst) {
+                    istart = i;
+                    thefirst--;
+                }
+            }
+            else if (x[i] > y[j])
+                break;
+            i++;
+        }
+    }
+}
+
+static inline X(sparse) * X(get_sparse_from_eigenvectors)(X(tb_eigen_FMM) * F1, X(triangular_banded) * A, X(triangular_banded) * B, FLT * D, int * p1, int * p2, int * p3, int * p4, int n, int s, int b, int idx) {
+    X(sparse) * S = X(malloc_sparse)(s, n-s, idx);
+    FLT * V = calloc(n, sizeof(FLT));
+    for (int l = 0; l < idx; l++) {
+        int j = p2[p4[l]]+s;
+        for (int i = 0; i < n; i++)
+            V[i] = 0;
+        V[j] = D[j];
+        FLT t, kt, d, kd, lam;
+        lam = X(get_triangular_banded_index)(A, j, j)/X(get_triangular_banded_index)(B, j, j);
+        for (int i = j-1; i >= 0; i--) {
+            t = kt = 0;
+            for (int k = i+1; k < MIN(i+b+1, n); k++) {
+                t += (X(get_triangular_banded_index)(A, i, k) - lam*X(get_triangular_banded_index)(B, i, k))*V[k];
+                kt += (Y(fabs)(X(get_triangular_banded_index)(A, i, k)) + Y(fabs)(lam*X(get_triangular_banded_index)(B, i, k)))*Y(fabs)(V[k]);
+            }
+            d = lam*X(get_triangular_banded_index)(B, i, i) - X(get_triangular_banded_index)(A, i, i);
+            kd = Y(fabs)(lam*X(get_triangular_banded_index)(B, i, i)) + Y(fabs)(X(get_triangular_banded_index)(A, i, i));
+            if (Y(fabs)(d) < 4*kd*Y(eps)() || Y(fabs)(t) < 4*kt*Y(eps)())
+                V[i] = 0;
+            else
+                V[i] = t/d;
+        }
+        X(bfsv)('N', F1, V);
+        S->p[l] = p1[p3[l]];
+        S->q[l] = p2[p4[l]];
+        S->v[l] = V[p1[p3[l]]];
+    }
+    free(V);
+    return S;
+}
+
 X(tb_eigen_FMM) * X(tb_eig_FMM)(X(triangular_banded) * A, X(triangular_banded) * B, FLT * D) {
     int n = A->n, b1 = A->b, b2 = B->b;
     int b = MAX(b1, b2);
@@ -599,9 +706,18 @@ X(tb_eigen_FMM) * X(tb_eig_FMM)(X(triangular_banded) * A, X(triangular_banded) *
             p2[i] = i;
         X(quicksort_1arg)(lambda+s, p2, 0, n-s-1, X(lt));
 
-        F->F0 = X(sample_hierarchicalmatrix)(X(cauchykernel), lambda, lambda, i, j, 'G');
+        int idx = X(count_intersections)(s, lambda, n-s, lambda+s, 16*Y(sqrt)(Y(eps)()));
+        int * p3 = malloc(idx*sizeof(int));
+        int * p4 = malloc(idx*sizeof(int));
+        X(produce_intersection_indices)(s, lambda, n-s, lambda+s, 16*Y(sqrt)(Y(eps)()), p3, p4);
+        X(sparse) * S = X(get_sparse_from_eigenvectors)(F->F1, A, B, D, p1, p2, p3, p4, n, s, b, idx);
+        free(p3);
+        free(p4);
+
+        F->F0 = X(sample_hierarchicalmatrix)(X(thresholded_cauchykernel), lambda, lambda, i, j, 'G');
         F->X = X;
         F->Y = Y;
+        F->S = S;
         F->t1 = calloc(s*FT_GET_MAX_THREADS(), sizeof(FLT));
         F->t2 = calloc((n-s)*FT_GET_MAX_THREADS(), sizeof(FLT));
         X(perm)('T', lambda, p1, s);
@@ -860,6 +976,7 @@ void X(bfmv)(char TRANS, X(tb_eigen_FMM) * F, FLT * x) {
         int s = n>>1, b = F->b;
         FLT * t1 = F->t1+s*FT_GET_THREAD_NUM(), * t2 = F->t2+(n-s)*FT_GET_THREAD_NUM();
         int * p1 = F->p1, * p2 = F->p2;
+        X(sparse) * S = F->S;
         if (TRANS == 'N') {
             // C(Λ₁, Λ₂) ∘ (-XYᵀ)
             for (int k = 0; k < b; k++) {
@@ -869,6 +986,8 @@ void X(bfmv)(char TRANS, X(tb_eigen_FMM) * F, FLT * x) {
                 for (int i = 0; i < s; i++)
                     x[p1[i]] += t1[i]*F->X[p1[i]+k*s];
             }
+            for (int l = 0; l < S->nnz; l++)
+                x[S->p[l]] += S->v[l]*x[S->q[l]+s];
             X(bfmv)(TRANS, F->F1, x);
             X(bfmv)(TRANS, F->F2, x+s);
         }
@@ -883,6 +1002,8 @@ void X(bfmv)(char TRANS, X(tb_eigen_FMM) * F, FLT * x) {
                 for (int i = 0; i < n-s; i++)
                     x[p2[i]+s] += t2[i]*F->Y[p2[i]+k*(n-s)];
             }
+            for (int l = 0; l < S->nnz; l++)
+                x[S->q[l]+s] += S->v[l]*x[S->p[l]];
         }
     }
 }
@@ -916,6 +1037,7 @@ void X(bfsv)(char TRANS, X(tb_eigen_FMM) * F, FLT * x) {
         int s = n>>1, b = F->b;
         FLT * t1 = F->t1+s*FT_GET_THREAD_NUM(), * t2 = F->t2+(n-s)*FT_GET_THREAD_NUM();
         int * p1 = F->p1, * p2 = F->p2;
+        X(sparse) * S = F->S;
         if (TRANS == 'N') {
             X(bfsv)(TRANS, F->F1, x);
             X(bfsv)(TRANS, F->F2, x+s);
@@ -927,6 +1049,8 @@ void X(bfsv)(char TRANS, X(tb_eigen_FMM) * F, FLT * x) {
                 for (int i = 0; i < s; i++)
                     x[p1[i]] += t1[i]*F->X[p1[i]+k*s];
             }
+            for (int l = 0; l < S->nnz; l++)
+                x[S->p[l]] -= S->v[l]*x[S->q[l]+s];
         }
         else if (TRANS == 'T') {
             // C(Λ₁, Λ₂) ∘ (-XYᵀ)
@@ -937,6 +1061,8 @@ void X(bfsv)(char TRANS, X(tb_eigen_FMM) * F, FLT * x) {
                 for (int i = 0; i < n-s; i++)
                     x[p2[i]+s] += t2[i]*F->Y[p2[i]+k*(n-s)];
             }
+            for (int l = 0; l < S->nnz; l++)
+                x[S->q[l]+s] -= S->v[l]*x[S->p[l]];
             X(bfsv)(TRANS, F->F1, x);
             X(bfsv)(TRANS, F->F2, x+s);
         }

src/banded_source.h

@@ -1,3 +1,12 @@
+typedef struct {
+    int * p;
+    int * q;
+    FLT * v;
+    int m;
+    int n;
+    int nnz;
+} X(sparse);
+
 typedef struct {
     FLT * data;
     int m;
@@ -23,6 +32,7 @@ struct X(tbstruct_FMM) {
     X(hierarchicalmatrix) * F0;
     X(tb_eigen_FMM) * F1;
     X(tb_eigen_FMM) * F2;
+    X(sparse) * S;
     FLT * V;
     FLT * X;
     FLT * Y;
@@ -47,6 +57,7 @@ struct X(tbstruct_ADI) {
     int b;
 };
 
+void X(destroy_sparse)(X(sparse) * A);
 void X(destroy_banded)(X(banded) * A);
 void X(destroy_triangular_banded)(X(triangular_banded) * A);
 void X(destroy_banded_qr)(X(banded_qr) * F);
@@ -56,6 +67,9 @@ void X(destroy_tb_eigen_ADI)(X(tb_eigen_ADI) * F);
 size_t X(summary_size_tb_eigen_FMM)(X(tb_eigen_FMM) * F);
 size_t X(summary_size_tb_eigen_ADI)(X(tb_eigen_ADI) * F);
 
+X(sparse) * X(malloc_sparse)(const int m, const int n, const int nnz);
+X(sparse) * X(calloc_sparse)(const int m, const int n, const int nnz);
+
 X(banded) * X(malloc_banded)(const int m, const int n, const int l, const int u);
 X(banded) * X(calloc_banded)(const int m, const int n, const int l, const int u);
 

src/drop_precision.c

@@ -129,6 +129,18 @@ X(symmetric_dpr1_eigen_FMM) * X(drop_precision_symmetric_dpr1_eigen_FMM)(X2(symm
     return F;
 }
 
+X(sparse) * X(drop_precision_sparse)(X2(sparse) * S2) {
+    X(sparse) * S = X(malloc_sparse)(S2->m, S2->n, S2->nnz);
+    for (int l = 0; l < S->nnz; l++) {
+        S->p[l] = S2->p[l];
+        S->q[l] = S2->q[l];
+        S->v[l] = S2->v[l];
+    }
+    return S;
+}
+
+static inline FLT X(thresholded_cauchykernel2)(FLT x, FLT y) {return X2(thresholded_cauchykernel)(x, y);}
+
 X(tb_eigen_FMM) * X(drop_precision_tb_eigen_FMM)(X2(tb_eigen_FMM) * F2) {
     int n = F2->n;
     X(tb_eigen_FMM) * F = malloc(sizeof(X(tb_eigen_FMM)));
@@ -155,11 +167,12 @@ X(tb_eigen_FMM) * X(drop_precision_tb_eigen_FMM)(X2(tb_eigen_FMM) * F2) {
             lambda[i] = F2->lambda[i];
         X(perm)('N', lambda, p1, s);
         X(perm)('N', lambda+s, p2, n-s);
-        F->F0 = X(sample_hierarchicalmatrix)(X(cauchykernel), lambda, lambda, (unitrange) {0, s}, (unitrange) {s, n}, 'G');
+        F->F0 = X(sample_hierarchicalmatrix)(X(thresholded_cauchykernel2), lambda, lambda, (unitrange) {0, s}, (unitrange) {s, n}, 'G');
         X(perm)('T', lambda, p1, s);
         X(perm)('T', lambda+s, p2, n-s);
         F->F1 = X(drop_precision_tb_eigen_FMM)(F2->F1);
         F->F2 = X(drop_precision_tb_eigen_FMM)(F2->F2);
+        F->S = X(drop_precision_sparse)(F2->S);
         F->X = malloc(s*b*sizeof(FLT));
         for (int i = 0; i < s*b; i++)
             F->X[i] = F2->X[i];

src/drop_precision.h

@@ -4,5 +4,6 @@ X(tdc_eigen_FMM) * X(drop_precision_tdc_eigen_FMM)(X2(tdc_eigen_FMM) * F2);
 X(symmetric_dpr1_eigen) * X(drop_precision_symmetric_dpr1_eigen)(X2(symmetric_dpr1_eigen) * F2);
 X(symmetric_dpr1_eigen_FMM) * X(drop_precision_symmetric_dpr1_eigen_FMM)(X2(symmetric_dpr1_eigen_FMM) * F2);
 
+X(sparse) * X(drop_precision_sparse)(X2(sparse) * S2);
 X(tb_eigen_FMM) * X(drop_precision_tb_eigen_FMM)(X2(tb_eigen_FMM) * F2);
 X(btb_eigen_FMM) * X(drop_precision_btb_eigen_FMM)(X2(btb_eigen_FMM) * F2);

src/hierarchical_source.c

@@ -937,3 +937,10 @@ FLT X(cauchykernel2)(FLT x, FLT ylo, FLT yhi) {return 1/(X(diff)(x, yhi) - ylo);
 FLT X(coulombkernel2)(FLT x, FLT ylo, FLT yhi) {return 1/((X(diff)(x, yhi) - ylo)*(X(diff)(x, yhi) - ylo));}
 FLT X(coulombprimekernel2)(FLT x, FLT ylo, FLT yhi) {return 1/(((X(diff)(x, yhi) - ylo)*(X(diff)(x, yhi) - ylo))*(X(diff)(x, yhi) - ylo));}
 FLT X(logkernel2)(FLT x, FLT ylo, FLT yhi) {return Y(log)(Y(fabs)(X(diff)(x, yhi) - ylo));}
+
+FLT X(thresholded_cauchykernel)(FLT x, FLT y) {
+    if (Y(fabs)(x - y) < 16*Y(sqrt)(Y(eps)())*MAX(Y(fabs)(x), Y(fabs)(y)))
+        return 0;
+    else
+        return 1/(x-y);
+}

src/hierarchical_source.h

@@ -114,3 +114,5 @@ FLT X(cauchykernel2)(FLT x, FLT ylo, FLT yhi);
 FLT X(coulombkernel2)(FLT x, FLT ylo, FLT yhi);
 FLT X(coulombprimekernel2)(FLT x, FLT ylo, FLT yhi);
 FLT X(logkernel2)(FLT x, FLT ylo, FLT yhi);
+
+FLT X(thresholded_cauchykernel)(FLT x, FLT y);