spbkp.c

Go to the documentation of this file.

#include <../../nrnconf.h>

/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**               Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/


/*
  Sparse matrix Bunch--Kaufman--Parlett factorisation and solve
  Radical revision started Thu 05th Nov 1992, 09:36:12 AM
  to use Karen George's suggestion of leaving the the row elements unordered
  Radical revision completed Mon 07th Dec 1992, 10:59:57 AM
*/

static  char    rcsid[] = "spbkp.c,v 1.1 1997/12/04 17:55:50 hines Exp";

#include    <stdio.h>
#include        "sparse2.h"
#include    <math.h>


#ifdef MALLOCDECL
#include <malloc.h>
#endif

#define alpha   0.6403882032022076 /* = (1+sqrt(17))/8 */


#define btos(x) ((x) ? "TRUE" : "FALSE")

/* assume no use of sqr() uses side-effects */
#define sqr(x)  ((x)*(x))

/* unord_get_idx -- returns index (encoded if entry not allocated)
    of the element of row r with column j
    -- uses linear search */
int unord_get_idx(r,j)
SPROW   *r;
int j;
{
    int     idx;
    row_elt *e;

    if ( ! r || ! r->elt )
    error(E_NULL,"unord_get_idx");
    for ( idx = 0, e = r->elt; idx < r->len; idx++, e++ )
    if ( e->col == j )
        break;
    if ( idx >= r->len )
    return -(r->len+2);
    else
    return idx;
}

/* unord_get_val -- returns value of the (i,j) entry of A
    -- same assumptions as unord_get_idx() */
double  unord_get_val(A,i,j)
SPMAT   *A;
int i, j;
{
    SPROW   *r;
    int     idx;

    if ( ! A )
    error(E_NULL,"unord_get_val");
    if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
    error(E_BOUNDS,"unord_get_val");

    r = &(A->row[i]);
    idx = unord_get_idx(r,j);
    if ( idx < 0 )
    return 0.0;
    else
    return r->elt[idx].val;
}

        
/* bkp_swap_elt -- swaps the (i,j) with the (k,l) entry of sparse matrix
    -- either or both of the entries may be unallocated */
static SPMAT    *bkp_swap_elt(A,i1,j1,idx1,i2,j2,idx2)
SPMAT   *A;
int i1, j1, idx1, i2, j2, idx2;
{
    int     tmp_row, tmp_idx;
    SPROW   *r1, *r2;
    row_elt *e1, *e2;
    Real    tmp;

    if ( ! A )
    error(E_NULL,"bkp_swap_elt");

    if ( i1 < 0 || j1 < 0 || i2 < 0 || j2 < 0 ||
     i1 >= A->m || j1 >= A->n || i2 >= A->m || j2 >= A->n )
    {
    error(E_BOUNDS,"bkp_swap_elt");
    }

    if ( i1 == i2 && j1 == j2 )
    return A;
    if ( idx1 < 0 && idx2 < 0 ) /* neither allocated */
    return A;

    r1 = &(A->row[i1]);     r2 = &(A->row[i2]);
    /* if ( idx1 >= r1->len || idx2 >= r2->len )
    error(E_BOUNDS,"bkp_swap_elt"); */
    if ( idx1 < 0 ) /* assume not allocated */
    {
    idx1 = r1->len;
    if ( idx1 >= r1->maxlen )
    {    tracecatch(sprow_xpd(r1,2*r1->maxlen+1,TYPE_SPMAT),
            "bkp_swap_elt");    }
    r1->len = idx1+1;
    r1->elt[idx1].col = j1;
    r1->elt[idx1].val = 0.0;
    /* now patch up column access path */
    tmp_row = -1;   tmp_idx = j1;
    chase_col(A,j1,&tmp_row,&tmp_idx,i1-1);

    if ( tmp_row < 0 )
    {
        r1->elt[idx1].nxt_row = A->start_row[j1];
        r1->elt[idx1].nxt_idx = A->start_idx[j1];
        A->start_row[j1] = i1;
        A->start_idx[j1] = idx1;
    }
    else
    {
        row_elt *tmp_e;

        tmp_e = &(A->row[tmp_row].elt[tmp_idx]);
        r1->elt[idx1].nxt_row = tmp_e->nxt_row;
        r1->elt[idx1].nxt_idx = tmp_e->nxt_idx;
        tmp_e->nxt_row = i1;
        tmp_e->nxt_idx = idx1;
    }
    }
    else if ( r1->elt[idx1].col != j1 )
    error(E_INTERN,"bkp_swap_elt");
    if ( idx2 < 0 )
    {
    idx2 = r2->len;
    if ( idx2 >= r2->maxlen )
    {    tracecatch(sprow_xpd(r2,2*r2->maxlen+1,TYPE_SPMAT),
            "bkp_swap_elt");    }

    r2->len = idx2+1;
    r2->elt[idx2].col = j2;
    r2->elt[idx2].val = 0.0;
    /* now patch up column access path */
    tmp_row = -1;   tmp_idx = j2;
    chase_col(A,j2,&tmp_row,&tmp_idx,i2-1);
    if ( tmp_row < 0 )
    {
        r2->elt[idx2].nxt_row = A->start_row[j2];
        r2->elt[idx2].nxt_idx = A->start_idx[j2];
        A->start_row[j2] = i2;
        A->start_idx[j2] = idx2;
    }
    else
    {
        row_elt *tmp_e;

        tmp_e = &(A->row[tmp_row].elt[tmp_idx]);
        r2->elt[idx2].nxt_row = tmp_e->nxt_row;
        r2->elt[idx2].nxt_idx = tmp_e->nxt_idx;
        tmp_e->nxt_row = i2;
        tmp_e->nxt_idx = idx2;
    }
    }
    else if ( r2->elt[idx2].col != j2 )
    error(E_INTERN,"bkp_swap_elt");

    e1 = &(r1->elt[idx1]);  e2 = &(r2->elt[idx2]);

    tmp = e1->val;
    e1->val = e2->val;
    e2->val = tmp;

    return A;
}

