home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / splufctr < prev    next >
Encoding:
Text File  |  1994-05-09  |  10.6 KB  |  408 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Stewart & 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 LU factorisation
  29.     See also: sparse.[ch] etc for details about sparse matrices
  30. */
  31.  
  32. #include    <stdio.h>
  33. #include    <math.h>
  34. #include        "sparse2.h"
  35.  
  36.  
  37.  
  38. /* Macro for speedup */
  39. /* #define    sprow_idx2(r,c,hint)    \
  40.    ( ( (hint) >= 0 && (r)->elt[hint].col == (c)) ? hint : sprow_idx((r),(c)) ) */
  41.  
  42.  
  43. /* spLUfactor -- sparse LU factorisation with pivoting
  44.     -- uses partial pivoting and Markowitz criterion
  45.             |a[p][k]| >= alpha * max_i |a[i][k]|
  46.     -- creates fill-in as needed
  47.     -- in situ factorisation */
  48. SPMAT    *spLUfactor(A,px,alpha)
  49. SPMAT    *A;
  50. PERM    *px;
  51. double    alpha;
  52. {
  53.     int    i, best_i, k, idx, len, best_len, m, n;
  54.     SPROW    *r, *r_piv, tmp_row;
  55.     static    SPROW    *merge = (SPROW *)NULL;
  56.     Real    max_val, tmp;
  57.     static VEC    *col_vals=VNULL;
  58.  
  59.     if ( ! A || ! px )
  60.         error(E_NULL,"spLUfctr");
  61.     if ( alpha <= 0.0 || alpha > 1.0 )
  62.         error(E_RANGE,"alpha in spLUfctr");
  63.     if ( px->size <= A->m )
  64.         px = px_resize(px,A->m);
  65.     px_ident(px);
  66.     col_vals = v_resize(col_vals,A->m);
  67.     MEM_STAT_REG(col_vals,TYPE_VEC);
  68.  
  69.     m = A->m;    n = A->n;
  70.     if ( ! A->flag_col )
  71.         sp_col_access(A);
  72.     if ( ! A->flag_diag )
  73.         sp_diag_access(A);
  74.     A->flag_col = A->flag_diag = FALSE;
  75.     if ( ! merge ) {
  76.        merge = sprow_get(20);
  77.        MEM_STAT_REG(merge,TYPE_SPROW);
  78.     }
  79.  
  80.     for ( k = 0; k < n; k++ )
  81.     {
  82.         /* find pivot row/elr ( idx = 0; idx < scan_row->dim; idx++ )
  83.         tmp_iv->ive[idx] = orig1_idx->ive[order->pe[idx]];
  84.         iv_copy(tmp_iv,orig1_idx);
  85.  
  86.         s_idx = 0;
  87.         old_col = -1;
  88.         for ( idx = 0; idx < scan_row->dim; idx++ )
  89.         {
  90.         if ( col_list->ive[idx] == old_col )
  91.         {
  92.             if ( scan_row->ive[idx] == i )
  93.             {
  94.             scan_row->ive[s_idx-1] = scan_row->ive[idx];
  95.             scan_idx->ive[s_idx-1] = scan_idx->ive[idx];
  96.             col_list->ive[s_idx-1] = col_list->ive[idx];
  97.             orig_idx->ive[s_idx-1] = orig_idx->ive[idx];
  98.             orig1_idx->ive[s_idx-1] = orig1_idx->ive[idx-1];
  99.             }
  100.             else if ( idx > 0 )
  101.             {
  102.             scan_row->ive[s_idx-1] = scan_row->ive[idx-1];
  103.             scan_idx->ive[s_idx-1] = scan_idx->ive[idx-1];
  104.             col_list->ive[s_idx-1] = col_list->ive[idx-1];
  105.             orig_idx->ive[s_idx-1] = orig_idx->ive[idx-1];
  106.             orig1_idx->ive[s_idx-1] = orig1_idx->ive[idx];
  107.             }
  108.         }
  109.         else
  110.         {
  111.             scan_row->ive[s_idx] = scan_row->ive[idx];
  112.             scan_idx->ive[s_idx] = scan_idx->ive[idx];
  113.             col_list->ive[s_idx] = col_list->ive[idx];
  114.             orig_idx->ive[s_idx] = orig_idx->ive[idx];
  115.             orig1_idx->ive[s_idx] = orig1_idx->ive[idx];
  116.             s_idx++;
  117.         }
  118.         old_col = col_list->ive[idx];
  119.         }
  120.         scan_row = iv_resize(scan_row,s_idx);
  121.         scan_idx = iv_resize(scan_idx,s_idx);
  122.         col_list = iv_resize(col_list,s_idx);
  123.         orig_idx = iv_resize(orig_idx,s_idx);
  124.         orig1_idx = iv_resize(orig1_idx,s_idx);
  125.  
  126.         /* for ( j = i+2; j < n; j++ )  { .... row operation .... } */
  127.         for ( s_idx = 0; s_idx < scan_row->dim; s_idx++ )
  128.         {
  129.         int    idx_piv, id if ( ! r ) 
  130.      return sprow_get(n);
  131.    
  132.    if (n == r->len)
  133.      return r;
  134.  
  135.    if ( ! r->elt )
  136.    {
  137.       r->elt = NEW_A((unsigned)n,row_elt);
  138.       if ( ! r->elt )
  139.     error(E_MEM,"sprow_resize");
  140.       else if (mem_info_is_on()) {
  141.      mem_bytes(type,0,n*sizeof(row_elt));
  142.       }
  143.       r->maxlen = r->len = n;
  144.       return r;
  145.    }
  146.  
  147.    if ( n <= r->maxlen )
  148.      r->len = n;
  149.    else
  150.    {
  151.       if (mem_info_is_on()) {
  152.      mem_bytes(type,r->maxlen*sizeof(row_elt),
  153.            n*sizeof(row_elt)); 
  154.       }
  155.       r->elt = RENEW(r->elt,n,row_elt);
  156.       if ( ! r->elt )
  157.     error(E_MEM,"sprow_resize");
  158.       r->maxlen = r->len = n;
  159.    }
  160.    
  161.    return r;
  162. }
  163.  
  164.  
  165. /* release a row of a matrix */
  166. int sprow_free(r)
  167. SPROW    *r;
  168. {
  169.    if ( ! r )
  170.      return -1;
  171.  
  172.    if (mem_info_is_on()) {
  173.       mem_bytes(TYPE_SPROW,sizeof(SPROW),0);
  174.       mem_numvar(TYPE_SPROW,-1);
  175.    }
  176.    
  177.    if ( r->elt )
  178.    {
  179.       if (mem_info_is_on()) {
  180.      mem_bytes(TYPE_SPROW,r->maxlen*sizeof(row_elt),0);
  181.       }
  182.       free((char *)r->elt);
  183.    }
  184.    free((char *)r);
  185.    return 0;
  186. }
  187.  
  188.  
  189. /* sprow_merge -- merges r1 and r2 into r_out
  190.    -- cannot be done in-situ
  191.    -- type must be SPMAT or SPROW depending on
  192.       whether r_out is a row of a SPMAT structure
  193.       or a SPROW variable
  194.    -- returns r_out */
  195. SPROW    *sprow_merge(r1,r2,r_out,type)
  196. SPROW    *r1, *r2, *r_out;
  197. int type;
  198. {
  199.    int    idx1, idx2, idx_out, len1, len2, len_out;
  200.    row_elt    *elt1, *elt2, *elt_out;
  201.    
  202.    if ( ! r1 || ! r2 )
  203.      error(E_NULL,"sprow_merge");
  204.    if ( ! r_out )
  205.      r_out = sprow_get(MINROWLEN);
  206.    if ( r1 == r_out || r2 == r_out )
  207.      error(E_INSITU,"sprow_merge");
  208.    
  209.    /* Initialise */
  210.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  211.    idx1 = idx2 = idx_out = 0;
  212.    elt1 = r1->elt;    elt2 = r2->elt;    elt_out = r_out->elt;
  213.    
  214.    while ( idx1 < len1 || idx2 < len2 )
  215.    {
  216.       if ( idx_out >= len_out )
  217.       {   /* r_out is too small */
  218.      r_out->len = idx_out;
  219.      r_out = sprow_xpd(r_out,0,type);
  220.      len_out = r_out->len;
  221.      elt_out = &(r_out->elt[idx_out]);
  222.       }
  223.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  224.       {
  225.      elt_out->col = elt1->col;
  226.      elt_out->val = elt1->val;
  227.      if ( elt1->col == elt2->col && idx2 < len2 )
  228.      {    elt2++;        idx2++;    }
  229.      elt1++;    idx1++;
  230.       }
  231.       else
  232.       {
  233.      elt_out->col = elt2->col;
  234.      elt_out->val = elt2->val;
  235.      elt2++;    idx2++;
  236.       }
  237.       elt_out++;    idx_out++;
  238.    }
  239.    r_out->len = idx_out;
  240.    
  241.    return r_out;
  242. }
  243.  
  244. /* sprow_copy -- copies r1 and r2 into r_out
  245.    -- cannot be done in-situ
  246.    -- type must be SPMAT or SPROW depending on
  247.       whether r_out is a row of a SPMAT structure
  248.       or a SPROW variable
  249.    -- returns r_out */
  250. SPROW    *sprow_copy(r1,r2,r_out,type)
  251. SPROW    *r1, *r2, *r_out;
  252. int type;
  253. {
  254.    int    idx1, idx2, idx_out, len1, len2, len_out;
  255.    row_elt    *elt1, *elt2, *elt_out;
  256.    
  257.    if ( ! r1 || ! r2 )
  258.      error(E_NULL,"sprow_copy");
  259.    if ( ! r_out )
  260.      r_out = sprow_get(MINROWLEN);
  261.    if ( r1 == r_out || r2 == r_out )
  262.      error(E_INSITU,"sprow_copy");
  263.    
  264.    /* Initialise */
  265.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  266.    idx1 = idx2 = idx_out = 0;
  267.    elt1 = r1->elt;    elt2 = r2->elt;    elt_out = r_out->elt;
  268.    
  269.    while ( idx1 < len1 || idx2 < len2 )
  270.    {
  271.       while ( idx_out >= len_out )
  272.       {   /* r_out is too small */
  273.      r_out->len = idx_out;
  274.      r_out = sprow_xpd(r_out,0,type);
  275.      len_out = r_out->maxlen;
  276.      elt_out = &(r_out->elt[idx_out]);
  277.       }
  278.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  279.       {
  280.      elt_out->col = elt1->col;
  281.      elt_out->val = elt1->val;
  282.      if ( elt1->col == elt2->col && idx2 < len2 )
  283.      {    elt2++;        idx2++;    }
  284.      elt1++;    idx1++;
  285.       }
  286.       else
  287.       {
  288.      elt_out->col = elt2->col;
  289.      elt_out->val = 0.0;
  290.      elt2++;    idx2++;
  291.       }
  292.       elt_out++;    idx_out++;
  293.    }
  294.    r_out->len = idx_out;
  295.    
  296.    return r_out;
  297. }
  298.  
  299. /* sprow_mltadd -- sets r_out <- r1 + alpha.r2
  300.    -- cannot be in situ
  301.    -- only for columns j0, j0+1, ...
  302.    -- type must be SPMAT or SPROW depending on
  303.       whether r_out is a row of a SPMAT structure
  304.       or a SPROW variable
  305.    -- returns r_out */
  306. SPROW    *sprow_mltadd(r1,r2,alpha,j0,r_out,type)
  307. SPROW    *r1, *r2, *r_out;
  308. double    alpha;
  309. int    j0, type;
  310. {
  311.    int    idx1, idx2, idx_out, len1, len2, len_out;
  312.    row_elt    *elt1, *elt2, *elt_out;
  313.    
  314.    if ( ! r1 || ! r2 )
  315.      error(E_NULL,"sprow_mltadd");
  316.    if ( r1 == r_out || r2 == r_out )
  317.      error(E_INSITU,"sprow_mltadd");
  318.    if ( j0 < 0 )
  319.      error(E_BOUNDS,"sprow_mltadd");
  320.    if ( ! r_out )
  321.      r_out = sprow_get(MINROWLEN);
  322.    
  323.    /* Initialise */
  324.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  325.    /* idx1 = idx2 = idx_out = 0; */
  326.    idx1    = sprow_idx(r1,j0);
  327.    idx2    = sprow_idx(r2,j0);
  328.    idx_out = sprow_idx(r_out,j0);
  329.    idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
  330.    idx2    = (idx2 < 0) ? -(idx2+2) : idx2;
  331.    idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
  332.    elt1    = &(r1->elt[idx1]);
  333.    elt2    = &(r2->elt[idx2]);
  334.    elt_out = &(r_out->elt[idx_out]);
  335.    
  336.    while ( idx1 < len1 || idx2 < len2 )
  337.    {
  338.       if ( idx_out >= len_out )
  339.       {   /* r_out is too small */
  340.      r_out->len = idx_out;
  341.      r_out = sprow_xpd(r_out,0,type);
  342.      len_out = r_out->maxlen;
  343.      elt_out = &(r_out->elt[idx_out]);
  344.       }
  345.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  346.       {
  347.      elt_out->col = elt1->col;
  348.      elt_out->val = elt1->val;
  349.      if ( idx2 < len2 && elt1->col == elt2->col )
  350.      {
  351.         elt_out->val += alpha*elt2->val;
  352.         elt2++;        idx2++;
  353.      }
  354.      elt1++;    idx1++;
  355.       }
  356.       else
  357.       {
  358.      elt_out->col = elt2->col;
  359.      elt_out->val = alpha*elt2->val;
  360.      elt2++;    idx2++;
  361.       }
  362.       elt_out++;    idx_out++;
  363.    }
  364.    r_out->len = idx_out;
  365.    
  366.    return r_out;
  367. }
  368.  
  369. /* sprow_add -- sets r_out <- r1 + r2
  370.    -- cannot be in situ
  371.    -- only for columns j0, j0+1, ...
  372.    -- type must be SPMAT or SPROW depending on
  373.       whether r_out is a row of a SPMAT structure
  374.       or a SPROW variable
  375.    -- returns r_out */
  376. SPROW    *sprow_add(r1,r2,j0,r_out,type)
  377. SPROW    *r1, *r2, *r_out;
  378. int    j0, type;
  379. {
  380.    int    idx1, idx2, idx_out, len1, len2, len_out;
  381.    row_elt    *elt1, *elt2, *elt_out;
  382.    
  383.    if ( ! r1 || ! r2 )
  384.      error(E_NULL,"sprow_add");
  385.    if ( r1 == r_out || r2 == r_out )
  386.      error(E_INSITU,"sprow_add");
  387.    if ( j0 < 0 )
  388.      error(E_BOUNDS,"sprow_add");
  389.    if ( ! r_out )
  390.      r_out = sprow_get(MINROWLEN);
  391.    
  392.    /* Initialise */
  393.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  394.    /* idx1 = idx2 = idx_out = 0; */
  395.    idx1    = sprow_idx(r1,j0);
  396.    idx2    = sprow_idx(r2,j0);
  397.    idx_out = sprow_idx(r_out,j0);
  398.    idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
  399.    idx2    = (idx2 < 0) ? -(idx2+2) : idx2;
  400.    idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
  401.    elt1    = &(r1->elt[idx1]);
  402.    elt2    = &(r2->elt[idx2]);
  403.    elt_out = &(r_out->elt[idx_out]);
  404.    
  405.    while ( idx1 < len1 || idx2 < len2 )
  406.    {
  407.       if ( idx_out >= len_out )
  408.       {   /* r_out is too small