add LU for coarsest level solver

example/amg.c

@@ -10,7 +10,7 @@ int main(void)
     char *mat_file = "A.dat";
     SX_RTN rtn;
     SX_INT prob = 3;
-    SX_INT ncx = 23, ncy = 13, ncz = 14;
+    SX_INT ncx = 130, ncy = 130, ncz = 14;
     SX_INT nglobal = 0;
 
     /* create distributed matrix */

include/amg-coarse-solver.h

@@ -10,7 +10,7 @@
 extern "C" {
 #endif
 
-void sx_amg_coarest_solve(SX_MAT *A, SX_VEC *b, SX_VEC *x, const SX_FLT ctol, const SX_INT verb);
+void sx_amg_coarest_solve(SX_AMG_COMP *cg, const SX_FLT ctol, const SX_INT verb);
 
 #ifdef __cplusplus
 }

include/amg-type-def.h

@@ -194,6 +194,9 @@ typedef struct SX_AMG_COMP_
 
     SX_VEC wp;                   /** cache work space */
 
+    SX_FLT *LU;
+    int *pvt;
+
 } SX_AMG_COMP;
 
 /** 
@@ -225,4 +228,32 @@ typedef struct SX_KRYLOV_
 
 } SX_KRYLOV;
 
+typedef enum
+{
+    ERROR_OPEN_FILE          = -10,  /**< fail to open a file */
+    ERROR_WRONG_FILE         = -11,  /**< input contains wrong format */
+    ERROR_INPUT_PAR          = -12,  /**< wrong input argument */
+    ERROR_MAT_SIZE           = -13,  /**< wrong problem size */
+    ERROR_MISC               = -14,  /**< other error */
+
+    ERROR_ALLOC_MEM          = -20,  /**< fail to allocate memory */
+    ERROR_DATA_STRUCTURE     = -21,  /**< problem with data structures */
+    ERROR_DATA_ZERODIAG      = -22,  /**< matrix has zero diagonal entries */
+    ERROR_DUMMY_VAR          = -23,  /**< unexpected input data */
+
+    ERROR_AMG_INTERP_TYPE    = -30,  /**< unknown interpolation type */
+    ERROR_AMG_SMOOTH_TYPE    = -31,  /**< unknown smoother type */
+    ERROR_AMG_COARSE_TYPE    = -32,  /**< unknown coarsening type */
+    ERROR_AMG_COARSEING      = -33,  /**< coarsening step failed to complete */
+
+    ERROR_SOLVER_STAG        = -42,  /**< solver stagnates */
+    ERROR_SOLVER_SOLSTAG     = -43,  /**< solver's solution is too small */
+    ERROR_SOLVER_TOLSMALL    = -44,  /**< solver's tolerance is too small */
+    ERROR_SOLVER_MAXIT       = -48,  /**< maximal iteration number exceeded */
+    ERROR_SOLVER_EXIT        = -49,  /**< solver does not quit successfully */
+
+    ERROR_UNKNOWN            = -99,  /**< an unknown error type */
+
+} SX_ERROR_CODE;
+
 #endif

include/utils.h

@@ -21,6 +21,15 @@ SX_FLT sx_get_mem(void);
 
 void sx_exit_on_errcode(const SX_INT status, const char *fctname);
 
+/* LU factorization */
+int sx_mat_dense_lu(int n, SX_FLT *A, int pvt[]);
+
+/* solves L*U*X = B  */
+void sx_mat_dense_sv(int n, SX_FLT *A, int pvt[], int m, SX_FLT *B);
+
+/* solves AX = B, returns 1 if successful and 0 if A is singular */
+int sx_mat_dense_solve(int n, int m, SX_FLT *A, SX_FLT *B);
+
 #ifdef __cplusplus
 }
 #endif

src/amg-coarse-solver.c

@@ -582,43 +582,91 @@ static SX_INT sx_solver_gmres(SX_KRYLOV *ks)
 }
 
 /**
- * \fn static void sx_amg_coarest_solve(SX_MAT *A, SX_VEC *b, SX_VEC *x,
- *                                      const SX_FLT ctol, const SX_INT verb)
+ * \fn static void sx_amg_coarest_solve(SX_AMG_COMP *cg, const SX_FLT ctol, const SX_INT verb)
  *
  * \brief Iterative on the coarset level
  *
- * \pars  A         pointer to matrix data
- * \pars  b         pointer to rhs data
- * \pars  x         pointer to sol data
- * \pars  ctol      tolerance for the coarsest level
+ * \pars  cg     pointer to matrix data
+ * \pars  ctol   tolerance for the coarsest level
  * \pars  verb   level of output
  *
  */