/* bkp_bump_col -- bumps row and idx to next entry in column j */
row_elt *bkp_bump_col(A, j, row, idx)
SPMAT   *A;
int j, *row, *idx;
{
    SPROW   *r;
    row_elt *e;

    if ( *row < 0 )
    {
    *row = A->start_row[j];
    *idx = A->start_idx[j];
    }
    else
    {
    r = &(A->row[*row]);
    e = &(r->elt[*idx]);
    if ( e->col != j )
        error(E_INTERN,"bkp_bump_col");
    *row = e->nxt_row;
    *idx = e->nxt_idx;
    }
    if ( *row < 0 )
    return (row_elt *)NULL;
    else
    return &(A->row[*row].elt[*idx]);
}

/* bkp_interchange -- swap rows/cols i and j (symmetric pivot)
    -- uses just the upper triangular part */
SPMAT   *bkp_interchange(A, i1, i2)
SPMAT   *A;
int i1, i2;
{
    int     tmp_row, tmp_idx;
    int     row1, row2, idx1, idx2, tmp_row1, tmp_idx1, tmp_row2, tmp_idx2;
    SPROW   *r1, *r2;
    row_elt *e1, *e2;
    IVEC    *done_list = IVNULL;

    if ( ! A )
    error(E_NULL,"bkp_interchange");
    if ( i1 < 0 || i1 >= A->n || i2 < 0 || i2 >= A->n )
    error(E_BOUNDS,"bkp_interchange");
    if ( A->m != A->n )
    error(E_SQUARE,"bkp_interchange");

    if ( i1 == i2 )
    return A;
    if ( i1 > i2 )
    {   tmp_idx = i1;   i1 = i2;    i2 = tmp_idx;   }

    done_list = iv_resize(done_list,A->n);
    for ( tmp_idx = 0; tmp_idx < A->n; tmp_idx++ )
    done_list->ive[tmp_idx] = FALSE;
    row1 = -1;      idx1 = i1;
    row2 = -1;      idx2 = i2;
    e1 = bkp_bump_col(A,i1,&row1,&idx1);
    e2 = bkp_bump_col(A,i2,&row2,&idx2);

    while ( (row1 >= 0 && row1 < i1) || (row2 >= 0 && row2 < i1) )
    /* Note: "row2 < i1" not "row2 < i2" as we must stop before the
       "knee bend" */
    {
    if ( row1 >= 0 && row1 < i1 && ( row1 < row2 || row2 < 0 ) )
    {
        tmp_row1 = row1;    tmp_idx1 = idx1;
        e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
        if ( ! done_list->ive[row1] )
        {
        if ( row1 == row2 )
            bkp_swap_elt(A,row1,i1,idx1,row1,i2,idx2);
        else
            bkp_swap_elt(A,row1,i1,idx1,row1,i2,-1);
        done_list->ive[row1] = TRUE;
        }
        row1 = tmp_row1;    idx1 = tmp_idx1;
    }
    else if ( row2 >= 0 && row2 < i1 && ( row2 < row1 || row1 < 0 ) )
    {
        tmp_row2 = row2;    tmp_idx2 = idx2;
        e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
        if ( ! done_list->ive[row2] )
        {
        if ( row1 == row2 )
            bkp_swap_elt(A,row2,i1,idx1,row2,i2,idx2);
        else
            bkp_swap_elt(A,row2,i1,-1,row2,i2,idx2);
        done_list->ive[row2] = TRUE;
        }
        row2 = tmp_row2;    idx2 = tmp_idx2;
    }
    else if ( row1 == row2 )
    {
        tmp_row1 = row1;    tmp_idx1 = idx1;
        e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
        tmp_row2 = row2;    tmp_idx2 = idx2;
        e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
        if ( ! done_list->ive[row1] )
        {
        bkp_swap_elt(A,row1,i1,idx1,row2,i2,idx2);
        done_list->ive[row1] = TRUE;
        }
        row1 = tmp_row1;    idx1 = tmp_idx1;
        row2 = tmp_row2;    idx2 = tmp_idx2;
    }
    }

    /* ensure we are **past** the first knee */
    while ( row2 >= 0 && row2 <= i1 )
    e2 = bkp_bump_col(A,i2,&row2,&idx2);

    /* at/after 1st "knee bend" */
    r1 = &(A->row[i1]);
    idx1 = 0;
    e1 = &(r1->elt[idx1]);
    while ( row2 >= 0 && row2 < i2 )
    {
    /* used for update of e2 at end of loop */
    tmp_row = row2; tmp_idx = idx2;
    if ( ! done_list->ive[row2] )
    {
        r2 = &(A->row[row2]);
        bkp_bump_col(A,i2,&tmp_row,&tmp_idx);
        done_list->ive[row2] = TRUE;
        tmp_idx1 = unord_get_idx(r1,row2);
        tracecatch(bkp_swap_elt(A,row2,i2,idx2,i1,row2,tmp_idx1),
               "bkp_interchange");
    }

    /* update e1 and e2 */
    row2 = tmp_row; idx2 = tmp_idx;
    e2 = ( row2 >= 0 ) ? &(A->row[row2].elt[idx2]) : (row_elt *)NULL;
    }

    idx1 = 0;
    e1 = r1->elt;
    while ( idx1 < r1->len )
    {
    if ( e1->col >= i2 || e1->col <= i1 )
    {
        idx1++;
        e1++;
        continue;
    }
    if ( ! done_list->ive[e1->col] )
    {
        tmp_idx2 = unord_get_idx(&(A->row[e1->col]),i2);
        tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,e1->col,i2,tmp_idx2),
               "bkp_interchange");
        done_list->ive[e1->col] = TRUE;
    }
    idx1++;
    e1++;
    }

    /* at/after 2nd "knee bend" */
    idx1 = 0;
    e1 = &(r1->elt[idx1]);
    r2 = &(A->row[i2]);
    idx2 = 0;
    e2 = &(r2->elt[idx2]);
    while ( idx1 < r1->len )
    {
    if ( e1->col <= i2 )
    {
        idx1++; e1++;
        continue;
    }
    if ( ! done_list->ive[e1->col] )
    {
        tmp_idx2 = unord_get_idx(r2,e1->col);
        tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,i2,e1->col,tmp_idx2),
               "bkp_interchange");
        done_list->ive[e1->col] = TRUE;
    }
    idx1++; e1++;
    }

    idx2 = 0;   e2 = r2->elt;
    while ( idx2 < r2->len )
    {
    if ( e2->col <= i2 )
    {
        idx2++; e2++;
        continue;
    }
    if ( ! done_list->ive[e2->col] )
    {
        tmp_idx1 = unord_get_idx(r1,e2->col);
        tracecatch(bkp_swap_elt(A,i2,e2->col,idx2,i1,e2->col,tmp_idx1),
               "bkp_interchange");
        done_list->ive[e2->col] = TRUE;
    }
    idx2++; e2++;
    }

    /* now interchange the digonal entries! */
    idx1 = unord_get_idx(&(A->row[i1]),i1);
    idx2 = unord_get_idx(&(A->row[i2]),i2);
    if ( idx1 >= 0 || idx2 >= 0 )
    {
    tracecatch(bkp_swap_elt(A,i1,i1,idx1,i2,i2,idx2),
           "bkp_interchange");
    }

    return A;
}


