diff --git a/Julia/basiclu.jl b/Julia/basiclu.jl
index 2715b3c..f15fb02 100644
--- a/Julia/basiclu.jl
+++ b/Julia/basiclu.jl
@@ -94,6 +94,7 @@ const BASICLU_MATRIX_INFNORM = 100
 const BASICLU_RESIDUAL_TEST = 112
 
 const BASICLU_MATRIX_NZ = 101
+const BASICLU_NBUCKETS = 113
 const BASICLU_RANK = 102
 const BASICLU_BUMP_SIZE = 103
 const BASICLU_BUMP_NZ = 104
@@ -223,31 +224,36 @@ end
 # factorize
 # =============================================================================
 """
-    factorize(this, B)
+    factorize(this, B, buckets=C_NULL)
 
 Load sparse matrix into BASICLU object and factorize it. @B must be square
-and have the same dimension for which @this was initialized.
+and have the same dimension for which @this was initialized. If @buckets is
+not C_NULL, then it must be of type Array{cint,1} and of length size(B,2).
 
 Return status code.
 """
-function factorize(this::BLU, B::spmatrix)
+function factorize(this::BLU, B::spmatrix, buckets=C_NULL)
     m = convert(cint, this.xstore[BASICLU_DIM])
     nrow, ncol = size(B)
     @assert nrow == ncol
     @assert nrow == m
+    if buckets != C_NULL
+        @assert typeof(buckets) == Array{cint,1}
+        @assert length(buckets) == m
+    end
     Bp = B.colptr-1
     Bi = B.rowval-1
     Bx = B.nzval                # don't need a copy
     c0ntinue = 0
     err = BASICLU_OK
     while true
-        err = ccall((:basiclu_factorize, "libbasiclu.so"), cint,
+        err = ccall((:basiclu_factorize_buckets, "libbasiclu.so"), cint,
                     (cint_ptr, cdbl_ptr,
                      cint_ptr, cdbl_ptr, cint_ptr, cdbl_ptr, cint_ptr, cdbl_ptr,
-                     cint_ptr, cint_ptr, cint_ptr, cdbl_ptr, cint),
+                     cint_ptr, cint_ptr, cint_ptr, cdbl_ptr, cint_ptr, cint),
                     this.istore, this.xstore,
                     this.Li, this.Lx, this.Ui, this.Ux, this.Wi, this.Wx,
-                    Bp, pointer(Bp, 2), Bi, Bx, c0ntinue)
+                    Bp, pointer(Bp, 2), Bi, Bx, buckets, c0ntinue)
         if err != BASICLU_REALLOCATE break; end
         realloc(this)
         c0ntinue = 1
diff --git a/Julia/test.jl b/Julia/test.jl
index 0cf0c08..d61c973 100644
--- a/Julia/test.jl
+++ b/Julia/test.jl
@@ -19,14 +19,18 @@ end
 # test_factorize
 # =============================================================================
 """
-    test_factorize(testdir, trans=false)
+    test_factorize(testdir, trans=false, nbuckets=0)
 
 For all *.mat files in @testdir read matrix B and factorize it. Monitor residual
 of factorization and forward/backward solves.
 
 @trans specifies if B or its transposed are factorized.
+@nbuckets specifies the number of column buckets.
+ nbuckets <= 0: no buckets used
+ nbuckets >= m: each column in its own bucket
+ 1 <= nbuckets <= m-1: column j in bucket mod(j-1, nbuckets)
 """
-function test_factorize(testdir::String, trans=false)
+function test_factorize(testdir::String, trans=false, nbuckets=0)
     files = readdir(testdir)
     for f in files
         if length(f) < 4 || f[end-3:end] != ".mat"
@@ -39,8 +43,16 @@ function test_factorize(testdir::String, trans=false)
             (B,Bt) = (Bt,B)
         end
         m = size(B,1)
+        buckets = C_NULL
+        if nbuckets > 0
+            buckets = Array{cint}(m)
+            nbuckets = convert(cint, min(nbuckets, m))
+            for j = 1:m
+                buckets[j] = mod(j-1, nbuckets)
+            end
+        end
         blu = initialize(m)
-        err = factorize(blu, B)
+        err = factorize(blu, B, buckets)
         if err != BASICLU_OK
             @printf(" failed (%d)\n", err)
             continue
diff --git a/doc/userguide.pdf b/doc/userguide.pdf
index 72cb069..1dbfe74 100644
Binary files a/doc/userguide.pdf and b/doc/userguide.pdf differ
diff --git a/doc/userguide.tex b/doc/userguide.tex
index 9591940..c0a4ea0 100644
--- a/doc/userguide.tex
+++ b/doc/userguide.tex
@@ -95,6 +95,12 @@ \subsection{basiclu\_factorize}
 \input{basiclu_factorize}
 }
 
+\newpage
+\subsection{basiclu\_factorize\_buckets}
+{\footnotesize
+\input{basiclu_factorize_buckets}
+}
+
 \newpage
 \subsection{basiclu\_get\_factors}
 {\footnotesize
@@ -154,6 +160,12 @@ \subsection{basiclu\_obj\_factorize}
 \input{basiclu_obj_factorize}
 }
 