-void sx_amg_coarest_solve(SX_MAT *A, SX_VEC *b, SX_VEC *x, const SX_FLT ctol, const SX_INT verb)
+void sx_amg_coarest_solve(SX_AMG_COMP *cg, const SX_FLT ctol, const SX_INT verb)
 {
+    SX_MAT *A = &cg->A;
+    SX_VEC *b = &cg->b;
+    SX_VEC *x = &cg->x;
     const SX_INT n = A->num_rows;
     const SX_INT maxit = SX_MAX(250, SX_MIN(n * n, 1000));
+    int use_krylov = 0;
 
-    SX_INT status;
-    SX_KRYLOV ks;
-
-    /* try cg first */
-    ks.A = A;
-    ks.b = b;
-    ks.u = x;
-    ks.tol = ctol;
-    ks.maxit = maxit;
-    ks.stop_type = 1;
-    ks.verb = 0;
-    status = sx_solver_cg(&ks);
-
-    /* try GMRES if cg fails */
-    if (status < 0) {
-        ks.restart = MAX_RESTART;
-        status = sx_solver_gmres(&ks);
+    /* init */
+    if (!use_krylov) {
+        if (cg->LU == NULL) {
+            SX_INT i, k, j, ibegin, iend;
+            int rtn;
+
+            cg->LU = sx_calloc(n * (size_t) n, sizeof(*cg->LU));
+            cg->pvt = sx_calloc(n, sizeof(*cg->pvt));
+
+            /* init LU */
+            for (i = 0; i < n; ++i) {
+                ibegin = A->Ap[i];
+                iend = A->Ap[i + 1];
+                for (k = ibegin; k < iend; ++k) {
+                    j = A->Aj[k];
+
+                    cg->LU[i * n + j] = A->Ax[k];
+                }
+            }
+
+            /* LU factorization */
+            rtn = sx_mat_dense_lu(n, cg->LU, cg->pvt);
+
+            /* failed */
+            if (!rtn) {
+                use_krylov = 1;
+
+                sx_free(cg->LU);
+                sx_free(cg->pvt);
+                cg->LU = NULL;
+                cg->pvt = NULL;
+            }
+        }
+    }
+
+    if (use_krylov) {
+        SX_INT status;
+        SX_KRYLOV ks;
+
+        /* try cg first */
+        ks.A = A;
+        ks.b = b;
+        ks.u = x;
+        ks.tol = ctol;
+        ks.maxit = maxit;
+        ks.stop_type = 1;
+        ks.verb = 0;
+        status = sx_solver_cg(&ks);
+
+        /* try GMRES if cg fails */
+        if (status < 0) {
+            ks.restart = MAX_RESTART;
+            status = sx_solver_gmres(&ks);
+        }
+
+        if (status < 0 && verb >= 2) {
+            sx_printf("### WARNING: Coarse level solver failed to converge!\n");
+        }
     }
+    else {
+        assert(cg->LU != NULL);
+        assert(cg->pvt != NULL);
+
+        /* init x */
+        memcpy(x->d, b->d, n * sizeof(*x->d));
 
-    if (status < 0 && verb >= 2) {
-        sx_printf("### WARNING: Coarse level solver failed to converge!\n");
+        /* solve */
+        sx_mat_dense_sv(n, cg->LU, cg->pvt, 1, x->d);
     }
 }

src/amg-cycle.c

@@ -59,7 +59,7 @@ void sx_amg_cycle(SX_AMG *mg)
     }
 
     // call the coarse space solver:
-    sx_amg_coarest_solve(&mg->cg[nl - 1].A, &mg->cg[nl - 1].b, &mg->cg[nl - 1].x, tol, verb);
+    sx_amg_coarest_solve(&mg->cg[nl - 1], tol, verb);
 
     // BackwardSweep:
     while (l > 0) {

src/amg-utils.c

@@ -58,6 +58,9 @@ void sx_amg_data_destroy(SX_AMG *mg)
 
         sx_vec_destroy(&mg->cg[i].wp);
         sx_ivec_destroy(&mg->cg[i].cfmark);
+
+        sx_free(mg->cg[i].LU);
+        sx_free(mg->cg[i].pvt);
     }
 
     sx_free(mg->cg);

src/internal.h

@@ -46,34 +46,6 @@ typedef enum
 
 } SX_STOP_TYPE;
 
