home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / sprow < prev    next >
Encoding:
Text File  |  1994-01-13  |  17.3 KB  |  715 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.   Sparse rows package
  28.   See also: sparse.h, matrix.h
  29.   */
  30.  
  31. #include    <stdio.h>
  32. #include    <math.h>
  33. #include        <stdlib.h>
  34. #include    "sparse.h"
  35.  
  36.  
  37. static char    rcsid[] = "$Id: sprow.c,v 1.1 1994/01/13 05:35:36 des Exp $";
  38.  
  39. #define    MINROWLEN    10
  40.  
  41.  
  42. /* sprow_dump - prints relevant information about the sparse row r */
  43.  
  44. void sprow_dump(fp,r)
  45. FILE *fp;
  46. SPROW *r;
  47. {
  48.    int  j_idx;
  49.    row_elt *elts;
  50.    
  51.    fprintf(fp,"SparseRow dump:\n");
  52.    if ( ! r )
  53.    {       fprintf(fp,"*** NULL row ***\n");   return; }
  54.    
  55.    fprintf(fp,"row: len = %d, maxlen = %d, diag idx = %d\n",
  56.        r->len,r->maxlen,r->diag);
  57.    fprintf(fp,"element list @ 0x%lx\n",(long)(r->elt));
  58.    if ( ! r->elt )
  59.    {
  60.       fprintf(fp,"*** NULL element list ***\n");
  61.       return;
  62.    }
  63.    elts = r->elt;
  64.    for ( j_idx = 0; j_idx < r->len; j_idx++, elts++ )
  65.      fprintf(fp,"Col: %d, Val: %g, nxt_row = %d, nxt_idx = %d\n",
  66.          elts->col,elts->val,elts->nxt_row,elts->nxt_idx);
  67.    fprintf(fp,"\n");
  68. }
  69.  
  70.  
  71. /* sprow_idx -- get index into row for a given column in a given row
  72.    -- return -1 on error
  73.    -- return -(idx+2) where idx is index to insertion point */
  74. int    sprow_idx(r,col)
  75. SPROW    *r;
  76. int    col;
  77. {
  78.    register int        lo, hi, mid;
  79.    int            tmp;
  80.    register row_elt    *r_elt;
  81.    
  82.    /*******************************************
  83.      if ( r == (SPROW *)NULL )
  84.      return -1;
  85.      if ( col < 0 )
  86.      return -1;
  87.      *******************************************/
  88.    
  89.    r_elt = r->elt;
  90.    if ( r->len <= 0 )
  91.      return -2;
  92.    
  93.    /* try the hint */
  94.    /* if ( hint >= 0 && hint < r->len && r_elt[hint].col == col )
  95.       return hint; */
  96.    
  97.    /* otherwise use binary search... */
  98.    /* code from K&R Ch. 6, p. 125 */
  99.    lo = 0;        hi = r->len - 1;    mid = lo;
  100.    while ( lo <= hi )
  101.    {
  102.       mid = (hi + lo)/2;
  103.       if ( (tmp=r_elt[mid].col-col) > 0 )
  104.     hi = mid-1;
  105.       else if ( tmp < 0 )
  106.     lo = mid+1;
  107.       else /* tmp == 0 */
  108.     return mid;
  109.    }
  110.    tmp = r_elt[mid].col - col;
  111.    
  112.    if ( tmp > 0 )
  113.      return -(mid+2);    /* insert at mid   */
  114.    else /* tmp < 0 */
  115.      return -(mid+3);    /* insert at mid+1 */
  116. }
  117.  
  118.  
  119. /* sprow_get -- gets, initialises and returns a SPROW structure
  120.    -- max. length is maxlen */
  121. SPROW    *sprow_get(maxlen)
  122. int    maxlen;
  123. {
  124.    SPROW    *r;
  125.    
  126.    if ( maxlen < 0 )
  127.      error(E_NEG,"sprow_get");
  128.  
  129.    r = NEW(SPROW);
  130.    if ( ! r )
  131.      error(E_MEM,"sprow_get");
  132.    else if (mem_info_is_on()) {
  133.       mem_bytes(TYPE_SPROW,0,sizeof(SPROW));
  134.       mem_numvar(TYPE_SPROW,1);
  135.    }
  136.    r->elt = NEW_A(maxlen,row_elt);
  137.    if ( ! r->elt )
  138.      error(E_MEM,"sprow_get");
  139.    else if (mem_info_is_on()) {
  140.       mem_bytes(TYPE_SPROW,0,maxlen*sizeof(row_elt));
  141.    }
  142.    r->len = 0;
  143.    r->maxlen = maxlen;
  144.    r->diag = -1;
  145.    
  146.    return r;
  147. }
  148.  
  149.  
  150. /* sprow_xpd -- expand row by means of realloc()
  151.    -- type must be TYPE_SPMAT if r is a row of a SPMAT structure,
  152.       otherwise it must be TYPE_SPROW
  153.    -- returns r */
  154. SPROW    *sprow_xpd(r,n,type)
  155. SPROW    *r;
  156. int    n,type;
  157. {
  158.    int    newlen;
  159.    
  160.    if ( ! r ) {
  161.      r = NEW(SPROW);
  162.      if (! r ) 
  163.        error(E_MEM,"sprow_xpd");
  164.      else if ( mem_info_is_on()) {
  165.     if (type != TYPE_SPMAT && type != TYPE_SPROW)
  166.       warning(WARN_WRONG_TYPE,"sprow_xpd");
  167.     mem_bytes(type,0,sizeof(SPROW));
  168.     if (type == TYPE_SPROW)
  169.       mem_numvar(type,1);
  170.      }
  171.    }
  172.  
  173.    if ( ! r->elt )
  174.    {
  175.       r->elt = NEW_A((unsigned)n,row_elt);
  176.       if ( ! r->elt )
  177.     error(E_MEM,"sprow_xpd");
  178.       else if (mem_info_is_on()) {
  179.      mem_bytes(type,0,n*sizeof(row_elt));
  180.       }
  181.       r->len = 0;
  182.       r->maxlen = n;
  183.       return r;
  184.    }
  185.    if ( n <= r->len )
  186.      newlen = max(2*r->len + 1,MINROWLEN);
  187.    else
  188.      newlen = n;
  189.    if ( newlen <= r->maxlen )
  190.    {
  191.       MEM_ZERO((char *)(&(r->elt[r->len])),
  192.            (newlen-r->len)*sizeof(row_elt));
  193.       r->len = newlen;
  194.    }
  195.    else
  196.    {
  197.       if (mem_info_is_on()) {
  198.      mem_bytes(type,r->maxlen*sizeof(row_elt),
  199.              newlen*sizeof(row_elt)); 
  200.       }
  201.       r->elt = RENEW(r->elt,newlen,row_elt);
  202.       if ( ! r->elt )
  203.     error(E_MEM,"sprow_xpd");
  204.       r->maxlen = newlen;
  205.       r->len = newlen;
  206.    }
  207.    
  208.    return r;
  209. }
  210.  
  211. /* sprow_resize -- resize a SPROW variable by means of realloc()
  212.    -- n is a new size
  213.    -- returns r */
  214. SPROW    *sprow_resize(r,n,type)
  215. SPROW    *r;
  216. int    n,type;
  217. {
  218.    if (n < 0)
  219.      error(E_NEG,"sprow_resize");
  220.  
  221.    if ( ! r ) 
  222.      return sprow_get(n);
  223.    
  224.    if (n == r->len)
  225.      return r;
  226.  
  227.    if ( ! r->elt )
  228.    {
  229.       r->elt = NEW_A((unsigned)n,row_elt);
  230.       if ( ! r->elt )
  231.     error(E_MEM,"sprow_resize");
  232.       else if (mem_info_is_on()) {
  233.      mem_bytes(type,0,n*sizeof(row_elt));
  234.       }
  235.       r->maxlen = r->len = n;
  236.       return r;
  237.    }
  238.  
  239.    if ( n <= r->maxlen )
  240.      r->len = n;
  241.    else
  242.    {
  243.       if (mem_info_is_on()) {
  244.      mem_bytes(type,r->maxlen*sizeof(row_elt),
  245.            n*sizeof(row_elt)); 
  246.       }
  247.       r->elt = RENEW(r->elt,n,row_elt);
  248.       if ( ! r->elt )
  249.     error(E_MEM,"sprow_resize");
  250.       r->maxlen = r->len = n;
  251.    }
  252.    
  253.    return r;
  254. }
  255.  
  256.  
  257. /* release a row of a matrix */
  258. int sprow_free(r)
  259. SPROW    *r;
  260. {
  261.    if ( ! r )
  262.      return -1;
  263.  
  264.    if (mem_info_is_on()) {
  265.       mem_bytes(TYPE_SPROW,sizeof(SPROW),0);
  266.       mem_numvar(TYPE_SPROW,-1);
  267.    }
  268.    
  269.    if ( r->elt )
  270.    {
  271.       if (mem_info_is_on()) {
  272.      mem_bytes(TYPE_SPROW,r->maxlen*sizeof(row_elt),0);
  273.       }
  274.       free((char *)r->elt);
  275.    }
  276.    free((char *)r);
  277.    return 0;
  278. }
  279.  
  280.  
  281. /* sprow_merge -- merges r1 and r2 into r_out
  282.    -- cannot be done in-situ
  283.    -- type must be SPMAT or SPROW depending on
  284.       whether r_out is a row of a SPMAT structure
  285.       or a SPROW variable
  286.    -- returns r_out */
  287. SPROW    *sprow_merge(r1,r2,r_out,type)
  288. SPROW    *r1, *r2, *r_out;
  289. int type;
  290. {
  291.    int    idx1, idx2, idx_out, len1, len2, len_out;
  292.    row_elt    *elt1, *elt2, *elt_out;
  293.    
  294.    if ( ! r1 || ! r2 )
  295.      error(E_NULL,"sprow_merge");
  296.    if ( ! r_out )
  297.      r_out = sprow_get(MINROWLEN);
  298.    if ( r1 == r_out || r2 == r_out )
  299.      error(E_INSITU,"sprow_merge");
  300.    
  301.    /* Initialise */
  302.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  303.    idx1 = idx2 = idx_out = 0;
  304.    elt1 = r1->elt;    elt2 = r2->elt;    elt_out = r_out->elt;
  305.    
  306.    while ( idx1 < len1 || idx2 < len2 )
  307.    {
  308.       if ( idx_out >= len_out )
  309.       {   /* r_out is too small */
  310.      r_out->len = idx_out;
  311.      r_out = sprow_xpd(r_out,0,type);
  312.      len_out = r_out->len;
  313.      elt_out = &(r_out->elt[idx_out]);
  314.       }
  315.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  316.       {
  317.      elt_out->col = elt1->col;
  318.      elt_out->val = elt1->val;
  319.      if ( elt1->col == elt2->col && idx2 < len2 )
  320.      {    elt2++;        idx2++;    }
  321.      elt1++;    idx1++;
  322.       }
  323.       else
  324.       {
  325.      elt_out->col = elt2->col;
  326.      elt_out->val = elt2->val;
  327.      elt2++;    idx2++;
  328.       }
  329.       elt_out++;    idx_out++;
  330.    }
  331.    r_out->len = idx_out;
  332.    
  333.    return r_out;
  334. }
  335.  
  336. /* sprow_copy -- copies r1 and r2 into r_out
  337.    -- cannot be done in-situ
  338.    -- type must be SPMAT or SPROW depending on
  339.       whether r_out is a row of a SPMAT structure
  340.       or a SPROW variable
  341.    -- returns r_out */
  342. SPROW    *sprow_copy(r1,r2,r_out,type)
  343. SPROW    *r1, *r2, *r_out;
  344. int type;
  345. {
  346.    int    idx1, idx2, idx_out, len1, len2, len_out;
  347.    row_elt    *elt1, *elt2, *elt_out;
  348.    
  349.    if ( ! r1 || ! r2 )
  350.      error(E_NULL,"sprow_copy");
  351.    if ( ! r_out )
  352.      r_out = sprow_get(MINROWLEN);
  353.    if ( r1 == r_out || r2 == r_out )
  354.      error(E_INSITU,"sprow_copy");
  355.    
  356.    /* Initialise */
  357.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  358.    idx1 = idx2 = idx_out = 0;
  359.    elt1 = r1->elt;    elt2 = r2->elt;    elt_out = r_out->elt;
  360.    
  361.    while ( idx1 < len1 || idx2 < len2 )
  362.    {
  363.       while ( idx_out >= len_out )
  364.       {   /* r_out is too small */
  365.      r_out->len = idx_out;
  366.      r_out = sprow_xpd(r_out,0,type);
  367.      len_out = r_out->maxlen;
  368.      elt_out = &(r_out->elt[idx_out]);
  369.       }
  370.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  371.       {
  372.      elt_out->col = elt1->col;
  373.      elt_out->val = elt1->val;
  374.      if ( elt1->col == elt2->col && idx2 < len2 )
  375.      {    elt2++;        idx2++;    }
  376.      elt1++;    idx1++;
  377.       }
  378.       else
  379.       {
  380.      elt_out->col = elt2->col;
  381.      elt_out->val = 0.0;
  382.      elt2++;    idx2++;
  383.       }
  384.       elt_out++;    idx_out++;
  385.    }
  386.    r_out->len = idx_out;
  387.    
  388.    return r_out;
  389. }
  390.  
  391. /* sprow_mltadd -- sets r_out <- r1 + alpha.r2
  392.    -- cannot be in situ
  393.    -- only for columns j0, j0+1, ...
  394.    -- type must be SPMAT or SPROW depending on
  395.       whether r_out is a row of a SPMAT structure
  396.       or a SPROW variable
  397.    -- returns r_out */
  398. SPROW    *sprow_mltadd(r1,r2,alpha,j0,r_out,type)
  399. SPROW    *r1, *r2, *r_out;
  400. double    alpha;
  401. int    j0, type;
  402. {
  403.    int    idx1, idx2, idx_out, len1, len2, len_out;
  404.    row_elt    *elt1, *elt2, *elt_out;
  405.    
  406.    if ( ! r1 || ! r2 )
  407.      error(E_NULL,"sprow_mltadd");
  408.    if ( r1 == r_out || r2 == r_out )
  409.      error(E_INSITU,"sprow_mltadd");
  410.    if ( j0 < 0 )
  411.      error(E_BOUNDS,"sprow_mltadd");
  412.    if ( ! r_out )
  413.      r_out = sprow_get(MINROWLEN);
  414.    
  415.    /* Initialise */
  416.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  417.    /* idx1 = idx2 = idx_out = 0; */
  418.    idx1    = sprow_idx(r1,j0);
  419.    idx2    = sprow_idx(r2,j0);
  420.    idx_out = sprow_idx(r_out,j0);
  421.    idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
  422.    idx2    = (idx2 < 0) ? -(idx2+2) : idx2;
  423.    idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
  424.    elt1    = &(r1->elt[idx1]);
  425.    elt2    = &(r2->elt[idx2]);
  426.    elt_out = &(r_out->elt[idx_out]);
  427.    
  428.    while ( idx1 < len1 || idx2 < len2 )
  429.    {
  430.       if ( idx_out >= len_out )
  431.       {   /* r_out is too small */
  432.      r_out->len = idx_out;
  433.      r_out = sprow_xpd(r_out,0,type);
  434.      len_out = r_out->maxlen;
  435.      elt_out = &(r_out->elt[idx_out]);
  436.       }
  437.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  438.       {
  439.      elt_out->col = elt1->col;
  440.      elt_out->val = elt1->val;
  441.      if ( idx2 < len2 && elt1->col == elt2->col )
  442.      {
  443.         elt_out->val += alpha*elt2->val;
  444.         elt2++;        idx2++;
  445.      }
  446.      elt1++;    idx1++;
  447.       }
  448.       else
  449.       {
  450.      elt_out->col = elt2->col;
  451.      elt_out->val = alpha*elt2->val;
  452.      elt2++;    idx2++;
  453.       }
  454.       elt_out++;    idx_out++;
  455.    }
  456.    r_out->len = idx_out;
  457.    
  458.    return r_out;
  459. }
  460.  
  461. /* sprow_add -- sets r_out <- r1 + r2
  462.    -- cannot be in situ
  463.    -- only for columns j0, j0+1, ...
  464.    -- type must be SPMAT or SPROW depending on
  465.       whether r_out is a row of a SPMAT structure
  466.       or a SPROW variable
  467.    -- returns r_out */
  468. SPROW    *sprow_add(r1,r2,j0,r_out,type)
  469. SPROW    *r1, *r2, *r_out;
  470. int    j0, type;
  471. {
  472.    int    idx1, idx2, idx_out, len1, len2, len_out;
  473.    row_elt    *elt1, *elt2, *elt_out;
  474.    
  475.    if ( ! r1 || ! r2 )
  476.      error(E_NULL,"sprow_add");
  477.    if ( r1 == r_out || r2 == r_out )
  478.      error(E_INSITU,"sprow_add");
  479.    if ( j0 < 0 )
  480.      error(E_BOUNDS,"sprow_add");
  481.    if ( ! r_out )
  482.      r_out = sprow_get(MINROWLEN);
  483.    
  484.    /* Initialise */
  485.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  486.    /* idx1 = idx2 = idx_out = 0; */
  487.    idx1    = sprow_idx(r1,j0);
  488.    idx2    = sprow_idx(r2,j0);
  489.    idx_out = sprow_idx(r_out,j0);
  490.    idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
  491.    idx2    = (idx2 < 0) ? -(idx2+2) : idx2;
  492.    idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
  493.    elt1    = &(r1->elt[idx1]);
  494.    elt2    = &(r2->elt[idx2]);
  495.    elt_out = &(r_out->elt[idx_out]);
  496.    
  497.    while ( idx1 < len1 || idx2 < len2 )
  498.    {
  499.       if ( idx_out >= len_out )
  500.       {   /* r_out is too small */
  501.      r_out->len = idx_out;
  502.      r_out = sprow_xpd(r_out,0,type);
  503.      len_out = r_out->maxlen;
  504.      elt_out = &(r_out->elt[idx_out]);
  505.       }
  506.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  507.       {
  508.      elt_out->col = elt1->col;
  509.      elt_out->val = elt1->val;
  510.      if ( idx2 < len2 && elt1->col == elt2->col )
  511.      {
  512.         elt_out->val += elt2->val;
  513.         elt2++;        idx2++;
  514.      }
  515.      elt1++;    idx1++;
  516.       }
  517.       else
  518.       {
  519.      elt_out->col = elt2->col;
  520.      elt_out->val = elt2->val;
  521.      elt2++;    idx2++;
  522.       }
  523.       elt_out++;    idx_out++;
  524.    }
  525.    r_out->len = idx_out;
  526.    
  527.    return r_out;
  528. }
  529.  
  530. /* sprow_sub -- sets r_out <- r1 - r2
  531.    -- cannot be in situ
  532.    -- only for columns j0, j0+1, ...
  533.    -- type must be SPMAT or SPROW depending on
  534.       whether r_out is a row of a SPMAT structure
  535.       or a SPROW variable
  536.    -- returns r_out */
  537. SPROW    *sprow_sub(r1,r2,j0,r_out,type)
  538. SPROW    *r1, *r2, *r_out;
  539. int    j0, type;
  540. {
  541.    int    idx1, idx2, idx_out, len1, len2, len_out;
  542.    row_elt    *elt1, *elt2, *elt_out;
  543.    
  544.    if ( ! r1 || ! r2 )
  545.      error(E_NULL,"sprow_sub");
  546.    if ( r1 == r_out || r2 == r_out )
  547.      error(E_INSITU,"sprow_sub");
  548.    if ( j0 < 0 )
  549.      error(E_BOUNDS,"sprow_sub");
  550.    if ( ! r_out )
  551.      r_out = sprow_get(MINROWLEN);
  552.    
  553.    /* Initialise */
  554.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  555.    /* idx1 = idx2 = idx_out = 0; */
  556.    idx1    = sprow_idx(r1,j0);
  557.    idx2    = sprow_idx(r2,j0);
  558.    idx_out = sprow_idx(r_out,j0);
  559.    idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
  560.    idx2    = (idx2 < 0) ? -(idx2+2) : idx2;
  561.    idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
  562.    elt1    = &(r1->elt[idx1]);
  563.    elt2    = &(r2->elt[idx2]);
  564.    elt_out = &(r_out->elt[idx_out]);
  565.    
  566.    while ( idx1 < len1 || idx2 < len2 )
  567.    {
  568.       if ( idx_out >= len_out )
  569.       {   /* r_out is too small */
  570.      r_out->len = idx_out;
  571.      r_out = sprow_xpd(r_out,0,type);
  572.      len_out = r_out->maxlen;
  573.      elt_out = &(r_out->elt[idx_out]);
  574.       }
  575.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  576.       {
  577.      elt_out->col = elt1->col;
  578.      elt_out->val = elt1->val;
  579.      if ( idx2 < len2 && elt1->col == elt2->col )
  580.      {
  581.         elt_out->val -= elt2->val;
  582.         elt2++;        idx2++;
  583.      }
  584.      elt1++;    idx1++;
  585.       }
  586.       else
  587.       {
  588.      elt_out->col = elt2->col;
  589.      elt_out->val = -elt2->val;
  590.      elt2++;    idx2++;
  591.       }
  592.       elt_out++;    idx_out++;
  593.    }
  594.    r_out->len = idx_out;
  595.    
  596.    return r_out;
  597. }
  598.  
  599.  
  600. /* sprow_smlt -- sets r_out <- alpha*r1 
  601.    -- can be in situ
  602.    -- only for columns j0, j0+1, ...
  603.    -- returns r_out */
  604. SPROW    *sprow_smlt(r1,alpha,j0,r_out,type)
  605. SPROW    *r1, *r_out;
  606. double    alpha;
  607. int    j0, type;
  608. {
  609.    int    idx1, idx_out, len1;
  610.    row_elt    *elt1, *elt_out;
  611.    
  612.    if ( ! r1 )
  613.      error(E_NULL,"sprow_smlt");
  614.    if ( j0 < 0 )
  615.      error(E_BOUNDS,"sprow_smlt");
  616.    if ( ! r_out )
  617.      r_out = sprow_get(MINROWLEN);
  618.    
  619.    /* Initialise */
  620.    len1 = r1->len;
  621.    idx1    = sprow_idx(r1,j0);
  622.    idx_out = sprow_idx(r_out,j0);
  623.    idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
  624.    idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
  625.    elt1    = &(r1->elt[idx1]);
  626.  
  627.    r_out = sprow_resize(r_out,idx_out+len1-idx1,type);  
  628.    elt_out = &(r_out->elt[idx_out]);
  629.  
  630.    for ( ; idx1 < len1; elt1++,elt_out++,idx1++,idx_out++ )
  631.    {
  632.       elt_out->col = elt1->col;
  633.       elt_out->val = alpha*elt1->val;
  634.    }
  635.  
  636.    r_out->len = idx_out;
  637.  
  638.    return r_out;
  639. }
  640.  
  641.   
  642. /* sprow_foutput -- print a representation of r on stream fp */
  643. void    sprow_foutput(fp,r)
  644. FILE    *fp;
  645. SPROW    *r;
  646. {
  647.    int    i, len;
  648.    row_elt    *e;
  649.    
  650.    if ( ! r )
  651.    {
  652.       fprintf(fp,"SparseRow: **** NULL ****\n");
  653.       return;
  654.    }
  655.    len = r->len;
  656.    fprintf(fp,"SparseRow: length: %d\n",len);
  657.    for ( i = 0, e = r->elt; i < len; i++, e++ )
  658.      fprintf(fp,"Column %d: %g, next row: %d, next index %d\n",
  659.          e->col, e->val, e->nxt_row, e->nxt_idx);
  660. }
  661.  
  662.  
  663. /* sprow_set_val -- sets the j-th column entry of the sparse row r
  664.    -- Note: destroys the usual column & row access paths */
  665. double  sprow_set_val(r,j,val)
  666. SPROW   *r;
  667. int     j;
  668. double  val;
  669. {
  670.    int  idx, idx2, new_len;
  671.    
  672.    if ( ! r )
  673.      error(E_NULL,"sprow_set_val");
  674.    
  675.    idx = sprow_idx(r,j);
  676.    if ( idx >= 0 )
  677.    {    r->elt[idx].val = val;  return val;     }
  678.    /* else */ if ( idx < -1 )
  679.    {
  680.       /* shift & insert new value */
  681.       idx = -(idx+2);   /* this is the intended insertion index */
  682.       if ( r->len >= r->maxlen )
  683.       {
  684.          r->len = r->maxlen;
  685.          new_len = max(2*r->maxlen+1,5);
  686.          if (mem_info_is_on()) {
  687.             mem_bytes(TYPE_SPROW,r->maxlen*sizeof(row_elt),
  688.                         new_len*sizeof(row_elt)); 
  689.          }
  690.          
  691.          r->elt = RENEW(r->elt,new_len,row_elt);
  692.          if ( ! r->elt )        /* can't allocate */
  693.            error(E_MEM,"sprow_set_val");
  694.          r->maxlen = 2*r->maxlen+1;
  695.       }
  696.       for ( idx2 = r->len-1; idx2 >= idx; idx2-- )
  697.         MEM_COPY((char *)(&(r->elt[idx2])),
  698.                  (char *)(&(r->elt[idx2+1])),sizeof(row_elt));
  699.       /************************************************************
  700.         if ( idx < r->len )
  701.         MEM_COPY((char *)(&(r->elt[idx])),(char *)(&(r->elt[idx+1])),
  702.         (r->len-idx)*sizeof(row_elt));
  703.         ************************************************************/
  704.       r->len++;
  705.       r->elt[idx].col = j;
  706.       r->elt[idx].nxt_row = -1;
  707.       r->elt[idx].nxt_idx = -1;
  708.       return r->elt[idx].val = val;
  709.    }
  710.    /* else -- idx == -1, error in index/matrix! */
  711.    return 0.0;
  712. }
  713.  
  714.  
  715.