+\newpage
+\subsection{basiclu\_obj\_factorize\_buckets}
+{\footnotesize
+\input{basiclu_obj_factorize_buckets}
+}
+
 \newpage
 \subsection{basiclu\_obj\_get\_factors}
 {\footnotesize
diff --git a/include/basiclu.h b/include/basiclu.h
index 9bc784e..146f865 100644
--- a/include/basiclu.h
+++ b/include/basiclu.h
@@ -111,6 +111,7 @@ typedef int64_t lu_int;
 #define BASICLU_RESIDUAL_TEST 111
 
 #define BASICLU_MATRIX_NZ 100
+#define BASICLU_NBUCKETS 112
 #define BASICLU_RANK 101
 #define BASICLU_BUMP_SIZE 102
 #define BASICLU_BUMP_NZ 103
@@ -130,6 +131,7 @@ typedef int64_t lu_int;
 
 #include "basiclu_initialize.h"
 #include "basiclu_factorize.h"
+#include "basiclu_factorize_buckets.h"
 #include "basiclu_get_factors.h"
 #include "basiclu_solve_dense.h"
 #include "basiclu_solve_sparse.h"
@@ -139,6 +141,7 @@ typedef int64_t lu_int;
 #include "basiclu_object.h"
 #include "basiclu_obj_initialize.h"
 #include "basiclu_obj_factorize.h"
+#include "basiclu_obj_factorize_buckets.h"
 #include "basiclu_obj_get_factors.h"
 #include "basiclu_obj_solve_dense.h"
 #include "basiclu_obj_solve_sparse.h"
diff --git a/include/basiclu_factorize_buckets.h b/include/basiclu_factorize_buckets.h
new file mode 100644
index 0000000..b2967af
--- /dev/null
+++ b/include/basiclu_factorize_buckets.h
@@ -0,0 +1,69 @@
+lu_int basiclu_factorize_buckets
+(
+    lu_int istore[],
+    double xstore[],
+    lu_int Li[],
+    double Lx[],
+    lu_int Ui[],
+    double Ux[],
+    lu_int Wi[],
+    double Wx[],
+    const lu_int Bbegin[],
+    const lu_int Bend[],
+    const lu_int Bi[],
+    const double Bx[],
+    const lu_int *buckets,
+    lu_int c0ntinue
+);
+
+/*
+Purpose:
+
+    Performs the same task as basiclu_factorize(), just that the column ordering
+    of the factorization is restricted to pivot on all columns in bucket 0,
+    before those in bucket 1, and so on. If buckets == NULL, then no restriction
+    is imposed and the function performs identically to basiclu_factorize().
+
+Return:
+
+    see basiclu_factorize().
+
+Arguments:
+
+    const lu_int *buckets
+
+        buckets[j] holds the bucket to which column j belongs. Buckets are
+        numbered from 0 to m-1. All columns in bucket 0 are chosen as pivot
+        column before those in bucket 1, and so on. When a column does not
+        contain an eligible pivot element (by means of the absolute pivot
+        tolerance), then it is ordered last in the pivot sequence, outside its
+        bucket. When buckets == NULL, then the function performs identically to
+        basiclu_factorize().
+
+    for the remaining arguments see basiclu_factorize().
+
+Parameters:
+
+    xstore[BASICLU_SEARCH_ROWS]
+
+        This parameter has no effect when buckets are given, in which case only
+        columns are searched for pivot elements.
+
+    for the remaining parameters see basiclu_factorize().
+
+Info:
+
+    xstore[BASICLU_STATUS]: status code.
+
+        BASICLU_ERROR_invalid_argument
+
+            The matrix is invalid (a column has a negative number of entries,
+            a row index is out of range, or a column has duplicate entries) or
+            buckets is invalid (buckets[j] < 0 or buckets[j] >= m for some j).
+
+        for the remaining status codes see basiclu_factorize().
+
+    xstore[BASICLU_NBUCKETS] number of (nonempty) buckets. 0 if buckets == NULL.
+
+    for the remaining info fields see basiclu_factorize().
+*/
diff --git a/include/basiclu_obj_factorize.h b/include/basiclu_obj_factorize.h
index 8852fc8..a2a710e 100644
--- a/include/basiclu_obj_factorize.h
+++ b/include/basiclu_obj_factorize.h
@@ -1,10 +1,10 @@
 lu_int basiclu_obj_factorize
 (
     struct basiclu_object *obj,
-    const lu_int *Bbegin,
-    const lu_int *Bend,
-    const lu_int *Bi,
-    const double *Bx
+    const lu_int Bbegin[],
+    const lu_int Bend[],
+    const lu_int Bi[],
+    const double Bx[]
 );
 
 /*
diff --git a/include/basiclu_obj_factorize_buckets.h b/include/basiclu_obj_factorize_buckets.h
new file mode 100644
index 0000000..37a4002
--- /dev/null
+++ b/include/basiclu_obj_factorize_buckets.h
@@ -0,0 +1,35 @@
+lu_int basiclu_obj_factorize_buckets
+(
+    struct basiclu_object *obj,
+    const lu_int Bbegin[],
+    const lu_int Bend[],
+    const lu_int Bi[],
+    const double Bx[],
+    const lu_int *buckets
+);
+
+/*
+Purpose:
+
+    Call basiclu_factorize_buckets() on a BASICLU object.
+
+Return:
+
+    BASICLU_ERROR_invalid_object
+
+        obj is NULL or initialized to a null object.
+
+    BASICLU_ERROR_out_of_memory
+
+        reallocation failed because of insufficient memory.
+
+    Other return codes are passed through from basiclu_factorize_buckets().
+
+Arguments:
+
+    struct basiclu_object *obj
+
+        Pointer to an initialized BASICLU object.
+
+    The other arguments are passed through to basiclu_factorize_buckets().
+*/
diff --git a/src/basiclu_factorize.c b/src/basiclu_factorize.c
index 28a83af..ac202bb 100644
--- a/src/basiclu_factorize.c
+++ b/src/basiclu_factorize.c
@@ -6,7 +6,6 @@
  */
 
 #include "lu_internal.h"
-#include "lu_timer.h"
 
 lu_int basiclu_factorize
 (
@@ -25,119 +24,6 @@ lu_int basiclu_factorize
     lu_int c0ntinue
 )
 {
-    struct lu this;
-    lu_int status;
-    double tic[2], elapsed, factor_cost;
-    lu_tic(tic);
-
-    status = lu_load(&this, istore, xstore, Li, Lx, Ui, Ux, Wi, Wx);
-    if (status != BASICLU_OK)
-        return status;
-
-    if (! (Li && Lx && Ui && Ux && Wi && Wx && Bbegin && Bend && Bi && Bx))
-    {
-        status = BASICLU_ERROR_argument_missing;
-        return lu_save(&this, istore, xstore, status);
-    }
-    if (!c0ntinue)
-    {
-        lu_reset(&this);
-        this.task = SINGLETONS;
-    }
-
-    /* continue factorization */
-    switch (this.task)
-    {
-    case SINGLETONS: goto singletons;
-    case SETUP_BUMP: goto setup_bump;
-    case FACTORIZE_BUMP: goto factorize_bump;
-    case BUILD_FACTORS: goto build_factors;
-    }
-    status = BASICLU_ERROR_invalid_call;
-    return lu_save(&this, istore, xstore, status);
-
-    /*
-     * Each of the following four parts of the factorization calls a routine
-     * lu_do_something() which may request reallocation. In this case return
-     * to the caller immediately, keeping the entry point in this.task.
-     */
-
-singletons:
-    this.task = SINGLETONS;
-    status = lu_singletons(&this, Bbegin, Bend, Bi, Bx);
-    if (status != BASICLU_OK)
-        goto return_to_caller;
-
-setup_bump:
-    this.task = SETUP_BUMP;
-    status = lu_setup_bump(&this, Bbegin, Bend, Bi, Bx);
-    if (status != BASICLU_OK)
-        goto return_to_caller;
-
-factorize_bump:
-    this.task = FACTORIZE_BUMP;
-    status = lu_factorize_bump(&this);
-    if (status != BASICLU_OK)
-        goto return_to_caller;
-
-build_factors:
-    this.task = BUILD_FACTORS;
-    status = lu_build_factors(&this);
-    if (status != BASICLU_OK)
-        goto return_to_caller;
-
-    /* factorization successfully finished */
-    this.task               = NO_TASK;
-    this.nupdate            =  0; /* make factorization valid */
-    this.ftran_for_update   = -1;
-    this.btran_for_update   = -1;
-    this.nfactorize++;
-
-    this.condestL = lu_condest(this.m, this.Lbegin, this.Lindex, this.Lvalue,
-                               NULL, this.p, 0, this.work1,
-                               &this.normL, &this.normestLinv);
-    this.condestU = lu_condest(this.m, this.Ubegin, this.Uindex, this.Uvalue,
-                               this.row_pivot, this.p, 1, this.work1,
-                               &this.normU, &this.normestUinv);
-
-    /* measure numerical stability of the factorization */
-    lu_residual_test(&this, Bbegin, Bend, Bi, Bx);
-
-    /*
-     * factor_cost is a deterministic measure of the factorization cost.
-     * The parameters have been adjusted such that (on my computer)
-     * 1e-6 * factor_cost =~ time_factorize.
-     *
-     * update_cost measures the accumulated cost of updates/solves compared
-     * to the last factorization. It is computed from
-     *
-     *   update_cost = update_cost_numer / update_cost_denom.
-     *
-     * update_cost_denom is fixed here.
-     * update_cost_numer is zero here and increased by solves/updates.
-     *
-     */
-    factor_cost =
-        0.04 * this.m +
-        0.07 * this.matrix_nz +
-        0.20 * this.bump_nz +
-        0.20 * this.nsearch_pivot +
-        0.008 * this.factor_flops;
-
-    this.update_cost_denom = factor_cost * 250;
-
-    #if 0
-    elapsed = this.time_factorize + lu_toc(tic);
-    printf(" 1e-6 * factor_cost / time_factorize: %.2f\n",
-           1e-6 * factor_cost / elapsed);
-    #endif
-
-    if (this.rank < this.m)
-        status = BASICLU_WARNING_singular_matrix;
-
-return_to_caller:
-    elapsed = lu_toc(tic);
-    this.time_factorize += elapsed;
-    this.time_factorize_total += elapsed;
-    return lu_save(&this, istore, xstore, status);
+    return basiclu_factorize_buckets(istore, xstore, Li, Lx, Ui, Ux, Wi, Wx,
+                                     Bbegin, Bend, Bi, Bx, NULL, c0ntinue);
 }
diff --git a/src/basiclu_factorize_buckets.c b/src/basiclu_factorize_buckets.c
new file mode 100644
index 0000000..716d08b
--- /dev/null
+++ b/src/basiclu_factorize_buckets.c
@@ -0,0 +1,144 @@
+/*
+ * basiclu_factorize_buckets.c
+ *
+ * Copyright (C) 2016-2019  ERGO-Code
+ *
+ */
+
+#include "lu_internal.h"
+#include "lu_timer.h"
+
+lu_int basiclu_factorize_buckets
+(
+    lu_int istore[],
+    double xstore[],
+    lu_int Li[],
+    double Lx[],
+    lu_int Ui[],
+    double Ux[],
+    lu_int Wi[],
+    double Wx[],
+    const lu_int Bbegin[],
+    const lu_int Bend[],
+    const lu_int Bi[],
+    const double Bx[],
+    const lu_int *buckets,
+    lu_int c0ntinue
+)
+{
+    struct lu this;
+    lu_int status;
+    double tic[2], elapsed, factor_cost;
+    lu_tic(tic);
+
+    status = lu_load(&this, istore, xstore, Li, Lx, Ui, Ux, Wi, Wx);
+    if (status != BASICLU_OK)
+        return status;
+
+    if (! (Li && Lx && Ui && Ux && Wi && Wx && Bbegin && Bend && Bi && Bx))
+    {
+        status = BASICLU_ERROR_argument_missing;
+        return lu_save(&this, istore, xstore, status);
+    }
+    if (!c0ntinue)
+    {
+        lu_reset(&this);
+        this.task = SINGLETONS;
+    }
+
+    /* continue factorization */
+    switch (this.task)
+    {
+    case SINGLETONS: goto singletons;
+    case SETUP_BUMP: goto setup_bump;
+    case FACTORIZE_BUMP: goto factorize_bump;
+    case BUILD_FACTORS: goto build_factors;
+    }
+    status = BASICLU_ERROR_invalid_call;
+    return lu_save(&this, istore, xstore, status);
+
+    /*
+     * Each of the following four parts of the factorization calls a routine
+     * lu_do_something() which may request reallocation. In this case return
+     * to the caller immediately, keeping the entry point in this.task.
+     */
+
+singletons:
+    this.task = SINGLETONS;
+    status = lu_singletons(&this, Bbegin, Bend, Bi, Bx, buckets);
+    if (status != BASICLU_OK)
+        goto return_to_caller;
+
+setup_bump:
+    this.task = SETUP_BUMP;
+    status = lu_setup_bump(&this, Bbegin, Bend, Bi, Bx);
+    if (status != BASICLU_OK)
+        goto return_to_caller;
+
+factorize_bump:
+    this.task = FACTORIZE_BUMP;
+    status = lu_factorize_bump(&this);
+    if (status != BASICLU_OK)
+        goto return_to_caller;
+
+build_factors:
+    this.task = BUILD_FACTORS;
+    status = lu_build_factors(&this);
+    if (status != BASICLU_OK)
+        goto return_to_caller;
+
+    /* factorization successfully finished */
+    this.task               = NO_TASK;
+    this.nupdate            =  0; /* make factorization valid */
+    this.ftran_for_update   = -1;
+    this.btran_for_update   = -1;
+    this.nfactorize++;
+
+    this.condestL = lu_condest(this.m, this.Lbegin, this.Lindex, this.Lvalue,
+                               NULL, this.p, 0, this.work1,
+                               &this.normL, &this.normestLinv);
+    this.condestU = lu_condest(this.m, this.Ubegin, this.Uindex, this.Uvalue,
+                               this.row_pivot, this.p, 1, this.work1,
+                               &this.normU, &this.normestUinv);
+
+    /* measure numerical stability of the factorization */
+    lu_residual_test(&this, Bbegin, Bend, Bi, Bx);
+
+    /*
+     * factor_cost is a deterministic measure of the factorization cost.
+     * The parameters have been adjusted such that (on my computer)
+     * 1e-6 * factor_cost =~ time_factorize.
+     *
+     * update_cost measures the accumulated cost of updates/solves compared
+     * to the last factorization. It is computed from
+     *
+     *   update_cost = update_cost_numer / update_cost_denom.
+     *
+     * update_cost_denom is fixed here.
+     * update_cost_numer is zero here and increased by solves/updates.
+     *
+     */
+    factor_cost =
+        0.04 * this.m +
+        0.07 * this.matrix_nz +
+        0.20 * this.bump_nz +
+        0.20 * this.nsearch_pivot +
+        0.008 * this.factor_flops;
+
+    this.update_cost_denom = factor_cost * 250;
+
+    #if 0
+    elapsed = this.time_factorize + lu_toc(tic);
+    printf(" 1e-6 * factor_cost / time_factorize: %.2f\n",
+           1e-6 * factor_cost / elapsed);
+    #endif
+
+    if (this.rank < this.m)
+        status = BASICLU_WARNING_singular_matrix;
+
+return_to_caller:
+    elapsed = lu_toc(tic);
+    this.time_factorize += elapsed;
+    this.time_factorize_total += elapsed;
+    return lu_save(&this, istore, xstore, status);
+}
diff --git a/src/basiclu_object.c b/src/basiclu_object.c
index 505f022..e650ba2 100644
--- a/src/basiclu_object.c
+++ b/src/basiclu_object.c
@@ -1,7 +1,7 @@
 /*
  * basiclu_object.c
  *
- * Copyright (C) 2016-2018  ERGO-Code
+ * Copyright (C) 2016-2019  ERGO-Code
  *
  */
 
@@ -200,27 +200,40 @@ void basiclu_obj_free(struct basiclu_object *obj)
 /*
  * basiclu_obj_factorize()
  */
-lu_int basiclu_obj_factorize(struct basiclu_object *obj, const lu_int *Bbegin,
-                             const lu_int *Bend, const lu_int *Bi,
-                             const double *Bx)
+lu_int basiclu_obj_factorize(struct basiclu_object *obj, const lu_int Bbegin[],
+                             const lu_int Bend[], const lu_int Bi[],
+                             const double Bx[])
+{
+    return basiclu_obj_factorize_buckets(obj, Bbegin, Bend, Bi, Bx, NULL);
+}
+
+/*
+ * basiclu_obj_factorize_buckets()
+ */
+lu_int basiclu_obj_factorize_buckets(struct basiclu_object *obj,
+                                     const lu_int Bbegin[], const lu_int Bend[],
+                                     const lu_int Bi[], const double Bx[],
+                                     const lu_int *buckets)
 {
     lu_int status;
 
     if (!isvalid(obj))
         return BASICLU_ERROR_invalid_object;
 
-    status = basiclu_factorize(obj->istore, obj->xstore, obj->Li, obj->Lx,
-                               obj->Ui, obj->Ux, obj->Wi, obj->Wx, Bbegin, Bend,
-                               Bi, Bx, 0);
+    status = basiclu_factorize_buckets(obj->istore, obj->xstore,
+                                       obj->Li, obj->Lx, obj->Ui, obj->Ux,
+                                       obj->Wi, obj->Wx,
+                                       Bbegin, Bend, Bi, Bx, buckets, 0);
 
     while (status == BASICLU_REALLOCATE)
     {
         status = lu_realloc_obj(obj);
         if (status != BASICLU_OK)
             break;
-        status = basiclu_factorize(obj->istore, obj->xstore, obj->Li, obj->Lx,
-                                   obj->Ui, obj->Ux, obj->Wi, obj->Wx, Bbegin,
-                                   Bend, Bi, Bx, 1);
+        status = basiclu_factorize_buckets(obj->istore, obj->xstore,
+                                           obj->Li, obj->Lx, obj->Ui, obj->Ux,
+                                           obj->Wi, obj->Wx,
+                                           Bbegin, Bend, Bi, Bx, buckets, 1);
     }
 
     return status;
diff --git a/src/lu_factorize_bump.c b/src/lu_factorize_bump.c
index 1c9e211..1643866 100644
--- a/src/lu_factorize_bump.c
+++ b/src/lu_factorize_bump.c
@@ -10,6 +10,8 @@
 #include "lu_internal.h"
 #include "lu_list.h"
 
+static void lu_next_bucket(struct lu *this);
+
 lu_int lu_factorize_bump(struct lu *this)
 {
     const lu_int m          = this->m;
@@ -21,6 +23,10 @@ lu_int lu_factorize_bump(struct lu *this)
 
     while (this->rank + this->rankdef < m)
     {
+        if (this->bucket_ptr == this->rank+this->rankdef+1)
+            lu_next_bucket(this);
+        assert(this->ncol_active > 0);
+
         /*
          * Find pivot element. Markowitz search need not be called if the
          * previous call to lu_pivot() returned for reallocation. In this case
@@ -36,6 +42,7 @@ lu_int lu_factorize_bump(struct lu *this)
             lu_list_remove(colcount_flink, colcount_blink, this->pivot_col);
             this->pivot_col = -1;
             this->rankdef++;
+            this->ncol_active--;
         }
         else
         {
@@ -50,7 +57,35 @@ lu_int lu_factorize_bump(struct lu *this)
             this->pivot_col = -1;
             this->pivot_row = -1;
             this->rank++;
+            this->ncol_active--;
         }
     }
     return status;
 }
+
+static void lu_next_bucket(struct lu *this)
+{
+    const lu_int m          = this->m;
+    lu_int *colcount_flink  = this->colcount_flink;
+    lu_int *colcount_blink  = this->colcount_blink;
+    const lu_int *Wbegin    = this->Wbegin;
+    const lu_int *Wend      = this->Wend;
+    const lu_int *Lbegin_p  = this->Lbegin_p;
+    lu_int bucket_ptr       = this->bucket_ptr;
+    lu_int j, nz;
+
+    assert(bucket_ptr >= 1 && bucket_ptr <= m);
+    assert(Lbegin_p[bucket_ptr] < 0);
+    assert(this->ncol_active == 0);
+
+    do {
+        j = Lbegin_p[bucket_ptr++];
+        if (j < 0) j = -j-1;    /* must unflip first index */
+        assert(this->qinv[j] < 0);
+        nz = Wend[j] - Wbegin[j];
+        lu_list_add(j, nz, colcount_flink, colcount_blink, m, &this->min_colnz);
+        this->ncol_active++;
+    } while (bucket_ptr <= m && Lbegin_p[bucket_ptr] >= 0);
+
+    this->bucket_ptr = bucket_ptr;
+}
diff --git a/src/lu_internal.c b/src/lu_internal.c
index 9aaa7d0..a243363 100644
--- a/src/lu_internal.c
+++ b/src/lu_internal.c
@@ -18,6 +18,8 @@
 #define BASICLU_RANKDEF 261
 #define BASICLU_MIN_COLNZ 262
 #define BASICLU_MIN_ROWNZ 263
+#define BASICLU_BUCKET_PTR 264
+#define BASICLU_NCOL_ACTIVE 265
 #define BASICLU_MARKER 266
 #define BASICLU_UPDATE_COST_NUMER 267
 #define BASICLU_UPDATE_COST_DENOM 268
@@ -102,6 +104,7 @@ lu_int lu_load(
     this->residual_test         = xstore[BASICLU_RESIDUAL_TEST];
 
     this->matrix_nz             = xstore[BASICLU_MATRIX_NZ];
+    this->nbuckets              = xstore[BASICLU_NBUCKETS];
     this->rank                  = xstore[BASICLU_RANK];
     this->bump_size             = xstore[BASICLU_BUMP_SIZE];
     this->bump_nz               = xstore[BASICLU_BUMP_NZ];
@@ -124,8 +127,10 @@ lu_int lu_load(
     this->marker                = xstore[BASICLU_MARKER];
     this->pivotlen              = xstore[BASICLU_PIVOTLEN];
     this->rankdef               = xstore[BASICLU_RANKDEF];
+    this->ncol_active           = xstore[BASICLU_NCOL_ACTIVE];
     this->min_colnz             = xstore[BASICLU_MIN_COLNZ];
     this->min_rownz             = xstore[BASICLU_MIN_ROWNZ];
+    this->bucket_ptr            = xstore[BASICLU_BUCKET_PTR];
 
     /* aliases to user arrays */
     this->Lindex = Li; this->Lvalue = Lx;
@@ -245,6 +250,7 @@ lu_int lu_save(
     xstore[BASICLU_RESIDUAL_TEST]           = this->residual_test;
 
     xstore[BASICLU_MATRIX_NZ]               = this->matrix_nz;
+    xstore[BASICLU_NBUCKETS]                = this->nbuckets;
     xstore[BASICLU_RANK]                    = this->rank;
     xstore[BASICLU_BUMP_SIZE]               = this->bump_size;
     xstore[BASICLU_BUMP_NZ]                 = this->bump_nz;
@@ -267,8 +273,10 @@ lu_int lu_save(
     xstore[BASICLU_MARKER]                  = this->marker;
     xstore[BASICLU_PIVOTLEN]                = this->pivotlen;
     xstore[BASICLU_RANKDEF]                 = this->rankdef;
+    xstore[BASICLU_NCOL_ACTIVE]             = this->ncol_active;
     xstore[BASICLU_MIN_COLNZ]               = this->min_colnz;
     xstore[BASICLU_MIN_ROWNZ]               = this->min_rownz;
+    xstore[BASICLU_BUCKET_PTR]              = this->bucket_ptr;
 
     return status;
 }
@@ -310,6 +318,7 @@ void lu_reset(struct lu *this)
     this->residual_test = 0;
 
     this->matrix_nz = 0;
+    this->nbuckets = 0;
     this->rank = 0;
     this->bump_size = 0;
     this->bump_nz = 0;
@@ -332,8 +341,10 @@ void lu_reset(struct lu *this)
     this->marker = 0;
     this->pivotlen = 0;
     this->rankdef = 0;
+    this->ncol_active = 0;
     this->min_colnz = 1;
     this->min_rownz = 1;
+    this->bucket_ptr = 0;
 
     /*
      * One past the final position in @Wend must hold the file size.
diff --git a/src/lu_internal.h b/src/lu_internal.h
index 46afb0d..cd23345 100644
--- a/src/lu_internal.h
+++ b/src/lu_internal.h
@@ -74,6 +74,7 @@ struct lu
     double residual_test;       /* computed by lu_residual_test() */
 
     lu_int matrix_nz;           /* nz in basis matrix when factorized */
+    lu_int nbuckets;            /* # column buckets; 0 if no buckets given */
     lu_int rank;                /* rank of basis matrix when factorized */
     lu_int bump_size;
     lu_int bump_nz;
@@ -95,10 +96,12 @@ struct lu
     lu_int btran_for_update;    /* >= 0 if BTRAN prepared for update */
     lu_int marker;              /* see @marked, below */
     lu_int pivotlen;            /* length of @pivotcol, @pivotrow; <= 2*m */
-    lu_int rankdef;             /* # columns removed from active submatrix
-                                   because maximum was 0 or < abstol */
+    lu_int rankdef;             /* # columns removed without pivot operation */
+    lu_int ncol_active;         /* # columns in lists */
     lu_int min_colnz;           /* colcount lists 1..min_colnz-1 are empty */
     lu_int min_rownz;           /* rowcount lists 1..min_rownz-1 are empty */
+    lu_int bucket_ptr;          /* 1 <= bucket_ptr <= m is the beginning of the
+                                   next bucket in Lbegin_p; 0 if no buckets */
 
     /* aliases to user arrays */
     lu_int *Lindex, *Uindex, *Windex;
@@ -170,7 +173,7 @@ lu_int lu_setup_bump(
 
 lu_int lu_singletons(
     struct lu *this, const lu_int *Bbegin, const lu_int *Bend, const lu_int *Bi,
-    const double *Bx);
+    const double *Bx, const lu_int *buckets);
 
 void lu_solve_dense(
     struct lu *this, const double *rhs, double *lhs, char trans);
diff --git a/src/lu_list.h b/src/lu_list.h
index 47329fc..468a586 100644
--- a/src/lu_list.h
+++ b/src/lu_list.h
@@ -37,6 +37,7 @@
  *  lu_list_add
  *  lu_list_remove
  *  lu_list_move
+ *  lu_list_move_if
  *  lu_list_swap
  *
  * The methods are defined in this header file as static inline. This header
@@ -117,6 +118,25 @@ static inline void lu_list_move(
 }
 
 
+/* ==========================================================================
+ * lu_list_move_if
+ *
+ * Remove element @elem from its list and add it to list @list. Do nothing if
+ * @elem is not in any list.
+ * ========================================================================== */
+
+static inline void lu_list_move_if(
+    lu_int elem, lu_int list, lu_int *flink, lu_int *blink, lu_int nelem,
+    lu_int *min_list)
+{
+    if (flink[elem] != elem)
+    {
+        lu_list_remove(flink, blink, elem);
+        lu_list_add(elem, list, flink, blink, nelem, min_list);
+    }
+}
+
+
 /* ==========================================================================
  * lu_list_swap
  *
diff --git a/src/lu_markowitz.c b/src/lu_markowitz.c
index 2f590aa..7a13f76 100644
--- a/src/lu_markowitz.c
+++ b/src/lu_markowitz.c
@@ -49,12 +49,12 @@ lu_int lu_markowitz(struct lu *this)
     const double abstol             = this->abstol;
     const double reltol             = this->reltol;
     const lu_int maxsearch          = this->maxsearch;
-    const lu_int search_rows        = this->search_rows;
+    const lu_int search_rows        = this->search_rows && this->bucket_ptr==0;
     const lu_int nz_start           = search_rows ?
         MIN(this->min_colnz, this->min_rownz) : this->min_colnz;
 
     lu_int i, j, pos, where, inext, nz, pivot_row, pivot_col;
-    lu_int nsearch, cheap, found, min_colnz, min_rownz;
+    lu_int ncol_searched, nrow_searched, cheap, found, min_colnz, min_rownz;
     double cmx, x, tol, tic[2];
 
     /* integers for Markowitz cost must be 64 bit to prevent overflow */
@@ -65,7 +65,8 @@ lu_int lu_markowitz(struct lu *this)
     pivot_row = -1;             /* row of best pivot so far */
     pivot_col = -1;             /* col of best pivot so far */
     MC = M*M;                   /* Markowitz cost of best pivot so far */
-    nsearch = 0;                /* count rows/columns searched */
+    ncol_searched = 0;          /* # columns searched for pivot */
+    nrow_searched = 0;          /* # rows searched for pivot */
     min_colnz = -1;             /* minimum col count in active submatrix */
     min_rownz = -1;             /* minimum row count in active submatrix */
     assert(nz_start >= 1);
@@ -82,9 +83,13 @@ lu_int lu_markowitz(struct lu *this)
 
     for (nz = nz_start; nz <= m; nz++)
     {
+        if (ncol_searched == this->ncol_active)
+            goto done;
+
         /* Search columns with nz nonzeros. */
         for (j = colcount_flink[m+nz]; j < m; j = colcount_flink[j])
         {
+            ncol_searched++;
             if (min_colnz == -1)
                 min_colnz = nz;
             assert(Wend[j] - Wbegin[j] == nz);
@@ -115,7 +120,7 @@ lu_int lu_markowitz(struct lu *this)
             }
             /* We have seen at least one eligible pivot in column j. */
             assert(MC < M*M);
-            if (++nsearch >= maxsearch)
+            if (ncol_searched+nrow_searched >= maxsearch)
                 goto done;
         }
         assert(j == m+nz);
@@ -126,6 +131,7 @@ lu_int lu_markowitz(struct lu *this)
         /* Search rows with nz nonzeros. */
         for (i = rowcount_flink[m+nz]; i < m; i = inext)
         {
+            nrow_searched++;
             if (min_rownz == -1)
                 min_rownz = nz;
             /* rowcount_flink[i] might be changed below, so keep a copy */
@@ -171,7 +177,7 @@ lu_int lu_markowitz(struct lu *this)
             else
             {
                 assert(MC < M*M);
-                if (++nsearch >= maxsearch)
+                if (ncol_searched+nrow_searched >= maxsearch)
                     goto done;
             }
         }
@@ -181,7 +187,7 @@ lu_int lu_markowitz(struct lu *this)
 done:
     this->pivot_row = pivot_row;
     this->pivot_col = pivot_col;
-    this->nsearch_pivot += nsearch;
+    this->nsearch_pivot += ncol_searched+nrow_searched;
     if (min_colnz >= 0)
         this->min_colnz = min_colnz;
     if (min_rownz >= 0)
diff --git a/src/lu_pivot.c b/src/lu_pivot.c
index eb264f4..9e11262 100644
--- a/src/lu_pivot.c
+++ b/src/lu_pivot.c
@@ -333,10 +333,11 @@ static lu_int lu_pivot_any(struct lu *this)
         assert(Windex[Wbegin[j]] == pivot_row);
         Wbegin[j]++;
 
-        /* Move column to new list and update min_colnz. */
+        /* Move column to new list and update min_colnz. If the column was not
+           in any list, then it belongs to another bucket and is kept out. */
         nz = Wend[j] - Wbegin[j];
-        lu_list_move(j, nz, colcount_flink, colcount_blink, m,
-                     &this->min_colnz);
+        lu_list_move_if(j, nz, colcount_flink, colcount_blink, m,
+                        &this->min_colnz);
 
         colmax[j] = cmx;
     }
@@ -640,10 +641,11 @@ static lu_int lu_pivot_small(struct lu *this)
         assert(Windex[Wbegin[j]] == pivot_row);
         Wbegin[j]++;
 
-        /* Move column to new list and update min_colnz. */
+        /* Move column to new list and update min_colnz. If the column was not
+           in any list, then it belongs to another bucket and is kept out. */
         nz = Wend[j] - Wbegin[j];
-        lu_list_move(j, nz, colcount_flink, colcount_blink, m,
-                     &this->min_colnz);
+        lu_list_move_if(j, nz, colcount_flink, colcount_blink, m,
+                        &this->min_colnz);
 
         colmax[j] = cmx;
     }
@@ -911,8 +913,8 @@ static lu_int lu_pivot_singleton_col(struct lu *this)
         Windex[where] = Windex[--Wend [j]];
         Wvalue[where] = Wvalue[Wend [j]];
         nz = Wend[j] - Wbegin[j];
-        lu_list_move(j, nz, colcount_flink, colcount_blink, m,
-                     &this->min_colnz);
+        lu_list_move_if(j, nz, colcount_flink, colcount_blink, m,
+                        &this->min_colnz);
         colmax[j] = cmx;
     }
     assert(found);
@@ -1083,8 +1085,8 @@ static lu_int lu_pivot_doubleton_col(struct lu *this)
 
                 /* Decrease column count. */
                 nz = end - Wbegin[j];
-                lu_list_move(j, nz, colcount_flink, colcount_blink, m,
-                             &this->min_colnz);
+                lu_list_move_if(j, nz, colcount_flink, colcount_blink, m,
+                                &this->min_colnz);
             }
         }
         else