-typedef enum
-{
-    ERROR_OPEN_FILE          = -10,  /**< fail to open a file */
-    ERROR_WRONG_FILE         = -11,  /**< input contains wrong format */
-    ERROR_INPUT_PAR          = -12,  /**< wrong input argument */
-    ERROR_MAT_SIZE           = -13,  /**< wrong problem size */
-    ERROR_MISC               = -14,  /**< other error */
-
-    ERROR_ALLOC_MEM          = -20,  /**< fail to allocate memory */
-    ERROR_DATA_STRUCTURE     = -21,  /**< problem with data structures */
-    ERROR_DATA_ZERODIAG      = -22,  /**< matrix has zero diagonal entries */
-    ERROR_DUMMY_VAR          = -23,  /**< unexpected input data */
-
-    ERROR_AMG_INTERP_TYPE    = -30,  /**< unknown interpolation type */
-    ERROR_AMG_SMOOTH_TYPE    = -31,  /**< unknown smoother type */
-    ERROR_AMG_COARSE_TYPE    = -32,  /**< unknown coarsening type */
-    ERROR_AMG_COARSEING      = -33,  /**< coarsening step failed to complete */
-
-    ERROR_SOLVER_STAG        = -42,  /**< solver stagnates */
-    ERROR_SOLVER_SOLSTAG     = -43,  /**< solver's solution is too small */
-    ERROR_SOLVER_TOLSMALL    = -44,  /**< solver's tolerance is too small */
-    ERROR_SOLVER_MAXIT       = -48,  /**< maximal iteration number exceeded */
-    ERROR_SOLVER_EXIT        = -49,  /**< solver does not quit successfully */
-
-    ERROR_UNKNOWN            = -99,  /**< an unknown error type */
-
-} SX_ERROR_CODE;
-
 #ifdef __cplusplus
 extern "C" {
 #endif

src/utils.c

@@ -342,3 +342,108 @@ void sx_iarray_cp(const SX_INT n, SX_INT *x, SX_INT *y)
 
     memcpy(y, x, n * sizeof(SX_INT));
 }
+
+/* LU factorization */
+int sx_mat_dense_lu(int n, SX_FLT *a0, int pvt[])
+{
+    int i, j, k;
+    SX_FLT d, r;
+
+#define a(i,j)	(a0[(i) * n + (j)])
+
+    for (i = 0; i < n; i++) {
+        d = fabs(a(i,i));
+        k = i;
+        for (j = i + 1; j < n; j++) {
+            if ((r = fabs(a(j,i))) > d) {
+                d = r;
+                k = j;
+            }
+        }
+
+        if (d == 0.0) {
+            return 0;
+        }
+
+        pvt[i] = k;
+        if (k != i) {
+            /* exchange row i and row k */
+            for (j = i; j < n; j++) {
+                d = a(i,j);
+                a(i,j) = a(k,j);
+                a(k,j) = d;
+            }
+        }
+
+        if ((d = a(i,i)) != 1.0) {
+            a(i,i) = d = 1.0 / d;
+            for (j = i + 1; j < n; j++) a(i,j) *= d;
+        }
+
+        for (j = i + 1; j < n; j++) {
+            if ((d = a(j,i)) == 0.0) continue;
+
+            for (k = i + 1; k < n; k++) a(j,k) -= d * a(i,k);
+        }
+    }
+
+#undef a
+
+    return 1;
+}
+
+/* solves L*U*X = B  */
+void sx_mat_dense_sv(int n, SX_FLT *a0, int pvt[], int m, SX_FLT *b0)
+{
+    int i, j, k;
+    SX_FLT d;
+
+#define a(i,j)	(a0[(i) * n + (j)])
+#define b(i,j)	(b0[(i) * m + (j)])
+
+    for (i = 0; i < n; i++) {
+        k = pvt[i];
+        if (k != i) {
+            /* exchange row i with row k */
+            for (j = 0; j < m; j++) {
+                d = b(i,j);
+                b(i,j) = b(k,j);
+                b(k,j) = d;
+            }
+        }
+
+        if ((d = a(i,i)) != 1.0) {
+            for (j = 0; j < m; j++)
+                b(i,j) *= d;
+        }
+
+        for (j = i + 1; j < n; j++) {
+            if ((d = a(j,i)) == 0.0) continue;
+
+            for (k = 0; k < m; k++) b(j,k) -= d * b(i,k);
+        }
+    }
+
+    for (i = n - 2; i >= 0; i--) {
+        for (j = i + 1; j < n; j++) {
+            d = a(i,j);
+            for (k = 0; k < m; k++) b(i,k) -= d * b(j,k);
+        }
+    }
+
+#undef a
+#undef b
+}
+
+/* solves AX = B, returns 1 if successful and 0 if A is singular */
+int sx_mat_dense_solve(int n, int m, SX_FLT *a0, SX_FLT *b0)
+{
+    int *pvt = malloc(n * sizeof(*pvt));
+
+    if (sx_mat_dense_lu(n, a0, pvt) == 0) return 0;
+
+    sx_mat_dense_sv(n, a0, pvt, m, b0);
+    free(pvt);
+
+    return 1;
+}