home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / spbkp < prev    next >
Encoding:
Text File  |  1994-01-13  |  34.8 KB  |  1,387 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*
  28.   Sparse matrix Bunch--Kaufman--Parlett factorisation and solve
  29.   Radical revision started Thu 05th Nov 1992, 09:36:12 AM
  30.   to use Karen George's suggestion of leaving the the row elements unordered
  31.   Radical revision completed Mon 07th Dec 1992, 10:59:57 AM
  32. */
  33.  
  34. static    char    rcsid[] = "$Id: spbkp.c,v 1.5 1994/01/13 05:44:35 des Exp $";
  35.  
  36. #include    <stdio.h>
  37. #include    <math.h>
  38. #include    "matrix.h"
  39. #include    "sparse.h"
  40. #include        "sparse2.h"
  41.  
  42.  
  43. #ifdef MALLOCDECL
  44. #include <malloc.h>
  45. #endif
  46.  
  47. #define alpha    0.6403882032022076 /* = (1+sqrt(17))/8 */
  48.  
  49.  
  50. #define    btos(x)    ((x) ? "TRUE" : "FALSE")
  51.  
  52. /* assume no use of sqr() uses side-effects */
  53. #define    sqr(x)    ((x)*(x))
  54.  
  55. /* unord_get_idx -- returns index (encoded if entry not allocated)
  56.     of the element of row r with column j
  57.     -- uses linear search */
  58. int    unord_get_idx(r,j)
  59. SPROW    *r;
  60. int    j;
  61. {
  62.     int        idx;
  63.     row_elt    *e;
  64.  
  65.     if ( ! r || ! r->elt )
  66.     error(E_NULL,"unord_get_idx");
  67.     for ( idx = 0, e = r->elt; idx < r->len; idx++, e++ )
  68.     if ( e->col == j )
  69.         break;
  70.     if ( idx >= r->len )
  71.     return -(r->len+2);
  72.     else
  73.     return idx;
  74. }
  75.  
  76. /* unord_get_val -- returns value of the (i,j) entry of A
  77.     -- same assumptions as unord_get_idx() */
  78. double    unord_get_val(A,i,j)
  79. SPMAT    *A;
  80. int    i, j;
  81. {
  82.     SPROW    *r;
  83.     int        idx;
  84.  
  85.     if ( ! A )
  86.     error(E_NULL,"unord_get_val");
  87.     if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
  88.     error(E_BOUNDS,"unord_get_val");
  89.  
  90.     r = &(A->row[i]);
  91.     idx = unord_get_idx(r,j);
  92.     if ( idx < 0 )
  93.     return 0.0;
  94.     else
  95.     return r->elt[idx].val;
  96. }
  97.  
  98.         
  99. /* bkp_swap_elt -- swaps the (i,j) with the (k,l) entry of sparse matrix
  100.     -- either or both of the entries may be unallocated */
  101. static SPMAT    *bkp_swap_elt(A,i1,j1,idx1,i2,j2,idx2)
  102. SPMAT    *A;
  103. int    i1, j1, idx1, i2, j2, idx2;
  104. {
  105.     int        tmp_row, tmp_idx;
  106.     SPROW    *r1, *r2;
  107.     row_elt    *e1, *e2;
  108.     Real    tmp;
  109.  
  110.     if ( ! A )
  111.     error(E_NULL,"bkp_swap_elt");
  112.  
  113.     if ( i1 < 0 || j1 < 0 || i2 < 0 || j2 < 0 ||
  114.      i1 >= A->m || j1 >= A->n || i2 >= A->m || j2 >= A->n )
  115.     {
  116.     error(E_BOUNDS,"bkp_swap_elt");
  117.     }
  118.  
  119.     if ( i1 == i2 && j1 == j2 )
  120.     return A;
  121.     if ( idx1 < 0 && idx2 < 0 )    /* neither allocated */
  122.     return A;
  123.  
  124.     r1 = &(A->row[i1]);        r2 = &(A->row[i2]);
  125.     /* if ( idx1 >= r1->len || idx2 >= r2->len )
  126.     error(E_BOUNDS,"bkp_swap_elt"); */
  127.     if ( idx1 < 0 )    /* assume not allocated */
  128.     {
  129.     idx1 = r1->len;
  130.     if ( idx1 >= r1->maxlen )
  131.     {    tracecatch(sprow_xpd(r1,2*r1->maxlen+1,TYPE_SPMAT),
  132.             "bkp_swap_elt");    }
  133.     r1->len = idx1+1;
  134.     r1->elt[idx1].col = j1;
  135.     r1->elt[idx1].val = 0.0;
  136.     /* now patch up column access path */
  137.     tmp_row = -1;    tmp_idx = j1;
  138.     chase_col(A,j1,&tmp_row,&tmp_idx,i1-1);
  139.  
  140.     if ( tmp_row < 0 )
  141.     {
  142.         r1->elt[idx1].nxt_row = A->start_row[j1];
  143.         r1->elt[idx1].nxt_idx = A->start_idx[j1];
  144.         A->start_row[j1] = i1;
  145.         A->start_idx[j1] = idx1;
  146.     }
  147.     else
  148.     {
  149.         row_elt    *tmp_e;
  150.  
  151.         tmp_e = &(A->row[tmp_row].elt[tmp_idx]);
  152.         r1->elt[idx1].nxt_row = tmp_e->nxt_row;
  153.         r1->elt[idx1].nxt_idx = tmp_e->nxt_idx;
  154.         tmp_e->nxt_row = i1;
  155.         tmp_e->nxt_idx = idx1;
  156.     }
  157.     }
  158.     else if ( r1->elt[idx1].col != j1 )
  159.     error(E_INTERN,"bkp_swap_elt");
  160.     if ( idx2 < 0 )
  161.     {
  162.     idx2 = r2->len;
  163.     if ( idx2 >= r2->maxlen )
  164.     {    tracecatch(sprow_xpd(r2,2*r2->maxlen+1,TYPE_SPMAT),
  165.             "bkp_swap_elt");    }
  166.  
  167.     r2->len = idx2+1;
  168.     r2->elt[idx2].col = j2;
  169.     r2->elt[idx2].val = 0.0;
  170.     /* now patch up column access path */
  171.     tmp_row = -1;    tmp_idx = j2;
  172.     chase_col(A,j2,&tmp_row,&tmp_idx,i2-1);
  173.     if ( tmp_row < 0 )
  174.     {
  175.         r2->elt[idx2].nxt_row = A->start_row[j2];
  176.         r2->elt[idx2].nxt_idx = A->start_idx[j2];
  177.         A->start_row[j2] = i2;
  178.         A->start_idx[j2] = idx2;
  179.     }
  180.     else
  181.     {
  182.         row_elt    *tmp_e;
  183.  
  184.         tmp_e = &(A->row[tmp_row].elt[tmp_idx]);
  185.         r2->elt[idx2].nxt_row = tmp_e->nxt_row;
  186.         r2->elt[idx2].nxt_idx = tmp_e->nxt_idx;
  187.         tmp_e->nxt_row = i2;
  188.         tmp_e->nxt_idx = idx2;
  189.     }
  190.     }
  191.     else if ( r2->elt[idx2].col != j2 )
  192.     error(E_INTERN,"bkp_swap_elt");
  193.  
  194.     e1 = &(r1->elt[idx1]);    e2 = &(r2->elt[idx2]);
  195.  
  196.     tmp = e1->val;
  197.     e1->val = e2->val;
  198.     e2->val = tmp;
  199.  
  200.     return A;
  201. }
  202.  
  203. /* bkp_bump_col -- bumps row and idx to next entry in column j */
  204. row_elt    *bkp_bump_col(A, j, row, idx)
  205. SPMAT    *A;
  206. int    j, *row, *idx;
  207. {
  208.     SPROW    *r;
  209.     row_elt    *e;
  210.  
  211.     if ( *row < 0 )
  212.     {
  213.     *row = A->start_row[j];
  214.     *idx = A->start_idx[j];
  215.     }
  216.     else
  217.     {
  218.     r = &(A->row[*row]);
  219.     e = &(r->elt[*idx]);
  220.     if ( e->col != j )
  221.         error(E_INTERN,"bkp_bump_col");
  222.     *row = e->nxt_row;
  223.     *idx = e->nxt_idx;
  224.     }
  225.     if ( *row < 0 )
  226.     return (row_elt *)NULL;
  227.     else
  228.     return &(A->row[*row].elt[*idx]);
  229. }
  230.  
  231. /* bkp_interchange -- swap rows/cols i and j (symmetric pivot)
  232.     -- uses just the upper triangular part */
  233. SPMAT    *bkp_interchange(A, i1, i2)
  234. SPMAT    *A;
  235. int    i1, i2;
  236. {
  237.     int        tmp_row, tmp_idx;
  238.     int        row1, row2, idx1, idx2, tmp_row1, tmp_idx1, tmp_row2, tmp_idx2;
  239.     SPROW    *r1, *r2;
  240.     row_elt    *e1, *e2;
  241.     IVEC    *done_list = IVNULL;
  242.  
  243.     if ( ! A )
  244.     error(E_NULL,"bkp_interchange");
  245.     if ( i1 < 0 || i1 >= A->n || i2 < 0 || i2 >= A->n )
  246.     error(E_BOUNDS,"bkp_interchange");
  247.     if ( A->m != A->n )
  248.     error(E_SQUARE,"bkp_interchange");
  249.  
  250.     if ( i1 == i2 )
  251.     return A;
  252.     if ( i1 > i2 )
  253.     {    tmp_idx = i1;    i1 = i2;    i2 = tmp_idx;    }
  254.  
  255.     done_list = iv_resize(done_list,A->n);
  256.     for ( tmp_idx = 0; tmp_idx < A->n; tmp_idx++ )
  257.     done_list->ive[tmp_idx] = FALSE;
  258.     row1 = -1;        idx1 = i1;
  259.     row2 = -1;        idx2 = i2;
  260.     e1 = bkp_bump_col(A,i1,&row1,&idx1);
  261.     e2 = bkp_bump_col(A,i2,&row2,&idx2);
  262.  
  263.     while ( (row1 >= 0 && row1 < i1) || (row2 >= 0 && row2 < i1) )
  264.     /* Note: "row2 < i1" not "row2 < i2" as we must stop before the
  265.        "knee bend" */
  266.     {
  267.     if ( row1 >= 0 && row1 < i1 && ( row1 < row2 || row2 < 0 ) )
  268.     {
  269.         tmp_row1 = row1;    tmp_idx1 = idx1;
  270.         e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
  271.         if ( ! done_list->ive[row1] )
  272.         {
  273.         if ( row1 == row2 )
  274.             bkp_swap_elt(A,row1,i1,idx1,row1,i2,idx2);
  275.         else
  276.             bkp_swap_elt(A,row1,i1,idx1,row1,i2,-1);
  277.         done_list->ive[row1] = TRUE;
  278.         }
  279.         row1 = tmp_row1;    idx1 = tmp_idx1;
  280.     }
  281.     else if ( row2 >= 0 && row2 < i1 && ( row2 < row1 || row1 < 0 ) )
  282.     {
  283.         tmp_row2 = row2;    tmp_idx2 = idx2;
  284.         e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
  285.         if ( ! done_list->ive[row2] )
  286.         {
  287.         if ( row1 == row2 )
  288.             bkp_swap_elt(A,row2,i1,idx1,row2,i2,idx2);
  289.         else
  290.             bkp_swap_elt(A,row2,i1,-1,row2,i2,idx2);
  291.         done_list->ive[row2] = TRUE;
  292.         }
  293.         row2 = tmp_row2;    idx2 = tmp_idx2;
  294.     }
  295.     else if ( row1 == row2 )
  296.     {
  297.         tmp_row1 = row1;    tmp_idx1 = idx1;
  298.         e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
  299.         tmp_row2 = row2;    tmp_idx2 = idx2;
  300.         e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
  301.         if ( ! done_list->ive[row1] )
  302.         {
  303.         bkp_swap_elt(A,row1,i1,idx1,row2,i2,idx2);
  304.         done_list->ive[row1] = TRUE;
  305.         }
  306.         row1 = tmp_row1;    idx1 = tmp_idx1;
  307.         row2 = tmp_row2;    idx2 = tmp_idx2;
  308.     }
  309.     }
  310.  
  311.     /* ensure we are **past** the first knee */
  312.     while ( row2 >= 0 && row2 <= i1 )
  313.     e2 = bkp_bump_col(A,i2,&row2,&idx2);
  314.  
  315.     /* at/after 1st "knee bend" */
  316.     r1 = &(A->row[i1]);
  317.     idx1 = 0;
  318.     e1 = &(r1->elt[idx1]);
  319.     while ( row2 >= 0 && row2 < i2 )
  320.     {
  321.     /* used for update of e2 at end of loop */
  322.     tmp_row = row2;    tmp_idx = idx2;
  323.     if ( ! done_list->ive[row2] )
  324.     {
  325.         r2 = &(A->row[row2]);
  326.         bkp_bump_col(A,i2,&tmp_row,&tmp_idx);
  327.         done_list->ive[row2] = TRUE;
  328.         tmp_idx1 = unord_get_idx(r1,row2);
  329.         tracecatch(bkp_swap_elt(A,row2,i2,idx2,i1,row2,tmp_idx1),
  330.                "bkp_interchange");
  331.     }
  332.  
  333.     /* update e1 and e2 */
  334.     row2 = tmp_row;    idx2 = tmp_idx;
  335.     e2 = ( row2 >= 0 ) ? &(A->row[row2].elt[idx2]) : (row_elt *)NULL;
  336.     }
  337.  
  338.     idx1 = 0;
  339.     e1 = r1->elt;
  340.     while ( idx1 < r1->len )
  341.     {
  342.     if ( e1->col >= i2 || e1->col <= i1 )
  343.     {
  344.         idx1++;
  345.         e1++;
  346.         continue;
  347.     }
  348.     if ( ! done_list->ive[e1->col] )
  349.     {
  350.         tmp_idx2 = unord_get_idx(&(A->row[e1->col]),i2);
  351.         tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,e1->col,i2,tmp_idx2),
  352.                "bkp_interchange");
  353.         done_list->ive[e1->col] = TRUE;
  354.     }
  355.     idx1++;
  356.     e1++;
  357.     }
  358.  
  359.     /* at/after 2nd "knee bend" */
  360.     idx1 = 0;
  361.     e1 = &(r1->elt[idx1]);
  362.     r2 = &(A->row[i2]);
  363.     idx2 = 0;
  364.     e2 = &(r2->elt[idx2]);
  365.     while ( idx1 < r1->len )
  366.     {
  367.     if ( e1->col <= i2 )
  368.     {
  369.         idx1++;    e1++;
  370.         continue;
  371.     }
  372.     if ( ! done_list->ive[e1->col] )
  373.     {
  374.         tmp_idx2 = unord_get_idx(r2,e1->col);
  375.         tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,i2,e1->col,tmp_idx2),
  376.                "bkp_interchange");
  377.         done_list->ive[e1->col] = TRUE;
  378.     }
  379.     idx1++;    e1++;
  380.     }
  381.  
  382.     idx2 = 0;    e2 = r2->elt;
  383.     while ( idx2 < r2->len )
  384.     {
  385.     if ( e2->col <= i2 )
  386.     {
  387.         idx2++;    e2++;
  388.         continue;
  389.     }
  390.     if ( ! done_list->ive[e2->col] )
  391.     {
  392.         tmp_idx1 = unord_get_idx(r1,e2->col);
  393.         tracecatch(bkp_swap_elt(A,i2,e2->col,idx2,i1,e2->col,tmp_idx1),
  394.                "bkp_interchange");
  395.         done_list->ive[e2->col] = TRUE;
  396.     }
  397.     idx2++;    e2++;
  398.     }
  399.  
  400.     /* now interchange the digonal entries! */
  401.     idx1 = unord_get_idx(&(A->row[i1]),i1);
  402.     idx2 = unord_get_idx(&(A->row[i2]),i2);
  403.     if ( idx1 >= 0 || idx2 >= 0 )
  404.     {
  405.     tracecatch(bkp_swap_elt(A,i1,i1,idx1,i2,i2,idx2),
  406.            "bkp_interchange");
  407.     }
  408.  
  409.     return A;
  410. }
  411.  
  412.  
  413. /* iv_min -- returns minimum of an integer vector
  414.    -- sets index to the position in iv if index != NULL */
  415. int    iv_min(iv,index)
  416. IVEC    *iv;
  417. int    *index;
  418. {
  419.     int        i, i_min, min_val, tmp;
  420.     
  421.     if ( ! iv ) 
  422.     error(E_NULL,"iv_min");
  423.     if ( iv->dim <= 0 )
  424.     error(E_SIZES,"iv_min");
  425.     i_min = 0;
  426.     min_val = iv->ive[0];
  427.     for ( i = 1; i < iv->dim; i++ )
  428.     {
  429.     tmp = iv->ive[i];
  430.     if ( tmp < min_val )
  431.     {
  432.         min_val = tmp;
  433.         i_min = i;
  434.     }
  435.     }
  436.     
  437.     if ( index != (int *)NULL )
  438.     *index = i_min;
  439.     
  440.     return min_val;
  441. }
  442.  
  443. /* max_row_col -- returns max { |A[j][k]| : k >= i, k != j, k != l } given j
  444.     using symmetry and only the upper triangular part of A */
  445. static double max_row_col(A,i,j,l)
  446. SPMAT    *A;
  447. int    i, j, l;
  448. {
  449.     int        row_num, idx;
  450.     SPROW    *r;
  451.     row_elt    *e;
  452.     Real    max_val, tmp;
  453.  
  454.     if ( ! A )
  455.     error(E_NULL,"max_row_col");
  456.     if ( i < 0 || i > A->n || j < 0 || j >= A->n )
  457.     error(E_BOUNDS,"max_row_col");
  458.  
  459.     max_val = 0.0;
  460.  
  461.     idx = unord_get_idx(&(A->row[i]),j);
  462.     if ( idx < 0 )
  463.     {
  464.     row_num = -1;    idx = j;
  465.     e = chase_past(A,j,&row_num,&idx,i);
  466.     }
  467.     else
  468.     {
  469.     row_num = i;
  470.     e = &(A->row[i].elt[idx]);
  471.     }
  472.     while ( row_num >= 0 && row_num < j )
  473.     {
  474.     if ( row_num != l )
  475.     {
  476.         tmp = fabs(e->val);
  477.         if ( tmp > max_val )
  478.         max_val = tmp;
  479.     }
  480.     e = bump_col(A,j,&row_num,&idx);
  481.     }
  482.     r = &(A->row[j]);
  483.     for ( idx = 0, e = r->elt; idx < r->len; idx++, e++ )
  484.     {
  485.     if ( e->col > j && e->col != l )
  486.     {
  487.         tmp = fabs(e->val);
  488.         if ( tmp > max_val )
  489.         max_val = tmp;
  490.     }
  491.     }
  492.  
  493.     return max_val;
  494. }
  495.  
  496. /* nonzeros -- counts non-zeros in A */
  497. static int    nonzeros(A)
  498. SPMAT    *A;
  499. {
  500.     int        cnt, i;
  501.  
  502.     if ( ! A )
  503.     return 0;
  504.     cnt = 0;
  505.     for ( i = 0; i < A->m; i++ )
  506.     cnt += A->row[i].len;
  507.  
  508.     return cnt;
  509. }
  510.  
  511. /* chk_col_access -- for spBKPfactor()
  512.     -- checks that column access path is OK */
  513. int    chk_col_access(A)
  514. SPMAT    *A;
  515. {
  516.     int        cnt_nz, j, row, idx;
  517.     SPROW    *r;
  518.     row_elt    *e;
  519.  
  520.     if ( ! A )
  521.     error(E_NULL,"chk_col_access");
  522.  
  523.     /* count nonzeros as we go down columns */
  524.     cnt_nz = 0;
  525.     for ( j = 0; j < A->n; j++ )
  526.     {
  527.     row = A->start_row[j];
  528.     idx = A->start_idx[j];
  529.     while ( row >= 0 )
  530.     {
  531.         if ( row >= A->m || idx < 0 )
  532.         return FALSE;
  533.         r = &(A->row[row]);
  534.         if ( idx >= r->len )
  535.         return FALSE;
  536.         e = &(r->elt[idx]);
  537.         if ( e->nxt_row >= 0 && e->nxt_row <= row )
  538.         return FALSE;
  539.         row = e->nxt_row;
  540.         idx = e->nxt_idx;
  541.         cnt_nz++;
  542.     }
  543.     }
  544.  
  545.     if ( cnt_nz != nonzeros(A) )
  546.     return FALSE;
  547.     else
  548.     return TRUE;
  549. }
  550.  
  551. /* col_cmp -- compare two columns -- for sorting rows using qsort() */
  552. static int    col_cmp(e1,e2)
  553. row_elt    *e1, *e2;
  554. {
  555.     return e1->col - e2->col;
  556. }
  557.  
  558. /* spBKPfactor -- sparse Bunch-Kaufman-Parlett factorisation of A in-situ
  559.    -- A is factored into the form P'AP = MDM' where 
  560.    P is a permutation matrix, M lower triangular and D is block
  561.    diagonal with blocks of size 1 or 2
  562.    -- P is stored in pivot; blocks[i]==i iff D[i][i] is a block */
  563. SPMAT    *spBKPfactor(A,pivot,blocks,tol)
  564. SPMAT    *A;
  565. PERM    *pivot, *blocks;
  566. double    tol;
  567. {
  568.     int        i, j, k, l, n, onebyone, r;
  569.     int        idx, idx1, idx_piv;
  570.     int        row_num;
  571.     int        best_deg, best_j, best_l, best_cost, mark_cost, deg, deg_j,
  572.             deg_l, ignore_deg;
  573.     int        list_idx, list_idx2, old_list_idx;
  574.     SPROW    *row, *r_piv, *r1_piv;
  575.     row_elt    *e, *e1;
  576.     Real    aii, aip1, aip1i;
  577.     Real    det, max_j, max_l, s, t;
  578.     static IVEC    *scan_row = IVNULL, *scan_idx = IVNULL, *col_list = IVNULL,
  579.         *tmp_iv = IVNULL;
  580.     static IVEC *deg_list = IVNULL;
  581.     static IVEC    *orig_idx = IVNULL, *orig1_idx = IVNULL;
  582.     static PERM    *order = PNULL;
  583.  
  584.     if ( ! A || ! pivot || ! blocks )
  585.     error(E_NULL,"spBKPfactor");
  586.     if ( A->m != A->n )
  587.     error(E_SQUARE,"spBKPfactor");
  588.     if ( A->m != pivot->size || pivot->size != blocks->size )
  589.     error(E_SIZES,"spBKPfactor");
  590.     if ( tol <= 0.0 || tol > 1.0 )
  591.     error(E_RANGE,"spBKPfactor");
  592.     
  593.     n = A->n;
  594.     
  595.     px_ident(pivot);    px_ident(blocks);
  596.     sp_col_access(A);    sp_diag_access(A);
  597.     ignore_deg = FALSE;
  598.  
  599.     deg_list = iv_resize(deg_list,n);
  600.     order = px_resize(order,n);
  601.     MEM_STAT_REG(deg_list,TYPE_IVEC);
  602.     MEM_STAT_REG(order,TYPE_PERM);
  603.  
  604.     scan_row = iv_resize(scan_row,5);
  605.     scan_idx = iv_resize(scan_idx,5);
  606.     col_list = iv_resize(col_list,5);
  607.     orig_idx = iv_resize(orig_idx,5);
  608.     orig_idx = iv_resize(orig1_idx,5);
  609.     orig_idx = iv_resize(tmp_iv,5);
  610.     MEM_STAT_REG(scan_row,TYPE_IVEC);
  611.     MEM_STAT_REG(scan_idx,TYPE_IVEC);
  612.     MEM_STAT_REG(col_list,TYPE_IVEC);
  613.     MEM_STAT_REG(orig_idx,TYPE_IVEC);
  614.     MEM_STAT_REG(orig1_idx,TYPE_IVEC);
  615.     MEM_STAT_REG(tmp_iv,TYPE_IVEC);
  616.  
  617.     for ( i = 0; i < n-1; i = onebyone ? i+1 : i+2 )
  618.     {
  619.     /* now we want to use a Markowitz-style selection rule for
  620.        determining which rows to swap and whether to use
  621.        1x1 or 2x2 pivoting */
  622.  
  623.     /* get list of degrees of nodes */
  624.     deg_list = iv_resize(deg_list,n-i);
  625.     if ( ! ignore_deg )
  626.         for ( j = i; j < n; j++ )
  627.         deg_list->ive[j-i] = 0;
  628.     else
  629.     {
  630.         for ( j = i; j < n; j++ )
  631.         deg_list->ive[j-i] = 1;
  632.         if ( i < n )
  633.         deg_list->ive[0] = 0;
  634.     }
  635.     order = px_resize(order,n-i);
  636.     px_ident(order);
  637.  
  638.     if ( ! ignore_deg )
  639.     {
  640.         for ( j = i; j < n; j++ )
  641.         {
  642.         /* idx = sprow_idx(&(A->row[j]),j+1); */
  643.         /* idx = fixindex(idx); */
  644.         idx = 0;
  645.         row = &(A->row[j]);
  646.         e = &(row->elt[idx]);
  647.         /* deg_list->ive[j-i] += row->len - idx; */
  648.         for ( ; idx < row->len; idx++, e++ )
  649.             if ( e->col >= i )
  650.             deg_list->ive[e->col - i]++;
  651.         }
  652.         /* now deg_list[k] == degree of node k+i */
  653.         
  654.         /* now sort them into increasing order */
  655.         iv_sort(deg_list,order);
  656.         /* now deg_list[idx] == degree of node i+order[idx] */
  657.     }
  658.  
  659.     /* now we can chase through the nodes in order of increasing
  660.        degree, picking out the ones that satisfy our stability
  661.        criterion */
  662.     list_idx = 0;    r = -1;
  663.     best_j = best_l = -1;
  664.     for ( deg = 0; deg <= n; deg++ )
  665.     {
  666.         Real    ajj, all, ajl;
  667.  
  668.         if ( list_idx >= deg_list->dim )
  669.         break;    /* That's all folks! */
  670.         old_list_idx = list_idx;
  671.         while ( list_idx < deg_list->dim &&
  672.             deg_list->ive[list_idx] <= deg )
  673.         {
  674.         j = i+order->pe[list_idx];
  675.         if ( j < i )
  676.             continue;
  677.         /* can we use row/col j for a 1 x 1 pivot? */
  678.         /* find max_j = max_{k>=i} {|A[k][j]|,|A[j][k]|} */
  679.         ajj = fabs(unord_get_val(A,j,j));
  680.         if ( ajj == 0.0 )
  681.         {
  682.             list_idx++;
  683.             continue;    /* can't use this for 1 x 1 pivot */
  684.         }
  685.  
  686.         max_j = max_row_col(A,i,j,-1);
  687.         if ( ajj >= tol/* *alpha */ *max_j )
  688.         {
  689.             onebyone = TRUE;
  690.             best_j = j;
  691.             best_deg = deg_list->ive[list_idx];
  692.             break;
  693.         }
  694.         list_idx++;
  695.         }
  696.         if ( best_j >= 0 )
  697.         break;
  698.         best_cost = 2*n;    /* > any possible Markowitz cost (bound) */
  699.         best_j = best_l = -1;
  700.         list_idx = old_list_idx;
  701.         while ( list_idx < deg_list->dim &&
  702.             deg_list->ive[list_idx] <= deg )
  703.         {
  704.         j = i+order->pe[list_idx];
  705.         ajj = fabs(unord_get_val(A,j,j));
  706.         for ( list_idx2 = 0; list_idx2 < list_idx; list_idx2++ )
  707.         {
  708.             deg_j = deg;
  709.             deg_l = deg_list->ive[list_idx2];
  710.             l = i+order->pe[list_idx2];
  711.             if ( l < i )
  712.             continue;
  713.             /* try using rows/cols (j,l) for a 2 x 2 pivot block */
  714.             all = fabs(unord_get_val(A,l,l));
  715.             ajl = ( j > l ) ? fabs(unord_get_val(A,l,j)) :
  716.                        fabs(unord_get_val(A,j,l));
  717.             det = fabs(ajj*all - ajl*ajl);
  718.             if ( det == 0.0 )
  719.             continue;
  720.             max_j = max_row_col(A,i,j,l);
  721.             max_l = max_row_col(A,i,l,j);
  722.             if ( tol*(all*max_j+ajl*max_l) < det &&
  723.              tol*(ajl*max_j+ajj*max_l) < det )
  724.             {
  725.             /* acceptably stable 2 x 2 pivot */
  726.             /* this is actually an overestimate of the
  727.                Markowitz cost for choosing (j,l) */
  728.             mark_cost = (ajj == 0.0) ?
  729.                 ((all == 0.0) ? deg_j+deg_l : deg_j+2*deg_l) :
  730.                 ((all == 0.0) ? 2*deg_j+deg_l :
  731.                  2*(deg_j+deg_l));
  732.             if ( mark_cost < best_cost )
  733.             {
  734.                 onebyone = FALSE;
  735.                 best_cost = mark_cost;
  736.                 best_j = j;
  737.                 best_l = l;
  738.                 best_deg = deg_j;
  739.             }
  740.             }
  741.         }
  742.         list_idx++;
  743.         }
  744.         if ( best_j >= 0 )
  745.         break;
  746.     }
  747.  
  748.     if ( best_deg > (int)floor(0.8*(n-i)) )
  749.         ignore_deg = TRUE;
  750.  
  751.     /* now do actual interchanges */
  752.     if ( best_j >= 0 && onebyone )
  753.     {
  754.         bkp_interchange(A,i,best_j);
  755.         px_transp(pivot,i,best_j);
  756.     }
  757.     else if ( best_j >= 0 && best_l >= 0 && ! onebyone )
  758.     {
  759.         if ( best_j == i || best_j == i+1 )
  760.         {
  761.         if ( best_l == i || best_l == i+1 )
  762.         {
  763.             /* no pivoting, but must update blocks permutation */
  764.             px_transp(blocks,i,i+1);
  765.             goto dopivot;
  766.         }
  767.         bkp_interchange(A,(best_j == i) ? i+1 : i,best_l);
  768.         px_transp(pivot,(best_j == i) ? i+1 : i,best_l);
  769.         }
  770.         else if ( best_l == i || best_l == i+1 )
  771.         {
  772.         bkp_interchange(A,(best_l == i) ? i+1 : i,best_j);
  773.         px_transp(pivot,(best_l == i) ? i+1 : i,best_j);
  774.         }
  775.         else /* best_j & best_l outside i, i+1 */
  776.         {
  777.         if ( i != best_j )
  778.         {
  779.             bkp_interchange(A,i,best_j);
  780.             px_transp(pivot,i,best_j);
  781.         }
  782.         if ( i+1 != best_l )
  783.         {
  784.             bkp_interchange(A,i+1,best_l);
  785.             px_transp(pivot,i+1,best_l);
  786.         }
  787.         }
  788.     }
  789.     else    /* can't pivot &/or nothing to pivot */
  790.         continue;
  791.  
  792.     /* update blocks permutation */
  793.     if ( ! onebyone )
  794.         px_transp(blocks,i,i+1);
  795.  
  796.     dopivot:
  797.     if ( onebyone )
  798.     {
  799.         int        idx_j, idx_k, s_idx, s_idx2;
  800.         row_elt    *e_ij, *e_ik;
  801.  
  802.         r_piv = &(A->row[i]);
  803.         idx_piv = unord_get_idx(r_piv,i);
  804.         /* if idx_piv < 0 then aii == 0 and no pivoting can be done;
  805.            -- this means that we should continue to the next iteration */
  806.         if ( idx_piv < 0 )
  807.         continue;
  808.         aii = r_piv->elt[idx_piv].val;
  809.         if ( aii == 0.0 )
  810.         continue;
  811.  
  812.         /* for ( j = i+1; j < n; j++ )  { ... pivot step ... } */
  813.         /* initialise scan_... etc for the 1 x 1 pivot */
  814.         scan_row = iv_resize(scan_row,r_piv->len);
  815.         scan_idx = iv_resize(scan_idx,r_piv->len);
  816.         col_list = iv_resize(col_list,r_piv->len);
  817.         orig_idx = iv_resize(orig_idx,r_piv->len);
  818.         row_num = i;    s_idx = idx = 0;
  819.         e = &(r_piv->elt[idx]);
  820.         for ( idx = 0; idx < r_piv->len; idx++, e++ )
  821.         {
  822.         if ( e->col < i )
  823.             continue;
  824.         scan_row->ive[s_idx] = i;
  825.         scan_idx->ive[s_idx] = idx;
  826.         orig_idx->ive[s_idx] = idx;
  827.         col_list->ive[s_idx] = e->col;
  828.         s_idx++;
  829.         }
  830.         scan_row = iv_resize(scan_row,s_idx);
  831.         scan_idx = iv_resize(scan_idx,s_idx);
  832.         col_list = iv_resize(col_list,s_idx);
  833.         orig_idx = iv_resize(orig_idx,s_idx);
  834.  
  835.         order = px_resize(order,scan_row->dim);
  836.         px_ident(order);
  837.         iv_sort(col_list,order);
  838.  
  839.         tmp_iv = iv_resize(tmp_iv,scan_row->dim);
  840.         for ( idx = 0; idx < order->size; idx++ )
  841.         tmp_iv->ive[idx] = scan_idx->ive[order->pe[idx]];
  842.         iv_copy(tmp_iv,scan_idx);
  843.         for ( idx = 0; idx < order->size; idx++ )
  844.         tmp_iv->ive[idx] = scan_row->ive[order->pe[idx]];
  845.         iv_copy(tmp_iv,scan_row);
  846.         for ( idx = 0; idx < scan_row->dim; idx++ )
  847.         tmp_iv->ive[idx] = orig_idx->ive[order->pe[idx]];
  848.         iv_copy(tmp_iv,orig_idx);
  849.  
  850.         /* now do actual pivot */
  851.         /* for ( j = i+1; j < n-1; j++ ) .... */
  852.  
  853.         for ( s_idx = 0; s_idx < scan_row->dim; s_idx++ )
  854.         {
  855.         idx_j = orig_idx->ive[s_idx];
  856.         if ( idx_j < 0 )
  857.             error(E_INTERN,"spBKPfactor");
  858.         e_ij = &(r_piv->elt[idx_j]);
  859.         j = e_ij->col;
  860.         if ( j < i+1 )
  861.             continue;
  862.         scan_to(A,scan_row,scan_idx,col_list,j);
  863.  
  864.         /* compute multiplier */
  865.         t = e_ij->val / aii;
  866.  
  867.         /* for ( k = j; k < n; k++ ) { .... update A[j][k] .... } */
  868.         /* this is the row in which pivoting is done */
  869.         row = &(A->row[j]);
  870.         for ( s_idx2 = s_idx; s_idx2 < scan_row->dim; s_idx2++ )
  871.         {
  872.             idx_k = orig_idx->ive[s_idx2];
  873.             e_ik = &(r_piv->elt[idx_k]);
  874.             k = e_ik->col;
  875.             /* k >= j since col_list has been sorted */
  876.  
  877.             if ( scan_row->ive[s_idx2] == j )
  878.             {    /* no fill-in -- can be done directly */
  879.             idx = scan_idx->ive[s_idx2];
  880.             /* idx = sprow_idx2(row,k,idx); */
  881.             row->elt[idx].val -= t*e_ik->val;
  882.             }
  883.             else
  884.             {    /* fill-in -- insert entry & patch column */
  885.             int    old_row, old_idx;
  886.             row_elt    *old_e, *new_e;
  887.  
  888.             old_row = scan_row->ive[s_idx2];
  889.             old_idx = scan_idx->ive[s_idx2];
  890.             /* old_idx = sprow_idx2(&(A->row[old_row]),k,old_idx); */
  891.  
  892.             if ( old_idx < 0 )
  893.                 error(E_INTERN,"spBKPfactor");
  894.             /* idx = sprow_idx(row,k); */
  895.             /* idx = fixindex(idx); */
  896.             idx = row->len;
  897.  
  898.             /* sprow_set_val(row,k,-t*e_ik->val); */
  899.             if ( row->len >= row->maxlen )
  900.             { tracecatch(sprow_xpd(row,2*row->maxlen+1,TYPE_SPMAT),
  901.                      "spBKPfactor");        }
  902.  
  903.             row->len = idx+1;
  904.  
  905.             new_e = &(row->elt[idx]);
  906.             new_e->val = -t*e_ik->val;
  907.             new_e->col = k;
  908.  
  909.             old_e = &(A->row[old_row].elt[old_idx]);
  910.             new_e->nxt_row = old_e->nxt_row;
  911.             new_e->nxt_idx = old_e->nxt_idx;
  912.             old_e->nxt_row = j;
  913.             old_e->nxt_idx = idx;
  914.             }
  915.         }
  916.         e_ij->val = t;
  917.         }
  918.     }
  919.     else /* onebyone == FALSE */
  920.     {    /* do 2 x 2 pivot */
  921.         int    idx_k, idx1_k, s_idx, s_idx2;
  922.         int    old_col;
  923.         row_elt    *e_tmp;
  924.  
  925.         r_piv = &(A->row[i]);
  926.         idx_piv = unord_get_idx(r_piv,i);
  927.         aii = aip1i = 0.0;
  928.         e_tmp = r_piv->elt;
  929.         for ( idx_piv = 0; idx_piv < r_piv->len; idx_piv++, e_tmp++ )
  930.         if ( e_tmp->col == i )
  931.             aii = e_tmp->val;
  932.             else if ( e_tmp->col == i+1 )
  933.             aip1i = e_tmp->val;
  934.  
  935.         r1_piv = &(A->row[i+1]);
  936.         e_tmp = r1_piv->elt;
  937.         aip1 = unord_get_val(A,i+1,i+1);
  938.         det = aii*aip1 - aip1i*aip1i;    /* Must have det < 0 */
  939.         if ( aii == 0.0 && aip1i == 0.0 )
  940.         {
  941.         /* error(E_RANGE,"spBKPfactor"); */
  942.         onebyone = TRUE;
  943.         continue;    /* cannot pivot */
  944.         }
  945.  
  946.         if ( det == 0.0 )
  947.         {
  948.         if ( aii != 0.0 )
  949.             error(E_RANGE,"spBKPfactor");
  950.         onebyone = TRUE;
  951.         continue;    /* cannot pivot */
  952.         }
  953.         aip1i = aip1i/det;
  954.         aii = aii/det;
  955.         aip1 = aip1/det;
  956.         
  957.         /* initialise scan_... etc for the 2 x 2 pivot */
  958.         s_idx = r_piv->len + r1_piv->len;
  959.         scan_row = iv_resize(scan_row,s_idx);
  960.         scan_idx = iv_resize(scan_idx,s_idx);
  961.         col_list = iv_resize(col_list,s_idx);
  962.         orig_idx = iv_resize(orig_idx,s_idx);
  963.         orig1_idx = iv_resize(orig1_idx,s_idx);
  964.  
  965.         e = r_piv->elt;
  966.         for ( idx = 0; idx < r_piv->len; idx++, e++ )
  967.         {
  968.         scan_row->ive[idx] = i;
  969.         scan_idx->ive[idx] = idx;
  970.         col_list->ive[idx] = e->col;
  971.         orig_idx->ive[idx] = idx;
  972.         orig1_idx->ive[idx] = -1;
  973.         }
  974.         e = r_piv->elt;
  975.         e1 = r1_piv->elt;
  976.         for ( idx = 0; idx < r1_piv->len; idx++, e1++ )
  977.         {
  978.         scan_row->ive[idx+r_piv->len] = i+1;
  979.         scan_idx->ive[idx+r_piv->len] = idx;
  980.         col_list->ive[idx+r_piv->len] = e1->col;
  981.         orig_idx->ive[idx+r_piv->len] = -1;
  982.         orig1_idx->ive[idx+r_piv->len] = idx;
  983.         }
  984.  
  985.         e1 = r1_piv->elt;
  986.         order = px_resize(order,scan_row->dim);
  987.         px_ident(order);
  988.         iv_sort(col_list,order);
  989.         tmp_iv = iv_resize(tmp_iv,scan_row->dim);
  990.         for ( idx = 0; idx < order->size; idx++ )
  991.         tmp_iv->ive[idx] = scan_idx->ive[order->pe[idx]];
  992.         iv_copy(tmp_iv,scan_idx);
  993.         for ( idx = 0; idx < order->size; idx++ )
  994.         tmp_iv->ive[idx] = scan_row->ive[order->pe[idx]];
  995.         iv_copy(tmp_iv,scan_row);
  996.         for ( idx = 0; idx < scan_row->dim; idx++ )
  997.         tmp_iv->ive[idx] = orig_idx->ive[order->pe[idx]];
  998.         iv_copy(tmp_iv,orig_idx);
  999.         for ( idx = 0; idx < scan_row->dim; idx++ )
  1000.         tmp_iv->ive[idx] = orig1_idx->ive[order->pe[idx]];
  1001.         iv_copy(tmp_iv,orig1_idx);
  1002.  
  1003.         s_idx = 0;
  1004.         old_col = -1;
  1005.         for ( idx = 0; idx < scan_row->dim; idx++ )
  1006.         {
  1007.         if ( col_list->ive[idx] == old_col )
  1008.         {
  1009.             if ( scan_row->ive[idx] == i )
  1010.             {
  1011.             scan_row->ive[s_idx-1] = scan_row->ive[idx];
  1012.             scan_idx->ive[s_idx-1] = scan_idx->ive[idx];
  1013.             col_list->ive[s_idx-1] = col_list->ive[idx];
  1014.             orig_idx->ive[s_idx-1] = orig_idx->ive[idx];
  1015.             orig1_idx->ive[s_idx-1] = orig1_idx->ive[idx-1];
  1016.             }
  1017.             else if ( idx > 0 )
  1018.             {
  1019.             scan_row->ive[s_idx-1] = scan_row->ive[idx-1];
  1020.             scan_idx->ive[s_idx-1] = scan_idx->ive[idx-1];
  1021.             col_list->ive[s_idx-1] = col_list->ive[idx-1];
  1022.             orig_idx->ive[s_idx-1] = orig_idx->ive[idx-1];
  1023.             orig1_idx->ive[s_idx-1] = orig1_idx->ive[idx];
  1024.             }
  1025.         }
  1026.         else
  1027.         {
  1028.             scan_row->ive[s_idx] = scan_row->ive[idx];
  1029.             scan_idx->ive[s_idx] = scan_idx->ive[idx];
  1030.             col_list->ive[s_idx] = col_list->ive[idx];
  1031.             orig_idx->ive[s_idx] = orig_idx->ive[idx];
  1032.             orig1_idx->ive[s_idx] = orig1_idx->ive[idx];
  1033.             s_idx++;
  1034.         }
  1035.         old_col = col_list->ive[idx];
  1036.         }
  1037.         scan_row = iv_resize(scan_row,s_idx);
  1038.         scan_idx = iv_resize(scan_idx,s_idx);
  1039.         col_list = iv_resize(col_list,s_idx);
  1040.         orig_idx = iv_resize(orig_idx,s_idx);
  1041.         orig1_idx = iv_resize(orig1_idx,s_idx);
  1042.  
  1043.         /* for ( j = i+2; j < n; j++ )  { .... row operation .... } */
  1044.         for ( s_idx = 0; s_idx < scan_row->dim; s_idx++ )
  1045.         {
  1046.         int    idx_piv, idx1_piv;
  1047.         Real    aip1j, aij, aik, aip1k;
  1048.         row_elt    *e_ik, *e_ip1k;
  1049.  
  1050.         j = col_list->ive[s_idx];
  1051.         if ( j < i+2 )
  1052.             continue;
  1053.         tracecatch(scan_to(A,scan_row,scan_idx,col_list,j),
  1054.                "spBKPfactor");
  1055.  
  1056.         idx_piv = orig_idx->ive[s_idx];
  1057.         aij = ( idx_piv < 0 ) ? 0.0 : r_piv->elt[idx_piv].val;
  1058.         /* aij = ( s_idx < r_piv->len ) ? r_piv->elt[s_idx].val :
  1059.             0.0; */
  1060.         /* aij   = sp_get_val(A,i,j); */
  1061.         idx1_piv = orig1_idx->ive[s_idx];
  1062.         aip1j = ( idx1_piv < 0 ) ? 0.0 : r1_piv->elt[idx1_piv].val;
  1063.         /* aip1j = ( s_idx < r_piv->len ) ? 0.0 :
  1064.             r1_piv->elt[s_idx-r_piv->len].val; */
  1065.         /* aip1j = sp_get_val(A,i+1,j); */
  1066.         s = - aip1i*aip1j + aip1*aij;
  1067.         t = - aip1i*aij + aii*aip1j;
  1068.  
  1069.         /* for ( k = j; k < n; k++ )  { .... update entry .... } */
  1070.         row = &(A->row[j]);
  1071.         /* set idx_k and idx1_k indices */
  1072.         s_idx2 = s_idx;
  1073.         k = col_list->ive[s_idx2];
  1074.         idx_k = orig_idx->ive[s_idx2];
  1075.         idx1_k = orig1_idx->ive[s_idx2];
  1076.  
  1077.         while ( s_idx2 < scan_row->dim )
  1078.         {
  1079.             k = col_list->ive[s_idx2];
  1080.             idx_k = orig_idx->ive[s_idx2];
  1081.             idx1_k = orig1_idx->ive[s_idx2];
  1082.             e_ik = ( idx_k < 0 ) ? (row_elt *)NULL :
  1083.             &(r_piv->elt[idx_k]);
  1084.             e_ip1k = ( idx1_k < 0 ) ? (row_elt *)NULL :
  1085.             &(r1_piv->elt[idx1_k]);
  1086.             aik = ( idx_k >= 0 ) ? e_ik->val : 0.0;
  1087.             aip1k = ( idx1_k >= 0 ) ? e_ip1k->val : 0.0;
  1088.             if ( scan_row->ive[s_idx2] == j )
  1089.             {    /* no fill-in */
  1090.             row = &(A->row[j]);
  1091.             /* idx = sprow_idx(row,k); */
  1092.             idx = scan_idx->ive[s_idx2];
  1093.             if ( idx < 0 )
  1094.                 error(E_INTERN,"spBKPfactor");
  1095.             row->elt[idx].val -= s*aik + t*aip1k;
  1096.             }
  1097.             else
  1098.             {    /* fill-in -- insert entry & patch column */
  1099.             Real    tmp;
  1100.             int    old_row, old_idx;
  1101.             row_elt    *old_e, *new_e;
  1102.  
  1103.             tmp = - s*aik - t*aip1k;
  1104.             if ( tmp != 0.0 )
  1105.             {
  1106.                 row = &(A->row[j]);
  1107.                 old_row = scan_row->ive[s_idx2];
  1108.                 old_idx = scan_idx->ive[s_idx2];
  1109.  
  1110.                 idx = row->len;
  1111.                 if ( row->len >= row->maxlen )
  1112.                 {  tracecatch(sprow_xpd(row,2*row->maxlen+1,
  1113.                             TYPE_SPMAT),
  1114.                        "spBKPfactor");        }
  1115.  
  1116.                 row->len = idx + 1;
  1117.                 /* idx = sprow_idx(row,k); */
  1118.                 new_e = &(row->elt[idx]);
  1119.                 new_e->val = tmp;
  1120.                 new_e->col = k;
  1121.  
  1122.                 if ( old_row < 0 )
  1123.                 error(E_INTERN,"spBKPfactor");
  1124.                 /* old_idx = sprow_idx2(&(A->row[old_row]),
  1125.                           k,old_idx); */
  1126.                 old_e = &(A->row[old_row].elt[old_idx]);
  1127.                 new_e->nxt_row = old_e->nxt_row;
  1128.                 new_e->nxt_idx = old_e->nxt_idx;
  1129.                 old_e->nxt_row = j;
  1130.                 old_e->nxt_idx = idx;
  1131.             }
  1132.             }
  1133.  
  1134.             /* update idx_k, idx1_k, s_idx2 etc */
  1135.             s_idx2++;
  1136.         }
  1137.  
  1138.         /* store multipliers -- may involve fill-in (!) */
  1139.         /* idx = sprow_idx(r_piv,j); */
  1140.         idx = orig_idx->ive[s_idx];
  1141.         if ( idx >= 0 )
  1142.         {
  1143.             r_piv->elt[idx].val = s;
  1144.         }
  1145.         else if ( s != 0.0 )
  1146.         {
  1147.             int        old_row, old_idx;
  1148.             row_elt    *new_e, *old_e;
  1149.  
  1150.             old_row = -1;    old_idx = j;
  1151.  
  1152.             if ( i > 0 )
  1153.             {
  1154.             tracecatch(chase_col(A,j,&old_row,&old_idx,i-1),
  1155.                    "spBKPfactor");
  1156.             }
  1157.             /* sprow_set_val(r_piv,j,s); */
  1158.             idx = r_piv->len;
  1159.             if ( r_piv->len >= r_piv->maxlen )
  1160.             {    tracecatch(sprow_xpd(r_piv,2*r_piv->maxlen+1,
  1161.                          TYPE_SPMAT),
  1162.                    "spBKPfactor");            }
  1163.  
  1164.             r_piv->len = idx + 1;
  1165.             /* idx = sprow_idx(r_piv,j); */
  1166.             /* if ( idx < 0 )
  1167.             error(E_INTERN,"spBKPfactor"); */
  1168.             new_e = &(r_piv->elt[idx]);
  1169.             new_e->val = s;
  1170.             new_e->col = j;
  1171.             if ( old_row < 0 )
  1172.             {
  1173.             new_e->nxt_row = A->start_row[j];
  1174.             new_e->nxt_idx = A->start_idx[j];
  1175.             A->start_row[j] = i;
  1176.             A->start_idx[j] = idx;
  1177.             }
  1178.             else
  1179.             {
  1180.             /* old_idx = sprow_idx2(&(A->row[old_row]),j,old_idx);*/
  1181.             if ( old_idx < 0 )
  1182.                 error(E_INTERN,"spBKPfactor");
  1183.             old_e = &(A->row[old_row].elt[old_idx]);
  1184.             new_e->nxt_row = old_e->nxt_row;
  1185.             new_e->nxt_idx = old_e->nxt_idx;
  1186.             old_e->nxt_row = i;
  1187.             old_e->nxt_idx = idx;
  1188.             }
  1189.         }
  1190.         /* idx1 = sprow_idx(r1_piv,j); */
  1191.         idx1 = orig1_idx->ive[s_idx];
  1192.         if ( idx1 >= 0 )
  1193.         {
  1194.             r1_piv->elt[idx1].val = t;
  1195.         }
  1196.         else if ( t != 0.0 )
  1197.         {
  1198.             int        old_row, old_idx;
  1199.             row_elt    *new_e, *old_e;
  1200.  
  1201.             old_row = -1;    old_idx = j;
  1202.             tracecatch(chase_col(A,j,&old_row,&old_idx,i),
  1203.                    "spBKPfactor");
  1204.             /* sprow_set_val(r1_piv,j,t); */
  1205.             idx1 = r1_piv->len;
  1206.             if ( r1_piv->len >= r1_piv->maxlen )
  1207.             {    tracecatch(sprow_xpd(r1_piv,2*r1_piv->maxlen+1,
  1208.                          TYPE_SPMAT),
  1209.                    "spBKPfactor");            }
  1210.  
  1211.             r1_piv->len = idx1 + 1;
  1212.             /* idx1 = sprow_idx(r1_piv,j); */
  1213.             /* if ( idx < 0 )
  1214.             error(E_INTERN,"spBKPfactor"); */
  1215.             new_e = &(r1_piv->elt[idx1]);
  1216.             new_e->val = t;
  1217.             new_e->col = j;
  1218.             if ( idx1 < 0 )
  1219.             error(E_INTERN,"spBKPfactor");
  1220.             new_e = &(r1_piv->elt[idx1]);
  1221.             if ( old_row < 0 )
  1222.             {
  1223.             new_e->nxt_row = A->start_row[j];
  1224.             new_e->nxt_idx = A->start_idx[j];
  1225.             A->start_row[j] = i+1;
  1226.             A->start_idx[j] = idx1;
  1227.             }
  1228.             else
  1229.             {
  1230.             old_idx = sprow_idx2(&(A->row[old_row]),j,old_idx);
  1231.             if ( old_idx < 0 )
  1232.                 error(E_INTERN,"spBKPfactor");
  1233.             old_e = &(A->row[old_row].elt[old_idx]);
  1234.             new_e->nxt_row = old_e->nxt_row;
  1235.             new_e->nxt_idx = old_e->nxt_idx;
  1236.             old_e->nxt_row = i+1;
  1237.             old_e->nxt_idx = idx1;
  1238.             }
  1239.         }
  1240.         }
  1241.     }
  1242.     }
  1243.  
  1244.     /* now sort the rows arrays */
  1245.     for ( i = 0; i < A->m; i++ )
  1246.     qsort(A->row[i].elt,A->row[i].len,sizeof(row_elt),(int(*)())col_cmp);
  1247.     A->flag_col = A->flag_diag = FALSE;
  1248.  
  1249.     return A;
  1250. }
  1251.  
  1252. /* spBKPsolve -- solves A.x = b where A has been factored a la BKPfactor()
  1253.    -- returns x, which is created if NULL */
  1254. VEC    *spBKPsolve(A,pivot,block,b,x)
  1255. SPMAT    *A;
  1256. PERM    *pivot, *block;
  1257. VEC    *b, *x;
  1258. {
  1259.     static VEC    *tmp=VNULL;    /* dummy storage needed */
  1260.     int        i /* , j */, n, onebyone;
  1261.     int        row_num, idx;
  1262.     Real    a11, a12, a22, b1, b2, det, sum, *tmp_ve, tmp_diag;
  1263.     SPROW    *r;
  1264.     row_elt    *e;
  1265.     
  1266.     if ( ! A || ! pivot || ! block || ! b )
  1267.     error(E_NULL,"spBKPsolve");
  1268.     if ( A->m != A->n )
  1269.     error(E_SQUARE,"spBKPsolve");
  1270.     n = A->n;
  1271.     if ( b->dim != n || pivot->size != n || block->size != n )
  1272.     error(E_SIZES,"spBKPsolve");
  1273.     x = v_resize(x,n);
  1274.     tmp = v_resize(tmp,n);
  1275.     MEM_STAT_REG(tmp,TYPE_VEC);
  1276.     
  1277.     tmp_ve = tmp->ve;
  1278.  
  1279.     if ( ! A->flag_col )
  1280.     sp_col_access(A);
  1281.  
  1282.     px_vec(pivot,b,tmp);
  1283.     /* printf("# BKPsolve: effect of pivot: tmp =\n");    v_output(tmp); */
  1284.  
  1285.     /* solve for lower triangular part */
  1286.     for ( i = 0; i < n; i++ )
  1287.     {
  1288.     sum = tmp_ve[i];
  1289.     if ( block->pe[i] < i )
  1290.     {
  1291.         /* for ( j = 0; j < i-1; j++ )
  1292.           sum -= A_me[j][i]*tmp_ve[j]; */
  1293.         row_num = -1;    idx = i;
  1294.         e = bump_col(A,i,&row_num,&idx);
  1295.         while ( row_num >= 0 && row_num < i-1 )
  1296.         {
  1297.         sum -= e->val*tmp_ve[row_num];
  1298.         e = bump_col(A,i,&row_num,&idx);
  1299.         }
  1300.     }
  1301.     else
  1302.     {
  1303.         /* for ( j = 0; j < i; j++ )
  1304.               sum -= A_me[j][i]*tmp_ve[j]; */
  1305.         row_num = -1; idx = i;
  1306.         e = bump_col(A,i,&row_num,&idx);
  1307.         while ( row_num >= 0 && row_num < i )
  1308.         {
  1309.         sum -= e->val*tmp_ve[row_num];
  1310.         e = bump_col(A,i,&row_num,&idx);
  1311.         }
  1312.     }
  1313.     tmp_ve[i] = sum;
  1314.     }
  1315.  
  1316.     /* printf("# BKPsolve: solving L part: tmp =\n");    v_output(tmp); */
  1317.     /* solve for diagonal part */
  1318.     for ( i = 0; i < n; i = onebyone ? i+1 : i+2 )
  1319.     {
  1320.     onebyone = ( block->pe[i] == i );
  1321.     if ( onebyone )
  1322.     {
  1323.         /* tmp_ve[i] /= A_me[i][i]; */
  1324.         tmp_diag = sp_get_val(A,i,i);
  1325.         if ( tmp_diag == 0.0 )
  1326.         error(E_SING,"spBKPsolve");
  1327.         tmp_ve[i] /= tmp_diag;
  1328.     }
  1329.     else
  1330.     {
  1331.         a11 = sp_get_val(A,i,i);
  1332.         a22 = sp_get_val(A,i+1,i+1);
  1333.         a12 = sp_get_val(A,i,i+1);
  1334.         b1 = tmp_ve[i];
  1335.         b2 = tmp_ve[i+1];
  1336.         det = a11*a22-a12*a12;    /* < 0 : see BKPfactor() */
  1337.         if ( det == 0.0 )
  1338.         error(E_SING,"BKPsolve");
  1339.         det = 1/det;
  1340.         tmp_ve[i]   = det*(a22*b1-a12*b2);
  1341.         tmp_ve[i+1] = det*(a11*b2-a12*b1);
  1342.     }
  1343.     }
  1344.  
  1345.     /* printf("# BKPsolve: solving D part: tmp =\n");    v_output(tmp); */
  1346.     /* solve for transpose of lower triangular part */
  1347.     for ( i = n-2; i >= 0; i-- )
  1348.     {
  1349.     sum = tmp_ve[i];
  1350.     if ( block->pe[i] > i )
  1351.     {
  1352.         /* onebyone is false */
  1353.         /* for ( j = i+2; j < n; j++ )
  1354.           sum -= A_me[i][j]*tmp_ve[j]; */
  1355.         if ( i+2 >= n )
  1356.         continue;
  1357.         r = &(A->row[i]);
  1358.         idx = sprow_idx(r,i+2);
  1359.         idx = fixindex(idx);
  1360.         e = &(r->elt[idx]);
  1361.         for ( ; idx < r->len; idx++, e++ )
  1362.         sum -= e->val*tmp_ve[e->col];
  1363.     }
  1364.     else /* onebyone */
  1365.     {
  1366.         /* for ( j = i+1; j < n; j++ )
  1367.           sum -= A_me[i][j]*tmp_ve[j]; */
  1368.         r = &(A->row[i]);
  1369.         idx = sprow_idx(r,i+1);
  1370.         idx = fixindex(idx);
  1371.         e = &(r->elt[idx]);
  1372.         for ( ; idx < r->len; idx++, e++ )
  1373.         sum -= e->val*tmp_ve[e->col];
  1374.     }
  1375.     tmp_ve[i] = sum;
  1376.     }
  1377.  
  1378.     /* printf("# BKPsolve: solving L^T part: tmp =\n");v_output(tmp); */
  1379.     /* and do final permutation */
  1380.     x = pxinv_vec(pivot,tmp,x);
  1381.  
  1382.     return x;
  1383. }
  1384.  
  1385.  
  1386.  
  1387.