Commits

Anonymous committed 082c241

implemented mzd_kernel_left_pluq to compute kernels via PLUQ

  • Participants
  • Parent commits 8260bc1

Comments (0)

Files changed (3)

src/packedmatrix.c

 
 #define SAFECHAR (int)(RADIX+RADIX/3)
 
+void mzd_copy_row_weird_to_even(mzd_t* B, size_t i, const mzd_t* A, size_t j);
+
 mzd_t *mzd_init(size_t r, size_t c) {
   mzd_t *A;
   size_t i,j;
       if (N->nrows < P->nrows || N->ncols < P->ncols)
 	m4ri_die("mzd_copy: Target matrix is too small.");
     }
-    for(size_t i=0; i<P->nrows; i++) {
-      mzd_copy_row(N, i, P, i);
+    if(N->offset == P->offset) {
+      for(size_t i=0; i<P->nrows; i++) {
+        mzd_copy_row(N, i, P, i);
+      }
+    } else if(N->offset == 0) {
+      for(size_t i=0; i<P->nrows; i++) {
+        mzd_copy_row_weird_to_even(N, i, P, i);
+      }
+    } else {
+      m4ri_die("mzd_copy: completely unaligned copy not implemented yet.");
     }
   }
 /*     size_t i, j, p_truerow, n_truerow; */
   return !status;
 }
 
+void mzd_copy_row_weird_to_even(mzd_t* B, size_t i, const mzd_t* A, size_t j) {
+  assert(B->offset == 0);
+  assert(B->ncols >= A->ncols);
+
+  word *b = B->rows[j];
+  size_t c;
+
+  size_t rest = A->ncols%RADIX;
+
+  for(c = 0; c+RADIX <= A->ncols; c+=RADIX) {
+    b[c/RADIX] = mzd_read_bits(A, i, c, RADIX);
+  }
+  if (rest) {
+    const word temp = mzd_read_bits(A, i, c, rest);
+    b[c/RADIX] &= LEFT_BITMASK(RADIX-rest);
+    b[c/RADIX] |= temp<<(RADIX-rest);
+  }
+}
+
 void mzd_copy_row(mzd_t* B, size_t i, const mzd_t* A, size_t j) {
   assert(B->offset == A->offset);
   assert(B->ncols >= A->ncols);
   /* P L U Q B5 = B1 */
 }
 
-
 void _mzd_solve_left (mzd_t *A, mzd_t *B, const int cutoff, const int inconsistency_check) {
   /**
    *  B is modified in place 
   mzp_free(Q);
 }
 
+mzd_t *mzd_kernel_left_pluq(mzd_t *A) {
+  mzp_t *P = mzp_init(A->nrows);
+  mzp_t *Q = mzp_init(A->ncols);
+
+  size_t r = mzd_pluq(A, P, Q, 0);
+
+  if (r == A->ncols) {
+    mzp_free(P);
+    mzp_free(Q);
+    return NULL;
+  }
+
+  mzd_t *U = mzd_init_window(A, 0, 0, r, r);
+  mzd_t *B = mzd_init_window(A, 0, r, r, A->ncols);
+
+  mzd_trsm_upper_left(U, B, 0);
+
+  mzd_t *R = mzd_init(A->ncols, A->ncols - r);
+  mzd_t *RU = mzd_init_window(R, 0, 0, r, R->ncols);
+  mzd_copy(RU, B);
+  for(size_t i=0;i<R->ncols;i++) {
+    mzd_write_bit(R, r+i, i, 1);
+  }
+  mzd_apply_p_left_trans(R, Q);
+  mzp_free(P);
+  mzp_free(Q);
+  mzd_free_window(RU);
+  mzd_free_window(U);
+  mzd_free_window(B);
+  return R;
+}
  */
 void _mzd_solve_left(mzd_t *A, mzd_t *B, const int cutoff, const int inconsistency_check);
 
+
+/**
+ * \brief Solve X for A X = 0.
+ *
+ * If r is the rank of the nr x nc matrix A, return the nc x (nc-r)
+ * matrix X such that A*X == 0 and that the columns of X are linearly
+ * independent.
+ *
+ * \param A Matrix.
+ *
+ * \wordoffset
+ *
+ * \sa mzd_pluq()
+ *
+ * \return X
+ */
+
+mzd_t *mzd_kernel_left_pluq(mzd_t *A);
+
 #endif // SOLVE_H