@@ -1113,8 +1115,8 @@ static lu_int lu_pivot_doubleton_col(struct lu *this)
 
             /* Decrease column count. */
             nz = Wend[j] - Wbegin[j];
-            lu_list_move(j, nz, colcount_flink, colcount_blink, m,
-                         &this->min_colnz);
+            lu_list_move_if(j, nz, colcount_flink, colcount_blink, m,
+                            &this->min_colnz);
         }
         colmax[j] = cmx;
     }
@@ -1247,5 +1249,5 @@ static void lu_remove_col(struct lu *this, lu_int j)
     /* Remove column j from column file. */
     colmax[j] = 0.0;
     Wend[j] = cbeg;
-    lu_list_move(j, 0, colcount_flink, colcount_blink, m, &this->min_colnz);
+    lu_list_move_if(j, 0, colcount_flink, colcount_blink, m, &this->min_colnz);
 }
diff --git a/src/lu_setup_bump.c b/src/lu_setup_bump.c
index 55d2ee4..6b77e60 100644
--- a/src/lu_setup_bump.c
+++ b/src/lu_setup_bump.c
@@ -63,6 +63,7 @@ lu_int lu_setup_bump(
     const double abstol     = this->abstol;
     const lu_int pad        = this->pad;
     const double stretch    = this->stretch;
+    const lu_int have_buckets = this->bucket_ptr != 0;
     lu_int *colcount_flink  = this->colcount_flink;
     lu_int *colcount_blink  = this->colcount_blink;
     lu_int *rowcount_flink  = this->rowcount_flink;
@@ -110,6 +111,7 @@ lu_int lu_setup_bump(
      * Build columnwise storage. Build row counts in iwork0.
      */
     lu_list_init(colcount_flink, colcount_blink, m, m+2, &min_colnz);
+    this->ncol_active = 0;
     put = 0;
     for (j = 0; j < m; j++)
     {
@@ -129,14 +131,13 @@ lu_int lu_setup_bump(
         {
             /* Leave column of active submatrix empty. */
             colmax[j] = 0.0;
-            lu_list_add(j, 0, colcount_flink, colcount_blink, m, &min_colnz);
             bump_nz -= cnz;
+            cnz = 0;
         }
         else
         {
             /* Copy column into active submatrix. */
             colmax[j] = cmx;
-            lu_list_add(j, cnz, colcount_flink, colcount_blink, m, &min_colnz);
             Wbegin[j] = put;
             for (pos = Bbegin[j]; pos < Bend[j]; pos++)
             {
@@ -152,6 +153,11 @@ lu_int lu_setup_bump(
             /* reappend line to list end */
             lu_list_move(j, 0, Wflink, Wblink, 2*m, NULL);
         }
+        if (!have_buckets)
+        {
+            lu_list_add(j, cnz, colcount_flink, colcount_blink, m, &min_colnz);
+            this->ncol_active++;
+        }
     }
 
     /*
diff --git a/src/lu_singletons.c b/src/lu_singletons.c
index c6aae2b..bad5cbd 100644
--- a/src/lu_singletons.c
+++ b/src/lu_singletons.c
@@ -140,7 +140,7 @@ static lu_int singleton_rows
 
 lu_int lu_singletons(
     struct lu *this, const lu_int *Bbegin, const lu_int *Bend, const lu_int *Bi,
-    const double *Bx)
+    const double *Bx, const lu_int *buckets)
 {
     const lu_int m      = this->m;
     const lu_int Lmem   = this->Lmem;
@@ -228,7 +228,7 @@ lu_int lu_singletons(
             i = Bi[pos];
             put = iwork1[i]++;
             Bti[put] = j;
-            Btx[put] = Bx [pos];
+            Btx[put] = Bx[pos];
             if (put > Btp[i] && Bti[put-1] == j)
                 ok = 0;
         }
@@ -236,20 +236,66 @@ lu_int lu_singletons(
     if (!ok)
         return BASICLU_ERROR_invalid_argument;
 
-    /* ---------------- */
-    /* Pivot singletons */
-    /* ---------------- */
+    /* If buckets are given, check that all of them are in range [0,m-1].
+     * Then sort column indices 0..m-1 by order of nondecreasing bucket value
+     * in Lbegin_p[1..m]. Use qinv as workspace for bucket sort, taking O(m)
+     * time. If column index j is the first one in its bucket, store it flipped.
+     */
+    if (buckets)
+    {
+        memset(qinv, 0, m*sizeof(lu_int));
+        ok = 1;
+        for (j = 0; j < m; j++)
+        {
+            if (buckets[j] < 0 || buckets[j] >= m)
+                ok = 0;
+            else
+                qinv[buckets[j]]++;
+        }
+        if (!ok)
+            return BASICLU_ERROR_invalid_argument;
+        for (i = 0; i < m-1; i++)
+            qinv[i+1] += qinv[i];
+        assert(qinv[m-1] == m);
+        for (j = 0; j < m; j++)
+            Lbegin_p[qinv[buckets[j]]--] = j;
+#ifndef NDEBUG
+        for (pos = 2; pos <= m; pos++)
+        {
+            assert(buckets[Lbegin_p[pos]] >= buckets[Lbegin_p[pos-1]]);
+        }
+#endif
+        this->nbuckets = 0;
+        for (pos = 1; pos <= m; pos++)
+        {
+            if (pos == 1 || buckets[Lbegin_p[pos]] != buckets[Lbegin_p[pos-1]])
+            {
+                Lbegin_p[pos] = -Lbegin_p[pos]-1;
+                this->nbuckets++;
+            }
+        }
+        this->bucket_ptr = 1;
+    }
 
     /* No pivot rows or pivot columns so far. */
     for (i = 0; i < m; i++)
         pinv[i] = -1;
     for (j = 0; j < m; j++)
         qinv[j] = -1;
+    Lbegin_p[0] = Ubegin[0] = 0;
+    this->matrix_nz = Bnz;
+
+    /* If buckets are given, then do not eliminate singletons. */
+    if (buckets)
+        return BASICLU_OK;
+
+    /* ---------------- */
+    /* Pivot singletons */
+    /* ---------------- */
 
+    rank = 0;
     if (nzbias >= 0)            /* put more in U */
     {
-        Lbegin_p[0] = Ubegin[0] = rank = 0;
-
         rank = singleton_cols(m, Bbegin, Bend, Bi, Bx, Btp, Bti, Btx,
                               Ubegin, Uindex, Uvalue, Lbegin_p, Lindex, Lvalue,
                               col_pivot, pinv, qinv, iwork1, iwork2, rank,
@@ -262,8 +308,6 @@ lu_int lu_singletons(
     }
     else                        /* put more in L */
     {
-        Lbegin_p[0] = Ubegin[0] = rank = 0;
-
         rank = singleton_rows(m, Bbegin, Bend, Bi, Bx, Btp, Bti, Btx,
                               Ubegin, Uindex, Uvalue, Lbegin_p, Lindex, Lvalue,
                               col_pivot, pinv, qinv, iwork1, iwork2, rank,
@@ -283,7 +327,6 @@ lu_int lu_singletons(
         if (qinv[j] < 0)
             qinv[j] = -1;
 
-    this->matrix_nz = Bnz;
     this->rank = rank;
     this->time_singletons = lu_toc(tic);
     return BASICLU_OK;