Commits

Anonymous committed 81c73f6

fixed a bug which escaped me for the last check in because I didnt check with cutoff=64

Comments (0)

Files changed (5)

 
 size_t _mzd_pluq(mzd_t *A, mzp_t * P, mzp_t * Q, const int cutoff) {
   size_t r  = _mzd_lqup(A,P,Q,cutoff);
-  mzd_apply_p_right_tri(A, Q, r);
+  mzd_apply_p_right_tri(A, Q);
   return r;
 }
 
 
     /* Compressing L */
 
+    mzp_t *shift = mzp_init(A->ncols);
     if (r1 < n1){
       for (i=r1,j=n1;i<r1+r2;i++,j++){
-	mzd_col_swap_in_rows (A, i, j, i, A->nrows);
+	mzd_col_swap_in_rows(A, i, j, i, r1+r2);
+        shift->values[i] = j;
       }
     }
+    mzd_apply_p_right_trans_even_capped(A, shift, r1+r2, 0);
+    mzp_free(shift);
 
     mzp_free_window(Q2);
     mzp_free_window(P2);

src/permutation.c

 /**
  * Implements both apply_p_right and apply_p_right_trans.
  */
-void _mzd_apply_p_right_even(mzd_t *A, mzp_t *P, size_t start_row, int notrans) {
+void _mzd_apply_p_right_even(mzd_t *A, mzp_t *P, size_t start_row, size_t start_col, int notrans) {
   assert(A->offset = 0);
   const size_t length = MIN(P->length,A->ncols);
   const size_t width = A->width;
     permutation[i] = i;
 
   if (!notrans) {
-    for(size_t i=0; i<length; i++) {
+    for(size_t i=start_col; i<length; i++) {
       size_t t = permutation[i];
       permutation[i] = permutation[P->values[i]];
       permutation[P->values[i]] = t;
     }
   } else {
-    for(size_t i=0; i<length; i++) {
+    for(size_t i=start_col; i<length; i++) {
       size_t t = permutation[length-i-1];
       permutation[length-i-1] = permutation[P->values[length-i-1]];
       permutation[P->values[length-i-1]] = t;
     _mzd_apply_p_right_trans(A,P);
     return;
   }
-  _mzd_apply_p_right_even(A, P, 0, 0); 
+  _mzd_apply_p_right_even(A, P, 0, 0, 0); 
 }
 
 void mzd_apply_p_right(mzd_t *A, mzp_t *P) {
     _mzd_apply_p_right(A,P);
     return;
   }
-  _mzd_apply_p_right_even(A, P, 0, 1); 
+  _mzd_apply_p_right_even(A, P, 0, 0, 1); 
 }
 
-void mzd_apply_p_right_trans_even_capped(mzd_t *A, mzp_t *P, size_t start_row) {
+void mzd_apply_p_right_trans_even_capped(mzd_t *A, mzp_t *P, size_t start_row, size_t start_col) {
   assert(A->offset == 0);
   if(!A->nrows)
     return;
-  _mzd_apply_p_right_even(A, P, start_row, 0); 
+  _mzd_apply_p_right_even(A, P, start_row, start_col, 0); 
 }
 
-void mzd_apply_p_right_even_capped(mzd_t *A, mzp_t *P, size_t start_row) {
+void mzd_apply_p_right_even_capped(mzd_t *A, mzp_t *P, size_t start_row, size_t start_col) {
   assert(A->offset == 0);
   if(!A->nrows)
     return;
-  _mzd_apply_p_right_even(A, P, start_row, 1); 
+  _mzd_apply_p_right_even(A, P, start_row, start_col, 1); 
 }
 
 
 }
 
 #if 0
-void  mzd_apply_p_right_tri (mzd_t * A, mzp_t * P){
+void  mzd_apply_p_right_tri (mzd_t * A, mzp_t * P, size_t rank){
   assert(P->length==A->ncols);
   for (size_t i =0 ; i<P->length; ++i) {
     assert(P->values[i] >= i);
   }
 }
 #else
-void  mzd_apply_p_right_tri(mzd_t *A, mzp_t *P, size_t rank) {
+void mzd_apply_p_right_tri(mzd_t *A, mzp_t *P) {
   assert(P->length==A->ncols);
   const size_t step_size = MAX((CPU_L1_CACHE>>2)/A->width,1);
 
   for(size_t r=0; r<A->nrows; r+=step_size) {
     const size_t row_bound = MIN(r+step_size, A->nrows);
-    for (size_t i =0 ; i<rank; ++i) {
+    for (size_t i =0 ; i<A->ncols; ++i) {
       assert(P->values[i] >= i);
       if (P->values[i] > i) {
         mzd_col_swap_in_rows(A, i, P->values[i], r, MIN(row_bound,i));

src/permutation.h

  * \wordoffset
  */
 
-void mzd_apply_p_right_even_capped(mzd_t *A, mzp_t *P, size_t start_row);
+void mzd_apply_p_right_even_capped(mzd_t *A, mzp_t *P, size_t start_row, size_t start_col);
 
 /**
  * Apply the permutation P^T to A from the right starting at start_row.
  * \wordoffset
  */
 
-void mzd_apply_p_right_trans_even_capped(mzd_t *A, mzp_t *P, size_t start_row);
+void mzd_apply_p_right_trans_even_capped(mzd_t *A, mzp_t *P, size_t start_row, size_t start_col);
 
 /**
  * Apply the mzp_t P to A from the right but transpose P before.
  *
  * \param A Matrix.
  * \param P Permutation.
- * \param rank a column where to stop.
  */
-void  mzd_apply_p_right_tri(mzd_t * A, mzp_t * Q, size_t rank);
+void  mzd_apply_p_right_tri(mzd_t * A, mzp_t * Q);
 
 /* /\** */
 /*  * Rotate zero columns to the end. */
     }
   }
   mzp_t *Qbar = mzp_init_window(Q,0,curr_row);
-  mzd_apply_p_right_trans_even_capped(A, Qbar, curr_row);
+  mzd_apply_p_right_trans_even_capped(A, Qbar, curr_row, 0);
   mzp_free_window(Qbar);
 /*   printf("\nAfter"); */
 /*   mzd_print(A); */

testsuite/test_lqup.c

     if(status)
       break;
   }
-  
   if (status)
     printf(" ... FAILED\n");
   else
 
   if (status) {
     printf("\n");
-    mzd_print(Acopy);
     printf(" ... FAILED\n");
   }  else
     printf (" ... passed\n");