Commits

Anonymous committed 048a94f

fixed some minor bug in TRSM

Comments (0)

Files changed (6)

 
     size_t i, j;
     size_t n1 = (((ncols - 1) / RADIX + 1) >> 1) * RADIX;
-    packedmatrix *A0  = mzd_init_window (A,  0,  0, nrows,    n1);
-    packedmatrix *A1  = mzd_init_window (A,  0, n1, nrows, ncols);
+    packedmatrix *A0  = mzd_init_window(A,  0,  0, nrows,    n1);
+    packedmatrix *A1  = mzd_init_window(A,  0, n1, nrows, ncols);
 
     size_t r1, r2;
     /* First recursive call */
-    r1 = _mzd_pluq (A0, P, Q, cutoff);
-
-    /*
-     *           r1           n1
+    r1 = _mzd_pluq(A0, P, Q, cutoff);
+    
+    /*           r1           n1
      *   ------------------------------------------
      *   | A00    |           | A01               |
      *   |        |           |                   |
      *   | A01    |           | A11               |
      *   |        |           |                   |
      *   ------------------------------------------
-     *
      */
 
-    packedmatrix *A00  = mzd_init_window (A,   0, 0, r1, r1);
-    packedmatrix *A10  = mzd_init_window (A,  r1, 0, nrows, r1);
-    packedmatrix *A01  = mzd_init_window (A,  0, n1, r1, ncols);
-    packedmatrix *A11  = mzd_init_window (A,  r1, n1, nrows, ncols);
+    packedmatrix *A00 = mzd_init_window(A,  0,  0, r1, r1);
+    packedmatrix *A10 = mzd_init_window(A, r1,  0, nrows, r1);
+    packedmatrix *A01 = mzd_init_window(A,  0, n1, r1, ncols);
+    packedmatrix *A11 = mzd_init_window(A, r1, n1, nrows, ncols);
 
     if (r1) {
       /* Computation of the Schur complement */
     }
 
     /* Second recursive call */
-    permutation * P2 = mzp_init_window(P, r1, nrows);
-    permutation * Q2 = mzp_init_window(Q, n1, ncols);
+    permutation *P2 = mzp_init_window(P, r1, nrows);
+    permutation *Q2 = mzp_init_window(Q, n1, ncols);
 
     r2 = _mzd_pluq(A11, P2, Q2, cutoff);
 
  * Crossover point for PLUQ factorization.
  */
 
-#define PLUQ_CUTOFF 512
+#define PLUQ_CUTOFF 64
 
 /**
  * \brief PLUQ matrix decomposition.
           Q->values[curr_pos] = j;
           mzd_row_swap(A, curr_pos, i);
           mzd_col_swap(A, curr_pos, j);
-          for(size_t l = i+1; l<A->nrows; l++) {
+          for(size_t l = i+1; l<A->nrows; l++)
             if(mzd_read_bit(A, l, curr_pos))
               mzd_row_add_offset(A, l, curr_pos, curr_pos + 1);
-          }
+          
           curr_pos++;
         } else {
           break;
  */
 void _mzd_trsm_upper_left_weird (packedmatrix *U, packedmatrix *B, const int cutoff);
 
