diff --git a/list.txt b/list.txt
new file mode 100644
index 0000000000..6f85b69219
--- /dev/null
+++ b/list.txt
@@ -0,0 +1,2318 @@
+./goblint --enable warn.debug --enable dbg.regression  --set ana.activated[+] affeq --set sem.int.signed_overflow assume_none  tests/regression/77-lin2vareq/20-function_call2.c
+Before assign_var_parallel m:
+⊤
+After assign_var_parallel multi_t:
+⊤
+Matrix in Domain m:
+
+Switched Matrix in Domain switched_m:
+
+Before normalizing we have m:
+After normalizing we have m:
+Before add_empty_columns m:
+indices: 0,0,0,
+Occ_cols is:
+(0,3)
+After add_empty_columns m:
+
+Before assign_var_parallel m:
+⊤
+After assign_var_parallel multi_t:
+⊤
+After affineEq_vec m:
+[1 0 0 0 ]
+
+Before Matrix.rref_vec x:
+ b:
+[-1 1 0 0 ]
+
+Before rref_vec we have m:
+v: [-1 1 0 0 ]
+
+Before normalizing we have m:
+[-1 1 0 0 ]
+After normalizing we have m:
+[1 -1 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 -1 0 0 ]
+
+After affineEq_vec m:
+[0 1 0 0 ]
+
+Before reduce_col 2 of m:
+[1 -1 0 0 ]
+
+After reduce_col 2 of m:
+[1 -1 0 0 ]
+
+Before Matrix.rref_vec x:
+[1 -1 0 0 ]
+ b:
+[0 -1 1 0 ]
+
+Before rref_vec we have m:
+[1 -1 0 0 ]
+v: [0 -1 1 0 ]
+
+Before normalizing we have m:
+[1 -1 0 0 ]
+[0 -1 1 0 ]
+After normalizing we have m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+
+Before add_empty_columns m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+indices: 0,1,2,
+Occ_cols is:
+(0,1)(1,1)(2,1)
+After add_empty_columns m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 1 0 0 0 0 0 ]
+
+Before reduce_col 0 of m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After reduce_col 0 of m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+ b:
+[1 -1 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+v: [1 -1 0 0 0 0 0 ]
+
+Before normalizing we have m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[1 -1 0 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 1 0 0 0 ]
+
+Before reduce_col 2 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After reduce_col 2 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+ b:
+[0 0 1 -1 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+v: [0 0 1 -1 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 1 -1 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 1 0 ]
+
+Before reduce_col 4 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After reduce_col 4 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+ b:
+[0 0 0 0 1 -1 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+v: [0 0 0 0 1 -1 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+
+Before dim_remove m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+Before reduce_col 1 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 1 -1 0 ]
+
+Before reduce_col 5 of m:
+[1 0 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 1 -1 0 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 0 0 0 ]
+
+After dim_remove m':
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before add_empty_columns m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+indices: 0,1,2,
+Occ_cols is:
+(0,1)(1,1)(2,1)
+After add_empty_columns m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+Before assign_var_parallel m:
+[|x#322#arg-z#324#arg=0; y#323#arg-z#324#arg=0|]
+Before add_empty_columns m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+indices: 0,0,0,
+Occ_cols is:
+(0,3)
+After add_empty_columns m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 1 0 0 0 0 0 ]
+
+Before reduce_col 0 of m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After reduce_col 0 of m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+ b:
+[1 0 0 0 -1 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+v: [1 0 0 0 -1 0 0 0 0 0 ]
+
+Before normalizing we have m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+[1 0 0 0 -1 0 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 0 1 0 0 0 ]
+
+Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+ b:
+[0 1 0 0 0 0 -1 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+v: [0 1 0 0 0 0 -1 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+[0 1 0 0 0 0 -1 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 0 0 0 1 0 ]
+
+Before reduce_col 2 of m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After reduce_col 2 of m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+ b:
+[0 0 1 0 0 0 0 0 -1 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+v: [0 0 1 0 0 0 0 0 -1 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After assign_var_parallel multi_t:
+[|0'-z#324#arg=0; 1'-z#324#arg=0; 2'-z#324#arg=0; x#322#arg-z#324#arg=0; y#323#arg-z#324#arg=0|]
+Matrix in Domain m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before set_col m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Switched Matrix in Domain switched_m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before dim_remove m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After dim_remove m':
+[1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [1 0 -1 0 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[1 0 -1 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [-1 0 1 0 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[-1 0 1 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [0 0 -1 0 1 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 -1 0 1 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [0 0 1 0 -1 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 1 0 -1 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [1 0 0 0 -1 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[1 0 0 0 -1 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [-1 0 0 0 1 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[-1 0 0 0 1 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before add_empty_columns m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+indices: 0,
+Occ_cols is:
+(0,1)
+After add_empty_columns m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 0 0 8 ]
+
+Before reduce_col 0 of m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 0 of m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+Before Matrix.rref_vec x:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+ b:
+[1 0 0 0 0 0 0 8 ]
+
+Before rref_vec we have m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+v: [1 0 0 0 0 0 0 8 ]
+
+Before normalizing we have m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+[1 0 0 0 0 0 0 8 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+Before dim_remove m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+Before reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+
+After dim_remove m':
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+
+Before add_empty_columns m:
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+indices: 2,3,4,
+Occ_cols is:
+(2,1)(3,1)(4,1)
+After add_empty_columns m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 -1 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After affineEq_vec m:
+[0 0 1 0 0 0 0 0 ]
+
+Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 -1 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 0 -1 0 0 ]
+ b:
+[0 1 -1 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 0 -1 0 0 ]
+v: [0 1 -1 0 0 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 0 -1 0 0 ]
+[0 1 -1 0 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 -1 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 -1 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+forget_vars m:
+[|#ret-8=0; x#322#arg-x#327=0; y#323#arg-z#324#arg=0|]Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 -1 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 1 0 0 0 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 8 ]
+ b:
+[0 0 0 1 -1 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 8 ]
+v: [0 0 0 1 -1 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 -1 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 -1 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 -1 0 0 0 ]
+
+forget_vars m:
+[|#ret-8=0; y#323#arg-y#328=0|]Before reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 -1 0 0 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 0 1 0 ]
+
+Before reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 8 ]
+ b:
+[0 0 0 0 0 1 -1 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 8 ]
+v: [0 0 0 0 0 1 -1 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 1 -1 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 1 -1 0 ]
+
+forget_vars m:
+[|#ret-8=0; z#324#arg-z#329=0|]Before reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+
+Before dim_remove m:
+[1 0 0 0 0 0 0 8 ]
+Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+
+Before reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 0 0 8 ]
+
+After dim_remove m':
+[1 0 0 0 8 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 0 0 8 ]
+
+Before add_empty_columns m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+indices: 0,
+Occ_cols is:
+(0,1)
+After add_empty_columns m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before rref_matrix m1 m2
+m1: [0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+m2: [1 0 0 0 8 ]
+
+Before normalizing we have m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+[1 0 0 0 8 ]
+After normalizing we have m:
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+After rref_matrix m, before removing zero rows:
+ [1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+meet a: [|x#327-z#329=0; y#328-z#329=0|] b: [|#ret-8=0|] -> [|#ret-8=0; x#327-z#329=0; y#328-z#329=0|] 
+Before dim_remove m:
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+Before reduce_col 0 of m:
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After reduce_col 0 of m:
+[0 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before del_cols cols_length=1 
+m:
+[0 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After dim_remove m':
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+Before del_cols cols_length=1 
+m:
+[0 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before add_empty_columns m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+indices: 0,
+Occ_cols is:
+(0,1)
+After add_empty_columns m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 ]
+
+Before reduce_col 0 of m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After reduce_col 0 of m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before Matrix.rref_vec x:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+ b:
+[1 0 0 0 0 ]
+
+Before rref_vec we have m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+v: [1 0 0 0 0 ]
+
+Before normalizing we have m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+[1 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before dim_remove m:
+[1 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+Before reduce_col 1 of m:
+[1 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 1 -1 0 ]
+
+Before reduce_col 2 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 1 -1 0 ]
+
+After reduce_col 2 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+After dim_remove m':
+[1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+Before add_empty_columns m:
+indices: 0,0,0,
+Occ_cols is:
+(0,3)
+After add_empty_columns m:
+
+Before assign_var_parallel m:
+⊤
+After assign_var_parallel multi_t:
+⊤
+After affineEq_vec m:
+[1 0 0 0 ]
+
+Before Matrix.rref_vec x:
+ b:
+[-1 1 0 0 ]
+
+Before rref_vec we have m:
+v: [-1 1 0 0 ]
+
+Before normalizing we have m:
+[-1 1 0 0 ]
+After normalizing we have m:
+[1 -1 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 -1 0 0 ]
+
+After affineEq_vec m:
+[0 1 0 0 ]
+
+Before reduce_col 2 of m:
+[1 -1 0 0 ]
+
+After reduce_col 2 of m:
+[1 -1 0 0 ]
+
+Before Matrix.rref_vec x:
+[1 -1 0 0 ]
+ b:
+[0 -1 1 0 ]
+
+Before rref_vec we have m:
+[1 -1 0 0 ]
+v: [0 -1 1 0 ]
+
+Before normalizing we have m:
+[1 -1 0 0 ]
+[0 -1 1 0 ]
+After normalizing we have m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+
+Before add_empty_columns m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+indices: 0,1,2,
+Occ_cols is:
+(0,1)(1,1)(2,1)
+After add_empty_columns m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 1 0 0 0 0 0 ]
+
+Before reduce_col 0 of m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After reduce_col 0 of m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+ b:
+[1 -1 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+v: [1 -1 0 0 0 0 0 ]
+
+Before normalizing we have m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[1 -1 0 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 1 0 0 0 ]
+
+Before reduce_col 2 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After reduce_col 2 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+ b:
+[0 0 1 -1 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+v: [0 0 1 -1 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 1 -1 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 1 0 ]
+
+Before reduce_col 4 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+After reduce_col 4 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+ b:
+[0 0 0 0 1 -1 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+v: [0 0 0 0 1 -1 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+
+Before dim_remove m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+Before reduce_col 1 of m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 1 -1 0 ]
+
+Before reduce_col 5 of m:
+[1 0 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 1 -1 0 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 0 0 0 ]
+
+After dim_remove m':
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 1 0 -1 0 0 ]
+[0 0 0 0 0 0 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before add_empty_columns m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+indices: 0,1,2,
+Occ_cols is:
+(0,1)(1,1)(2,1)
+After add_empty_columns m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+
+Before assign_var_parallel m:
+[|x#322#arg-z#324#arg=0; y#323#arg-z#324#arg=0|]
+Before add_empty_columns m:
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+indices: 0,0,0,
+Occ_cols is:
+(0,3)
+After add_empty_columns m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 1 0 0 0 0 0 ]
+
+Before reduce_col 0 of m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After reduce_col 0 of m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+ b:
+[1 0 0 0 -1 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+v: [1 0 0 0 -1 0 0 0 0 0 ]
+
+Before normalizing we have m:
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+[1 0 0 0 -1 0 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 0 1 0 0 0 ]
+
+Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+ b:
+[0 1 0 0 0 0 -1 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+v: [0 1 0 0 0 0 -1 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+[0 1 0 0 0 0 -1 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 0 0 0 1 0 ]
+
+Before reduce_col 2 of m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After reduce_col 2 of m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+ b:
+[0 0 1 0 0 0 0 0 -1 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+v: [0 0 1 0 0 0 0 0 -1 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After assign_var_parallel multi_t:
+[|0'-z#324#arg=0; 1'-z#324#arg=0; 2'-z#324#arg=0; x#322#arg-z#324#arg=0; y#323#arg-z#324#arg=0|]
+Matrix in Domain m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before set_col m:
+[1 0 0 0 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 0 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After set_col m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Switched Matrix in Domain switched_m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before dim_remove m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+After dim_remove m':
+[1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 1 0 0 -1 0 ]
+[0 0 1 0 0 0 0 1 -1 0 ]
+[0 0 0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 1 0 -1 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [1 0 -1 0 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[1 0 -1 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [-1 0 1 0 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[-1 0 1 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [0 0 -1 0 1 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 -1 0 1 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [0 0 1 0 -1 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 1 0 -1 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [1 0 0 0 -1 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[1 0 0 0 -1 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+v: [-1 0 0 0 1 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[-1 0 0 0 1 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+[0 0 0 0 0 0 0 ]
+
+Before add_empty_columns m:
+[1 0 0 0 0 -1 0 ]
+[0 1 0 0 0 -1 0 ]
+[0 0 1 0 0 -1 0 ]
+[0 0 0 1 0 -1 0 ]
+[0 0 0 0 1 -1 0 ]
+indices: 0,
+Occ_cols is:
+(0,1)
+After add_empty_columns m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 0 0 8 ]
+
+Before reduce_col 0 of m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 0 of m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+Before Matrix.rref_vec x:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+ b:
+[1 0 0 0 0 0 0 8 ]
+
+Before rref_vec we have m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+v: [1 0 0 0 0 0 0 8 ]
+
+Before normalizing we have m:
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+[1 0 0 0 0 0 0 8 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+Before dim_remove m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 0 -1 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 1 0 0 -1 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+Before reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+
+After dim_remove m':
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 1 0 0 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 0 1 0 -1 0 ]
+[0 0 0 0 0 0 0 0 ]
+
+Before add_empty_columns m:
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+indices: 2,3,4,
+Occ_cols is:
+(2,1)(3,1)(4,1)
+After add_empty_columns m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 -1 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After affineEq_vec m:
+[0 0 1 0 0 0 0 0 ]
+
+Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 0 0 0 -1 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 0 -1 0 0 ]
+ b:
+[0 1 -1 0 0 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 0 -1 0 0 ]
+v: [0 1 -1 0 0 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 0 -1 0 0 ]
+[0 1 -1 0 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 -1 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 -1 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+forget_vars m:
+[|#ret-8=0; x#322#arg-x#327=0; y#323#arg-z#324#arg=0|]Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 1 -1 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 1 0 0 0 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 0 -1 0 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 8 ]
+ b:
+[0 0 0 1 -1 0 0 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 8 ]
+v: [0 0 0 1 -1 0 0 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 -1 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 -1 0 0 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 -1 0 0 0 ]
+
+forget_vars m:
+[|#ret-8=0; y#323#arg-y#328=0|]Before reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 1 -1 0 0 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 0 1 0 ]
+
+Before reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+
+Before Matrix.rref_vec x:
+[1 0 0 0 0 0 0 8 ]
+ b:
+[0 0 0 0 0 1 -1 0 ]
+
+Before rref_vec we have m:
+[1 0 0 0 0 0 0 8 ]
+v: [0 0 0 0 0 1 -1 0 ]
+
+Before normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 1 -1 0 ]
+After normalizing we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 1 -1 0 ]
+
+forget_vars m:
+[|#ret-8=0; z#324#arg-z#329=0|]Before reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 1 -1 0 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+[0 0 0 0 0 0 0 0 ]
+
+Before dim_remove m:
+[1 0 0 0 0 0 0 8 ]
+Before reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 0 0 8 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 0 0 8 ]
+
+Before reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+
+After reduce_col 5 of m:
+[1 0 0 0 0 0 0 8 ]
+
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 0 0 8 ]
+
+After dim_remove m':
+[1 0 0 0 8 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 0 0 8 ]
+
+Before add_empty_columns m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+indices: 0,
+Occ_cols is:
+(0,1)
+After add_empty_columns m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before rref_matrix m1 m2
+m1: [0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+m2: [1 0 0 0 8 ]
+
+Before normalizing we have m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+[1 0 0 0 8 ]
+After normalizing we have m:
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+After rref_matrix m, before removing zero rows:
+ [1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+meet a: [|x#327-z#329=0; y#328-z#329=0|] b: [|#ret-8=0|] -> [|#ret-8=0; x#327-z#329=0; y#328-z#329=0|] 
+Before dim_remove m:
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+Before reduce_col 0 of m:
+[1 0 0 0 8 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After reduce_col 0 of m:
+[0 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before del_cols cols_length=1 
+m:
+[0 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After dim_remove m':
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+Before del_cols cols_length=1 
+m:
+[0 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before add_empty_columns m:
+[1 0 -1 0 ]
+[0 1 -1 0 ]
+indices: 0,
+Occ_cols is:
+(0,1)
+After add_empty_columns m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After affineEq_vec m:
+[0 0 0 0 0 ]
+
+Before reduce_col 0 of m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After reduce_col 0 of m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before Matrix.rref_vec x:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+ b:
+[1 0 0 0 0 ]
+
+Before rref_vec we have m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+v: [1 0 0 0 0 ]
+
+Before normalizing we have m:
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+[1 0 0 0 0 ]
+After normalizing we have m:
+[1 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+After rref_vec, before removing zero rows, we have m:
+[1 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+Before dim_remove m:
+[1 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+Before reduce_col 1 of m:
+[1 0 0 0 0 ]
+[0 1 0 -1 0 ]
+[0 0 1 -1 0 ]
+
+After reduce_col 1 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 1 -1 0 ]
+
+Before reduce_col 2 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 1 -1 0 ]
+
+After reduce_col 2 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+Before reduce_col 3 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+After reduce_col 3 of m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+After dim_remove m':
+[1 0 ]
+Before del_cols cols_length=3 
+m:
+[1 0 0 0 0 ]
+[0 0 0 0 0 ]
+[0 0 0 0 0 ]
+
+[Info][Deadcode] Logical lines of code (LLoC) summary:
+  live: 10
+  dead: 0
+  total lines: 10
+See result/index.xml
diff --git a/regtest.sh b/regtest.sh
index 488dd0bab4..f263b2b968 100755
--- a/regtest.sh
+++ b/regtest.sh
@@ -14,7 +14,7 @@ if [[ $OSTYPE == 'darwin'* ]]; then
   grep="ggrep"
 fi
 params="`$grep -oP "PARAM: \K.*" $file`"
-cmd="./goblint --enable warn.debug --enable dbg.regression --html $params ${@:3} $file" # -v
+cmd="./goblint --enable warn.debug --enable dbg.regression  $params ${@:3} $file -v" # -v --html 
 echo "$cmd"
 eval $cmd
 echo "See result/index.xml"
diff --git a/src/analyses/apron/affineEqualityAnalysis.apron.ml b/src/analyses/apron/affineEqualityAnalysis.apron.ml
index d4a1e5be2e..be94d389d1 100644
--- a/src/analyses/apron/affineEqualityAnalysis.apron.ml
+++ b/src/analyses/apron/affineEqualityAnalysis.apron.ml
@@ -4,15 +4,22 @@
 
 open Analyses
 
+open ArrayVector
+open ArrayMatrix
+
+open SparseVector
+open ListMatrix
+
 include RelationAnalysis
 
 let spec_module: (module MCPSpec) Lazy.t =
   lazy (
-    let module AD = AffineEqualityDomain.D2 (VectorMatrix.ArrayVector) (VectorMatrix.ArrayMatrix) in
+    let module AD = AffineEqualityDomain.D2 (SparseVector) (ListMatrix) in
+    let module AD_A = AffineEqualityDomain.D2 (ArrayVector) (ArrayMatrix) in (* TODO: Remove this! Just to suppress warning *)
     let module Priv = (val RelationPriv.get_priv ()) in
     let module Spec =
     struct
-      include SpecFunctor (Priv) (AD) (RelationPrecCompareUtil.DummyUtil)
+      include SpecFunctor (Priv) (AD_A) (RelationPrecCompareUtil.DummyUtil)
       let name () = "affeq"
     end
     in
diff --git a/src/cdomains/affineEquality/abstractMatrix.ml b/src/cdomains/affineEquality/abstractMatrix.ml
new file mode 100644
index 0000000000..64181c1a2b
--- /dev/null
+++ b/src/cdomains/affineEquality/abstractMatrix.ml
@@ -0,0 +1,11 @@
+open RatOps
+open AbstractVector
+open Matrix
+
+(** Some functions inside have the suffix _with, which means that the function has side effects. *)
+module type AbstractMatrix =
+  functor (A: RatOps) (V: AbstractVector) ->
+  sig
+    include Matrix with type vec := V(A).t and type num := A.t
+  end
+
diff --git a/src/cdomains/affineEquality/abstractVector.ml b/src/cdomains/affineEquality/abstractVector.ml
new file mode 100644
index 0000000000..e5d7cd40e6
--- /dev/null
+++ b/src/cdomains/affineEquality/abstractVector.ml
@@ -0,0 +1,10 @@
+open RatOps
+open Vector
+
+(** Some functions inside have the suffix _with, which means that the function has side effects. *)
+module type AbstractVector =
+  functor (A: RatOps) ->
+  sig
+    include Vector with type num:= A.t
+  end
+
diff --git a/src/cdomains/vectorMatrix.ml b/src/cdomains/affineEquality/arrayImplementation/arrayMatrix.ml
similarity index 55%
rename from src/cdomains/vectorMatrix.ml
rename to src/cdomains/affineEquality/arrayImplementation/arrayMatrix.ml
index a3712274ab..83907f29dd 100644
--- a/src/cdomains/vectorMatrix.ml
+++ b/src/cdomains/affineEquality/arrayImplementation/arrayMatrix.ml
@@ -1,280 +1,10 @@
-(** OCaml implementations of vectors and matrices. *)
+open AbstractVector
+open RatOps
+open ConvenienceOps
+open AbstractMatrix
 
 open Batteries
 module Array = Batteries.Array
-module M = Messages
-
-(* let timing_wrap = Timing.wrap *)
-(* Disable timing of VectorMatrix and AffineEqualityDomain.
-   This is cleaner than a timing functor because the timed functions also call each other. *)
-let timing_wrap _ f x = f x
-
-(** Abstracts the functions of the Mpqf module for rationals from Apron that implements multi-precision rationals.
-    One could later exchange "Mpqf" with a different module that provides the functions specified by this interface. *)
-module type RatOps =
-sig
-  type t [@@deriving eq, ord, hash]
-  val add : t -> t -> t
-  val sub : t -> t -> t
-  val mul : t -> t -> t
-  val div : t -> t -> t
-  val neg : t -> t
-  val abs : t -> t
-  val to_string:  t -> string
-  val of_int: int -> t
-  val zero: t
-  val one: t
-  val get_den: t -> Z.t
-  val get_num: t -> Z.t
-end
-
-(** It provides more readable infix operators for the functions of RatOps.
-    It is designed to be included by modules that make use of RatOps's functions. *)
-module ConvenienceOps (A: RatOps) =
-struct
-  let ( *: ) = A.mul
-  let (+:) = A.add
-  let (-:) = A.sub
-  let (/:) = A.div
-  let (=:) x y = A.equal x y
-  let (<>:) x y = not (A.equal x y)
-  let (<:) x y = A.compare x y < 0
-  let (>:) x y = A.compare x y > 0
-  let (<=:) x y = A.compare x y <= 0
-  let (>=:) x y = A.compare x y >= 0
-  let of_int x = A.of_int x
-end
-
-(** High-level abstraction of a vector. *)
-module type Vector =
-sig
-  type num
-  type t [@@deriving eq, ord, hash]
-
-  val show: t -> string
-
-  val keep_vals: t -> int ->  t
-
-  val remove_val: t -> int ->  t
-
-  val set_val: t -> int -> num ->  t
-
-  val set_val_with: t -> int -> num -> unit
-
-  val insert_val: int -> num ->  t ->  t
-
-  val apply_with_c: (num -> num -> num) -> num ->  t ->  t
-
-  val apply_with_c_with: (num -> num -> num) -> num -> t -> unit
-
-  val zero_vec: int -> t
-
-  val nth: t -> int -> num
-
-  val length: t -> int
-
-  val map2: (num -> num -> num) -> t -> t -> t
-
-  val map2_with: (num -> num -> num) -> t -> t -> unit
-
-  val findi: (num -> bool) ->  t -> int
-
-  val map: (num -> num) -> t -> t
-
-  val map_with: (num -> num) -> t -> unit
-
-  val compare_length_with: t -> int -> int
-
-  val of_list: num list -> t
-
-  val to_list: t -> num list
-
-  val filteri: (int -> num -> bool) -> t -> t
-
-  val append: t -> t -> t
-
-  val exists: (num -> bool) -> t -> bool
-
-  val rev: t -> t
-
-  val rev_with: t -> unit
-
-  val map2i: (int -> num -> num -> num) -> t -> t -> t
-
-  val map2i_with: (int -> num -> num -> num) -> t -> t -> unit
-
-  val mapi: (int -> num -> num) -> t -> t
-
-  val mapi_with: (int -> num -> num) -> t -> unit
-
-  val find2i: (num -> num -> bool) -> t -> t -> int
-
-  val to_array: t -> num array
-
-  val of_array: num array -> t
-
-  val copy: t -> t
-end
-
-(** Some functions inside have the suffix _with, which means that the function has side effects. *)
-module type AbstractVector =
-  functor (A: RatOps) ->
-  sig
-    include Vector with type num:= A.t
-  end
-
-
-(** High-level abstraction of a matrix. *)
-module type Matrix =
-sig
-  type num
-  type vec
-  type t [@@deriving eq, ord, hash]
-
-  val empty: unit -> t (* TODO: needs unit? *)
-
-  val is_empty: t -> bool
-
-  val show: t -> string
-
-  val add_empty_columns: t -> int array -> t
-
-  val append_row: t -> vec -> t
-
-  val get_row: t -> int -> vec
-
-  val del_col: t -> int -> t
-
-  val del_cols: t -> int array -> t
-
-  val remove_row: t -> int -> t
-
-  val get_col: t -> int -> vec
-
-  val append_matrices: t -> t -> t
-
-  val num_rows: t -> int
-
-  val num_cols: t -> int
-
-  val reduce_col: t -> int -> t
-
-  val reduce_col_with: t -> int -> unit
-
-  val normalize: t -> t Option.t (*Gauss-Jordan Elimination to get matrix in reduced row echelon form (rref) + deletion of zero rows. None matrix has no solution*)
-
-  val normalize_with: t -> bool
-
-  val rref_vec_with: t -> vec -> t Option.t
-
-  val rref_matrix_with: t -> t -> t Option.t
-
-  val find_opt: (vec -> bool) -> t -> vec option
-
-  val map2: (vec -> num -> vec) -> t -> vec -> t
-
-  val map2_with: (vec -> num -> vec) -> t -> vec -> unit
-
-  val map2i: (int -> vec-> num -> vec) -> t -> vec -> t
-
-  val map2i_with: (int -> vec -> num -> vec) -> t -> vec -> unit
-
-  val set_col: t -> vec -> int -> t
-
-  val set_col_with: t -> vec -> int -> t
-
-  val init_with_vec: vec -> t
-
-  val remove_zero_rows: t -> t
-
-  val is_covered_by: t -> t -> bool
-
-  val copy: t -> t
-
-end
-
-(** Some functions inside have the suffix _with, which means that the function has side effects. *)
-module type AbstractMatrix =
-  functor (A: RatOps) (V: AbstractVector) ->
-  sig
-    include Matrix with type vec := V(A).t and type num := A.t
-  end
-
-
-(** Array-based vector implementation. *)
-module ArrayVector: AbstractVector =
-  functor (A: RatOps) ->
-  struct
-    include ConvenienceOps (A)
-    include Array
-    type t = A.t array [@@deriving eq, ord, hash]
-
-    let show t =
-      let t = Array.to_list t in
-      let rec list_str l =
-        match l with
-        | [] -> "]"
-        | x :: xs -> (A.to_string x) ^" "^(list_str xs)
-      in
-      "["^list_str t^"\n"
-
-    let keep_vals v n =
-      if n >= Array.length v then v else
-        Array.filteri (fun i x -> i < n) v (* TODO: take? *)
-
-    let compare_length_with v len =
-      Int.compare (Array.length v) len
-
-    let remove_val v n =
-      if n >= Array.length v then failwith "n outside of Array range" else
-        Array.init (Array.length v - 1) (fun i -> if i < n then Array.get v i else Array.get v (i + 1)) (* TODO: remove_at? *)
-
-    let set_val_with v n new_val =
-      if n >= Array.length v then failwith "n outside of Array range" else
-        Array.set v n new_val
-
-    let set_val v n new_val =
-      let copy = copy v in
-      set_val_with copy n new_val; copy
-
-    let insert_val n new_val v =
-      if n > Array.length v then failwith "n too large" else
-        Array.init (Array.length v + 1) (fun i -> if i < n then Array.get v i else if i = n then new_val else Array.get v (i -1)) (* insert? *)
-
-    let apply_with_c f c v =
-      Array.map (fun x -> f x c) v
-
-    let zero_vec n = Array.make n A.zero
-
-    let nth = Array.get
-
-    let map2i f v1 v2 =
-      let f' i = uncurry (f i) in
-      Array.mapi f' (Array.combine v1 v2) (* TODO: iter2i? *)
-
-    let map2i_with f v1 v2 = Array.iter2i (fun i x y -> v1.(i) <- f i x y) v1 v2
-
-    let find2i f v1 v2 =
-      Array.findi (uncurry f) (Array.combine v1 v2) (* TODO: iter2i? *)
-
-    let to_array v = v
-
-    let of_array v = v
-
-    let apply_with_c_with f c v = Array.modify (fun x -> f x c) v
-
-    let rev_with v = Array.rev_in_place v
-
-    let map_with f v = Array.modify f v
-
-    let map2_with f v1 v2 = Array.iter2i (fun i x y -> v1.(i) <- f x y) v1 v2
-
-    let copy v = Array.copy v
-
-    let mapi_with f v = Array.iteri (fun i x -> v.(i) <- f i x) v
-  end
-
-open Batteries.Array
 
 (** Array-based matrix implementation.
     It provides a normalization function to reduce a matrix into reduced row echelon form.
@@ -306,13 +36,14 @@ module ArrayMatrix: AbstractMatrix =
       let cp = Array.make_matrix (num_rows m) (num_cols m) A.zero in
       Array.iteri (fun i x -> Array.blit x 0 cp.(i) 0 (num_cols m)) m; cp
 
-    let copy m = timing_wrap "copy" (copy) m
+    let copy m = Timing.wrap "copy" (copy) m
 
     let add_empty_columns m cols =
+      Array.modifyi (+) cols;
       let nnc = Array.length cols in
       if is_empty m || nnc = 0 then m else
         let nr, nc = num_rows m, num_cols m in
-        let m' = make_matrix nr (nc + nnc) A.zero in
+        let m' = Array.make_matrix nr (nc + nnc) A.zero in
         for i = 0 to nr - 1 do
           let offset = ref 0 in
           for j = 0 to nc - 1 do
@@ -322,11 +53,11 @@ module ArrayMatrix: AbstractMatrix =
         done;
         m'
 
-    let add_empty_columns m cols = timing_wrap "add_empty_cols" (add_empty_columns m) cols
+    let add_empty_columns m cols = Timing.wrap "add_empty_cols" (add_empty_columns m) cols
 
     let append_row m row  =
       let size = num_rows m in
-      let new_matrix = make_matrix (size + 1) (num_cols m) A.zero in
+      let new_matrix = Array.make_matrix (size + 1) (num_cols m) A.zero in
       for i = 0 to size - 1 do
         new_matrix.(i) <- m.(i)
       done;
@@ -337,7 +68,7 @@ module ArrayMatrix: AbstractMatrix =
       V.of_array m.(n)
 
     let remove_row m n =
-      let new_matrix = make_matrix (num_rows m - 1) (num_cols m) A.zero in
+      let new_matrix = Array.make_matrix (num_rows m - 1) (num_cols m) A.zero in
       if not @@ is_empty new_matrix then
         if n = 0 then
           Array.blit m 1 new_matrix 0 (num_rows m - 1)
@@ -349,14 +80,15 @@ module ArrayMatrix: AbstractMatrix =
     let get_col m n =
       V.of_array @@ Array.init (Array.length m) (fun i -> m.(i).(n))
 
-    let get_col m n = timing_wrap "get_col" (get_col m) n
+    let get_col m n = Timing.wrap "get_col" (get_col m) n
 
     let set_col_with m new_col n =
       for i = 0 to num_rows m - 1 do
         m.(i).(n) <- V.nth new_col i
-      done; m
+      done;
+      m
 
-    let set_col_with m new_col n = timing_wrap "set_col" (set_col_with m new_col) n
+    let set_col_with m new_col n = Timing.wrap "set_col" (set_col_with m new_col) n
 
     let set_col m new_col n =
       let copy = copy m in
@@ -365,7 +97,7 @@ module ArrayMatrix: AbstractMatrix =
     let append_matrices m1 m2  =
       Array.append m1 m2
 
-    let equal m1 m2 = timing_wrap "equal" (equal m1) m2
+    let equal m1 m2 = Timing.wrap "equal" (equal m1) m2
 
     let reduce_col_with m j =
       if not @@ is_empty m then
@@ -383,7 +115,7 @@ module ArrayMatrix: AbstractMatrix =
          done;
          if !r >= 0 then Array.fill m.(!r) 0 (num_cols m) A.zero)
 
-    let reduce_col_with m j  = timing_wrap "reduce_col_with" (reduce_col_with m) j
+    let reduce_col_with m j  = Timing.wrap "reduce_col_with" (reduce_col_with m) j
     let reduce_col m j =
       let copy = copy m in
       reduce_col_with copy j;
@@ -412,12 +144,7 @@ module ArrayMatrix: AbstractMatrix =
           done;
           m'
 
-    let del_cols m cols = timing_wrap "del_cols" (del_cols m) cols
-
-    let map2i f m v =
-      let f' x (i,y) = V.to_array @@ f i (V.of_array x) y in
-      let range_array = Array.init (V.length v) Fun.id in
-      Array.map2 f' m (Array.combine range_array (V.to_array v))
+    let del_cols m cols = Timing.wrap "del_cols" (del_cols m) cols
 
     let remove_zero_rows m =
       Array.filter (fun x -> Array.exists (fun y -> y <>: A.zero) x) m
@@ -436,19 +163,19 @@ module ArrayMatrix: AbstractMatrix =
         for i = 0 to num_rows-1 do
           let exception Found in
           try (
-            for j = i to num_cols -2 do
+            for j = i to num_cols -2 do (* Find pivot *)
               for k = i to num_rows -1 do
                 if m.(k).(j) <>: A.zero then
                   (
                     if k <> i then swap_rows k i;
                     let piv = m.(i).(j) in
-                    Array.iteri(fun j' x -> m.(i).(j') <- x /: piv) m.(i);
-                    for l = 0 to num_rows-1 do
+                    Array.iteri(fun j' x -> m.(i).(j') <- x /: piv) m.(i); (* Normalize pivot *)
+                    for l = 0 to num_rows-1 do (* Subtract from each row *)
                       if l <> i && m.(l).(j) <>: A.zero then (
                         let is_only_zero = ref true in
                         let m_lj = m.(l).(j) in
                         for k = 0 to num_cols - 2 do
-                          m.(l).(k) <- m.(l).(k) -: m.(i).(k) *: m_lj /: m.(i).(j);
+                          m.(l).(k) <- m.(l).(k) -: m.(i).(k) *: m_lj /: m.(i).(j); (* Subtraction *)
                           if m.(l).(k) <>: A.zero then is_only_zero := false;
                         done;
                         let k_end = num_cols - 1 in
@@ -466,7 +193,7 @@ module ArrayMatrix: AbstractMatrix =
         true)
       with Unsolvable -> false
 
-    let rref_with m = timing_wrap "rref_with" rref_with m
+    let rref_with m = Timing.wrap "rref_with" rref_with m
 
     let init_with_vec v =
       let new_matrix = Array.make_matrix 1 (V.length v) A.zero in
@@ -480,11 +207,12 @@ module ArrayMatrix: AbstractMatrix =
           Array.iteri (fun j' x ->  m.(i).(j') <- x -: beta *: v.(j')) m.(i)
       done
 
+
     let get_pivot_positions m =
       let pivot_elements = Array.make (num_rows m) 0
       in Array.iteri (fun i x -> pivot_elements.(i) <- Array.findi (fun z -> z =: A.one) x) m; pivot_elements
 
-    let rref_vec m pivot_positions v =
+    let rref_vec_helper m pivot_positions v =
       let insert = ref (-1) in
       for j = 0 to Array.length v -2 do
         if v.(j) <>: A.zero then
@@ -508,23 +236,41 @@ module ArrayMatrix: AbstractMatrix =
         else (Array.blit m 0 new_m 0 i; new_m.(i) <- v; Array.blit m i new_m (i + 1) (Array.length m - j));
         Some new_m
 
+    let normalize_with m =
+      rref_with m
+
+    let normalize_with m = Timing.wrap "normalize_with" normalize_with m
 
+    let normalize m =
+      let copy = copy m in
+      if normalize_with copy then
+        Some copy
+      else
+        None
     let rref_vec_with m v =
       (*This function yields the same result as appending vector v to m and normalizing it afterwards would. However, it is usually faster than performing those ops manually.*)
       (*m must be in rref form and contain the same num of cols as v*)
       (*If m is empty then v is simply normalized and returned*)
-      let v = V.to_array v in
-      if is_empty m then
+      (*let v = V.to_array v in
+        if is_empty m then
         match Array.findi (fun x -> x <>: A.zero) v with
         | exception Not_found -> None
         | i -> if i = Array.length v - 1 then None else
             let v_i = v.(i) in
             Array.iteri (fun j x -> v.(j) <- x /: v_i) v; Some (init_with_vec @@ V.of_array v)
-      else
+        else
         let pivot_elements = get_pivot_positions m in
-        rref_vec m pivot_elements v
+        rref_vec_helper m pivot_elements v*)
+      normalize @@ append_row m v
 
-    let rref_vec_with m v = timing_wrap "rref_vec_with" (rref_vec_with m) v
+    let rref_vec_with m v = Timing.wrap "rref_vec_with" (rref_vec_with m) v
+
+    let rref_vec m v = (* !! There was another rref_vec function that has been renamed to rref_vec_helper !!*)
+      let m' = copy m in
+      let v' = V.copy v in 
+      match rref_vec_with m' v' with
+      | Some res -> Some (remove_zero_rows res)
+      | None -> None
 
     let rref_matrix_with m1 m2 =
       (*Similar to rref_vec_with but takes two matrices instead.*)
@@ -535,7 +281,7 @@ module ArrayMatrix: AbstractMatrix =
       try (
         for i = 0 to num_rows s_m - 1 do
           let pivot_elements = get_pivot_positions !b in
-          let res = rref_vec !b pivot_elements s_m.(i) in
+          let res = rref_vec_helper !b pivot_elements s_m.(i) in
           match res with
           | None -> raise Unsolvable
           | Some res -> b := res
@@ -544,19 +290,14 @@ module ArrayMatrix: AbstractMatrix =
       )
       with Unsolvable -> None
 
-    let rref_matrix_with m1 m2 = timing_wrap "rref_matrix_with" (rref_matrix_with m1) m2
-
-    let normalize_with m =
-      rref_with m
+    let rref_matrix_with m1 m2 = Timing.wrap "rref_matrix_with" (rref_matrix_with m1) m2
 
-    let normalize_with m = timing_wrap "normalize_with" normalize_with m
-
-    let normalize m =
-      let copy = copy m in
-      if normalize_with copy then
-        Some copy
-      else
-        None
+    let rref_matrix m1 m2 = 
+      let m1' = copy m1 in
+      let m2' = copy m2 in 
+      match rref_matrix_with m1' m2' with
+      | Some m -> Some m
+      | None -> None
 
     let is_covered_by m1 m2 =
       (*Performs a partial rref reduction to check if concatenating both matrices and afterwards normalizing them would yield a matrix <> m2 *)
@@ -565,6 +306,7 @@ module ArrayMatrix: AbstractMatrix =
         let p2 = lazy (get_pivot_positions m2) in
         try (
           for i = 0 to num_rows m1 - 1 do
+            (* check if there are rows in m1 and m2 that aren't equal *)
             if Array.exists2 (<>:) m1.(i) m2.(i) then
               let m1_i = Array.copy m1.(i) in
               for j = 0 to Array.length m1_i - 2 do
@@ -581,13 +323,11 @@ module ArrayMatrix: AbstractMatrix =
         )
         with Stdlib.Exit -> false;;
 
-    let is_covered_by m1 m2 = timing_wrap "is_covered_by" (is_covered_by m1) m2
+    let is_covered_by m1 m2 = Timing.wrap "is_covered_by" (is_covered_by m1) m2
 
     let find_opt f m =
       let f' x = f (V.of_array x) in Option.map V.of_array (Array.find_opt f' m)
 
-    let map2 f m v =
-      let f' x y = V.to_array @@ f (V.of_array x) y in Array.map2 f' m (V.to_array v)
 
     let map2_with f m v =
       if num_rows m = V.length v then
@@ -597,7 +337,12 @@ module ArrayMatrix: AbstractMatrix =
           m.(i) <- V.to_array @@ f (V.of_array m.(i)) (V.nth v i)
         done
 
-    let map2_with f m v = timing_wrap "map2_with" (map2_with f m) v
+    let map2_with f m v = Timing.wrap "map2_with" (map2_with f m) v
+
+    let map2 f m v =
+      let m' = copy m in
+      map2_with f m' v;
+      m'
 
     let map2i_with f m v =
       if num_rows m = V.length v then
@@ -607,5 +352,14 @@ module ArrayMatrix: AbstractMatrix =
           m.(i) <- V.to_array @@ f i (V.of_array m.(i)) (V.nth v i)
         done
 
-    let map2i_with f m v = timing_wrap "map2i_with" (map2i_with f m) v
-  end
+    let map2i_with f m v = Timing.wrap "map2i_with" (map2i_with f m) v    
+
+    let map2i f m v =
+      let m' = copy m in
+      map2i_with f m' v;
+      m'
+
+    let swap_rows m j k =
+      failwith "TODO"
+
+  end
\ No newline at end of file
diff --git a/src/cdomains/affineEquality/arrayImplementation/arrayVector.ml b/src/cdomains/affineEquality/arrayImplementation/arrayVector.ml
new file mode 100644
index 0000000000..e90a5afbce
--- /dev/null
+++ b/src/cdomains/affineEquality/arrayImplementation/arrayVector.ml
@@ -0,0 +1,145 @@
+open AbstractVector
+open RatOps
+open ConvenienceOps
+
+open Batteries
+module Array = Batteries.Array
+
+
+(** Array-based vector implementation. *)
+module ArrayVector: AbstractVector =
+  functor (A: RatOps) ->
+  struct
+    include ConvenienceOps (A)
+    include Array
+    type t = A.t array [@@deriving eq, ord, hash]
+
+    let show t =
+      let t = Array.to_list t in
+      let rec list_str l =
+        match l with
+        | [] -> "]"
+        | x :: xs -> (A.to_string x) ^" "^(list_str xs)
+      in
+      "["^list_str t^"\n"
+
+    let keep_vals v n =
+      if n >= Array.length v then v else
+        Array.filteri (fun i x -> i < n) v (* TODO: take? *)
+
+    let compare_length_with v len =
+      Int.compare (Array.length v) len
+
+    let remove_nth v n =
+      if n >= Array.length v then failwith "n outside of Array range" else
+        Array.init (Array.length v - 1) (fun i -> if i < n then Array.get v i else Array.get v (i + 1)) (* TODO: remove_at? *)
+
+    let remove_at_indices v idx = failwith "TODO"
+
+    let insert_zero_at_indices v idx count = failwith "TODO"
+
+    let set_nth_with v n new_val =
+      if n >= Array.length v then failwith "n outside of Array range" else
+        Array.set v n new_val
+
+    let set_nth v n new_val =
+      let copy = copy v in
+      set_nth_with copy n new_val; copy
+
+    let insert_val_at n new_val v =
+      if n > Array.length v then failwith "n too large" else
+        Array.init (Array.length v + 1) (fun i -> if i < n then Array.get v i else if i = n then new_val else Array.get v (i -1)) (* insert? *)
+
+    let apply_with_c f c v =
+      Array.map (fun x -> f x c) v
+
+    let zero_vec n = Array.make n A.zero
+
+    let is_zero_vec v = not (Array.exists (fun x -> x <>: A.zero) v)
+
+    let is_const_vec v =
+      compare_length_with (filteri (fun i x -> (*Inefficient*)
+        compare_length_with v (i + 1) > 0 && x <>: A.zero) v) 1 = 0
+
+    let nth = Array.get
+
+    let map2i f v1 v2 =
+      let f' i = uncurry (f i) in
+      Array.mapi f' (Array.combine v1 v2) (* TODO: iter2i? *)
+
+    let find2i f v1 v2 =
+      Array.findi (uncurry f) (Array.combine v1 v2) (* TODO: iter2i? *)
+
+    let to_array v = v
+
+    let of_array v = v
+
+    let rev v = Array.rev v
+
+    let map f v = Array.map f v
+
+    let map2_with f v1 v2 = Array.iter2i (fun i x y -> v1.(i) <- f x y) v1 v2
+
+    let map2 f v1 v2 =
+      let copy_v1 = copy v1 in
+      map2_with f copy_v1 v2; copy_v1
+
+    let copy v = Array.copy v
+
+    let mapi_with f v = Array.iteri (fun i x -> v.(i) <- f i x) v
+
+    let mapi f v =
+      let copy = copy v in
+      mapi_with f copy; copy
+
+    let of_sparse_list col_count ls =
+      let vec = Array.make col_count A.zero in
+      List.iter (fun (idx, value) -> vec.(idx) <- value) ls;
+      vec
+
+    let to_sparse_list v = 
+      let rec aux idx acc =
+        if idx < 0 then acc
+        else
+          let value = v.(idx) in
+          let acc = if value <> A.zero then (idx, value):: acc else acc in
+          aux (idx - 1) acc
+      in aux (length v - 1) []
+
+
+    let find_opt f v =
+      failwith "Never implemented!"
+
+    let map_f_preserves_zero f v =
+      map f v
+
+    let map2_f_preserves_zero f v1 v2 =
+      map2 f v1 v2
+
+    let mapi_f_preserves_zero f v =
+      mapi f v
+
+    let map2i_f_preserves_zero f v v' =
+      map2i f v v'
+
+    let fold_left_f_preserves_zero f acc v =
+      failwith "Never implemented!"
+
+    let fold_left2_f_preserves_zero f acc v v' =
+      failwith "Never implemented!"
+
+    let findi_val_opt f v =
+      failwith "Never implemented!"
+
+    let exists2 f v1 v1 =
+      failwith "Never implemented!"
+
+    let starting_from_nth n v =
+      failwith "Never implemented!"
+
+    let find_first_non_zero v =
+      failwith "Never implemented!"
+
+    let apply_with_c_f_preserves_zero f c v =
+      apply_with_c f c v
+  end
\ No newline at end of file
diff --git a/src/cdomains/affineEquality/convenienceOps.ml b/src/cdomains/affineEquality/convenienceOps.ml
new file mode 100644
index 0000000000..39aeeeb5e8
--- /dev/null
+++ b/src/cdomains/affineEquality/convenienceOps.ml
@@ -0,0 +1,18 @@
+open RatOps
+
+(** It provides more readable infix operators for the functions of RatOps.
+    It is designed to be included by modules that make use of RatOps's functions. *)
+module ConvenienceOps (A: RatOps) =
+struct
+  let ( *: ) = A.mul
+  let (+:) = A.add
+  let (-:) = A.sub
+  let (/:) = A.div
+  let (=:) x y = A.equal x y
+  let (<>:) x y = not (A.equal x y)
+  let (<:) x y = A.compare x y < 0
+  let (>:) x y = A.compare x y > 0
+  let (<=:) x y = A.compare x y <= 0
+  let (>=:) x y = A.compare x y >= 0
+  let of_int x = A.of_int x
+end
diff --git a/src/cdomains/affineEquality/matrix.ml b/src/cdomains/affineEquality/matrix.ml
new file mode 100644
index 0000000000..b3433ee899
--- /dev/null
+++ b/src/cdomains/affineEquality/matrix.ml
@@ -0,0 +1,60 @@
+(** High-level abstraction of a matrix. *)
+module type Matrix =
+sig
+  type num
+  type vec
+  type t [@@deriving eq, ord, hash]
+
+  val show: t -> string
+
+  val copy: t -> t
+
+  val empty: unit -> t (* TODO: needs unit? *)
+
+  val is_empty: t -> bool
+
+  val num_rows: t -> int
+
+  val num_cols: t -> int
+
+  val init_with_vec: vec -> t
+
+  val append_row: t -> vec -> t
+
+  val get_row: t -> int -> vec
+
+  val remove_row: t -> int -> t
+
+  val remove_zero_rows: t -> t
+
+  val swap_rows: t -> int -> int -> t
+
+  val add_empty_columns: t -> int array -> t
+
+  val get_col: t -> int -> vec
+
+  val set_col: t -> vec -> int -> t
+
+  val del_col: t -> int -> t
+
+  val del_cols: t -> int array -> t
+
+  val map2: (vec -> num -> vec) -> t -> vec -> t
+
+  val map2i: (int -> vec-> num -> vec) -> t -> vec -> t
+
+  val find_opt: (vec -> bool) -> t -> vec option
+
+  val append_matrices: t -> t -> t
+
+  val reduce_col: t -> int -> t
+
+  val normalize: t -> t Option.t (*Gauss-Jordan Elimination to get matrix in reduced row echelon form (rref) + deletion of zero rows. None matrix has no solution*)
+
+  val rref_vec: t -> vec -> t Option.t
+
+  val rref_matrix: t -> t -> t Option.t
+
+  val is_covered_by: t -> t -> bool
+
+end
\ No newline at end of file
diff --git a/src/cdomains/affineEquality/ratOps.ml b/src/cdomains/affineEquality/ratOps.ml
new file mode 100644
index 0000000000..c8d8017693
--- /dev/null
+++ b/src/cdomains/affineEquality/ratOps.ml
@@ -0,0 +1,18 @@
+(** Abstracts the functions of the Mpqf module for rationals from Apron that implements multi-precision rationals.
+    One could later exchange "Mpqf" with a different module that provides the functions specified by this interface. *)
+module type RatOps =
+sig
+  type t [@@deriving eq, ord, hash]
+  val add : t -> t -> t
+  val sub : t -> t -> t
+  val mul : t -> t -> t
+  val div : t -> t -> t
+  val neg : t -> t
+  val abs : t -> t
+  val to_string:  t -> string
+  val of_int: int -> t
+  val zero: t
+  val one: t
+  val get_den: t -> Z.t
+  val get_num: t -> Z.t
+end
\ No newline at end of file
diff --git a/src/cdomains/affineEquality/sparseImplementation/listMatrix.ml b/src/cdomains/affineEquality/sparseImplementation/listMatrix.ml
new file mode 100644
index 0000000000..e8bc108d91
--- /dev/null
+++ b/src/cdomains/affineEquality/sparseImplementation/listMatrix.ml
@@ -0,0 +1,328 @@
+open AbstractMatrix
+open AbstractVector
+open RatOps
+open ConvenienceOps
+
+open BatList
+module List = BatList
+
+(** Sparse matrix implementation.
+    It provides a normalization function to reduce a matrix into reduced row echelon form.
+    Operations exploit that the input matrix/matrices are in reduced row echelon form already. *)
+module ListMatrix: AbstractMatrix =
+  functor (A: RatOps) (V: AbstractVector) ->
+  struct
+    include ConvenienceOps(A)
+    module V = V(A)
+
+    type t = V.t list (*List of rows*)
+    [@@deriving eq, ord, hash]
+
+    let show x =
+      List.fold_left (^) "" (List.map (fun x -> (V.show x)) x)
+
+    let empty () = []
+
+    let num_rows = List.length
+
+    let is_empty m =
+      num_rows m = 0
+    (*This should be different if the implimentation is sound*)
+    (*m.column_count = 0*)
+
+    let num_cols m = if m = [] then 0 else V.length (hd m)
+
+    let copy m = m
+    (* Lists are immutable, so this should suffice? A.t is mutuable currently, but is treated like its not in ArrayMatrix*)
+
+    let copy m =
+      Timing.wrap "copy" (copy) m
+
+    (* This only works if Array.modifyi has been removed from dim_add *)
+    let add_empty_columns (m : t) (cols : int array) : t =    
+      let cols = Array.to_list cols in 
+      let sorted_cols = List.sort Stdlib.compare cols in
+      let rec count_sorted_occ acc cols last count =
+        match cols with
+        | [] -> if count > 0 then (last, count) :: acc else acc
+        | x :: xs when x = last -> count_sorted_occ acc xs x (count + 1)
+        | x :: xs -> let acc = if count > 0 then (last, count) :: acc else acc in
+          count_sorted_occ acc xs x 1      
+      in
+      let occ_cols = List.rev @@ count_sorted_occ [] sorted_cols 0 0 in
+      List.map (fun row -> V.insert_zero_at_indices row occ_cols (List.length cols)) m
+
+    let add_empty_columns m cols =
+      Timing.wrap "add_empty_cols" (add_empty_columns m) cols
+
+    let append_row m row  =
+      m @ [row]
+
+    let get_row m n =
+      List.nth m n
+
+    let remove_row m n =
+      List.remove_at n m
+
+    let get_col m n =
+      (*let () = Printf.printf "get_col %i of m:\n%s\n%s\n" n (show m) (V.show (V.of_list @@ List.map (fun row -> V.nth row n) m)) in*)
+      V.of_list @@ List.map (fun row -> V.nth row n) m (* builds full col including zeros, maybe use sparselist instead? *)
+
+    let get_col m n =
+      Timing.wrap "get_col" (get_col m) n
+
+    let set_col m new_col n = 
+      (* List.mapi (fun row_idx row -> V.set_nth row n (V.nth new_col row_idx)) m *)
+      List.map2 (fun row value -> V.set_nth row n value) m (V.to_list new_col)
+
+    let append_matrices m1 m2  = (* keeps dimensions of first matrix, what if dimensions differ?*)
+      m1 @ m2
+
+    let equal m1 m2 = Timing.wrap "equal" (equal m1) m2
+
+    let div_row (row : V.t) (pivot : A.t) : V.t =
+      V.map_f_preserves_zero (fun a -> a /: pivot) row
+
+    let swap_rows m j k =
+      List.mapi (fun i row -> if i = j then List.nth m k else if i = k then List.nth m j else row) m
+
+    let sub_scaled_row row1 row2 s =
+      V.map2_f_preserves_zero (fun x y -> x -: (s *: y)) row1 row2
+
+    (* Reduces the jth column with the last row that has a non-zero element in this column. *)
+    let reduce_col m j = 
+      if is_empty m then m 
+      else
+        let rec find_pivot idx entries = (* Finds non-zero element in column j and returns pair of row idx and the pivot value *)
+          match entries with
+          | [] -> None
+          | row :: rest -> let value = V.nth row j in 
+            if value =: A.zero then find_pivot (idx - 1) rest else Some (idx, value)
+        in
+        match (find_pivot (num_rows m - 1) (List.rev m)) with
+        | None -> m (* column is already filled with zeroes *)
+        | Some (row_idx, pivot) -> 
+          let pivot_row = List.nth m row_idx in (* use the pivot row to reduce all entries in column j to zero, then "delete" the pivot row *)
+          List.mapi (fun idx row ->
+              if idx = row_idx then 
+                V.zero_vec (num_cols m)
+              else
+                let row_value = V.nth row j in 
+                if row_value = A.zero then row
+                else (let s = row_value /: pivot in
+                      sub_scaled_row row pivot_row s)            
+            ) m
+
+    let del_col m j =
+      if num_cols m = 1 then empty () 
+      else 
+        List.map (fun row -> V.remove_nth row j) m
+
+    let del_cols m cols =
+      let cols = Array.to_list cols in (* TODO: Is it possible to use list for Apron dimchange? *)
+      let sorted_cols = List.sort_uniq Stdlib.compare cols in (* Apron Docs:  Repetitions are meaningless (and are not correct specification) *)
+      if (List.length sorted_cols) = num_cols m then empty() 
+      else
+        List.map (fun row -> V.remove_at_indices row sorted_cols) m
+
+    let del_cols m cols = Timing.wrap "del_cols" (del_cols m) cols
+
+    let map2i f m v =
+      let rec map2i_min i acc m v =
+        match m, v with
+        | [], _  -> List.rev acc
+        | row :: rs, [] -> List.rev_append (row :: acc) rs
+        | row :: rs, value :: vs -> map2i_min (i + 1) (f i row value :: acc) rs vs
+      in  
+      map2i_min 0 [] m (V.to_list v)
+
+    let remove_zero_rows m =
+      List.filter (fun row -> not (V.is_zero_vec row)) m
+
+    let init_with_vec v =
+      [v]
+
+    let normalize m =
+      let col_count = num_cols m in
+      let dec_mat_2D (m : t) (row_idx : int) (col_idx : int) : t = 
+        List.filteri_map (fun i row -> if i < row_idx then None else Some (V.starting_from_nth col_idx row)) m
+      in
+      let dec_mat_2D m row_idx col_idx = Timing.wrap "dec_mat_2D" (dec_mat_2D m row_idx) col_idx in
+      (* Function for finding first pivot in an extracted part of the matrix (row_idx and col_idx indicate which part of the original matrix) *)
+      (* The last column represents the constant in the affeq *)
+      let find_first_pivot m' row_idx col_idx =
+        if col_idx >= col_count then None else 
+          let max_piv_col_idx = num_cols m' - 2 in (* col at num_cols - 1 is the constant of the affeq *)
+          (* Finding pivot by extracting the minimum column index of the first non zero value of each row*)
+          let (piv_row, piv_col, piv_val) = List.fold_lefti (fun (cur_row, cur_col, cur_val) i row  -> 
+              let row_first_non_zero = V.find_first_non_zero row in
+              match row_first_non_zero with
+              | None -> (cur_row, cur_col, cur_val)
+              | Some (idx, value) -> (* let () = Printf.printf "We found first non-zero at index %i in row %i\n" idx i in *)
+                if idx < cur_col then (i, idx, value) else (cur_row, cur_col, cur_val)
+            ) (num_rows m', max_piv_col_idx + 1, A.zero) m' (* Initializing with max, so num_cols m indicates that pivot is not found *)
+          in
+          if piv_col = (max_piv_col_idx + 1) then None else Some (row_idx + piv_row, col_idx + piv_col, piv_val)
+      in
+      let find_first_pivot m' row_idx col_idx = Timing.wrap "find_first_pivot" (find_first_pivot m' row_idx) col_idx in
+      let affeq_rows_are_valid m = (* Check if the semantics of an rref-affeq matrix are correct *)
+        let col_count = num_cols m in
+        let row_is_valid row = (* TODO: Vector findi_opt *)
+          match V.find_first_non_zero row with
+          | Some (idx, _) -> if idx < col_count - 1 then true else false (* If all cofactors of the affeq are zero, but the constant is non-zero, the row is invalid *)
+          | None -> true (* Full zero row is valid *)
+        in
+        List.for_all row_is_valid m in
+      let rec main_loop m m' row_idx col_idx = 
+        if col_idx >= (col_count - 1) then m  (* In this case the whole bottom of the matrix starting from row_index is Zero, so it is normalized *)
+        else
+          match find_first_pivot m' row_idx col_idx with
+          | None -> m (* No pivot found means already normalized*)
+          | Some (piv_row_idx, piv_col_idx, piv_val) -> (
+              (* let () = Printf.printf "The current matrix is: \n%s and the pivot is (%i, %i, %s)\n" (show m) piv_row_idx piv_col_idx (A.to_string piv_val) in *)
+              let m = if piv_row_idx <> row_idx then swap_rows m row_idx piv_row_idx else m in
+              let normalized_m = List.mapi (fun idx row -> if idx = row_idx then div_row row piv_val else row) m in
+              let piv_row = (List.nth normalized_m row_idx) in
+              let subtracted_m = List.mapi (fun idx row -> if idx <> row_idx then 
+                                               let scale = V.nth row piv_col_idx in
+                                               sub_scaled_row row piv_row scale else row) normalized_m in
+              let m' = dec_mat_2D subtracted_m (row_idx + 1) (piv_col_idx + 1) in
+              main_loop subtracted_m m' (row_idx + 1) (piv_col_idx + 1)) (* We start at piv_col_idx + 1 because every other col before that is zero at the bottom*)
+      in 
+      let m' = main_loop m m 0 0 in
+      if affeq_rows_are_valid m' then Some m' else None (* TODO: We can check this for each row, using the helper function row_is_invalid *)
+
+
+    (* This function return a tuple of row index and pivot position (column) in m *)
+    (* TODO: maybe we could use a Hashmap instead of a list? *)
+    let get_pivot_positions (m : t) : (int * int) list =
+      List.rev @@ List.fold_lefti (
+        fun acc i row -> match V.find_first_non_zero row with
+          | None -> acc
+          | Some (pivot_col, _) -> (i, pivot_col) :: acc
+      ) [] m
+
+    let assert_rref m = 
+      let pivot_l = get_pivot_positions m in
+      let rec validate m i =
+        match m with 
+        | [] -> ()
+        | v::vs when (V.is_zero_vec v) -> 
+          if List.exists (fun v -> not @@ V.is_zero_vec v) vs 
+          then raise (Invalid_argument "Matrix not in rref: zero row!")
+          else ()
+        | v::vs -> 
+          let rec validate_vec pl =  
+            match pivot_l with
+            | [] -> true
+            | (pr, pc)::ps -> 
+              let target = if pr <> i then A.zero else A.one in 
+              if V.nth v pc <>: target then false else validate_vec ps 
+          in if validate_vec pivot_l then validate vs (i+1) else raise (Invalid_argument "Matrix not in rref: pivot column not empty!")
+      in validate m 0
+
+    (* TODO: Remove this! Just to suppress warning *)
+    let () = assert_rref (empty ())
+
+    (* Sets the jth column to zero by subtracting multiples of v *)
+    let reduce_col_with_vec m j v = 
+      let pivot_element = V.nth v j in
+      if pivot_element = A.zero then m
+      else List.mapi (fun idx row ->
+          let row_value = V.nth row j in 
+          if row_value =: A.zero then row
+          else (let s = row_value /: pivot_element in
+                V.map2_f_preserves_zero (fun x y -> x -: (s *: y)) row v)
+        ) m
+
+    (* Inserts the vector v with pivot at piv_idx at the correct position in m. m has to be in rref form and v is only <>: A.zero on piv_idx or idx not included in piv_positions *)
+    let insert_v_according_to_piv m v piv_idx pivot_positions = 
+      let reduced_m = reduce_col_with_vec m piv_idx v in
+      match List.find_opt (fun (row_idx, piv_col) -> piv_col > piv_idx) pivot_positions with
+      | None -> append_row reduced_m v
+      | Some (row_idx, _) ->
+        let (before, after) = List.split_at row_idx reduced_m in (* TODO: Optimize *)
+        before @ (v :: after)
+
+    (* This function yields the same result as appending v to m, normalizing and removing zero rows would. *)
+    (* m must be in rref form and must contain the same number of columns as v *)
+    let rref_vec m v =
+      if is_empty m then (* In this case, v is normalized and returned *)
+        match V.find_first_non_zero v with 
+        | None -> None
+        | Some (_, value) -> 
+          let normalized_v = V.map_f_preserves_zero (fun x -> x /: value) v in
+          Some (init_with_vec normalized_v)
+      else (* We try to normalize v and check if a contradiction arises. If not, we insert v at the appropriate place in m (depending on the pivot) *)
+        let pivot_positions = get_pivot_positions m in
+        let v_after_elim = List.fold_left (
+            fun acc (row_idx, pivot_position) ->
+              let v_at_piv = V.nth acc pivot_position in 
+              if v_at_piv =: A.zero then 
+                acc
+              else
+                let piv_row = List.nth m row_idx in
+                sub_scaled_row acc piv_row v_at_piv
+          ) v pivot_positions
+        in 
+        match V.find_first_non_zero v_after_elim with (* now we check for contradictions and finally insert v *)
+        | None -> Some m (* v is zero vector and was therefore already covered by m *)
+        | Some (idx, value) -> 
+          if idx = (num_cols m - 1) then 
+            None
+          else
+            let normalized_v = V.map_f_preserves_zero (fun x -> x /: value) v_after_elim in
+            Some (insert_v_according_to_piv m normalized_v idx pivot_positions)
+
+    let rref_vec m v = Timing.wrap "rref_vec" (rref_vec m) v
+
+    (* This should yield the same result as appending m2 to m1, normalizing and removing zero rows. However, it is usually faster. *)
+    (* Both input matrices are assumed to be in rref form *)
+    let rref_matrix (m1 : t) (m2 : t) =
+      let big_m, small_m = if num_rows m1 > num_rows m2 then m1, m2 else m2, m1 in
+      fst @@ List.fold_while (fun acc _ -> Option.is_some acc) 
+        (fun acc_big_m small -> rref_vec (Option.get acc_big_m) small ) (Some big_m) small_m (* TODO: pivot_positions are recalculated at each step, but since they need to be adjusted after each step it might not make sense to keep track of them here.*)
+
+    let rref_matrix m1 m2 = Timing.wrap "rref_matrix" (rref_matrix m1) m2
+
+    (* Performs a partial rref reduction to check if concatenating both matrices and afterwards normalizing them would yield a matrix <> m2 *)
+    (* Both input matrices must be in rref form! *)
+    let is_covered_by m1 m2 =
+      let rec is_linearly_independent_rref v m = 
+        let pivot_opt = V.find_first_non_zero v in
+        match pivot_opt with
+        | None -> false (* When we found no pivot, the vector is already A.zero. *)
+        | Some (pivot_id, pivot) ->
+          let m' = List.drop_while (fun v2 -> 
+              match V.find_first_non_zero v2 with
+              | None -> true  (* In this case, m2 only has zero rows after that *)
+              | Some (idx', _) -> idx' < pivot_id
+            ) m in 
+          match m' with
+          | [] -> not @@ V.is_zero_vec v
+          | x::xs ->
+            let new_v = V.map2_f_preserves_zero (fun v1 v2 -> v1 -: (pivot *: v2)) v x in
+            is_linearly_independent_rref new_v m'
+      in
+      if num_rows m1 > num_rows m2 then false else
+        let rec is_covered_by_helper m1 m2 = 
+          match m1 with 
+          | [] -> true
+          | v1::vs1 -> 
+            let first_non_zero = V.find_first_non_zero v1 in
+            match first_non_zero with
+            | None -> true  (* vs1 must also be zero-vectors because of rref *)
+            | Some (idx, _) -> 
+              let linearly_indep = is_linearly_independent_rref v1 m2 in
+              if linearly_indep then false else is_covered_by_helper vs1 m2
+        in is_covered_by_helper m1 m2
+
+    let is_covered_by m1 m2 = Timing.wrap "is_covered_by" (is_covered_by m1) m2
+
+    let find_opt f m =
+      List.find_opt f m
+
+    let map2 f m v =
+      let vector_length = V.length v in
+      List.mapi (fun index row -> if index < vector_length then f row (V.nth v index) else row ) m
+  end
\ No newline at end of file
diff --git a/src/cdomains/affineEquality/sparseImplementation/sparseVector.ml b/src/cdomains/affineEquality/sparseImplementation/sparseVector.ml
new file mode 100644
index 0000000000..a3d499f922
--- /dev/null
+++ b/src/cdomains/affineEquality/sparseImplementation/sparseVector.ml
@@ -0,0 +1,314 @@
+open AbstractVector
+open RatOps
+open ConvenienceOps
+
+open BatList
+open Batteries
+
+module List = BatList
+module Array = BatArray
+
+module SparseVector: AbstractVector =
+  functor (A: RatOps) ->
+  struct
+    include ConvenienceOps (A)
+
+    type t = {
+      entries: (int * A.t) list ;
+      len: int
+    }[@@deriving eq, ord, hash]
+
+    let copy v = v 
+
+    let of_list l = 
+      let entries' = List.rev @@ List.fold_lefti (fun acc i x -> if x <> A.zero then (i, x) :: acc else acc) [] l
+      in {entries = entries'; len = List.length l}
+
+    let of_array a =
+      let len' = Array.length a in
+      let entries' = List.rev @@ Array.fold_lefti (fun acc i x -> if x <> A.zero then (i, x) :: acc else acc ) [] a in
+      {entries = entries'; len = len'}
+
+    let of_sparse_list col_count ls =
+      {entries = ls; len = col_count}
+
+    let to_list v = 
+      let[@tail_mod_cons] rec extend_zero_aux i v' =
+        match v', i with
+        | (xi,xv)::xs, _ -> if xi = i then xv::(extend_zero_aux (i+1) xs) else A.zero ::(extend_zero_aux (i+1) v')
+        | [], j when i < v.len -> A.zero :: (extend_zero_aux (i+1) v')
+        | [], _ -> []
+      in 
+      (extend_zero_aux 0 v.entries)
+
+    let to_array v = 
+      let vec = Array.make v.len A.zero in 
+      List.iter (fun (idx, value) -> vec.(idx) <- value) v.entries;
+      vec
+
+    let to_sparse_list v = 
+      v.entries
+
+    let show v = 
+      let t = to_list v in 
+      let rec list_str l =
+        match l with
+        | [] -> "]"
+        | x :: xs -> (A.to_string x) ^ " " ^ (list_str xs)
+      in
+      "["^list_str t^"\n"
+
+    let length v =
+      v.len
+    let compare_length_with v n = 
+      Int.compare v.len n
+
+    let zero_vec n = 
+      {entries = []; len = n}
+
+    let is_zero_vec v = (v.entries = [])
+
+    (* Checks if exactly one variable is part of the affine equality, i.e. whether there is exactly one non-zero entry, excluding the constant at the end *)
+    let is_const_vec v = 
+      match v.entries with 
+      | [] -> false
+      | (idx, _) :: (const_idx , _) :: [] when const_idx = (v.len - 1) -> true
+      | (idx, _)::[] when idx <> v.len -1 -> true 
+      | _ -> false
+
+    let nth v n = (* Note: This exception HAS TO BE THROWN! It is expected by the domain *)
+      if n >= v.len then raise (Invalid_argument "Cannot access vector element (out of bounds)")
+      else
+        let rec nth v = match v with (* List.assoc would also work, but we use the fact that v is sorted *)
+          | [] -> A.zero
+          | (col_idx, value) :: xs when col_idx > n -> A.zero
+          | (col_idx, value) :: xs when col_idx = n -> value
+          | (col_idx, value) :: xs -> nth xs 
+        in nth v.entries
+
+    let set_nth v n num = (* TODO: Optimize! *)
+      if n >= v.len then failwith "Out of bounds" 
+      else
+        let rec set_nth_helper vec acc =
+          match vec with
+          | [] -> if num <>: A.zero then List.rev ((n, num) :: acc) else List.rev acc
+          | (idx, value) :: xs when n = idx -> if num <>: A.zero then List.rev_append ((idx, num) :: acc) xs else List.rev_append acc xs
+          | (idx, value) :: xs when n < idx -> if num <>: A.zero then List.rev_append ((n, num) :: acc) vec else List.rev_append acc vec
+          | x :: xs -> set_nth_helper xs (x :: acc)
+        in
+        {entries = set_nth_helper v.entries []; len = v.len}
+
+    let insert_val_at n new_val v = 
+      if n > v.len then failwith "n too large" else (* Does this happen? Otherwise we can omit this comparison, here right now to be semantically equivalent *)
+        let entries' = List.rev @@ List.fold_left (fun acc (idx, value) -> 
+            if idx < n then (idx, value) :: acc 
+            else (if new_val <>: A.zero then (idx + 1, value) :: (n, new_val) :: acc else (idx + 1, value) :: acc)
+          ) [] v.entries in
+        {entries = entries'; len = v.len + 1}
+
+    (* Note that idx is assumed to be sorted. *)
+    let insert_zero_at_indices v idx count = 
+      let rec add_indices_helper vec idx added_count acc = 
+        match vec, idx with 
+        | [], _ -> List.rev acc (* inserting at the end only means changing the dimension *)
+        | ((col_idx, value) :: xs), [] -> add_indices_helper xs [] added_count ((col_idx + added_count, value) :: acc)
+        | ((col_idx, value) :: xs), ((i, count) :: ys) when i = col_idx -> add_indices_helper vec ys (added_count + count) acc
+        | ((col_idx, value) :: xs), ((i, count) :: ys) when i < col_idx -> add_indices_helper vec ys (added_count + count) acc
+        | ((col_idx, value) :: xs), ((i, count) :: ys) -> add_indices_helper xs idx added_count ((col_idx + added_count, value) :: acc)
+      in
+      {entries = add_indices_helper v.entries idx 0 []; len = v.len + count}
+
+    let remove_nth v n =
+      if n >= v.len then failwith "Out of bounds"
+      else
+        let new_entries = List.filter_map (fun (col_idx, value) ->
+            if col_idx = n then None
+            else if col_idx > n 
+            then Some (col_idx - 1, value)
+            else Some (col_idx, value)
+          ) v.entries in
+        {entries = new_entries; len = v.len - 1}
+
+    (* Note: It is assumed and necessary here that idx is sorted!!! *)
+    let remove_at_indices v idx = 
+      let rec remove_indices_helper vec idx deleted_count acc = 
+        match vec, idx with 
+        | [], [] -> List.rev acc
+        | [], (y :: ys) when deleted_count >= v.len || y >= v.len -> failwith "remove at indices: no more columns to delete"
+        | [], (y :: ys) -> remove_indices_helper [] ys (deleted_count + 1) acc (* Removing zero (also in next iteration, else failwith ) *)
+        | ((col_idx, value) :: xs), [] -> remove_indices_helper xs [] deleted_count ((col_idx - deleted_count, value) :: acc)
+        | ((col_idx, value) :: xs), (y :: ys) when y = col_idx -> remove_indices_helper xs ys (deleted_count + 1) acc
+        | ((col_idx, value) :: xs), (y :: ys) when y < col_idx -> remove_indices_helper vec ys (deleted_count + 1) acc
+        | ((col_idx, value) :: xs), (y :: ys) -> remove_indices_helper xs idx deleted_count ((col_idx - deleted_count, value) :: acc)
+      in
+      {entries = remove_indices_helper v.entries idx 0 []; len = v.len - List.length idx}
+
+    let keep_vals v n = 
+      if n >= v.len then v else
+        {entries = List.take_while (fun (idx, _) -> idx < n) v.entries; len=n}
+
+    let starting_from_nth n v =
+      let entries' = List.filter_map (fun (idx, value) -> if idx < n then None else Some (idx - n, value)) v.entries in
+      {entries = entries'; len = v.len - n}
+
+    let find_opt f v = (* TODO: Do we need this? And optimize!!!*)
+      List.find_opt f (to_list v)
+
+    let findi f v = 
+      if f A.zero then  
+        fst @@ List.findi (fun i (idx, value) -> if idx > i then true else f value) v.entries (* Here fst is the iteration variable i, not the tuple idx *)
+      else
+        fst @@ List.find (fun (idx, value) -> f value) v.entries (* Here fst is the idx contained in the found tuple *)
+
+    let find2i f v v' = (* TODO: optimize! *)
+      fst @@ List.findi (fun _ (val1, val2) -> (uncurry f) (val1, val2)) (List.combine (to_list v) (to_list v'))
+
+    (* Returns optional of (index * value) where f value evaluated to true *)
+    let findi_val_opt f v =
+      let rec find_zero_or_val vec last_col_idx =
+        let f0 = f A.zero in
+        match vec with
+        | [] -> if f0 && v.len <> last_col_idx + 1 then Some (last_col_idx + 1, A.zero) else None
+        | (idx, value) :: xs -> 
+          if f0  && idx <> last_col_idx + 1 then Some (last_col_idx + 1, A.zero) 
+          else if f value then Some (idx, value) 
+          else find_zero_or_val xs idx
+      in find_zero_or_val v.entries (-1)
+
+    let findi_val_opt f v = Timing.wrap "findi_val_opt" (findi_val_opt f) v
+
+    let find_first_non_zero v =
+      if v.entries = [] then None
+      else Some (List.hd v.entries)
+
+    let find_first_non_zero v = Timing.wrap "find_first_non_zero" (find_first_non_zero) v
+
+    let exists f v  = 
+      let c = v.len in
+      let rec exists_aux at f v =
+        match v with
+        | [] -> if at = 0 then false else f A.zero
+        | (xi, xv)::xs -> if f xv then true else exists_aux (at - 1) f xs
+      in (exists_aux c f v.entries)
+
+    let exists2 f v1 v2 = (* TODO: optimize! *)
+      List.exists2 f (to_list v1) (to_list v2)
+
+    let filteri f v = (* TODO: optimize! *)
+      of_list (List.filteri f (to_list v))
+
+    let map f v = 
+      let f_zero = f A.zero in
+      let rec map2_helper acc vec i =
+        match vec with
+        | [] when i >= v.len || f_zero =: A.zero -> List.rev acc
+        | []  -> map2_helper ((i, f_zero) :: acc) [] (i  + 1)
+        | (idx, value) :: xs when idx = i -> let new_val = f value in if new_val <>: A.zero then map2_helper ((idx, new_val) :: acc) xs (i + 1) else map2_helper acc xs (i + 1)
+        | (idx, _) :: xs when idx > i -> if f_zero <>: A.zero then map2_helper ((i, f_zero) :: acc) vec (i + 1) else map2_helper acc vec (i + 1)
+        | (_, _) :: _ -> failwith "This should not happen"
+      in
+      {entries = map2_helper [] v.entries 0; len = v.len}
+
+    let map_f_preserves_zero f v = (* map for functions f such that f 0 = 0 since f won't be applied to zero values. See also map *)
+      let entries' = List.filter_map (
+          fun (idx, value) -> let new_val = f value in 
+            if new_val = A.zero then None else Some (idx, new_val)) v.entries in 
+      {entries = entries'; len = v.len}
+
+    let map_f_preserves_zero f v = Timing.wrap "map_f_preserves_zero" (map_f_preserves_zero f) v
+
+    let mapi f v  = (* TODO: optimize! *)
+      of_list (List.mapi f (to_list v))
+
+    let mapi_f_preserves_zero f v =
+      let entries' = List.filter_map (
+          fun (idx, value) -> let new_val = f idx value in 
+            if new_val = A.zero then None else Some (idx, new_val)) v.entries in 
+      {entries = entries'; len = v.len}
+
+    let map2 f v v' = (* TODO: Optimize! *)
+      if v.len <> v'.len then failwith "Unequal vector length" else 
+        of_list (List.map2 f (to_list v) (to_list v'))
+
+    (* map for functions f such that f 0 0 = 0 since f won't be applied to if both values are zero. See also map *)
+    let map2_f_preserves_zero f v1 v2 =
+      let f_rem_zero acc idx e1 e2 =
+        let r = f e1 e2 in 
+        if r =: A.zero then acc else (idx, r)::acc
+      in  
+      let rec aux acc v1 v2 =
+        match v1, v2 with 
+        | [], [] -> acc 
+        | [], (yidx, yval)::ys -> aux (f_rem_zero acc yidx A.zero yval) [] ys
+        | (xidx, xval)::xs, [] -> aux (f_rem_zero acc xidx xval A.zero) xs []
+        | (xidx, xval)::xs, (yidx, yval)::ys -> 
+          match xidx - yidx with
+          | d when d < 0 -> aux (f_rem_zero acc xidx xval A.zero) xs v2
+          | d when d > 0 -> aux (f_rem_zero acc yidx A.zero yval) v1 ys
+          | _            -> aux (f_rem_zero acc xidx xval yval) xs ys
+      in
+      if v1.len <> v2.len then raise (Invalid_argument "Unequal lengths") else 
+        {entries = List.rev (aux [] v1.entries v2.entries); len = v1.len}
+
+    let map2_f_preserves_zero f v1 v2 = Timing.wrap "map2_f_preserves_zero" (map2_f_preserves_zero f v1) v2
+
+    let map2i f v1 v2 = (* TODO: optimize! *)
+      if v1.len <> v2.len then raise (Invalid_argument "Unequal lengths") else
+        (*of_list (List.map2i f (to_list v) (to_list v'))*)
+        let f_rem_zero idx acc e1 e2 =
+          let r = f idx e1 e2 in 
+          if r =: A.zero then acc else (idx, r)::acc
+        in  
+        let rec aux acc vec1 vec2 i =
+          match vec1, vec2 with 
+          | [], [] when i = v1.len -> acc 
+          | [], [] -> aux (f_rem_zero i acc A.zero A.zero) [] [] (i + 1)
+          | [], (yidx, yval)::ys when i = yidx -> aux (f_rem_zero i acc A.zero yval) [] ys (i + 1)
+          | (xidx, xval)::xs, [] when i = xidx -> aux (f_rem_zero i acc xval A.zero) xs [] (i + 1)
+          | [], (_, _)::_ | (_, _)::_, [] -> aux (f_rem_zero i acc A.zero A.zero) vec1 vec2 (i + 1) (* When one vec is not zero_vec, but has implicit zeroes at front *)
+          | (xidx, xval)::xs, (yidx, yval)::ys -> 
+            if xidx <> i && yidx <> i then aux (f_rem_zero i acc A.zero A.zero) vec1 vec2 (i+1) (* When both vectors have implicit zeroes at front *)
+            else
+              match xidx - yidx with (* Here at least one of the idx is i, which is the smaller one *)
+              | d when d < 0 -> aux (f_rem_zero i acc xval A.zero) xs vec2 (i + 1)
+              | d when d > 0 -> aux (f_rem_zero i acc A.zero yval) vec1 ys (i + 1)
+              | _            -> aux (f_rem_zero i acc xval yval) xs ys (i + 1)
+        in
+        {entries = List.rev (aux [] v1.entries v2.entries 0); len = v1.len}
+
+    let map2i_f_preserves_zero f v v' = failwith "TODO"
+
+    let fold_left_f_preserves_zero f acc v =
+      List.fold_left (fun acc (_, value) -> f acc value) acc v.entries
+
+    let fold_left2_f_preserves_zero f acc v v' =  
+      let rec aux acc v1 v2 =
+        match v1, v2 with 
+        | [], [] -> acc 
+        | [], (yidx, yval)::ys -> aux (f acc A.zero yval) [] ys
+        | (xidx, xval)::xs, [] -> aux (f acc xval A.zero) xs []
+        | (xidx, xval)::xs, (yidx, yval)::ys -> 
+          match xidx - yidx with
+          | d when d < 0 -> aux (f acc xval A.zero) xs v2
+          | d when d > 0 -> aux (f acc A.zero yval) v1 ys
+          | _            -> aux (f acc xval yval) xs ys
+      in
+      if v.len <> v'.len then raise (Invalid_argument "Unequal lengths") else 
+        (aux acc v.entries v'.entries)
+
+    let apply_with_c f c v  = (* TODO: optimize! *)
+      of_list @@ List.map (fun value -> f value c) (to_list v)
+
+    let apply_with_c_f_preserves_zero f c v =
+      let entries' = List.filter_map (fun (idx, value) -> let new_val = f value c in if new_val =: A.zero then None else Some (idx, new_val)) v.entries in
+      {entries = entries'; len = v.len}
+
+    let rev v = 
+      let entries' = List.rev @@ List.map (fun (idx, value) -> (v.len - 1 - idx, value)) v.entries in 
+      {entries = entries'; len = v.len}
+
+    let append v v' = 
+      let entries' = v.entries @ List.map (fun (idx, value) -> (idx + v.len), value) v'.entries in
+      {entries = entries'; len = v.len + v'.len}
+
+  end 
diff --git a/src/cdomains/affineEquality/vector.ml b/src/cdomains/affineEquality/vector.ml
new file mode 100644
index 0000000000..799df35c0e
--- /dev/null
+++ b/src/cdomains/affineEquality/vector.ml
@@ -0,0 +1,94 @@
+(** High-level abstraction of a vector. *)
+module type Vector =
+sig
+  type num
+  type t [@@deriving eq, ord, hash]
+
+  val show: t -> string
+
+  val copy: t -> t
+
+  val of_list: num list -> t
+
+  val of_array: num array -> t
+
+  val of_sparse_list: int -> (int * num) list -> t
+
+  val to_list: t -> num list
+
+  val to_array: t -> num array
+
+  val to_sparse_list: t -> (int * num) list
+
+  val length: t -> int
+
+  val compare_length_with: t -> int -> int
+
+  val zero_vec: int -> t
+
+  val is_zero_vec: t -> bool
+
+  val is_const_vec: t -> bool
+
+  val nth: t -> int -> num
+
+  val set_nth: t -> int -> num ->  t
+
+  val insert_val_at: int -> num ->  t ->  t
+
+  val insert_zero_at_indices: t -> (int * int) list -> int -> t
+
+  val remove_nth: t -> int ->  t
+
+  val remove_at_indices: t -> int list -> t
+
+  val keep_vals: t -> int ->  t
+
+  (* Returns the part of the vector starting from index n*)
+  val starting_from_nth : int -> t -> t
+
+  val find_opt: (num -> bool) -> t -> num Option.t
+
+  val findi: (num -> bool) ->  t -> int
+
+  val find2i: (num -> num -> bool) -> t -> t -> int
+
+  (* Returns optional tuple of position and value which was found*)
+  val findi_val_opt: (num -> bool) ->  t -> (int * num) Option.t
+
+  val find_first_non_zero : t -> (int * num) option
+
+  val exists: (num -> bool) -> t -> bool
+
+  val exists2: (num -> num -> bool) -> t -> t -> bool
+
+  val filteri: (int -> num -> bool) -> t -> t
+
+  val map: (num -> num) -> t -> t
+
+  val map_f_preserves_zero: (num -> num) -> t -> t
+
+  val mapi: (int -> num -> num) -> t -> t
+
+  val mapi_f_preserves_zero: (int -> num -> num) -> t -> t
+
+  val map2: (num -> num -> num) -> t -> t -> t
+
+  val map2_f_preserves_zero: (num -> num -> num) -> t ->  t -> t
+
+  val map2i: (int -> num -> num -> num) -> t -> t -> t
+
+  val map2i_f_preserves_zero: (int -> num -> num -> num) -> t -> t -> t
+
+  val fold_left_f_preserves_zero: ('acc -> num -> 'acc) -> 'acc -> t -> 'acc
+
+  val fold_left2_f_preserves_zero: ('acc -> num -> num -> 'acc) -> 'acc -> t -> t -> 'acc
+
+  val apply_with_c: (num -> num -> num) -> num ->  t ->  t
+
+  val apply_with_c_f_preserves_zero: (num -> num -> num) -> num ->  t ->  t
+
+  val rev: t -> t
+
+  val append: t -> t -> t
+end
\ No newline at end of file
diff --git a/src/cdomains/apron/affineEqualityDomain.apron.ml b/src/cdomains/apron/affineEqualityDomain.apron.ml
index 06268130b2..2e9676382f 100644
--- a/src/cdomains/apron/affineEqualityDomain.apron.ml
+++ b/src/cdomains/apron/affineEqualityDomain.apron.ml
@@ -8,20 +8,27 @@
 
 open GoblintCil
 open Pretty
+
 module M = Messages
 open GobApron
-open VectorMatrix
+
+open ConvenienceOps
+open AbstractVector
+open AbstractMatrix
+
+open Batteries
+module Array = Batteries.Array
 
 module Mpqf = SharedFunctions.Mpqf
 
+
 module AffineEqualityMatrix (Vec: AbstractVector) (Mx: AbstractMatrix) =
 struct
   include Mx(Mpqf) (Vec)
   let dim_add (ch: Apron.Dim.change) m =
-    Array.modifyi (+) ch.dim;
-    add_empty_columns m ch.dim
+    add_empty_columns m ch.dim 
 
-  let dim_add ch m = timing_wrap "dim add" (dim_add ch) m
+  let dim_add ch m = Timing.wrap "dim add" (dim_add ch) m
 
 
   let dim_remove (ch: Apron.Dim.change) m ~del =
@@ -29,10 +36,10 @@ struct
       m
     else (
       Array.modifyi (+) ch.dim;
-      let m' = if not del then let m = copy m in Array.fold_left (fun y x -> reduce_col_with y x; y) m ch.dim else m in
+      let m' = if not del then Array.fold_left (fun y x -> reduce_col y x) m ch.dim else m in
       remove_zero_rows @@ del_cols m' ch.dim)
 
-  let dim_remove ch m ~del = VectorMatrix.timing_wrap "dim remove" (fun del -> dim_remove ch m ~del:del) del
+  let dim_remove ch m ~del = Timing.wrap "dim remove" (fun del -> dim_remove ch m ~del:del) del
 end
 
 (** It defines the type t of the affine equality domain (a struct that contains an optional matrix and an apron environment) and provides the functions needed for handling variables (which are defined by RelationDomain.D2) such as add_vars remove_vars.
@@ -59,9 +66,8 @@ struct
     let open Apron.Texpr1 in
     let exception NotLinear in
     let zero_vec = Vector.zero_vec @@ Environment.size t.env + 1 in
-    let neg v = Vector.map_with Mpqf.neg v; v in
-    let is_const_vec v = Vector.compare_length_with (Vector.filteri (fun i x -> (*Inefficient*)
-        Vector.compare_length_with v (i + 1) > 0 && x <>: Mpqf.zero) v) 1 = 0
+    let neg v = Vector.map_f_preserves_zero Mpqf.neg v in
+    let is_const_vec = Vector.is_const_vec 
     in
     let rec convert_texpr = function
       (*If x is a constant, replace it with its const. val. immediately*)
@@ -72,36 +78,35 @@ struct
           | Scalar Mpqf x -> x
           | Scalar Mpfrf x -> Mpfr.to_mpq x
         in
-        Vector.set_val zero_vec ((Vector.length zero_vec) - 1) (of_union x)
+        Vector.set_nth zero_vec ((Vector.length zero_vec) - 1) (of_union x)
       | Var x ->
-        let zero_vec_cp = Vector.copy zero_vec in
-        let entry_only v = Vector.set_val_with v (Environment.dim_of_var t.env x) Mpqf.one; v in
+        let entry_only v = Vector.set_nth v (Environment.dim_of_var t.env x) Mpqf.one in
         begin match t.d with
           | Some m ->
             let row = Matrix.find_opt (fun r -> Vector.nth r (Environment.dim_of_var t.env x) =: Mpqf.one) m in
             begin match row with
               | Some v when is_const_vec v ->
-                Vector.set_val_with zero_vec_cp ((Vector.length zero_vec) - 1) (Vector.nth v (Vector.length v - 1)); zero_vec_cp
-              | _ -> entry_only zero_vec_cp
+                Vector.set_nth zero_vec ((Vector.length zero_vec) - 1) (Vector.nth v (Vector.length v - 1))
+              | _ -> entry_only zero_vec
             end
-          | None -> entry_only zero_vec_cp end
+          | None -> entry_only zero_vec end
       | Unop (Neg, e, _, _) -> neg @@ convert_texpr e
       | Unop (Cast, e, _, _) -> convert_texpr e (*Ignore since casts in apron are used for floating point nums and rounding in contrast to CIL casts*)
       | Unop (Sqrt, e, _, _) -> raise NotLinear
       | Binop (Add, e1, e2, _, _) ->
         let v1 = convert_texpr e1 in
         let v2 = convert_texpr e2 in
-        Vector.map2_with (+:) v1 v2; v1
+        Vector.map2_f_preserves_zero (+:) v1 v2
       | Binop (Sub, e1, e2, _, _) ->
         let v1 = convert_texpr e1 in
         let v2 = convert_texpr e2 in
-        Vector.map2_with (+:) v1 (neg @@ v2); v1
+        Vector.map2_f_preserves_zero (+:) v1 (neg @@ v2)
       | Binop (Mul, e1, e2, _, _) ->
         let v1 = convert_texpr e1 in
         let v2 = convert_texpr e2 in
         begin match to_constant_opt v1, to_constant_opt v2 with
-          | _, Some c -> Vector.apply_with_c_with ( *:) c v1; v1
-          | Some c, _ -> Vector.apply_with_c_with ( *:) c v2; v2
+          | _, Some c -> Vector.apply_with_c_f_preserves_zero ( *:) c v1
+          | Some c, _ -> Vector.apply_with_c_f_preserves_zero ( *:) c v2
           | _, _ -> raise NotLinear
         end
       | Binop _ -> raise NotLinear
@@ -110,7 +115,7 @@ struct
       Some (convert_texpr texp)
     with NotLinear -> None
 
-  let get_coeff_vec t texp = timing_wrap "coeff_vec" (get_coeff_vec t) texp
+  let get_coeff_vec t texp = Timing.wrap "coeff_vec" (get_coeff_vec t) texp
 end
 
 (** As it is specifically used for the new affine equality domain, it can only provide bounds if the expression contains known constants only and in that case, min and max are the same. *)
@@ -137,7 +142,7 @@ struct
     res
 
 
-  let bound_texpr d texpr1 = timing_wrap "bounds calculation" (bound_texpr d) texpr1
+  let bound_texpr d texpr1 = Timing.wrap "bounds calculation" (bound_texpr d) texpr1
 end
 
 module D(Vc: AbstractVector) (Mx: AbstractMatrix) =
@@ -192,7 +197,7 @@ struct
       in
       let res = (String.concat "" @@ Array.to_list @@ Array.map dim_to_str vars)
                 ^ (const_to_str arr.(Array.length arr - 1)) ^ "=0" in
-      if String.starts_with res ~prefix:"+" then
+      if String.starts_with res "+" then
         Str.string_after res 1
       else
         res
@@ -240,11 +245,9 @@ struct
       else if is_top_env t2 then
         {d = Some (dim_add (Environment.dimchange t1.env sup_env) m1); env = sup_env}
       else
-        let rref_matr = Matrix.rref_matrix_with (Matrix.copy m1) (Matrix.copy m2) in
-        if Option.is_none rref_matr then
-          bot ()
-        else
-          {d = rref_matr; env = sup_env}
+        match Matrix.rref_matrix m1 m2 with
+        | None -> bot ()
+        | rref_matr -> {d = rref_matr; env = sup_env}
 
 
   let meet t1 t2 =
@@ -252,7 +255,7 @@ struct
     if M.tracing then M.tracel "meet" "meet a: %s b: %s -> %s " (show t1) (show t2) (show res) ;
     res
 
-  let meet t1 t2 = timing_wrap "meet" (meet t1) t2
+  let meet t1 t2 = Timing.wrap "meet" (meet t1) t2
 
   let leq t1 t2 =
     let env_comp = Environment.cmp t1.env t2.env in (* Apron's Environment.cmp has defined return values. *)
@@ -271,7 +274,7 @@ struct
       let m1' = if env_comp = 0 then m1 else dim_add (Environment.dimchange t1.env t2.env) m1 in
       Matrix.is_covered_by m2 m1'
 
-  let leq a b = timing_wrap "leq" (leq a) b
+  let leq a b = Timing.wrap "leq" (leq a) b
 
   let leq t1 t2 =
     let res = leq t1 t2 in
@@ -283,29 +286,26 @@ struct
       if s >= Matrix.num_cols a then a else
         let case_two a r col_b =
           let a_r = Matrix.get_row a r in
-          Matrix.map2i_with (fun i x y -> if i < r then
-                                Vector.map2_with (fun u j -> u +: y *: j) x a_r; x) a col_b;
+          let a = Matrix.map2i (fun i x y -> if i < r then Vector.map2_f_preserves_zero (fun u j -> u +: y *: j) x a_r else x) a col_b; in
           Matrix.remove_row a r
         in
         let case_three a b col_a col_b max =
-          let col_a, col_b = Vector.copy col_a, Vector.copy col_b in
           let col_a, col_b = Vector.keep_vals col_a max, Vector.keep_vals col_b max in
           if Vector.equal col_a col_b then
             (a, b, max)
           else
             (
-              Vector.rev_with col_a;
-              Vector.rev_with col_b;
+              let col_a = Vector.rev col_a in
+              let col_b = Vector.rev col_b in
               let i = Vector.find2i (<>:) col_a col_b in
               let (x, y) = Vector.nth col_a i, Vector.nth col_b i in
               let r, diff = Vector.length col_a - (i + 1), x -: y  in
               let a_r, b_r = Matrix.get_row a r, Matrix.get_row b r in
-              Vector.map2_with (-:) col_a col_b;
-              Vector.rev_with col_a;
+              let col_a = Vector.map2_f_preserves_zero (-:) col_a col_b in
+              let col_a = Vector.rev col_a in
               let multiply_by_t m t =
-                Matrix.map2i_with (fun i' x c -> if i' <= max then (let beta = c /: diff in
-                                                                    Vector.map2_with (fun u j -> u -: (beta *: j)) x t); x) m col_a;
-                m
+                Matrix.map2i (fun i' x c -> if i' <= max then (let beta = c /: diff in Vector.map2_f_preserves_zero (fun u j -> u -: (beta *: j)) x t) else x) m col_a;
+
               in
               Matrix.remove_row (multiply_by_t a a_r) r, Matrix.remove_row (multiply_by_t b b_r) r, (max - 1)
             )
@@ -336,11 +336,11 @@ struct
         let sup_env = Environment.lce a.env b.env in
         let mod_x = dim_add (Environment.dimchange a.env sup_env) x in
         let mod_y = dim_add (Environment.dimchange b.env sup_env) y in
-        {d = Some (lin_disjunc 0 0 (Matrix.copy mod_x) (Matrix.copy mod_y)); env = sup_env}
+        {d = Some (lin_disjunc 0 0 mod_x mod_y); env = sup_env}
       | x, y when Matrix.equal x y -> {d = Some x; env = a.env}
-      | x, y  -> {d = Some(lin_disjunc 0 0 (Matrix.copy x) (Matrix.copy y)); env = a.env}
+      | x, y  -> {d = Some(lin_disjunc 0 0 x y); env = a.env}
 
-  let join a b = timing_wrap "join" (join a) b
+  let join a b = Timing.wrap "join" (join a) b
 
   let join a b =
     let res = join a b in
@@ -358,65 +358,65 @@ struct
   let pretty_diff () (x, y) =
     dprintf "%s: %a not leq %a" (name ()) pretty x pretty y
 
-  let remove_rels_with_var x var env inplace =
-    let j0 = Environment.dim_of_var env var in
-    if inplace then
-      (Matrix.reduce_col_with x j0; x)
-    else
-      Matrix.reduce_col x j0
+  let remove_rels_with_var x var env =
+    let j0 = Environment.dim_of_var env var in Matrix.reduce_col x j0
 
-  let remove_rels_with_var x var env inplace = timing_wrap "remove_rels_with_var" (remove_rels_with_var x var env) inplace
+  let remove_rels_with_var x var env = Timing.wrap "remove_rels_with_var" remove_rels_with_var x var env 
 
   let forget_vars t vars =
     if is_bot t || is_top_env t || vars = [] then
       t
     else
       let m = Option.get t.d in
-      let rem_from m = List.fold_left (fun m' x -> remove_rels_with_var m' x t.env true) m vars in
-      {d = Some (Matrix.remove_zero_rows @@ rem_from (Matrix.copy m)); env = t.env}
+      let rem_from m = List.fold_left (fun m' x -> remove_rels_with_var m' x t.env) m vars in
+      {d = Some (Matrix.remove_zero_rows @@ rem_from  m); env = t.env}
 
   let forget_vars t vars =
     let res = forget_vars t vars in
     if M.tracing then M.tracel "ops" "forget_vars %s -> %s" (show t) (show res);
     res
 
-  let forget_vars t vars = timing_wrap "forget_vars" (forget_vars t) vars
+  let forget_vars t vars = Timing.wrap "forget_vars" (forget_vars t) vars
 
   let assign_texpr (t: VarManagement(Vc)(Mx).t) var texp =
     let assign_invertible_rels x var b env =
       let j0 = Environment.dim_of_var env var in
       let a_j0 = Matrix.get_col x j0  in (*Corresponds to Axj0*)
       let b0 = Vector.nth b j0 in
-      Vector.apply_with_c_with (/:) b0 a_j0; (*Corresponds to Axj0/Bj0*)
-      let recalc_entries m rd_a = Matrix.map2_with (fun x y -> Vector.map2i_with (fun j z d ->
+      let a_j0 = Vector.apply_with_c_f_preserves_zero (/:) b0 a_j0 in (*Corresponds to Axj0/Bj0*)
+      let recalc_entries m rd_a = Matrix.map2 (fun x y -> Vector.map2i (fun j z d ->
           if j = j0 then y
           else if Vector.compare_length_with b (j + 1) > 0 then z -: y *: d
-          else z +: y *: d) x b; x) m rd_a
+          else z +: y *: d) x b) m rd_a
       in
-      recalc_entries x a_j0;
-      if Matrix.normalize_with x then {d = Some x; env = env} else bot ()
+      let x = recalc_entries x a_j0 in
+      match Matrix.normalize x with 
+      | None -> bot ()
+      | some_normalized_matrix -> {d = some_normalized_matrix; env = env}
     in
-    let assign_invertible_rels x var b env = timing_wrap "assign_invertible" (assign_invertible_rels x var b) env in
+    let assign_invertible_rels x var b env = Timing.wrap "assign_invertible" (assign_invertible_rels x var b) env in
     let assign_uninvertible_rel x var b env =
       let b_length = Vector.length b in
-      Vector.mapi_with (fun i z -> if i < b_length - 1 then Mpqf.neg z else z) b;
-      Vector.set_val_with b (Environment.dim_of_var env var) Mpqf.one;
-      let opt_m = Matrix.rref_vec_with x b in
-      if Option.is_none opt_m then bot () else
-        {d = opt_m; env = env}
+      let b = Vector.mapi_f_preserves_zero (fun i z -> if i < b_length - 1 then Mpqf.neg z else z) b in
+      let b = Vector.set_nth b (Environment.dim_of_var env var) Mpqf.one in
+      match Matrix.rref_vec x b with
+      | None -> bot ()
+      | some_matrix -> {d = some_matrix; env = env}
     in
-    (* let assign_uninvertible_rel x var b env = timing_wrap "assign_uninvertible" (assign_uninvertible_rel x var b) env in *)
+    (* let assign_uninvertible_rel x var b env = Timing.wrap "assign_uninvertible" (assign_uninvertible_rel x var b) env in *)
     let is_invertible v = Vector.nth v @@ Environment.dim_of_var t.env var <>: Mpqf.zero
-    in let affineEq_vec = get_coeff_vec t texp
-    in if is_bot t then t else let m = Option.get t.d in
+    in let affineEq_vec = get_coeff_vec t texp in
+    if is_bot t then t else let m = Option.get t.d in
       match affineEq_vec with
-      | Some v when is_top_env t -> if is_invertible v then t else assign_uninvertible_rel m var v t.env
-      | Some v -> if is_invertible v then let t' = assign_invertible_rels (Matrix.copy m) var v t.env in {d = t'.d; env = t'.env}
-        else let new_m = Matrix.remove_zero_rows @@ remove_rels_with_var m var t.env false
+      | Some v when is_top_env t -> 
+        if is_invertible v then t else assign_uninvertible_rel m var v t.env
+      | Some v -> 
+        if is_invertible v then let t' = assign_invertible_rels m var v t.env in {d = t'.d; env = t'.env}
+        else let new_m = Matrix.remove_zero_rows @@ remove_rels_with_var m var t.env 
           in assign_uninvertible_rel new_m var v t.env
-      | None -> {d = Some (Matrix.remove_zero_rows @@ remove_rels_with_var m var t.env false); env = t.env}
+      | None -> {d = Some (Matrix.remove_zero_rows @@ remove_rels_with_var m var t.env); env = t.env}
 
-  let assign_texpr t var texp = timing_wrap "assign_texpr" (assign_texpr t var) texp
+  let assign_texpr t var texp = Timing.wrap "assign_texpr" (assign_texpr t var) texp
 
   let assign_exp ask (t: VarManagement(Vc)(Mx).t) var exp (no_ov: bool Lazy.t) =
     let t = if not @@ Environment.mem_var t.env var then add_vars t [var] else t in
@@ -449,20 +449,18 @@ struct
     let t_primed = add_vars t primed_vars in
     let multi_t = List.fold_left2 (fun t' v_prime (_,v') -> assign_var t' v_prime v') t_primed primed_vars vv's in
     match multi_t.d with
-    | Some m when not @@ is_top_env multi_t ->
+    | Some m when not @@ is_top_env multi_t -> 
       let replace_col m x y =
         let dim_x, dim_y = Environment.dim_of_var multi_t.env x, Environment.dim_of_var multi_t.env y in
         let col_x = Matrix.get_col m dim_x in
-        Matrix.set_col_with m col_x dim_y
+        Matrix.set_col m col_x dim_y
       in
-      let m_cp = Matrix.copy m in
-      let switched_m = List.fold_left2 replace_col m_cp primed_vars assigned_vars in
+      let switched_m = List.fold_left2 replace_col m primed_vars assigned_vars in
       let res = drop_vars {d = Some switched_m; env = multi_t.env} primed_vars ~del:true in
       let x = Option.get res.d in
-      if Matrix.normalize_with x then
-        {d = Some x; env = res.env}
-      else
-        bot ()
+      (match Matrix.normalize x with 
+       | None -> bot ()
+       | some_matrix -> {d = some_matrix; env = res.env})
     | _ -> t
 
   let assign_var_parallel t vv's =
@@ -470,7 +468,7 @@ struct
     if M.tracing then M.tracel "ops" "assign_var parallel: %s -> %s " (show t) (show res);
     res
 
-  let assign_var_parallel t vv's = timing_wrap "var_parallel" (assign_var_parallel t) vv's
+  let assign_var_parallel t vv's = Timing.wrap "var_parallel" (assign_var_parallel t) vv's
 
   let assign_var_parallel_with t vv's =
     let t' = assign_var_parallel t vv's in
@@ -500,7 +498,7 @@ struct
     if M.tracing then M.tracel "ops" "Substitute_expr t: \n %s \n var: %a \n exp: %a \n -> \n %s" (show t) Var.pretty var d_exp exp (show res);
     res
 
-  let substitute_exp ask t var exp no_ov = timing_wrap "substitution" (substitute_exp ask t var exp) no_ov
+  let substitute_exp ask t var exp no_ov = Timing.wrap "substitution" (substitute_exp ask t var exp) no_ov
 
   (** Assert a constraint expression.
 
@@ -514,13 +512,13 @@ struct
     let meet_vec e =
       (* Flip the sign of the const. val in coeff vec *)
       let coeff = Vector.nth e (Vector.length e - 1) in
-      Vector.set_val_with e (Vector.length e - 1) (Mpqf.neg coeff);
+      let e = Vector.set_nth e (Vector.length e - 1) (Mpqf.neg coeff) in
       if is_bot t then
         bot ()
       else
-        let opt_m = Matrix.rref_vec_with (Matrix.copy @@ Option.get t.d) e in
-        if Option.is_none opt_m then bot () else {d = opt_m; env = t.env}
-
+        match Matrix.rref_vec (Option.get t.d) e with
+        | None -> bot ()
+        | some_matrix -> {d = some_matrix; env = t.env}
     in
     match get_coeff_vec t (Texpr1.to_expr @@ Tcons1.get_texpr1 tcons) with
     | Some v ->
@@ -545,7 +543,7 @@ struct
       end
     | None -> t
 
-  let meet_tcons t tcons expr = timing_wrap "meet_tcons" (meet_tcons t tcons) expr
+  let meet_tcons t tcons expr = Timing.wrap "meet_tcons" (meet_tcons t tcons) expr
 
   let unify a b =
     meet a b
@@ -561,7 +559,7 @@ struct
     | tcons1 -> meet_tcons ask d tcons1 e
     | exception Convert.Unsupported_CilExp _ -> d
 
-  let assert_constraint ask d e negate no_ov = timing_wrap "assert_constraint" (assert_constraint ask d e negate) no_ov
+  let assert_constraint ask d e negate no_ov = Timing.wrap "assert_constraint" (assert_constraint ask d e negate) no_ov
 
   let relift t = t
 
diff --git a/src/cdomains/apron/linearTwoVarEqualityDomain.apron.ml b/src/cdomains/apron/linearTwoVarEqualityDomain.apron.ml
index c1ca3661a5..72f24610c6 100644
--- a/src/cdomains/apron/linearTwoVarEqualityDomain.apron.ml
+++ b/src/cdomains/apron/linearTwoVarEqualityDomain.apron.ml
@@ -12,7 +12,8 @@ open GoblintCil
 open Pretty
 module M = Messages
 open GobApron
-open VectorMatrix
+
+open ConvenienceOps
 
 module Mpqf = SharedFunctions.Mpqf
 
@@ -173,7 +174,7 @@ module EqualitiesConjunction = struct
   let dim_add (ch: Apron.Dim.change) m =
     modify_variables_in_domain m ch.dim (+)
 
-  let dim_add ch m = timing_wrap "dim add" (dim_add ch) m
+  let dim_add ch m = Timing.wrap "dim add" (dim_add ch) m
 
   let dim_remove (ch: Apron.Dim.change) m =
     if Array.length ch.dim = 0 || is_empty m then
@@ -184,7 +185,7 @@ module EqualitiesConjunction = struct
       let m' = Array.fold_lefti (fun y i x -> forget_variable y (x)) m cpy in  (* clear m' from relations concerning ch.dim *)
       modify_variables_in_domain m' cpy (-))
 
-  let dim_remove ch m = VectorMatrix.timing_wrap "dim remove" (fun m -> dim_remove ch m) m
+  let dim_remove ch m = Timing.wrap "dim remove" (fun m -> dim_remove ch m) m
 
   let dim_remove ch m ~del = let res = dim_remove ch m in if M.tracing then
       M.tracel "dim_remove" "dim remove at positions [%s] in { %s } -> { %s }"
@@ -343,7 +344,7 @@ struct
           | [(coeff,var,divi)] -> Some (Rhs.canonicalize (Some (Z.mul divisor coeff,var), Z.mul constant divi,Z.mul divisor divi))
           |_ -> None))
 
-  let simplify_to_ref_and_offset t texp = timing_wrap "coeff_vec" (simplify_to_ref_and_offset t) texp
+  let simplify_to_ref_and_offset t texp = Timing.wrap "coeff_vec" (simplify_to_ref_and_offset t) texp
 
   let assign_const t var const divi = match t.d with
     | None -> t
@@ -365,7 +366,7 @@ struct
          Some res, Some res)
       | _ -> None, None
 
-  let bound_texpr d texpr1 = timing_wrap "bounds calculation" (bound_texpr d) texpr1
+  let bound_texpr d texpr1 = Timing.wrap "bounds calculation" (bound_texpr d) texpr1
 end
 
 module D =
@@ -456,7 +457,7 @@ struct
     if M.tracing then M.tracel "meet" "meet a: %s\n U  \n b: %s \n -> %s" (show t1) (show t2) (show res) ;
     res
 
-  let meet t1 t2 = timing_wrap "meet" (meet t1) t2
+  let meet t1 t2 = Timing.wrap "meet" (meet t1) t2
 
   let leq t1 t2 =
     let env_comp = Environment.cmp t1.env t2.env in (* Apron's Environment.cmp has defined return values. *)
@@ -475,7 +476,7 @@ struct
       let m1' = if env_comp = 0 then m1 else VarManagement.dim_add (Environment.dimchange t1.env t2.env) m1 in
       EConj.IntMap.for_all (implies m1') (snd m2) (* even on sparse m2, it suffices to check the non-trivial equalities, still present in sparse m2 *)
 
-  let leq a b = timing_wrap "leq" (leq a) b
+  let leq a b = Timing.wrap "leq" (leq a) b
 
   let leq t1 t2 =
     let res = leq t1 t2 in
@@ -551,7 +552,7 @@ struct
     | Some x, Some y when EConj.equal x y -> {d = Some x; env = a.env}
     | Some x, Some y  -> {d = join_d x y; env = a.env}
 
-  let join a b = timing_wrap "join" (join a) b
+  let join a b = Timing.wrap "join" (join a) b
 
   let join a b =
     let res = join a b in
@@ -593,7 +594,7 @@ struct
     if M.tracing then M.tracel "ops" "forget_vars %s -> %s" (show t) (show res);
     res
 
-  let forget_vars t vars = timing_wrap "forget_vars" (forget_vars t) vars
+  let forget_vars t vars = Timing.wrap "forget_vars" (forget_vars t) vars
 
   (** implemented as described on page 10 in the paper about Fast Interprocedural Linear Two-Variable Equalities in the Section "Abstract Effect of Statements"
       This makes a copy of the data structure, it doesn't change it in-place. *)
@@ -617,7 +618,7 @@ struct
       end
     | None -> bot_env
 
-  let assign_texpr t var texp = timing_wrap "assign_texpr" (assign_texpr t var) texp
+  let assign_texpr t var texp = Timing.wrap "assign_texpr" (assign_texpr t var) texp
 
   (* no_ov -> no overflow
      if it's true then there is no overflow
@@ -664,7 +665,7 @@ struct
     if M.tracing then M.tracel "ops" "assign_var parallel: %s -> %s" (show t) (show res);
     res
 
-  let assign_var_parallel t vv's = timing_wrap "var_parallel" (assign_var_parallel t) vv's
+  let assign_var_parallel t vv's = Timing.wrap "var_parallel" (assign_var_parallel t) vv's
 
   let assign_var_parallel_with t vv's =
     (* TODO: If we are angling for more performance, this might be a good place ot try. `assign_var_parallel_with` is used whenever a function is entered (body),
@@ -696,7 +697,7 @@ struct
     if M.tracing then M.tracel "ops" "Substitute_expr t: \n %s \n var: %a \n exp: %a \n -> \n %s" (show t) Var.pretty var d_exp exp (show res);
     res
 
-  let substitute_exp ask t var exp no_ov = timing_wrap "substitution" (substitute_exp ask t var exp) no_ov
+  let substitute_exp ask t var exp no_ov = Timing.wrap "substitution" (substitute_exp ask t var exp) no_ov
 
 
   (** Assert a constraint expression.
@@ -751,7 +752,7 @@ struct
     if M.tracing then M.tracel "meet_tcons" "meet_tcons with expr: %a no_ov:%b" d_exp original_expr (Lazy.force no_ov);
     meet_tcons ask t tcons original_expr no_ov
 
-  let meet_tcons t tcons expr = timing_wrap "meet_tcons" (meet_tcons t tcons) expr
+  let meet_tcons t tcons expr = Timing.wrap "meet_tcons" (meet_tcons t tcons) expr
 
   let unify a b =
     meet a b
@@ -777,7 +778,7 @@ struct
     | tcons1 -> meet_tcons ask d tcons1 e no_ov
     | exception Convert.Unsupported_CilExp _ -> d
 
-  let assert_constraint ask d e negate no_ov = timing_wrap "assert_constraint" (assert_constraint ask d e negate) no_ov
+  let assert_constraint ask d e negate no_ov = Timing.wrap "assert_constraint" (assert_constraint ask d e negate) no_ov
 
   let relift t = t
 
diff --git a/src/cdomains/apron/sharedFunctions.apron.ml b/src/cdomains/apron/sharedFunctions.apron.ml
index 86b5f2770f..bda88a3f62 100644
--- a/src/cdomains/apron/sharedFunctions.apron.ml
+++ b/src/cdomains/apron/sharedFunctions.apron.ml
@@ -4,6 +4,7 @@ open Batteries
 
 open GobApron
 
+
 module M = Messages
 
 
@@ -402,7 +403,7 @@ struct
         in
         {d = Some (if add then RelDomain.dim_add dim_change m else RelDomain.dim_remove dim_change m ~del:del); env = new_env}
 
-  let change_d t new_env ~add ~del = VectorMatrix.timing_wrap "dimension change" (fun del -> change_d t new_env ~add:add ~del:del) del
+  let change_d t new_env ~add ~del = Timing.wrap "dimension change" (fun del -> change_d t new_env ~add:add ~del:del) del
 
   let vars x = Environment.ivars_only x.env
 
@@ -411,18 +412,18 @@ struct
     let env' = Environment.add_vars t.env vars in
     change_d t env' ~add:true ~del:false
 
-  let add_vars t vars = VectorMatrix.timing_wrap "add_vars" (add_vars t) vars
+  let add_vars t vars = Timing.wrap "add_vars" (add_vars t) vars
 
   let drop_vars t vars ~del =
     let t = copy t in
     let env' = Environment.remove_vars t.env vars in
     change_d t env' ~add:false ~del:del
 
-  let drop_vars t vars = VectorMatrix.timing_wrap "drop_vars" (drop_vars t) vars
+  let drop_vars t vars = Timing.wrap "drop_vars" (drop_vars t) vars
 
   let remove_vars t vars = drop_vars t vars ~del:false
 
-  let remove_vars t vars = VectorMatrix.timing_wrap "remove_vars" (remove_vars t) vars
+  let remove_vars t vars = Timing.wrap "remove_vars" (remove_vars t) vars
 
   let remove_vars_with t vars =
     let t' = remove_vars t vars in
@@ -433,7 +434,7 @@ struct
     let env' = Environment.remove_filter t.env f in
     change_d t env' ~add:false ~del:false
 
-  let remove_filter t f = VectorMatrix.timing_wrap "remove_filter" (remove_filter t) f
+  let remove_filter t f = Timing.wrap "remove_filter" (remove_filter t) f
 
   let remove_filter_with t f =
     let t' = remove_filter t f in
@@ -445,14 +446,14 @@ struct
     let env' = Environment.keep_filter t.env f in
     change_d t env' ~add:false ~del:false
 
-  let keep_filter t f = VectorMatrix.timing_wrap "keep_filter" (keep_filter t) f
+  let keep_filter t f = Timing.wrap "keep_filter" (keep_filter t) f
 
   let keep_vars t vs =
     let t = copy t in
     let env' = Environment.keep_vars t.env vs in
     change_d t env' ~add:false ~del:false
 
-  let keep_vars t vs = VectorMatrix.timing_wrap "keep_vars" (keep_vars t) vs
+  let keep_vars t vs = Timing.wrap "keep_vars" (keep_vars t) vs
 
   let mem_var t var = Environment.mem_var t.env var
 
diff --git a/src/goblint_lib.ml b/src/goblint_lib.ml
index 91f9837419..0bf0be6bfb 100644
--- a/src/goblint_lib.ml
+++ b/src/goblint_lib.ml
@@ -430,7 +430,18 @@ module BaseInvariant = BaseInvariant
 module CommonPriv = CommonPriv
 module WideningThresholds = WideningThresholds
 
-module VectorMatrix = VectorMatrix
+(* There might be a more elegant solution. *)
+module Vector = Vector
+module Matrix = Matrix
+module AbstractVector = AbstractVector 
+module AbstractMatrix = AbstractMatrix
+module ArrayVector = ArrayVector
+module ArrayMatrix = ArrayMatrix
+module SparseVector = SparseVector
+module ListMatrix = ListMatrix
+module RatOps = RatOps
+module ConvenienceOps = ConvenienceOps
+
 module SharedFunctions = SharedFunctions
 module GobApron = GobApron
 
diff --git a/tests/regression/63-affeq/21-function_call.c b/tests/regression/63-affeq/21-function_call.c
new file mode 100644
index 0000000000..723b75e14a
--- /dev/null
+++ b/tests/regression/63-affeq/21-function_call.c
@@ -0,0 +1,21 @@
+// SKIP PARAM: --set ana.activated[+] affeq --set sem.int.signed_overflow assume_none
+// This test was added from 77-lin2vareq because it found mistakes in dim_add that weren't detected by the other tests
+#include <goblint.h>
+
+int check_equal(int x, int y, int z) {
+  __goblint_check(x == y);
+  __goblint_check(z == y);
+  __goblint_check(x == z);
+  return 8;
+}
+
+int main(void) {
+  int x, y, z;
+
+  y = x;
+  z = y;
+
+  check_equal(x, y, z);
+
+  return 0;
+}
diff --git a/tests/unit/cdomains/affineEqualityDomain/sparseImplementation/sparseMatrixImplementationTest.ml b/tests/unit/cdomains/affineEqualityDomain/sparseImplementation/sparseMatrixImplementationTest.ml
new file mode 100644
index 0000000000..b562760a45
--- /dev/null
+++ b/tests/unit/cdomains/affineEqualityDomain/sparseImplementation/sparseMatrixImplementationTest.ml
@@ -0,0 +1,370 @@
+(*
+  To run this, type `dune runtest tests/unit/`. 
+*)
+
+open Goblint_lib
+open OUnit2
+open SparseVector
+open ListMatrix
+open ArrayVector
+open ArrayMatrix
+open ConvenienceOps
+
+module D = SharedFunctions.Mpqf
+module Vector = SparseVector (D)
+module Matrix = ListMatrix (D) (SparseVector)
+include ConvenienceOps(D)
+
+(** Shorthands for common functions. *)
+let int x = D.of_int x
+
+let frac numerator denominator = D.of_frac numerator denominator
+
+(** Shorthands for common functions. *)
+let float x = D.of_float x
+
+let make_matrix_of_2d_list l =
+  List.fold_left
+    (fun acc row -> Matrix.append_row acc (Vector.of_list row))
+    (Matrix.empty ())
+    l
+
+(** This function runs the equality assertion with the solution after normalizing the matrix. *)
+let normalize_and_assert (matrix : Matrix.t) (solution : Matrix.t) =
+  let get_dimensions m = (Matrix.num_rows m, Matrix.num_cols m) in
+  let do_dimensions_match dim1 dim2 =
+    match (dim1, dim2) with (r1, c1), (r2, c2) -> r1 == r2 && c1 == c2
+  in
+  let matrix_dim = get_dimensions matrix in
+  let solution_dim = get_dimensions solution in
+  if not (do_dimensions_match solution_dim matrix_dim) then
+    failwith
+      "The matrix to normalize and the solution have different dimensions!"
+  else
+    match Matrix.normalize matrix with
+    | None -> assert_failure "The normalization returned None."
+    | Some reduced_matrix -> assert_equal reduced_matrix solution
+
+(**
+   Example from a [Youtube video](https://www.youtube.com/watch?v=TYs4h-AoqyQ)
+   but extended by a solution vector b = [0;0;25]^T.
+*)
+let standard_normalize _ =
+  let width = 5 in
+  let any_matrix =
+    Matrix.append_row
+      (Matrix.append_row
+         (Matrix.append_row (Matrix.empty ())
+            (Vector.of_sparse_list width
+               [ (0, int 2); (2, int 4); (3, int (-1)) ]))
+         (Vector.of_sparse_list width [ (0, int 2); (1, int 3); (2, int 4) ]))
+      (Vector.of_sparse_list width
+         [ (0, int 1); (1, int 1); (2, int 2); (3, int 4); (4, int 25) ])
+  in
+  let normalized_matrix =
+    Matrix.append_row
+      (Matrix.append_row
+         (Matrix.append_row (Matrix.empty ())
+            (Vector.of_sparse_list width [ (0, int 1); (2, int 2); (4, int 3) ]))
+         (Vector.of_sparse_list width [ (1, int 1); (4, int (-2)) ]))
+      (Vector.of_sparse_list width [ (3, int 1); (4, int 6) ])
+  in
+  normalize_and_assert any_matrix normalized_matrix
+
+(**
+   Normalization just sorts the matrix, as the rows are already reduced. 
+*)
+let does_just_sort _ =
+  let width = 7 in
+  let chaotic_matrix =
+    Matrix.append_row
+      (Matrix.append_row
+         (Matrix.append_row
+            (Matrix.append_row (Matrix.empty ())
+               (Vector.of_sparse_list width [ (2, int 1) ]))
+            (Vector.of_sparse_list width [ (5, int 1) ]))
+         (Vector.of_sparse_list width [ (0, int 1); (3, int 1) ]))
+      (Vector.of_sparse_list width [ (1, int 1); (4, int (-3)) ])
+  in
+
+  let sorted_matrix =
+    Matrix.append_row
+      (Matrix.append_row
+         (Matrix.append_row
+            (Matrix.append_row (Matrix.empty ())
+               (Vector.of_sparse_list width [ (0, int 1); (3, int 1) ]))
+            (Vector.of_sparse_list width [ (1, int 1); (4, int (-3)) ]))
+         (Vector.of_sparse_list width [ (2, int 1) ]))
+      (Vector.of_sparse_list width [ (5, int 1) ])
+  in
+  normalize_and_assert chaotic_matrix sorted_matrix
+
+(**
+   Normalization should eliminate both linearly dependent rows.
+*)
+let does_eliminate_dependent_rows _ =
+  let width = 3 in
+  let linearly_dependent_rows =
+    Matrix.append_row
+      (Matrix.append_row
+         (Matrix.append_row (Matrix.empty ())
+            (Vector.of_sparse_list width [ (0, int 1); (2, int 4) ]))
+         (Vector.of_sparse_list width [ (0, int 2); (2, int 8) ]))
+      (Vector.of_sparse_list width [ (0, int 3); (2, int 12) ])
+  in
+
+  let minimized_matrix =
+    Matrix.append_row
+      (Matrix.append_row
+         (Matrix.append_row (Matrix.empty ())
+            (Vector.of_sparse_list width [ (0, int 1); (2, int 4) ]))
+         (Vector.zero_vec width))
+      (Vector.zero_vec width)
+  in
+  normalize_and_assert linearly_dependent_rows minimized_matrix
+
+let does_handle_floats _ =
+  let width = 3 in
+  let any_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width [ (0, float 5.); (2, float (1. /. 2.)) ]))
+      (Vector.of_sparse_list width
+         [ (0, float (1. /. 4.)); (1, float 23.); (2, float 2.) ])
+  in
+
+  let normalized_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width [ (0, float 1.); (2, frac 1 10) ]))
+      (Vector.of_sparse_list width [ (1, float 1.); (2, frac 79 920) ])
+  in
+  normalize_and_assert any_matrix normalized_matrix
+
+let does_handle_fractions _ =
+  let width = 3 in
+  let any_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width [ (0, frac 5 1); (2, frac 1 2) ]))
+      (Vector.of_sparse_list width
+         [ (0, frac 1 4); (1, frac 23 1); (2, frac 2 1) ])
+  in
+
+  let normalized_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width [ (0, frac 1 1); (2, frac 1 10) ]))
+      (Vector.of_sparse_list width [ (1, frac 1 1); (2, frac 79 920) ])
+  in
+  normalize_and_assert any_matrix normalized_matrix
+
+let does_negate_negative _ =
+  let width = 5 in
+  let negative_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width [ (0, int (-1)); (3, int (-3)) ]))
+      (Vector.of_sparse_list width [ (2, int (-1)); (4, int (-5)) ])
+  in
+
+  let negated_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width [ (0, int 1); (3, int 3) ]))
+      (Vector.of_sparse_list width [ (2, int 1); (4, int 5) ])
+  in
+  normalize_and_assert negative_matrix negated_matrix
+
+(**
+   Normalization is idempotent.
+*)
+let does_not_change_normalized_matrix _ =
+  let width = 5 in
+  let already_normalized =
+    Matrix.append_row
+      (Matrix.append_row
+         (Matrix.append_row (Matrix.empty ())
+            (Vector.of_sparse_list width [ (0, int 1); (3, int 1) ]))
+         (Vector.of_sparse_list width [ (1, int 1); (3, int 3) ]))
+      (Vector.of_sparse_list width [ (2, int 1); (3, int 1) ])
+  in
+  normalize_and_assert already_normalized already_normalized
+
+let is_covered_by_simple _ =
+  let m1 = Matrix.init_with_vec (Vector.of_list [int 1; int 1; int 2; int 10]) in
+  let m2 = Matrix.append_row 
+      (Matrix.append_row 
+         (Matrix.empty ()) 
+         (Vector.of_list [int 1; int 1; int 0; int 6])) 
+      (Vector.of_list [int 0; int 0; int 1; int 2]) in
+  let result = Matrix.is_covered_by m1 m2 in
+  assert_bool "Matrix m1 is covered by m2, but was false" result
+
+let is_covered_by_vector_first_row _ =
+  let m1 = Matrix.init_with_vec (Vector.of_list [int 1; int 2; int 0; int 7]) in
+  let m2 = Matrix.append_row 
+      (Matrix.append_row 
+         (Matrix.empty ()) 
+         (Vector.of_list [int 1; int 2; int 0; int 7])) 
+      (Vector.of_list [int 0; int 0; int 1; int 2]) in
+  let result = Matrix.is_covered_by m1 m2 in
+  assert_bool "Matrix m1 is covered by m2, but was false" result
+
+let is_zero_vec_covered _ = 
+  let m1 = Matrix.init_with_vec (Vector.zero_vec 4) in
+  let m2 = Matrix.append_row 
+      (Matrix.append_row 
+         (Matrix.empty ()) 
+         (Vector.of_list [int 1; int 2; int 0; int 7])) 
+      (Vector.of_list [int 0; int 0; int 1; int 2]) in
+  let result = Matrix.is_covered_by m1 m2 in
+  assert_bool "Matrix m1 is covered by m2, but was false" result
+
+let is_not_covered _ =
+  let m1 = Matrix.init_with_vec (Vector.of_list [int 1; int 1; int 2; int 10]) in
+  let m2 = Matrix.append_row 
+      (Matrix.append_row 
+         (Matrix.empty ()) 
+         (Vector.of_list [int 1; int 1; int 0; int 6])) 
+      (Vector.of_list [int 0; int 0; int 1; int 3]) in
+  let result = Matrix.is_covered_by m2 m1 in
+  assert_bool "Matrix m1 is not covered by m2, but was true" (not result)
+
+let is_covered_big _ =
+  let m1 = make_matrix_of_2d_list @@
+    [[int 1; int 0; int 0; int 0; int 0; int (-1); int 0];
+     [int 0; int 1; int 0; int 0; int 0; int (-2); int 0];
+     [int 0; int 0; int 1;  (frac (-1) 3); frac 1 3; int 0; frac 1 3]] in
+
+  let m2 = make_matrix_of_2d_list @@
+    [[int 1; int 0; int 0; int 0; int 0; int 0; int 0];
+     [int 0; int 1; int 0; int 0; int 0; int 0; int 0];
+     [int 0; int 0; int 1; frac (-1) 3; frac 1 3; int 0; frac 1 3];
+     [int 0; int 0; int 0; int 0; int 0; int 1; int 0]] in
+
+  let result = Matrix.is_covered_by m1 m2 in
+  assert_bool "Matrix m1 is covered by m2, but was false" (result)
+
+let is_covered_big2 _ = 
+  let m1 = make_matrix_of_2d_list @@
+    [[int 1; int 0; int 0; int 0; int 0; int 1; int 0]
+    ] in
+
+  let m2 = make_matrix_of_2d_list @@
+    [[int 1; int 0; int 0; int 0; int 0; int 0; int 0];
+     [int 0; int 1; int 0; int 0; int 0; int 0; int 0];
+     [int 0; int 0; int 0; int 0; int 0; int 1; int 0]] in
+
+  let result = Matrix.is_covered_by m1 m2 in
+  assert_bool "Matrix m1 is covered by m2, but was false" (result)
+(**
+   Normalization works on an empty matrix.
+*)
+let normalize_empty _ =
+  let width = 3 in
+  let empty_matrix =
+    Matrix.append_row
+      (Matrix.append_row
+         (Matrix.append_row (Matrix.empty ()) (Vector.zero_vec width))
+         (Vector.zero_vec width))
+      (Vector.zero_vec width)
+  in
+  normalize_and_assert empty_matrix empty_matrix
+
+let normalize_two_columns _ =
+  let width = 2 in
+  let two_col_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width [ (0, int 3); (1, int 2) ]))
+      (Vector.of_sparse_list width [ (0, int 9); (1, int 6) ])
+  in
+  let normalized_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width
+            [ (0, int 1); (1, D.div (int 2) (int 3)) ]))
+      (Vector.of_sparse_list width [])
+  in
+  normalize_and_assert two_col_matrix normalized_matrix
+
+let int_domain_to_rational _ =
+  let width = 3 in
+  let int_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width [ (0, int 3); (2, int 1) ]))
+      (Vector.of_sparse_list width [ (0, int 2); (1, int 1); (2, int 1) ])
+  in
+  let normalized_matrix =
+    Matrix.append_row
+      (Matrix.append_row (Matrix.empty ())
+         (Vector.of_sparse_list width
+            [ (0, int 1); (2, D.div (int 1) (int 3)) ]))
+      (Vector.of_sparse_list width [ (1, int 1); (2, D.div (int 1) (int 3)) ])
+  in
+  normalize_and_assert int_matrix normalized_matrix
+
+
+let vectorMap2i _ =
+  let v1 = Vector.of_list [int 0; int 1; int 0; int 2; int 3; int 0; int 4; int 0; int 1] in
+  let v2 = Vector.of_list [int 4; int 0; int 0; int 0; int 5; int 6; int 0; int 0; int 2] in
+  let result = Vector.map2i (fun i x y -> (int i) *: (x +: y)) v1 v2 in
+  let expected = Vector.of_list [int 0; int 1; int 0; int 6; int 32; int 30; int 24; int 0; int 24] in
+  assert_equal expected result
+
+
+let vectorMap2i_empty _ = 
+  let v1 = Vector.of_list [] in
+  let v2 = Vector.of_list [] in
+  let result = Vector.map2i (fun i x y -> (int i) *: (x +: y)) v1 v2 in
+  let expected = Vector.of_list [] in
+  assert_equal expected result
+
+let vectorMap2i_one_zero _ = 
+    let v1 = Vector.of_list [int 0; int 0; int 0; int 0] in
+    let v2 = Vector.of_list [int 1; int 2; int 3; int 4] in
+    let result = Vector.map2i (fun i x y -> (int i) *: (x +: y)) v1 v2 in
+    let expected = Vector.of_list [int 0; int 2; int 6; int 12] in
+    assert_equal expected result
+
+let vectorMap _ =
+  let v1 = Vector.of_list [int 0; int 1; int 2; int 0; int 0; int 3; int 4; int 0; int 0; int 5] in
+  let result = Vector.map (fun x -> x +: int 1) v1 in
+  let expected = Vector.of_list [int 1; int 2; int 3; int 1; int 1; int 4; int 5; int 1; int 1; int 6] in
+  assert_equal expected result
+
+let vectorMap_zero_preserving_normal _ =
+  let v1 = Vector.of_list [int 0; int 1; int 2; int 0; int 0; int 4; int 5; int 0; int 0;] in
+  let result = Vector.map (fun x -> x *: x) v1 in
+  let expected = Vector.of_list [int 0; int 1; int 4; int 0; int 0; int 16; int 25; int 0; int 0;] in
+  assert_equal expected result
+
+  let tests =
+    "SparseMatrixImplementationTest"
+    >::: [
+      "can solve a standard normalization" >:: standard_normalize;
+      "does sort already reduzed" >:: does_just_sort;
+      "does eliminate dependent rows" >:: does_eliminate_dependent_rows;
+      "can handle float domain" >:: does_handle_floats;
+      "can handle fraction domain" >:: does_handle_fractions;
+      "does negate negative matrix" >:: does_negate_negative;
+      "does not change already normalized matrix" >:: does_not_change_normalized_matrix;
+      "m1 is covered by m2" >:: is_covered_by_simple;
+      "m1 is covered by m2 with vector in first row" >:: is_covered_by_vector_first_row;
+      "zero vector is covered by m2" >:: is_zero_vec_covered;
+      "m1 is not covered by m2" >:: is_not_covered;
+      "m1 is covered by m2 with big matrix" >:: is_covered_big;
+      "does not change an empty matrix" >:: normalize_empty;
+      "can correctly normalize a two column matrix" >:: normalize_two_columns;
+      "can handle a rational solution" >:: int_domain_to_rational;
+      "m1 is covered by m2 with big matrix2" >:: is_covered_big2;
+      "map2i two vectors" >:: vectorMap2i;
+      "map2i two empty vectors" >:: vectorMap2i_empty;
+      "map2i one zero vector" >:: vectorMap2i_one_zero;
+      "map one vector" >:: vectorMap;
+      "map zero preserving normal" >:: vectorMap_zero_preserving_normal;
+    ]
+
+let () = run_test_tt_main tests
diff --git a/tests/unit/mainTest.ml b/tests/unit/mainTest.ml
index 4f071ea25c..9527679fe7 100644
--- a/tests/unit/mainTest.ml
+++ b/tests/unit/mainTest.ml
@@ -7,14 +7,19 @@ let all_tests =
     MapDomainTest.test ();
     SolverTest.test ();
     LvalTest.test ();
+    SparseMatrixImplementationTest.tests ;
     CompilationDatabaseTest.tests;
     LibraryDslTest.tests;
     CilfacadeTest.tests;
     (* etc *)
     "domaintest" >::: QCheck_ounit.to_ounit2_test_list Maindomaintest.all_testsuite;
     IntOpsTest.tests;
+    (* SparseMatrixImplementationTest.tests; *) (* Uncomment this to add the sparse matrix tests to all tests *)
   ]
 
+let subset_tests = "" >::: [SparseMatrixImplementationTest.tests]
+
 let () =
   print_string "\027[0;1munit: \027[0;0;00m";
-  run_test_tt_main all_tests
+  run_test_tt_main subset_tests (* Remove this and uncomment the line below to run all tests.*)
+(* run_test_tt_main all_tests *)