home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / spswap.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  8KB  |  305 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 swap and permutation routines
  29.     Modified Mon 09th Nov 1992, 08:50:54 PM
  30.     to use Karen George's suggestion to use unordered rows
  31. */
  32.  
  33. static    char    rcsid[] = "$Id: spswap.c,v 1.3 1994/01/13 05:44:43 des Exp $";
  34.  
  35. #include    <stdio.h>
  36. #include    <math.h>
  37. #include    "matrix.h"
  38. #include    "sparse.h"
  39. #include        "sparse2.h"
  40.  
  41.  
  42. #define    btos(x)    ((x) ? "TRUE" : "FALSE")
  43.  
  44. /* scan_to -- updates scan (int) vectors to point to the last row in each
  45.     column with row # <= max_row, if any */
  46. void    scan_to(A, scan_row, scan_idx, col_list, max_row)
  47. SPMAT    *A;
  48. IVEC    *scan_row, *scan_idx, *col_list;
  49. int    max_row;
  50. {
  51.     int        col, idx, j_idx, row_num;
  52.     SPROW    *r;
  53.     row_elt    *e;
  54.  
  55.     if ( ! A || ! scan_row || ! scan_idx || ! col_list )
  56.     error(E_NULL,"scan_to");
  57.     if ( scan_row->dim != scan_idx->dim || scan_idx->dim != col_list->dim )
  58.     error(E_SIZES,"scan_to");
  59.  
  60.     if ( max_row < 0 )
  61.     return;
  62.  
  63.     if ( ! A->flag_col )
  64.     sp_col_access(A);
  65.  
  66.     for ( j_idx = 0; j_idx < scan_row->dim; j_idx++ )
  67.     {
  68.     row_num = scan_row->ive[j_idx];
  69.     idx = scan_idx->ive[j_idx];
  70.     col = col_list->ive[j_idx];
  71.  
  72.     if ( col < 0 || col >= A->n )
  73.         error(E_BOUNDS,"scan_to");
  74.     if ( row_num < 0 )
  75.     {
  76.         idx = col;
  77.         continue;
  78.     }
  79.     r = &(A->row[row_num]);
  80.     if ( idx < 0 )
  81.         error(E_INTERN,"scan_to");
  82.     e = &(r->elt[idx]);
  83.     if ( e->col != col )
  84.         error(E_INTERN,"scan_to");
  85.     if ( idx < 0 )
  86.     {
  87.         printf("scan_to: row_num = %d, idx = %d, col = %d\n",
  88.            row_num, idx, col);
  89.         error(E_INTERN,"scan_to");
  90.     }
  91.     /* if ( e->nxt_row <= max_row )
  92.         chase_col(A, col, &row_num, &idx, max_row); */
  93.     while ( e->nxt_row >= 0 && e->nxt_row <= max_row )
  94.     {
  95.         row_num = e->nxt_row;
  96.         idx = e->nxt_idx;
  97.         e = &(A->row[row_num].elt[idx]);
  98.     }
  99.         
  100.     /* printf("scan_to: computed j_idx = %d, row_num = %d, idx = %d\n",
  101.            j_idx, row_num, idx); */
  102.     scan_row->ive[j_idx] = row_num;
  103.     scan_idx->ive[j_idx] = idx;
  104.     }
  105. }
  106.  
  107. /* patch_col -- patches column access paths for fill-in */
  108. void patch_col(A, col, old_row, old_idx, row_num, idx)
  109. SPMAT    *A;
  110. int    col, old_row, old_idx, row_num, idx;
  111. {
  112.     SPROW    *r;
  113.     row_elt    *e;
  114.     
  115.     if ( old_row >= 0 )
  116.     {
  117.     r = &(A->row[old_row]);
  118.     old_idx = sprow_idx2(r,col,old_idx);
  119.     e = &(r->elt[old_idx]);
  120.     e->nxt_row = row_num;
  121.     e->nxt_idx = idx;
  122.     }
  123.     else
  124.     {
  125.     A->start_row[col] = row_num;
  126.     A->start_idx[col] = idx;
  127.     }
  128. }
  129.  
  130. /* chase_col -- chases column access path in column col, starting with
  131.    row_num and idx, to find last row # in this column <= max_row
  132.    -- row_num is returned; idx is also set by this routine
  133.    -- assumes that the column access paths (possibly without the
  134.    nxt_idx fields) are set up */
  135. row_elt *chase_col(A, col, row_num, idx, max_row)
  136. SPMAT    *A;
  137. int    col, *row_num, *idx, max_row;
  138. {
  139.     int        old_idx, old_row, tmp_idx, tmp_row;
  140.     SPROW    *r;
  141.     row_elt    *e;
  142.     
  143.     if ( col < 0 || col >= A->n )
  144.     error(E_BOUNDS,"chase_col");
  145.     tmp_row = *row_num;
  146.     if ( tmp_row < 0 )
  147.     {
  148.     if ( A->start_row[col] > max_row )
  149.     {
  150.         tmp_row = -1;
  151.         tmp_idx = col;
  152.         return (row_elt *)NULL;
  153.     }
  154.     else
  155.     {
  156.         tmp_row = A->start_row[col];
  157.         tmp_idx = A->start_idx[col];
  158.     }
  159.     }
  160.     else
  161.     tmp_idx = *idx;
  162.     
  163.     old_row = tmp_row;
  164.     old_idx = tmp_idx;
  165.     while ( tmp_row >= 0 && tmp_row < max_row )
  166.     {
  167.     r = &(A->row[tmp_row]);
  168.     /* tmp_idx = sprow_idx2(r,col,tmp_idx); */
  169.     if ( tmp_idx < 0 || tmp_idx >= r->len ||
  170.          r->elt[tmp_idx].col != col )
  171.     {
  172. #ifdef DEBUG
  173.         printf("chase_col:error: col = %d, row # = %d, idx = %d\n",
  174.            col, tmp_row, tmp_idx);
  175.         printf("chase_col:error: old_row = %d, old_idx = %d\n",
  176.            old_row, old_idx);
  177.         printf("chase_col:error: A =\n");
  178.         sp_dump(stdout,A);
  179. #endif
  180.         error(E_INTERN,"chase_col");
  181.     }
  182.     e = &(r->elt[tmp_idx]);
  183.     old_row = tmp_row;
  184.     old_idx = tmp_idx;
  185.     tmp_row = e->nxt_row;
  186.     tmp_idx = e->nxt_idx;
  187.     }
  188.     if ( old_row > max_row )
  189.     {
  190.     old_row = -1;
  191.     old_idx = col;
  192.     e = (row_elt *)NULL;
  193.     }
  194.     else if ( tmp_row <= max_row && tmp_row >= 0 )
  195.     {
  196.     old_row = tmp_row;
  197.     old_idx = tmp_idx;
  198.     }
  199.  
  200.     *row_num = old_row;
  201.     if ( old_row >= 0 )
  202.     *idx = old_idx;
  203.     else
  204.     *idx = col;
  205.  
  206.     return e;
  207. }
  208.  
  209. /* chase_past -- as for chase_col except that we want the first
  210.     row whose row # >= min_row; -1 indicates no such row */
  211. row_elt *chase_past(A, col, row_num, idx, min_row)
  212. SPMAT    *A;
  213. int    col, *row_num, *idx, min_row;
  214. {
  215.     SPROW    *r;
  216.     row_elt    *e;
  217.     int        tmp_idx, tmp_row;
  218.  
  219.     tmp_row = *row_num;
  220.     tmp_idx = *idx;
  221.     chase_col(A,col,&tmp_row,&tmp_idx,min_row);
  222.     if ( tmp_row < 0 )    /* use A->start_row[..] etc. */
  223.     {
  224.     if ( A->start_row[col] < 0 )
  225.         tmp_row = -1;
  226.     else
  227.     {
  228.         tmp_row = A->start_row[col];
  229.         tmp_idx = A->start_idx[col];
  230.     }
  231.     }
  232.     else if ( tmp_row < min_row )
  233.     {
  234.     r = &(A->row[tmp_row]);
  235.     if ( tmp_idx < 0 || tmp_idx >= r->len ||
  236.          r->elt[tmp_idx].col != col )
  237.         error(E_INTERN,"chase_past");
  238.     tmp_row = r->elt[tmp_idx].nxt_row;
  239.     tmp_idx = r->elt[tmp_idx].nxt_idx;
  240.     }
  241.  
  242.     *row_num = tmp_row;
  243.     *idx = tmp_idx;
  244.     if ( tmp_row < 0 )
  245.     e = (row_elt *)NULL;
  246.     else
  247.     {
  248.     if ( tmp_idx < 0 || tmp_idx >= A->row[tmp_row].len ||
  249.          A->row[tmp_row].elt[tmp_idx].col != col )
  250.         error(E_INTERN,"bump_col");
  251.     e = &(A->row[tmp_row].elt[tmp_idx]);
  252.     }
  253.  
  254.     return e;
  255. }
  256.  
  257. /* bump_col -- move along to next nonzero entry in column col after row_num
  258.     -- update row_num and idx */
  259. row_elt *bump_col(A, col, row_num, idx)
  260. SPMAT    *A;
  261. int    col, *row_num, *idx;
  262. {
  263.     SPROW    *r;
  264.     row_elt    *e;
  265.     int        tmp_row, tmp_idx;
  266.  
  267.     tmp_row = *row_num;
  268.     tmp_idx = *idx;
  269.     /* printf("bump_col: col = %d, row# = %d, idx = %d\n",
  270.        col, *row_num, *idx); */
  271.     if ( tmp_row < 0 )
  272.     {
  273.     tmp_row = A->start_row[col];
  274.     tmp_idx = A->start_idx[col];
  275.     }
  276.     else
  277.     {
  278.     r = &(A->row[tmp_row]);
  279.     if ( tmp_idx < 0 || tmp_idx >= r->len ||
  280.          r->elt[tmp_idx].col != col )
  281.         error(E_INTERN,"bump_col");
  282.     e = &(r->elt[tmp_idx]);
  283.     tmp_row = e->nxt_row;
  284.     tmp_idx = e->nxt_idx;
  285.     }
  286.     if ( tmp_row < 0 )
  287.     {
  288.     e = (row_elt *)NULL;
  289.     tmp_idx = col;
  290.     }
  291.     else
  292.     {
  293.     if ( tmp_idx < 0 || tmp_idx >= A->row[tmp_row].len ||
  294.          A->row[tmp_row].elt[tmp_idx].col != col )
  295.         error(E_INTERN,"bump_col");
  296.     e = &(A->row[tmp_row].elt[tmp_idx]);
  297.     }
  298.     *row_num = tmp_row;
  299.     *idx = tmp_idx;
  300.  
  301.     return e;
  302. }
  303.  
  304.  
  305.