/* iv_min -- returns minimum of an integer vector
   -- sets index to the position in iv if index != NULL */
int iv_min(iv,index)
IVEC    *iv;
int *index;
{
    int     i, i_min, min_val, tmp;
    
    if ( ! iv ) 
    error(E_NULL,"iv_min");
    if ( iv->dim <= 0 )
    error(E_SIZES,"iv_min");
    i_min = 0;
    min_val = iv->ive[0];
    for ( i = 1; i < iv->dim; i++ )
    {
    tmp = iv->ive[i];
    if ( tmp < min_val )
    {
        min_val = tmp;
        i_min = i;
    }
    }
    
    if ( index != (int *)NULL )
    *index = i_min;
    
    return min_val;
}

/* max_row_col -- returns max { |A[j][k]| : k >= i, k != j, k != l } given j
    using symmetry and only the upper triangular part of A */
static double max_row_col(A,i,j,l)
SPMAT   *A;
int i, j, l;
{
    int     row_num, idx;
    SPROW   *r;
    row_elt *e;
    Real    max_val, tmp;

    if ( ! A )
    error(E_NULL,"max_row_col");
    if ( i < 0 || i > A->n || j < 0 || j >= A->n )
    error(E_BOUNDS,"max_row_col");

    max_val = 0.0;

    idx = unord_get_idx(&(A->row[i]),j);
    if ( idx < 0 )
    {
    row_num = -1;   idx = j;
    e = chase_past(A,j,&row_num,&idx,i);
    }
    else
    {
    row_num = i;
    e = &(A->row[i].elt[idx]);
    }
    while ( row_num >= 0 && row_num < j )
    {
    if ( row_num != l )
    {
        tmp = fabs(e->val);
        if ( tmp > max_val )
        max_val = tmp;
    }
    e = bump_col(A,j,&row_num,&idx);
    }
    r = &(A->row[j]);
    for ( idx = 0, e = r->elt; idx < r->len; idx++, e++ )
    {
    if ( e->col > j && e->col != l )
    {
        tmp = fabs(e->val);
        if ( tmp > max_val )
        max_val = tmp;
    }
    }

    return max_val;
}

/* nonzeros -- counts non-zeros in A */
static int  nonzeros(A)
SPMAT   *A;
{
    int     cnt, i;

    if ( ! A )
    return 0;
    cnt = 0;
    for ( i = 0; i < A->m; i++ )
    cnt += A->row[i].len;

    return cnt;
}

/* chk_col_access -- for spBKPfactor()
    -- checks that column access path is OK */
int chk_col_access(A)
SPMAT   *A;
{
    int     cnt_nz, j, row, idx;
    SPROW   *r;
    row_elt *e;

    if ( ! A )
    error(E_NULL,"chk_col_access");

    /* count nonzeros as we go down columns */
    cnt_nz = 0;
    for ( j = 0; j < A->n; j++ )
    {
    row = A->start_row[j];
    idx = A->start_idx[j];
    while ( row >= 0 )
    {
        if ( row >= A->m || idx < 0 )
        return FALSE;
        r = &(A->row[row]);
        if ( idx >= r->len )
        return FALSE;
        e = &(r->elt[idx]);
        if ( e->nxt_row >= 0 && e->nxt_row <= row )
        return FALSE;
        row = e->nxt_row;
        idx = e->nxt_idx;
        cnt_nz++;
    }
    }

    if ( cnt_nz != nonzeros(A) )
    return FALSE;
    else
    return TRUE;
}

/* col_cmp -- compare two columns -- for sorting rows using qsort() */
static int  col_cmp(e1,e2)
row_elt *e1, *e2;
{
    return e1->col - e2->col;
}

/* spBKPfactor -- sparse Bunch-Kaufman-Parlett factorisation of A in-situ
   -- A is factored into the form P'AP = MDM' where 
   P is a permutation matrix, M lower triangular and D is block
   diagonal with blocks of size 1 or 2
   -- P is stored in pivot; blocks[i]==i iff D[i][i] is a block */
