From: Matthew Fernandez <matthew.fernandez@gmail.com>
Date: Sat, 30 Jul 2022 00:30:23 +0000 (-0700)
Subject: sfdpgen SparseMatrix_solve: remove unused Jacobi matrix solving
X-Git-Tag: 5.0.1~22^2~9
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=05dab56d6568697188058fcfa191d24129548e9f;p=graphviz

sfdpgen SparseMatrix_solve: remove unused Jacobi matrix solving

This also removes the `method` parameter to this function, which was only ever
set to a single value. This squashes 8 compiler warnings and drops a lot of dead
code.
---

diff --git a/lib/sfdpgen/post_process.c b/lib/sfdpgen/post_process.c
index eb894074f..a16bd785e 100644
--- a/lib/sfdpgen/post_process.c
+++ b/lib/sfdpgen/post_process.c
@@ -704,8 +704,7 @@ double StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim,
     if (sm->scheme == SM_SCHEME_UNIFORM_STRESS){
       uniform_stress_solve(Lw, alpha, dim, x, y, sm->tol_cg, sm->maxit_cg);
     } else {
-      SparseMatrix_solve(Lw, dim, x, y,  sm->tol_cg, sm->maxit_cg, SOLVE_METHOD_CG, &flag);
-      //SparseMatrix_solve(Lw, dim, x, y,  sm->tol_cg, 1, SOLVE_METHOD_JACOBI, &flag);
+      SparseMatrix_solve(Lw, dim, x, y,  sm->tol_cg, sm->maxit_cg, &flag);
     }
 
     if (flag) goto RETURN;
diff --git a/lib/sfdpgen/sparse_solve.c b/lib/sfdpgen/sparse_solve.c
index 05be86d40..f5fefd071 100644
--- a/lib/sfdpgen/sparse_solve.c
+++ b/lib/sfdpgen/sparse_solve.c
@@ -236,77 +236,17 @@ double cg(Operator Ax, Operator precond, int n, int dim, double *x0, double *rhs
   return res;
 }
 
-static double* jacobi(SparseMatrix A, int dim, double *x0, double *rhs, int maxit) {
-  /* maxit iteration of jacobi */
-  double *x, *y, *b, sum, diag, *a;
-  int k, i, j, n = A->n, *ia, *ja, iter;
-  x = MALLOC(sizeof(double)*n);
-  y = MALLOC(sizeof(double)*n);
-  b = MALLOC(sizeof(double)*n);
-  assert(A->type == MATRIX_TYPE_REAL);
-  ia = A->ia;
-  ja = A->ja;
-  a = A->a;
-
-  for (k = 0; k < dim; k++){
-    for (i = 0; i < n; i++) {
-      x[i] = x0[i*dim+k];
-      b[i] = rhs[i*dim+k];
-    }
-
-    for (iter = 0; iter < maxit; iter++){
-      for (i = 0; i < n; i++){
-	sum = 0; diag = 0;
-	for (j = ia[i]; j < ia[i+1]; j++){
-	  if (ja[j] != i){
-	    sum += a[j]*x[ja[j]];
-	  } else {
-	    diag = a[j];
-	  }
-	}
-	if (sum == 0) fprintf(stderr,"neighb=%d\n",ia[i+1]-ia[i]);
-	assert(diag != 0);
-	y[i] = (b[i] - sum)/diag;
-
-      }
-      memcpy(x, y, sizeof(double)*n);
-    }
-
-    for (i = 0; i < n; i++) {
-      rhs[i*dim+k] = x[i];
-    }
-  }
-
-
-  free(x);
-  free(y);
-  free(b);
-  return rhs;
-}
-
-double SparseMatrix_solve(SparseMatrix A, int dim, double *x0, double *rhs, double tol, int maxit, int method, int *flag){
+double SparseMatrix_solve(SparseMatrix A, int dim, double *x0, double *rhs, double tol, int maxit, int *flag){
   Operator Ax, precond;
   int n = A->m;
   double res = 0;
   *flag = 0;
 
-  switch (method){
-  case SOLVE_METHOD_CG:
-    Ax =  Operator_matmul_new(A);
-    precond = Operator_diag_precon_new(A);
-    res = cg(Ax, precond, n, dim, x0, rhs, tol, maxit);
-    Operator_matmul_delete(Ax);
-    Operator_diag_precon_delete(precond);
-    break;
-  case SOLVE_METHOD_JACOBI:{
-    jacobi(A, dim, x0, rhs, maxit);
-    break;
-  }
-  default:
-    assert(0);
-    break;
-
-  }
+  Ax =  Operator_matmul_new(A);
+  precond = Operator_diag_precon_new(A);
+  res = cg(Ax, precond, n, dim, x0, rhs, tol, maxit);
+  Operator_matmul_delete(Ax);
+  Operator_diag_precon_delete(precond);
   return res;
 }
 
diff --git a/lib/sfdpgen/sparse_solve.h b/lib/sfdpgen/sparse_solve.h
index 69e743645..39827129e 100644
--- a/lib/sfdpgen/sparse_solve.h
+++ b/lib/sfdpgen/sparse_solve.h
@@ -12,8 +12,6 @@
 
 #include <sparse/SparseMatrix.h>
 
-enum {SOLVE_METHOD_CG, SOLVE_METHOD_JACOBI};
-
 typedef struct Operator_struct *Operator;
 
 struct Operator_struct {
@@ -23,7 +21,7 @@ struct Operator_struct {
 
 double cg(Operator Ax, Operator precond, int n, int dim, double *x0, double *rhs, double tol, int maxit);
 
-double SparseMatrix_solve(SparseMatrix A, int dim, double *x0, double *rhs, double tol, int maxit, int method, int *flag);
+double SparseMatrix_solve(SparseMatrix A, int dim, double *x0, double *rhs, double tol, int maxit, int *flag);
 
 Operator Operator_uniform_stress_matmul(SparseMatrix A, double alpha);