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 eb894074f8b28a31f2a4ea04d4b0f10e49c187b7..a16bd785e4a85230ee99ddefdf484462ee86c5e2 100644 (file)
--- 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 05be86d40a6d6775db90f36287e1f6010cc39556..f5fefd07139311c4ed47addad027fae4e3306d06 100644 (file)
--- 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 69e7436453e73f360801fb7ad7b8623d0276475c..39827129e42780d612a44fe3b6a553bfabb1e56a 100644 (file)
--- 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);