SPMAT   *spBKPfactor(A,pivot,blocks,tol)
SPMAT   *A;
PERM    *pivot, *blocks;
double  tol;
{
    int     i, j, k, l, n, onebyone=0, r;
    int     idx, idx1, idx_piv;
    int     row_num;
    int     best_deg=0, best_j, best_l, best_cost, mark_cost, deg, deg_j,
            deg_l, ignore_deg;
    int     list_idx, list_idx2, old_list_idx;
    SPROW   *row, *r_piv, *r1_piv;
    row_elt *e, *e1;
    Real    aii, aip1, aip1i;
    Real    det, max_j, max_l, s, t;
    static IVEC *scan_row = IVNULL, *scan_idx = IVNULL, *col_list = IVNULL,
        *tmp_iv = IVNULL;
    static IVEC *deg_list = IVNULL;
    static IVEC *orig_idx = IVNULL, *orig1_idx = IVNULL;
    static PERM *order = PNULL;

    if ( ! A || ! pivot || ! blocks )
    error(E_NULL,"spBKPfactor");
    if ( A->m != A->n )
    error(E_SQUARE,"spBKPfactor");
    if ( A->m != pivot->size || pivot->size != blocks->size )
    error(E_SIZES,"spBKPfactor");
    if ( tol <= 0.0 || tol > 1.0 )
    error(E_RANGE,"spBKPfactor");
    
    n = A->n;
    
    px_ident(pivot);    px_ident(blocks);
    sp_col_access(A);   sp_diag_access(A);
    ignore_deg = FALSE;

    deg_list = iv_resize(deg_list,n);
    order = px_resize(order,n);
    MEM_STAT_REG(deg_list,TYPE_IVEC);
    MEM_STAT_REG(order,TYPE_PERM);

    scan_row = iv_resize(scan_row,5);
    scan_idx = iv_resize(scan_idx,5);
    col_list = iv_resize(col_list,5);
    orig_idx = iv_resize(orig_idx,5);
    orig_idx = iv_resize(orig1_idx,5);
    orig_idx = iv_resize(tmp_iv,5);
    MEM_STAT_REG(scan_row,TYPE_IVEC);
    MEM_STAT_REG(scan_idx,TYPE_IVEC);
    MEM_STAT_REG(col_list,TYPE_IVEC);
    MEM_STAT_REG(orig_idx,TYPE_IVEC);
    MEM_STAT_REG(orig1_idx,TYPE_IVEC);
    MEM_STAT_REG(tmp_iv,TYPE_IVEC);

    for ( i = 0; i < n-1; i = onebyone ? i+1 : i+2 )
    {
    /* now we want to use a Markowitz-style selection rule for
       determining which rows to swap and whether to use
       1x1 or 2x2 pivoting */

    /* get list of degrees of nodes */
    deg_list = iv_resize(deg_list,n-i);
    if ( ! ignore_deg )
        for ( j = i; j < n; j++ )
        deg_list->ive[j-i] = 0;
    else
    {
        for ( j = i; j < n; j++ )
        deg_list->ive[j-i] = 1;
        if ( i < n )
        deg_list->ive[0] = 0;
    }
    order = px_resize(order,n-i);
    px_ident(order);

    if ( ! ignore_deg )
    {
        for ( j = i; j < n; j++ )
        {
        /* idx = sprow_idx(&(A->row[j]),j+1); */
        /* idx = fixindex(idx); */
        idx = 0;
        row = &(A->row[j]);
        e = &(row->elt[idx]);
        /* deg_list->ive[j-i] += row->len - idx; */
        for ( ; idx < row->len; idx++, e++ )
            if ( e->col >= i )
            deg_list->ive[e->col - i]++;
        }
        /* now deg_list[k] == degree of node k+i */
        
        /* now sort them into increasing order */
        iv_sort(deg_list,order);
        /* now deg_list[idx] == degree of node i+order[idx] */
    }

    /* now we can chase through the nodes in order of increasing
       degree, picking out the ones that satisfy our stability
       criterion */
    list_idx = 0;   r = -1;
    best_j = best_l = -1;
    for ( deg = 0; deg <= n; deg++ )
    {
        Real    ajj, all, ajl;

        if ( list_idx >= deg_list->dim )
        break;  /* That's all folks! */
        old_list_idx = list_idx;
        while ( list_idx < deg_list->dim &&
            deg_list->ive[list_idx] <= deg )
        {
        j = i+order->pe[list_idx];
        if ( j < i )
            continue;
        /* can we use row/col j for a 1 x 1 pivot? */
        /* find max_j = max_{k>=i} {|A[k][j]|,|A[j][k]|} */
        ajj = fabs(unord_get_val(A,j,j));
        if ( ajj == 0.0 )
        {
            list_idx++;
            continue;   /* can't use this for 1 x 1 pivot */
        }

        max_j = max_row_col(A,i,j,-1);
        if ( ajj >= tol/* *alpha */ *max_j )
        {
            onebyone = TRUE;
            best_j = j;
            best_deg = deg_list->ive[list_idx];
            break;
        }
        list_idx++;
        }
        if ( best_j >= 0 )
        break;
        best_cost = 2*n;    /* > any possible Markowitz cost (bound) */
        best_j = best_l = -1;
        list_idx = old_list_idx;
        while ( list_idx < deg_list->dim &&
            deg_list->ive[list_idx] <= deg )
        {
        j = i+order->pe[list_idx];
        ajj = fabs(unord_get_val(A,j,j));
        for ( list_idx2 = 0; list_idx2 < list_idx; list_idx2++ )
        {
            deg_j = deg;
            deg_l = deg_list->ive[list_idx2];
            l = i+order->pe[list_idx2];
            if ( l < i )
            continue;
            /* try using rows/cols (j,l) for a 2 x 2 pivot block */
            all = fabs(unord_get_val(A,l,l));
            ajl = ( j > l ) ? fabs(unord_get_val(A,l,j)) :
                       fabs(unord_get_val(A,j,l));
            det = fabs(ajj*all - ajl*ajl);
            if ( det == 0.0 )
            continue;
            max_j = max_row_col(A,i,j,l);
            max_l = max_row_col(A,i,l,j);
            if ( tol*(all*max_j+ajl*max_l) < det &&
             tol*(ajl*max_j+ajj*max_l) < det )
            {
            /* acceptably stable 2 x 2 pivot */
            /* this is actually an overestimate of the
               Markowitz cost for choosing (j,l) */
            mark_cost = (ajj == 0.0) ?
                ((all == 0.0) ? deg_j+deg_l : deg_j+2*deg_l) :
                ((all == 0.0) ? 2*deg_j+deg_l :
                 2*(deg_j+deg_l));
            if ( mark_cost < best_cost )
            {
                onebyone = FALSE;
                best_cost = mark_cost;
                best_j = j;
                best_l = l;
                best_deg = deg_j;
            }
            }
        }
        list_idx++;
        }
        if ( best_j >= 0 )
        break;
    }

    if ( best_deg > (int)floor(0.8*(n-i)) )
        ignore_deg = TRUE;

    /* now do actual interchanges */
    if ( best_j >= 0 && onebyone )
    {
        bkp_interchange(A,i,best_j);
        px_transp(pivot,i,best_j);
    }
    else if ( best_j >= 0 && best_l >= 0 && ! onebyone )
    {
        if ( best_j == i || best_j == i+1 )
        {
        if ( best_l == i || best_l == i+1 )
        {
            /* no pivoting, but must update blocks permutation */
            px_transp(blocks,i,i+1);
            goto dopivot;
        }
        bkp_interchange(A,(best_j == i) ? i+1 : i,best_l);
        px_transp(pivot,(best_j == i) ? i+1 : i,best_l);
        }
        else if ( best_l == i || best_l == i+1 )
        {
        bkp_interchange(A,(best_l == i) ? i+1 : i,best_j);
        px_transp(pivot,(best_l == i) ? i+1 : i,best_j);
        }
        else /* best_j & best_l outside i, i+1 */
        {
        if ( i != best_j )
        {
            bkp_interchange(A,i,best_j);
            px_transp(pivot,i,best_j);
        }
        if ( i+1 != best_l )
        {
            bkp_interchange(A,i+1,best_l);
            px_transp(pivot,i+1,best_l);
        }
        }
    }
    else    /* can't pivot &/or nothing to pivot */
        continue;

    /* update blocks permutation */
    if ( ! onebyone )
        px_transp(blocks,i,i+1);

    dopivot:
    if ( onebyone )
    {
        int     idx_j, idx_k, s_idx, s_idx2;
        row_elt *e_ij, *e_ik;

        r_piv = &(A->row[i]);
        idx_piv = unord_get_idx(r_piv,i);
        /* if idx_piv < 0 then aii == 0 and no pivoting can be done;
           -- this means that we should continue to the next iteration */
        if ( idx_piv < 0 )
        continue;
        aii = r_piv->elt[idx_piv].val;
        if ( aii == 0.0 )
        continue;

        /* for ( j = i+1; j < n; j++ )  { ... pivot step ... } */
        /* initialise scan_... etc for the 1 x 1 pivot */
        scan_row = iv_resize(scan_row,r_piv->len);
        scan_idx = iv_resize(scan_idx,r_piv->len);
        col_list = iv_resize(col_list,r_piv->len);
        orig_idx = iv_resize(orig_idx,r_piv->len);
        row_num = i;    s_idx = idx = 0;
        e = &(r_piv->elt[idx]);
        for ( idx = 0; idx < r_piv->len; idx++, e++ )
        {
        if ( e->col < i )
            continue;
        scan_row->ive[s_idx] = i;
        scan_idx->ive[s_idx] = idx;
        orig_idx->ive[s_idx] = idx;
        col_list->ive[s_idx] = e->col;
        s_idx++;
        }
        scan_row = iv_resize(scan_row,s_idx);
        scan_idx = iv_resize(scan_idx,s_idx);
        col_list = iv_resize(col_list,s_idx);
        orig_idx = iv_resize(orig_idx,s_idx);

        order = px_resize(order,scan_row->dim);
        px_ident(order);
        iv_sort(col_list,order);

        tmp_iv = iv_resize(tmp_iv,scan_row->dim);
        for ( idx = 0; idx < order->size; idx++ )
        tmp_iv->ive[idx] = scan_idx->ive[order->pe[idx]];
        iv_copy(tmp_iv,scan_idx);
        for ( idx = 0; idx < order->size; idx++ )
        tmp_iv->ive[idx] = scan_row->ive[order->pe[idx]];
        iv_copy(tmp_iv,scan_row);
        for ( idx = 0; idx < scan_row->dim; idx++ )
        tmp_iv->ive[idx] = orig_idx->ive[order->pe[idx]];
        iv_copy(tmp_iv,orig_idx);

        /* now do actual pivot */
        /* for ( j = i+1; j < n-1; j++ ) .... */

        for ( s_idx = 0; s_idx < scan_row->dim; s_idx++ )
        {
        idx_j = orig_idx->ive[s_idx];
        if ( idx_j < 0 )
            error(E_INTERN,"spBKPfactor");
        e_ij = &(r_piv->elt[idx_j]);
        j = e_ij->col;
        if ( j < i+1 )
            continue;
        scan_to(A,scan_row,scan_idx,col_list,j);

        /* compute multiplier */
        t = e_ij->val / aii;

        /* for ( k = j; k < n; k++ ) { .... update A[j][k] .... } */
        /* this is the row in which pivoting is done */
        row = &(A->row[j]);
        for ( s_idx2 = s_idx; s_idx2 < scan_row->dim; s_idx2++ )
        {
            idx_k = orig_idx->ive[s_idx2];
            e_ik = &(r_piv->elt[idx_k]);
            k = e_ik->col;
            /* k >= j since col_list has been sorted */

            if ( scan_row->ive[s_idx2] == j )
            {   /* no fill-in -- can be done directly */
            idx = scan_idx->ive[s_idx2];
            /* idx = sprow_idx2(row,k,idx); */
            row->elt[idx].val -= t*e_ik->val;
            }
            else
            {   /* fill-in -- insert entry & patch column */
            int old_row, old_idx;
            row_elt *old_e, *new_e;

            old_row = scan_row->ive[s_idx2];
            old_idx = scan_idx->ive[s_idx2];
            /* old_idx = sprow_idx2(&(A->row[old_row]),k,old_idx); */

            if ( old_idx < 0 )
                error(E_INTERN,"spBKPfactor");
            /* idx = sprow_idx(row,k); */
            /* idx = fixindex(idx); */
            idx = row->len;

            /* sprow_set_val(row,k,-t*e_ik->val); */
            if ( row->len >= row->maxlen )
            { tracecatch(sprow_xpd(row,2*row->maxlen+1,TYPE_SPMAT),
                     "spBKPfactor");        }

            row->len = idx+1;

            new_e = &(row->elt[idx]);
            new_e->val = -t*e_ik->val;
            new_e->col = k;

            old_e = &(A->row[old_row].elt[old_idx]);
            new_e->nxt_row = old_e->nxt_row;
            new_e->nxt_idx = old_e->nxt_idx;
            old_e->nxt_row = j;
            old_e->nxt_idx = idx;
            }
        }
        e_ij->val = t;
        }
    }
    else /* onebyone == FALSE */
    {   /* do 2 x 2 pivot */
        int idx_k, idx1_k, s_idx, s_idx2;
        int old_col;
        row_elt *e_tmp;

        r_piv = &(A->row[i]);
        idx_piv = unord_get_idx(r_piv,i);
        aii = aip1i = 0.0;
        e_tmp = r_piv->elt;
        for ( idx_piv = 0; idx_piv < r_piv->len; idx_piv++, e_tmp++ )
        if ( e_tmp->col == i )
            aii = e_tmp->val;
            else if ( e_tmp->col == i+1 )
            aip1i = e_tmp->val;

        r1_piv = &(A->row[i+1]);
        e_tmp = r1_piv->elt;
        aip1 = unord_get_val(A,i+1,i+1);
        det = aii*aip1 - aip1i*aip1i;   /* Must have det < 0 */
        if ( aii == 0.0 && aip1i == 0.0 )
        {
        /* error(E_RANGE,"spBKPfactor"); */
        onebyone = TRUE;
        continue;   /* cannot pivot */
        }

        if ( det == 0.0 )
        {
        if ( aii != 0.0 )
            error(E_RANGE,"spBKPfactor");
        onebyone = TRUE;
        continue;   /* cannot pivot */
        }
        aip1i = aip1i/det;
        aii = aii/det;
        aip1 = aip1/det;
        
        /* initialise scan_... etc for the 2 x 2 pivot */
        s_idx = r_piv->len + r1_piv->len;
        scan_row = iv_resize(scan_row,s_idx);
        scan_idx = iv_resize(scan_idx,s_idx);
        col_list = iv_resize(col_list,s_idx);
        orig_idx = iv_resize(orig_idx,s_idx);
        orig1_idx = iv_resize(orig1_idx,s_idx);

        e = r_piv->elt;
        for ( idx = 0; idx < r_piv->len; idx++, e++ )
        {
        scan_row->ive[idx] = i;
        scan_idx->ive[idx] = idx;
        col_list->ive[idx] = e->col;
        orig_idx->ive[idx] = idx;
        orig1_idx->ive[idx] = -1;
        }
        e = r_piv->elt;
        e1 = r1_piv->elt;
        for ( idx = 0; idx < r1_piv->len; idx++, e1++ )
        {
        scan_row->ive[idx+r_piv->len] = i+1;
        scan_idx->ive[idx+r_piv->len] = idx;
        col_list->ive[idx+r_piv->len] = e1->col;
        orig_idx->ive[idx+r_piv->len] = -1;
        orig1_idx->ive[idx+r_piv->len] = idx;
        }

        e1 = r1_piv->elt;
        order = px_resize(order,scan_row->dim);
        px_ident(order);
        iv_sort(col_list,order);
        tmp_iv = iv_resize(tmp_iv,scan_row->dim);
        for ( idx = 0; idx < order->size; idx++ )
        tmp_iv->ive[idx] = scan_idx->ive[order->pe[idx]];
        iv_copy(tmp_iv,scan_idx);
        for ( idx = 0; idx < order->size; idx++ )
        tmp_iv->ive[idx] = scan_row->ive[order->pe[idx]];
        iv_copy(tmp_iv,scan_row);
        for ( idx = 0; idx < scan_row->dim; idx++ )
        tmp_iv->ive[idx] = orig_idx->ive[order->pe[idx]];
        iv_copy(tmp_iv,orig_idx);
        for ( idx = 0; idx < scan_row->dim; idx++ )
        tmp_iv->ive[idx] = orig1_idx->ive[order->pe[idx]];
        iv_copy(tmp_iv,orig1_idx);

        s_idx = 0;
        old_col = -1;
        for ( idx = 0; idx < scan_row->dim; idx++ )
        {
        if ( col_list->ive[idx] == old_col )
        {
            if ( scan_row->ive[idx] == i )
            {
            scan_row->ive[s_idx-1] = scan_row->ive[idx];
            scan_idx->ive[s_idx-1] = scan_idx->ive[idx];
            col_list->ive[s_idx-1] = col_list->ive[idx];
            orig_idx->ive[s_idx-1] = orig_idx->ive[idx];
            orig1_idx->ive[s_idx-1] = orig1_idx->ive[idx-1];
            }
            else if ( idx > 0 )
            {
            scan_row->ive[s_idx-1] = scan_row->ive[idx-1];
            scan_idx->ive[s_idx-1] = scan_idx->ive[idx-1];
            col_list->ive[s_idx-1] = col_list->ive[idx-1];
            orig_idx->ive[s_idx-1] = orig_idx->ive[idx-1];
            orig1_idx->ive[s_idx-1] = orig1_idx->ive[idx];
            }
        }
        else
        {
            scan_row->ive[s_idx] = scan_row->ive[idx];
            scan_idx->ive[s_idx] = scan_idx->ive[idx];
            col_list->ive[s_idx] = col_list->ive[idx];
            orig_idx->ive[s_idx] = orig_idx->ive[idx];
            orig1_idx->ive[s_idx] = orig1_idx->ive[idx];
            s_idx++;
        }
        old_col = col_list->ive[idx];
        }
        scan_row = iv_resize(scan_row,s_idx);
        scan_idx = iv_resize(scan_idx,s_idx);
        col_list = iv_resize(col_list,s_idx);
        orig_idx = iv_resize(orig_idx,s_idx);
        orig1_idx = iv_resize(orig1_idx,s_idx);

        /* for ( j = i+2; j < n; j++ )  { .... row operation .... } */
        for ( s_idx = 0; s_idx < scan_row->dim; s_idx++ )
        {
        int idx_piv, idx1_piv;
        Real    aip1j, aij, aik, aip1k;
        row_elt *e_ik, *e_ip1k;

        j = col_list->ive[s_idx];
        if ( j < i+2 )
            continue;
        tracecatch(scan_to(A,scan_row,scan_idx,col_list,j),
               "spBKPfactor");

        idx_piv = orig_idx->ive[s_idx];
        aij = ( idx_piv < 0 ) ? 0.0 : r_piv->elt[idx_piv].val;
        /* aij = ( s_idx < r_piv->len ) ? r_piv->elt[s_idx].val :
            0.0; */
        /* aij   = sp_get_val(A,i,j); */
        idx1_piv = orig1_idx->ive[s_idx];
        aip1j = ( idx1_piv < 0 ) ? 0.0 : r1_piv->elt[idx1_piv].val;
        /* aip1j = ( s_idx < r_piv->len ) ? 0.0 :
            r1_piv->elt[s_idx-r_piv->len].val; */
        /* aip1j = sp_get_val(A,i+1,j); */
        s = - aip1i*aip1j + aip1*aij;
        t = - aip1i*aij + aii*aip1j;

        /* for ( k = j; k < n; k++ )  { .... update entry .... } */
        row = &(A->row[j]);
        /* set idx_k and idx1_k indices */
        s_idx2 = s_idx;
        k = col_list->ive[s_idx2];
        idx_k = orig_idx->ive[s_idx2];
        idx1_k = orig1_idx->ive[s_idx2];

        while ( s_idx2 < scan_row->dim )
        {
            k = col_list->ive[s_idx2];
            idx_k = orig_idx->ive[s_idx2];
            idx1_k = orig1_idx->ive[s_idx2];
            e_ik = ( idx_k < 0 ) ? (row_elt *)NULL :
            &(r_piv->elt[idx_k]);
            e_ip1k = ( idx1_k < 0 ) ? (row_elt *)NULL :
            &(r1_piv->elt[idx1_k]);
            aik = ( idx_k >= 0 ) ? e_ik->val : 0.0;
            aip1k = ( idx1_k >= 0 ) ? e_ip1k->val : 0.0;
            if ( scan_row->ive[s_idx2] == j )
            {   /* no fill-in */
            row = &(A->row[j]);
            /* idx = sprow_idx(row,k); */
            idx = scan_idx->ive[s_idx2];
            if ( idx < 0 )
                error(E_INTERN,"spBKPfactor");
            row->elt[idx].val -= s*aik + t*aip1k;
            }
            else
            {   /* fill-in -- insert entry & patch column */
            Real    tmp;
            int old_row, old_idx;
            row_elt *old_e, *new_e;

            tmp = - s*aik - t*aip1k;
            if ( tmp != 0.0 )
            {
                row = &(A->row[j]);
                old_row = scan_row->ive[s_idx2];
                old_idx = scan_idx->ive[s_idx2];

                idx = row->len;
                if ( row->len >= row->maxlen )
                {  tracecatch(sprow_xpd(row,2*row->maxlen+1,
                            TYPE_SPMAT),
                       "spBKPfactor");      }

                row->len = idx + 1;
                /* idx = sprow_idx(row,k); */
                new_e = &(row->elt[idx]);
                new_e->val = tmp;
                new_e->col = k;

                if ( old_row < 0 )
                error(E_INTERN,"spBKPfactor");
                /* old_idx = sprow_idx2(&(A->row[old_row]),
                          k,old_idx); */
                old_e = &(A->row[old_row].elt[old_idx]);
                new_e->nxt_row = old_e->nxt_row;
                new_e->nxt_idx = old_e->nxt_idx;
                old_e->nxt_row = j;
                old_e->nxt_idx = idx;
            }
            }

            /* update idx_k, idx1_k, s_idx2 etc */
            s_idx2++;
        }

        /* store multipliers -- may involve fill-in (!) */
        /* idx = sprow_idx(r_piv,j); */
        idx = orig_idx->ive[s_idx];
        if ( idx >= 0 )
        {
            r_piv->elt[idx].val = s;
        }
        else if ( s != 0.0 )
        {
            int     old_row, old_idx;
            row_elt *new_e, *old_e;

            old_row = -1;   old_idx = j;

            if ( i > 0 )
            {
            tracecatch(chase_col(A,j,&old_row,&old_idx,i-1),
                   "spBKPfactor");
            }
            /* sprow_set_val(r_piv,j,s); */
            idx = r_piv->len;
            if ( r_piv->len >= r_piv->maxlen )
            {   tracecatch(sprow_xpd(r_piv,2*r_piv->maxlen+1,
                         TYPE_SPMAT),
                   "spBKPfactor");          }

            r_piv->len = idx + 1;
            /* idx = sprow_idx(r_piv,j); */
            /* if ( idx < 0 )
            error(E_INTERN,"spBKPfactor"); */
            new_e = &(r_piv->elt[idx]);
            new_e->val = s;
            new_e->col = j;
            if ( old_row < 0 )
            {
            new_e->nxt_row = A->start_row[j];
            new_e->nxt_idx = A->start_idx[j];
            A->start_row[j] = i;
            A->start_idx[j] = idx;
            }
            else
            {
            /* old_idx = sprow_idx2(&(A->row[old_row]),j,old_idx);*/
            if ( old_idx < 0 )
                error(E_INTERN,"spBKPfactor");
            old_e = &(A->row[old_row].elt[old_idx]);
            new_e->nxt_row = old_e->nxt_row;
            new_e->nxt_idx = old_e->nxt_idx;
            old_e->nxt_row = i;
            old_e->nxt_idx = idx;
            }
        }
        /* idx1 = sprow_idx(r1_piv,j); */
        idx1 = orig1_idx->ive[s_idx];
        if ( idx1 >= 0 )
        {
            r1_piv->elt[idx1].val = t;
        }
        else if ( t != 0.0 )
        {
            int     old_row, old_idx;
            row_elt *new_e, *old_e;

            old_row = -1;   old_idx = j;
            tracecatch(chase_col(A,j,&old_row,&old_idx,i),
                   "spBKPfactor");
            /* sprow_set_val(r1_piv,j,t); */
            idx1 = r1_piv->len;
            if ( r1_piv->len >= r1_piv->maxlen )
            {   tracecatch(sprow_xpd(r1_piv,2*r1_piv->maxlen+1,
                         TYPE_SPMAT),
                   "spBKPfactor");          }

            r1_piv->len = idx1 + 1;
            /* idx1 = sprow_idx(r1_piv,j); */
            /* if ( idx < 0 )
            error(E_INTERN,"spBKPfactor"); */
            new_e = &(r1_piv->elt[idx1]);
            new_e->val = t;
            new_e->col = j;
            if ( idx1 < 0 )
            error(E_INTERN,"spBKPfactor");
            new_e = &(r1_piv->elt[idx1]);
            if ( old_row < 0 )
            {
            new_e->nxt_row = A->start_row[j];
            new_e->nxt_idx = A->start_idx[j];
            A->start_row[j] = i+1;
            A->start_idx[j] = idx1;
            }
            else
            {
            old_idx = sprow_idx2(&(A->row[old_row]),j,old_idx);
            if ( old_idx < 0 )
                error(E_INTERN,"spBKPfactor");
            old_e = &(A->row[old_row].elt[old_idx]);
            new_e->nxt_row = old_e->nxt_row;
            new_e->nxt_idx = old_e->nxt_idx;
            old_e->nxt_row = i+1;
            old_e->nxt_idx = idx1;
            }
        }
        }
    }
    }

    /* now sort the rows arrays */
    for ( i = 0; i < A->m; i++ )
    qsort(A->row[i].elt,A->row[i].len,sizeof(row_elt),(int(*)())col_cmp);
    A->flag_col = A->flag_diag = FALSE;

    return A;
}

