From 7e3fe070129a133565f6f4fa8baf317814763914 Mon Sep 17 00:00:00 2001
From: "Jostein R. Natvig" <jostein.r.natvig@sintef.no>
Date: Fri, 6 Aug 2010 10:53:07 +0000
Subject: [PATCH] Add system assembly.

---
 hybsys.c                   | 184 +++++++++++++++++++++++++++++++++++++
 hybsys.h                   |  11 ++-
 mex_schur_comp_symm.c      |  55 ++++++++++-
 test_mex_schur_comp_symm.m |   2 +-
 4 files changed, 249 insertions(+), 3 deletions(-)

diff --git a/hybsys.c b/hybsys.c
index aef139bc..35258934 100644
--- a/hybsys.c
+++ b/hybsys.c
@@ -215,3 +215,187 @@ hybsys_compute_press_flux(int nc, const int *nconn, const int *conn,
         p2 += nconn[c] * nconn[c];
     }
 }
+
+
+/*
+ * Routines to assemble global matrix
+ *
+ */
+
+/* ---------------------------------------------------------------------- */
+static int *
+hybsys_build_ia(int nc, int nf, int *nconn, int *conn)
+/* ---------------------------------------------------------------------- */
+{
+   int    *ia = malloc((nf+1) * sizeof *ia);
+   
+   int i;
+   for(i=0; i<nf+1; ++i)
+   {
+      ia[i] = 0;
+   }
+   
+   /* 
+    *   Compute rowsizes 
+    */   
+   int c, pos = 0; 
+   for(c=0; c<nc; ++c)
+   {
+      int n = nconn[c];
+      for (i=pos; i<pos+n; ++i)
+      {
+         mxAssert(conn[i]<nf, "conn out of bounds");
+         ia[1+conn[i]] += n - 1;
+      }
+      pos += n;
+   }
+   
+   /*
+    *   cumulative sum... 
+    */
+   for(i=1; i<nf+1; ++i)
+   {
+      ia[i] = ia[i-1] + ia[i]+1;
+   }
+   return ia;
+}
+
+/* ---------------------------------------------------------------------- */
+static int*
+hybsys_build_ja(int nc, int nf, int *nconn, int *conn, 
+                              int *ia, int *work)
+/* ---------------------------------------------------------------------- */
+{
+   int *ja = malloc(ia[nf] * sizeof *ja);   
+   int  i,j;
+
+   /*
+    *   For each row, diagonal entries are positioned first.
+    */
+   for(i=0; i<nf; ++i)
+   {
+      work[i] = 1;
+      ja[ia[i]] = i;
+   }
+   
+   int c, pos = 0; 
+   for(c=0; c<nc; ++c)
+   {
+      int n = nconn[c];
+      for (i=pos; i<pos+n; ++i)
+      {
+         int fi = conn[i]; 
+         /* mxAssert(fi<nf, "fi out of bounds"); */
+
+         for (j=pos; j<i; ++j)
+         {
+            int fj = conn[j]; 
+            /* mxAssert(fj<nf, "fj out of bounds"); */
+
+            /*
+             *   No conditionals since off-diagonals entries are
+             *   visited only once.
+             */
+            ja[ia[fi] + work[fi]++] = fj;    
+            ja[ia[fj] + work[fj]++] = fi;
+         }
+      }
+      pos += n;
+   }
+   return ja;
+}
+/* ---------------------------------------------------------------------- */
+static double*
+hybsys_build_sa(int nc, int nf, int *nconn, int *conn, int *ia, 
+                double *S, double *R, int *work)
+/* ---------------------------------------------------------------------- */
+{
+   double *sa = malloc(ia[nf] * sizeof *sa);   
+   int     i,j;
+
+   /*
+    *   Clear diagonal and work array
+    */
+   for(i=0; i<nf; ++i)
+   {
+      work[i]   = 1;
+      sa[ia[i]] = 0;
+      R[i]      = 0;
+   }
+   
+   double *s = S;
+   double *r = R;
+
+   int c, pos = 0; 
+   for(c=0; c<nc; ++c)
+   {
+      int n = nconn[c];
+      for (i=pos; i<pos+n; ++i)
+      {
+         int fi = conn[i];
+         int ii = i-pos;
+
+         for (j=pos; j<i; ++j)
+         {
+            int fj = conn[j]; 
+            int jj = j-pos;
+            
+            /*
+             *   We can use assignment since off-diagonal entries are
+             *   visited only once.
+             */
+            sa[ia[fi] + work[fi]++] = s[ii + jj*n];
+            sa[ia[fj] + work[fj]++] = s[jj + ii*n];
+         }
+
+         /*
+          *   Diagonal entries are more than once.
+          */
+         sa[ia[fi]] += s[ii + ii*n];
+         R[fi]      += r[ii];
+      }
+
+      s   += n*n;
+      pos += n;
+   }
+
+   return sa;
+}   
+
+/* ---------------------------------------------------------------------- */
+void hybsys_build_matrix_structure(int nc, int nf, int *nconn, int *conn, 
+                                   int **ia, int **ja)
+/* ---------------------------------------------------------------------- */
+{
+   int *work = malloc(nf * sizeof *work);
+
+   *ia       = hybsys_build_ia(nc, nf, nconn, conn);
+   *ja       = hybsys_build_ja(nc, nf, nconn, conn, *ia, work);
+
+   free(work);
+}
+/* ---------------------------------------------------------------------- */
+void hybsys_assemble_global_system(int nc, int nf, int *nconn, int *conn, 
+                                   double *S, double *R, 
+                                   double **sa, int *ia)
+/* ---------------------------------------------------------------------- */
+{
+   int *work  = malloc(nf * sizeof *work);
+
+   *sa        = hybsys_build_sa(nc, nf, nconn, conn, ia, S, R, work); 
+
+   free(work);
+}
+
+/* ---------------------------------------------------------------------- */
+struct Sparse*
+hybsys_assemble(int nc, int nf, int *nconn, int *conn, double *S, double *R)
+/* ---------------------------------------------------------------------- */
+{
+   struct Sparse *A  = malloc(sizeof *A);
+   A->m    =  A->n   = nf;
+   hybsys_build_matrix_structure(nc, nf, nconn, conn, &A->ia, &A->ja);
+   hybsys_assemble_global_system(nc, nf, nconn, conn, S, R, &A->sa, A->ia);
+ 
+   return A;
+}
diff --git a/hybsys.h b/hybsys.h
index fa08500d..4f3da3c9 100644
--- a/hybsys.h
+++ b/hybsys.h
@@ -8,7 +8,14 @@ struct hybsys {
     double *S;                  /* system matrix in single cell */
     double *one;                /* ones(max_ncf, 1) */
 };