-/*
- * Variant where U and B start at an odd bit position
- * Assumes that U->ncols < 64
- */
 void _mzd_trsm_upper_left_even(packedmatrix *U, packedmatrix *B, const int cutoff);
 
 void mzd_trsm_upper_left(packedmatrix *U, packedmatrix *B, const int cutoff) {
   size_t nb = B->ncols;
   size_t Boffset = B->offset;
   size_t nbrest = (nb + Boffset) % RADIX;
-  if (nb + Boffset >= RADIX) {
+  if (nb + Boffset > RADIX) {
 
     // Large B
     word mask_begin = RIGHT_BITMASK(RADIX-B->offset);
   if (mb <= RADIX){
     /* base case */
 
-    if (nb + B->offset >= RADIX) {
+    if (nb + B->offset > RADIX) {
       // B is large
       word mask_begin = RIGHT_BITMASK(RADIX-B->offset);
       word mask_end = LEFT_BITMASK(nbrest);

testsuite/test_lqup.c

   mzd_randomize (L);
 
   size_t i,j;
-  for (i=0; i<m; ++i){
+  for (i=0; i<m && i<n; ++i){
     mzd_write_bit(U,i,i, 1);
     for (j=0; j<i;++j)
       mzd_write_bit(U,i,j, 0);

testsuite/test_trsm.c

 #include "m4ri/m4ri.h"
 
 
-int test_trsm_upper_right (int m, int n, int offset){
+int test_trsm_upper_right (int m, int n, int offset, const char* description){
+  printf("upper_right: %s  m: %4d n: %4d offset: %4d ... ",description, m, n, offset);
+
   packedmatrix* Ubase = mzd_init (2048,2048);
   packedmatrix* Bbase = mzd_init (2048,2048);
   mzd_randomize (Ubase);
   return status;
 }
 
-int test_trsm_lower_right (int m, int n, int offset){
+int test_trsm_lower_right (int m, int n, int offset, const char *description){
+  printf("lower_right: %s  m: %4d n: %4d offset: %4d ... ",description, m, n, offset);
   packedmatrix* Lbase = mzd_init (2048,2048);
   packedmatrix* Bbase = mzd_init (2048,2048);
   mzd_randomize (Lbase);
 
 
 
-int test_trsm_upper_left (int m, int n, int offsetU, int offsetB){
+int test_trsm_upper_left (int m, int n, int offsetU, int offsetB, const char *description) {
+  printf("upper_left: %s  m: %4d n: %4d offset: %4d ... ",description, m, n, offsetU);
   packedmatrix* Ubase = mzd_init (2048,2048);
   packedmatrix* Bbase = mzd_init (2048,2048);
   mzd_randomize (Ubase);
 int main(int argc, char **argv) {
   int status = 0;
 
-  printf("UpperRight: small, even placed ... ");
-  status += test_trsm_upper_right(57, 10, 0);
-  printf("UpperRight: large, even placed ... ");
-  status += test_trsm_upper_right(57, 150, 0);
-  printf("UpperRight: small, odd placed  ... ");
-  status += test_trsm_upper_right(57, 3, 4);
-  printf("UpperRight: medium, odd placed ... ");
-  status += test_trsm_upper_right(57, 4, 62);
-  printf("UpperRight: large, odd placed  ... ");
-  status += test_trsm_upper_right(57, 80, 60);
-  printf("UpperRight: larger, odd placed ... ");
-  status += test_trsm_upper_right(1577, 1802, 189);
+  status += test_trsm_upper_right(  57,   10,   0, "small, even placed");
+  status += test_trsm_upper_right(  57,  150,   0, "large, even placed");
+  status += test_trsm_upper_right(  57,    3,   4, " small, odd placed");
+  status += test_trsm_upper_right(  57,    4,  62, "medium, odd placed");
+  status += test_trsm_upper_right(  57,   80,  60, " large, odd placed");
+  status += test_trsm_upper_right(1577, 1802, 189, "larger, odd placed");
 
   printf("\n");
 
-  printf("LowerRight: small, even placed ... ");
-  status += test_trsm_lower_right(57, 10, 0);
-  printf("LowerRight: large, even placed ... ");
-  status += test_trsm_lower_right(57, 150, 0);
-  printf("LowerRight: small, odd placed  ... ");
-  status += test_trsm_lower_right(57, 3, 4);
-  printf("LowerRight: medium, odd placed ... ");
-  status += test_trsm_lower_right(57, 4, 62);
-  printf("LowerRight: large, odd placed  ... ");
-  status += test_trsm_lower_right(57, 80, 60);
-  printf("LowerRight: larger, odd placed ... ");
-  status += test_trsm_lower_right(1577, 1802, 189);
+  status += test_trsm_lower_right(  57,   10,  0,"small, even placed");
+  status += test_trsm_lower_right(  57,  150,  0,"large, even placed");
+  status += test_trsm_lower_right(  57,    3,  4," small, odd placed");
+  status += test_trsm_lower_right(  57,    4, 62,"medium, odd placed");
+  status += test_trsm_lower_right(  57,   80, 60," large, odd placed");
+  status += test_trsm_lower_right(1577, 1802,189,"larger, odd placed");
 
   printf("\n");
 
 
   printf("\n");
 
-  printf("UpperLeft: small U even, small B even ... ");
-  status += test_trsm_upper_left (10, 20, 0, 0);
-  printf("UpperLeft: small U even, large B even ... ");
-  status += test_trsm_upper_left (10, 80, 0, 0);
-  printf("UpperLeft: small U even, small B odd  ... ");
-  status += test_trsm_upper_left (10, 20, 0, 15);
-  printf("UpperLeft: small U even, large B odd  ... ");
-  status += test_trsm_upper_left (10, 80, 0, 15);
-  printf("UpperLeft: small U odd, small B even  ... ");
-  status += test_trsm_upper_left (10, 20, 15, 0);
-  printf("UpperLeft: small U odd, large B even  ... ");
-  status += test_trsm_upper_left (10, 80, 15, 0);
-  printf("UpperLeft: small U odd, small B odd   ... ");
-  status += test_trsm_upper_left (10, 20, 15, 20);
-  printf("UpperLeft: small U odd, large B odd   ... ");
-  status += test_trsm_upper_left (10, 80, 15, 20);
-  printf("UpperLeft: large U even, small B even ... ");
-  status += test_trsm_upper_left (70, 20, 0, 0);
-  printf("UpperLeft: large U even, large B even ... ");
-  status += test_trsm_upper_left (70, 80, 0, 0);
-  printf("UpperLeft: large U even, small B odd  ... ");
-  status += test_trsm_upper_left (70, 10, 0, 15);
-  printf("UpperLeft: large U even, large B odd  ... ");
-  status += test_trsm_upper_left (70, 80, 0, 15);
-  printf("UpperLeft: large U odd, small B even  ... ");
-  status += test_trsm_upper_left (70, 20, 15, 0);
-  printf("UpperLeft: large U odd, large B even  ... ");
-  status += test_trsm_upper_left (70, 80, 15, 0);
-  printf("UpperLeft: large U odd, small B odd   ... ");
-  status += test_trsm_upper_left (70, 20, 15, 20);
-  printf("UpperLeft: large U odd, large B odd   ... ");
-  status += test_trsm_upper_left (70, 80, 15, 20);
-  printf("UpperLeft: larger U odd, larger B odd ... ");
-  status += test_trsm_upper_left (770, 1600, 75, 89);
-  printf("UpperLeft: larger U odd, larger B odd ... ");
-  status += test_trsm_upper_left (1764, 1345, 198, 123);
+  status += test_trsm_upper_left(  10,  20,  0,  0,"small U even, small B even");
+  status += test_trsm_upper_left(  10,  80,  0,  0,"small U even, large B even");
+  status += test_trsm_upper_left(  10,  20,  0, 15," small U even, small B odd");
+  status += test_trsm_upper_left(  10,  80,  0, 15," small U even, large B odd");
+  status += test_trsm_upper_left(  10,  20, 15,  0," small U odd, small B even");
+  status += test_trsm_upper_left(  10,  80, 15,  0," small U odd, large B even");
+  status += test_trsm_upper_left(  10,  20, 15, 20,"  small U odd, small B odd");
+  status += test_trsm_upper_left(  10,  80, 15, 20,"  small U odd, large B odd");
+  status += test_trsm_upper_left(  70,  20,  0,  0,"large U even, small B even");
+  status += test_trsm_upper_left(  63,   1,  0,  0,"                          ");
+  status += test_trsm_upper_left(  70,  80,  0,  0,"large U even, large B even");
+  status += test_trsm_upper_left(  70,  10,  0, 15," large U even, small B odd");
+  status += test_trsm_upper_left(  70,  80,  0, 15," large U even, large B odd");
+  status += test_trsm_upper_left(  70,  20, 15,  0," large U odd, small B even");
+  status += test_trsm_upper_left(  70,  80, 15,  0," large U odd, large B even");
+  status += test_trsm_upper_left(  70,  20, 15, 20,"  large U odd, small B odd");
+  status += test_trsm_upper_left(  70,  80, 15, 20,"  large U odd, large B odd");
+  status += test_trsm_upper_left( 770,1600, 75, 89,"larger U odd, larger B odd");
+  status += test_trsm_upper_left(1764,1345,198,123,"larger U odd, larger B odd");
 
   if (!status) {
     printf("All tests passed.\n");