/* spBKPsolve -- solves A.x = b where A has been factored a la BKPfactor()
   -- returns x, which is created if NULL */
VEC *spBKPsolve(A,pivot,block,b,x)
SPMAT   *A;
PERM    *pivot, *block;
VEC *b, *x;
{
    static VEC  *tmp=VNULL; /* dummy storage needed */
    int     i /* , j */, n, onebyone;
    int     row_num, idx;
    Real    a11, a12, a22, b1, b2, det, sum, *tmp_ve, tmp_diag;
    SPROW   *r;
    row_elt *e;
    
    if ( ! A || ! pivot || ! block || ! b )
    error(E_NULL,"spBKPsolve");
    if ( A->m != A->n )
    error(E_SQUARE,"spBKPsolve");
    n = A->n;
    if ( b->dim != n || pivot->size != n || block->size != n )
    error(E_SIZES,"spBKPsolve");
    x = v_resize(x,n);
    tmp = v_resize(tmp,n);
    MEM_STAT_REG(tmp,TYPE_VEC);
    
    tmp_ve = tmp->ve;

    if ( ! A->flag_col )
    sp_col_access(A);

    px_vec(pivot,b,tmp);
    /* printf("# BKPsolve: effect of pivot: tmp =\n");  v_output(tmp); */

    /* solve for lower triangular part */
    for ( i = 0; i < n; i++ )
    {
    sum = tmp_ve[i];
    if ( block->pe[i] < i )
    {
        /* for ( j = 0; j < i-1; j++ )
          sum -= A_me[j][i]*tmp_ve[j]; */
        row_num = -1;   idx = i;
        e = bump_col(A,i,&row_num,&idx);
        while ( row_num >= 0 && row_num < i-1 )
        {
        sum -= e->val*tmp_ve[row_num];
        e = bump_col(A,i,&row_num,&idx);
        }
    }
    else
    {
        /* for ( j = 0; j < i; j++ )
              sum -= A_me[j][i]*tmp_ve[j]; */
        row_num = -1; idx = i;
        e = bump_col(A,i,&row_num,&idx);
        while ( row_num >= 0 && row_num < i )
        {
        sum -= e->val*tmp_ve[row_num];
        e = bump_col(A,i,&row_num,&idx);
        }
    }
    tmp_ve[i] = sum;
    }

    /* printf("# BKPsolve: solving L part: tmp =\n");   v_output(tmp); */
    /* solve for diagonal part */
    for ( i = 0; i < n; i = onebyone ? i+1 : i+2 )
    {
    onebyone = ( block->pe[i] == i );
    if ( onebyone )
    {
        /* tmp_ve[i] /= A_me[i][i]; */
        tmp_diag = sp_get_val(A,i,i);
        if ( tmp_diag == 0.0 )
        error(E_SING,"spBKPsolve");
        tmp_ve[i] /= tmp_diag;
    }
    else
    {
        a11 = sp_get_val(A,i,i);
        a22 = sp_get_val(A,i+1,i+1);
        a12 = sp_get_val(A,i,i+1);
        b1 = tmp_ve[i];
        b2 = tmp_ve[i+1];
        det = a11*a22-a12*a12;  /* < 0 : see BKPfactor() */
        if ( det == 0.0 )
        error(E_SING,"BKPsolve");
        det = 1/det;
        tmp_ve[i]   = det*(a22*b1-a12*b2);
        tmp_ve[i+1] = det*(a11*b2-a12*b1);
    }
    }

    /* printf("# BKPsolve: solving D part: tmp =\n");   v_output(tmp); */
    /* solve for transpose of lower triangular part */
    for ( i = n-2; i >= 0; i-- )
    {
    sum = tmp_ve[i];
    if ( block->pe[i] > i )
    {
        /* onebyone is false */
        /* for ( j = i+2; j < n; j++ )
          sum -= A_me[i][j]*tmp_ve[j]; */
        if ( i+2 >= n )
        continue;
        r = &(A->row[i]);
        idx = sprow_idx(r,i+2);
        idx = fixindex(idx);
        e = &(r->elt[idx]);
        for ( ; idx < r->len; idx++, e++ )
        sum -= e->val*tmp_ve[e->col];
    }
    else /* onebyone */
    {
        /* for ( j = i+1; j < n; j++ )
          sum -= A_me[i][j]*tmp_ve[j]; */
        r = &(A->row[i]);
        idx = sprow_idx(r,i+1);
        idx = fixindex(idx);
        e = &(r->elt[idx]);
        for ( ; idx < r->len; idx++, e++ )
        sum -= e->val*tmp_ve[e->col];
    }
    tmp_ve[i] = sum;
    }

    /* printf("# BKPsolve: solving L^T part: tmp =\n");v_output(tmp); */
    /* and do final permutation */
    x = pxinv_vec(pivot,tmp,x);

    return x;
}