-
+struct Sparse
+{
+   int     m;
+   int     n;
+   int    *ia;
+   int    *ja; 
+   double *sa;
+};
 struct hybsys *
 hybsys_allocate(int max_ncf, int nc, int ncf_tot);
 
@@ -38,4 +45,6 @@ hybsys_compute_press_flux(int nc, const int *nconn, const int *conn,
                           const double *pi, double *press, double *flux,
                           double *work, const int lwork);
 
+struct Sparse*
+hybsys_assemble(int nc, int nf, int *nconn, int *conn, double *S, double *R);
 #endif  /* HYBSYS_H_INCLUDED */
diff --git a/mex_schur_comp_symm.c b/mex_schur_comp_symm.c
index 8928f5a8..fbeccfca 100644
--- a/mex_schur_comp_symm.c
+++ b/mex_schur_comp_symm.c
@@ -14,7 +14,7 @@ verify_args(int nlhs, int nrhs, const mxArray *prhs[])
 {
     int ok;
 
-    ok =       (nlhs == 4) && (nrhs == 2);
+    ok =       (nlhs == 4) && (nrhs == 3);
     ok = ok && mxIsDouble(prhs[0]);
     ok = ok && (mxIsDouble(prhs[1]) || mxIsInt32(prhs[1]));
     ok = ok && (mxGetNumberOfElements(prhs[0]) >
@@ -124,6 +124,51 @@ get_nconn(const mxArray *M_nconn, int *nconn)
     }
 }
 
+/* ---------------------------------------------------------------------- */
+static void
+get_conn(const mxArray *M_conn, int *conn)
+/* ---------------------------------------------------------------------- */
+{
+    size_t nel, i;
+
+    int    *pi;
+    double *pd;
+
+    nel = mxGetNumberOfElements(M_conn);
+
+    if (mxIsDouble(M_conn)) {
+        pd = mxGetPr(M_conn);
+
+        for (i = 0; i < nel; i++) { conn[i] = pd[i] - 1; }
+    } else {
+        pi = mxGetData(M_conn);
+
+        for (i = 0; i < nel; i++) { conn[i] = pi[i] - 1; }
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+static int 
+get_number_of_faces(int nc, int *nconn, int *conn)
+/* ---------------------------------------------------------------------- */
+{
+   int N = 0;
+   int i;
+
+   for (i=0; i<nc; ++i)
+   {
+      N += nconn[i];
+   }
+
+   int nf=0; 
+   for(i=0; i<N; ++i)
+   {
+      nf = MAX(nf, conn[i]+1);
+   }
+   
+   return nf;
+}
+
 
 /*
  * [S, r, F, L] = mex_schur_comp_symm(BI, nconn)
@@ -173,7 +218,14 @@ mexFunction(int nlhs,       mxArray *plhs[],
             p1 += nconn[c];
             p2 += nconn[c] * nconn[c];
         }
+        
+        int *conn= mxMalloc(mxGetNumberOfElements(prhs[2])*sizeof *conn);
+        get_conn(prhs[2], conn);
 
+        int  nf          = get_number_of_faces(nc, nconn, conn);           
+        struct Sparse *A = hybsys_assemble(nc, nf, nconn, conn, ptr, sys->r);
+        free(A->ia);  free(A->ja); free(A->sa); free(A);
+        
         ptr = mxGetPr(plhs[1]);
         memcpy(ptr, sys->r, ncf_tot * sizeof *ptr);
 
@@ -185,5 +237,6 @@ mexFunction(int nlhs,       mxArray *plhs[],
 
         hybsys_free(sys);
         deallocate_aux_arrays(nconn, src, gflux);
+        mxFree(conn);
     }
 }
diff --git a/test_mex_schur_comp_symm.m b/test_mex_schur_comp_symm.m
index ef44ad79..e25941bd 100644
--- a/test_mex_schur_comp_symm.m
+++ b/test_mex_schur_comp_symm.m
@@ -5,7 +5,7 @@ BI = mex_ip_simple(G, rock);
 
 nconn = diff(G.cells.facePos);
 conn  = G.cells.faces(:, 1);
-[S, r, F, L] = mex_schur_comp_symm(BI, nconn);
+[S, r, F, L] = mex_schur_comp_symm(BI, nconn, conn);
 
 [i, j] = blockDiagIndex(nconn, nconn);