[Commit-gnuradio] [gnuradio] 14/39: fec: LDPC: Reducing complexity of encoder by adding back solve.

[git]

This is an automated email from the git hooks/post-receive script.

jcorgan pushed a commit to branch master
in repository gnuradio.

commit aa13153c895853c74fe2540271a72542ec34307d
Author: tracierenea <address@hidden>
Date:   Mon Jul 28 19:52:57 2014 -0500

    fec: LDPC: Reducing complexity of encoder by adding back solve.
    
    Exploit that T submatrix is upper triangular by using back
    substitution rather than standard matrix multiplication when finding
    parity part of codeword.
---
 gr-fec/include/gnuradio/fec/ldpc_R_U_mtrx.h |  5 +--
 gr-fec/lib/ldpc_R_U_encoder_impl.cc         | 62 ++++++++++++++++++++++++++---
 gr-fec/lib/ldpc_R_U_encoder_impl.h          |  2 +
 gr-fec/lib/ldpc_R_U_mtrx.cc                 | 36 +++++++++--------
 4 files changed, 80 insertions(+), 25 deletions(-)

diff --git a/gr-fec/include/gnuradio/fec/ldpc_R_U_mtrx.h 
b/gr-fec/include/gnuradio/fec/ldpc_R_U_mtrx.h
index 216cc64..7b36589 100644
--- a/gr-fec/include/gnuradio/fec/ldpc_R_U_mtrx.h
+++ b/gr-fec/include/gnuradio/fec/ldpc_R_U_mtrx.h
@@ -53,8 +53,7 @@ namespace gr {
         gsl_matrix_view d_B_view;
         gsl_matrix_view d_D_view;
         gsl_matrix_view d_E_view;
-        gsl_matrix d_inv_phi;
-        gsl_matrix *d_T_inverse_ptr;
+        gsl_matrix_view d_T_view;
         gsl_matrix *d_phi_inverse_ptr;
 
         // Read the matrix from a file in alist format
@@ -76,7 +75,7 @@ namespace gr {
         const gsl_matrix *B();
         const gsl_matrix *D();
         const gsl_matrix *E();
-        const gsl_matrix *T_inverse();
+        const gsl_matrix *T();
         const gsl_matrix *phi_inverse();
         // Subtract matrices using mod2 operation
         gsl_matrix *add_matrices_mod2(const gsl_matrix *,
diff --git a/gr-fec/lib/ldpc_R_U_encoder_impl.cc 
b/gr-fec/lib/ldpc_R_U_encoder_impl.cc
index 37ec192..61a60b5 100644
--- a/gr-fec/lib/ldpc_R_U_encoder_impl.cc
+++ b/gr-fec/lib/ldpc_R_U_encoder_impl.cc
@@ -79,6 +79,56 @@ namespace gr {
         return (d_H->n())/static_cast<double>(d_frame_size);
       }
 
+      gsl_matrix*
+      ldpc_R_U_encoder_impl::back_solve_mod2(const gsl_matrix *U,
+                                             const gsl_matrix *y)
+      {
+        // Exploit the fact that the matrix T is upper triangular and
+        // sparse. In the steps to find p1 and p2, back solve rather 
+        // than do matrix multiplication to reduce number of
+        // operations required. 
+
+        // Form is Ux = y where U is upper triangular and y is column
+        // vector. Solve for x.
+
+        // Allocate memory for the result
+        int num_rows   = (*U).size1;
+        int num_cols_U = (*U).size2;
+        gsl_matrix *x  = gsl_matrix_alloc(num_rows,1);
+
+        // Back solve
+        for (int i = num_rows-1; i >= 0; i--) {
+          // x[i] = y[i]
+          gsl_matrix_set(x, i, 0, gsl_matrix_get(y, i, 0));
+
+          int j;
+          for (j = i+1; j < num_cols_U; j++) {
+            int U_i_j = gsl_matrix_get(U, i, j);
+            int x_i   = gsl_matrix_get(x, i, 0);
+            int x_j   = gsl_matrix_get(x, j, 0);
+            int temp1 = (U_i_j * x_j) % 2;
+            int temp2 = (x_i + temp1) % 2;
+            gsl_matrix_set(x, i, 0, temp2);
+          }
+          // Perform x[i] /= U[i,i], GF(2) operations
+          int U_i_i = gsl_matrix_get(U, i, i);
+          int x_i   = gsl_matrix_get(x, i, 0);
+          if      (x_i==0 && U_i_i==1)
+            gsl_matrix_set(x, i, 0, 0);
+          else if (x_i==0 && U_i_i==0)
+            gsl_matrix_set(x, i, 0, 0);
+          else if (x_i==1 && U_i_i==1)
+            gsl_matrix_set(x, i, 0, 1);
+          else if (x_i==1 && U_i_i==0)
+            std::cout << "Error in "
+                      << " ldpc_R_U_encoder_impl::back_solve_mod2,"
+                      << " division not defined.\n";  
+          else
+            std::cout << "Error in ldpc_R_U_encoder_impl::back_solve_mod2\n";  
+        }
+        return x;
+      }
+
       void
       ldpc_R_U_encoder_impl::generic_work(void *inbuffer,
                                           void *outbuffer)
@@ -92,19 +142,21 @@ namespace gr {
           gsl_matrix_set(s, index, 0, value);
         }
 
-        // Solve for p2 and p1 (parity parts)
+        // Solve for p2 (parity part). By using back substitution,
+        // the overall complexity of determining p2 is O(n + g^2).
         gsl_matrix *temp1 = d_H->mult_matrices_mod2(d_H->B(), s);
-        gsl_matrix *temp2 = d_H->mult_matrices_mod2(
-                                         d_H->T_inverse(), temp1);
+        gsl_matrix *temp2 = back_solve_mod2(d_H->T(),temp1);
         gsl_matrix *temp3 = d_H->mult_matrices_mod2(d_H->E(), temp2);
         gsl_matrix *temp4 = d_H->mult_matrices_mod2(d_H->D(), s);
         gsl_matrix *temp5 = d_H->add_matrices_mod2(temp4, temp3);
         gsl_matrix *p2    = d_H->mult_matrices_mod2(
                                          d_H->phi_inverse(), temp5);
+
+        // Solve for p1 (parity part). By using back substitution,
+        // the overall complexity of determining p1 is O(n).
         gsl_matrix *temp6 = d_H->mult_matrices_mod2(d_H->A(), p2);
         gsl_matrix *temp7 = d_H->add_matrices_mod2(temp6, temp1);
-        gsl_matrix *p1    = d_H->mult_matrices_mod2(
-                                         d_H->T_inverse(), temp7);
+        gsl_matrix *p1 =  back_solve_mod2(d_H->T(),temp7);
 
         // Populate the codeword to be output
         unsigned int p1_length = (*p1).size1;
diff --git a/gr-fec/lib/ldpc_R_U_encoder_impl.h 
b/gr-fec/lib/ldpc_R_U_encoder_impl.h
index 6417380..a70f101 100644
--- a/gr-fec/lib/ldpc_R_U_encoder_impl.h
+++ b/gr-fec/lib/ldpc_R_U_encoder_impl.h
@@ -35,6 +35,8 @@ namespace gr {
         void generic_work(void *inbuffer, void *outbuffer);
         int get_output_size();
         int get_input_size();
+        gsl_matrix *back_solve_mod2(const gsl_matrix *,
+                                    const gsl_matrix *);
 
         // Number of bits in the information word
         unsigned int d_frame_size;
diff --git a/gr-fec/lib/ldpc_R_U_mtrx.cc b/gr-fec/lib/ldpc_R_U_mtrx.cc
index a1e9116..30a5f3b 100644
--- a/gr-fec/lib/ldpc_R_U_mtrx.cc
+++ b/gr-fec/lib/ldpc_R_U_mtrx.cc
@@ -96,10 +96,10 @@ namespace gr {
       }
 
       const gsl_matrix*
-      ldpc_R_U_mtrx::T_inverse()
+      ldpc_R_U_mtrx::T()
       {
-        const gsl_matrix *T_inverse_ptr = d_T_inverse_ptr;
-        return T_inverse_ptr; 
+        const gsl_matrix *T_ptr = &d_T_view.matrix;
+        return T_ptr; 
       }
 
       const gsl_matrix*
@@ -113,9 +113,9 @@ namespace gr {
       ldpc_R_U_mtrx::read_matrix_from_file(const std::string filename)
       {
         /* This function reads in an alist file and creates the
-          corresponding parity check matrix. The format of alist
-          files is described at:
-          http://www.inference.phy.cam.ac.uk/mackay/codes/alist.html
+           corresponding parity check matrix. The format of alist
+           files is described at:
+           http://www.inference.phy.cam.ac.uk/mackay/codes/alist.html
         */
         std::ifstream inputFile;
 
@@ -185,13 +185,11 @@ namespace gr {
         unsigned int t = d_num_rows - d_gap; 
 
         // T submatrix
-        gsl_matrix_view T_view = gsl_matrix_submatrix(d_H_ptr,
-                                                      0,
-                                                      0,
-                                                      t,
-                                                      t);
+        d_T_view = gsl_matrix_submatrix(d_H_ptr, 0, 0, t, t);
+
+        gsl_matrix *d_T_inverse_ptr;
         try {
-          d_T_inverse_ptr = calc_inverse_mod2(&T_view.matrix);
+          d_T_inverse_ptr = calc_inverse_mod2(&d_T_view.matrix);
         }
         catch (char const *exceptionString) {
           std::cout << "Error in set_parameters_for_encoding while "
@@ -229,7 +227,6 @@ namespace gr {
             gsl_matrix *inverse_phi = calc_inverse_mod2(phi);
               
             // At this point, an inverse was found. 
-            d_inv_phi = *inverse_phi;
             d_phi_inverse_ptr = inverse_phi;
 
           }
@@ -250,13 +247,19 @@ namespace gr {
         // D submatrix
         d_D_view = gsl_matrix_submatrix(d_H_ptr, t, t + d_gap, d_gap,
                                         d_n-d_gap-t);
+
+        // Free memory
+        gsl_matrix_free(temp1);
+        gsl_matrix_free(temp2);
+        gsl_matrix_free(phi);
+        gsl_matrix_free(d_T_inverse_ptr);
+
       }
 
       gsl_matrix*
       ldpc_R_U_mtrx::add_matrices_mod2(const gsl_matrix *matrix1, const 
gsl_matrix *matrix2)
       {
         // This function returns ((matrix1 + matrix2) % 2). 
-        // (same thing as ((matrix1 - matrix2) % 2)
 
         // Verify that matrix sizes are appropriate
         unsigned int matrix1_rows = (*matrix1).size1;
@@ -316,14 +319,14 @@ namespace gr {
           exit(1);
         }
 
-        // Allocate memory for the result
+        // Allocate memory for the result.
         gsl_matrix *result = gsl_matrix_alloc(a,d);
 
         // Perform matrix multiplication. This is not mod 2.
         gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, matrix1,
                         matrix2, 0.0, result);
 
-        // Take care of mod 2 manually
+        // Take care of mod 2 manually.
         unsigned int row_index, col_index;
         unsigned int rows = (*result).size1, 
         cols = (*result).size2;
@@ -430,7 +433,6 @@ namespace gr {
       {
         // Call the gsl_matrix_free function to free memory.
         gsl_matrix_free (d_H_ptr);
-        gsl_matrix_free (d_T_inverse_ptr);
         gsl_matrix_free (d_phi_inverse_ptr);
       }
     } /* namespace code */