home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / sprow < prev    next >
Encoding:
Text File  |  1994-01-13  |  17.3 KB  |  769 lines
  1.  sprow_xpd(r_out,0,type);
  2.      len_out = r_out->maxlen;
  3.      elt_out = &(r_out->elt[idx_out]);
  4.       }
  5.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  6.       {
  7.      elt_out->col = elt1->col;
  8.      elt_out->val = elt1->val;
  9.      if ( idx2 < len2 && elt1->col == elt2->col )
  10.      {
  11.         elt_out->val += elt2->val;
  12.         elt2++;        idx2++;
  13.      }
  14.      elt1++;    idx1++;
  15.       }
  16.       else
  17.       {
  18.      elt_out->col = elt2->col;
  19.      elt_out->val = elt2->val;
  20.      elt2++;    idx2++;
  21.       }
  22.       elt_out++;    idx_out++;
  23.    }
  24.    r_out->len = idx_out;
  25.    
  26.    return r_out;
  27. }
  28.  
  29. /* sprow_sub -- sets r_out <- r1 - r2
  30.    -- cannot be in situ
  31.    -- only for columns j0, j0+1, ...
  32.    -- type must be SPMAT or SPROW depending on
  33.       whether r_out is a row of a SPMAT structure
  34.       or a SPROW variable
  35.    -- returns r_out */
  36. SPROW    *sprow_sub(r1,r2,j0,r_out,type)
  37. SPROW    *r1, *r2, *r_out;
  38. int    j0, type;
  39. {
  40.    int    idx1, idx2, idx_out, len1, len2, len_out;
  41.    row_elt    *elt1, *elt2, *elt_out;
  42.    
  43.    if ( ! r1 || ! r2 )
  44.      error(E_NULL,"sprow_sub");
  45.    if ( r1 == r_out || r2 == r_out )
  46.      error(E_INSITU,"sprow_sub");
  47.    if ( j0 < 0 )
  48.      error(E_BOUNDS,"sprow_sub");
  49.    if ( ! r_out )
  50.      r_out = sprow_get(MINROWLEN);
  51.    
  52.    /* Initialise */
  53.    len1 = r1->len;    len2 = r2->len;    len_out = r_out->maxlen;
  54.    /* idx1 = idx2 = idx_out = 0; */
  55.    idx1    = sprow_idx(r1,j0);
  56.    idx2    = sprow_idx(r2,j0);
  57.    idx_out = sprow_idx(r_out,j0);
  58.    idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
  59.    idx2    = (idx2 < 0) ? -(idx2+2) : idx2;
  60.    idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
  61.    elt1    = &(r1->elt[idx1]);
  62.    elt2    = &(r2->elt[idx2]);
  63.    elt_out = &(r_out->elt[idx_out]);
  64.    
  65.    while ( idx1 < len1 || idx2 < len2 )
  66.    {
  67.       if ( idx_out >= len_out )
  68.       {   /* r_out is too small */
  69.      r_out->len = idx_out;
  70.      r_out = sprow_xpd(r_out,0,type);
  71.      len_out = r_out->maxlen;
  72.      elt_out = &(r_out->elt[idx_out]);
  73.       }
  74.       if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
  75.       {
  76.      elt_out->col = elt1->col;
  77.      elt_out->val = elt1->val;
  78.      if ( idx2 < len2 && elt1->col == elt2->col )
  79.      {
  80.         elt_out->val -= elt2->val;
  81.         elt2++;        idx2++;
  82.      }
  83.      elt1++;    idx1++;
  84.       }
  85.       else
  86.       {
  87.      elt_out->col = elt2->col;
  88.      elt_out->val = -elt2->val;
  89.      elt2++;    idx2++;
  90.       }
  91.       elt_out++;    idx_out++;
  92.    }
  93.    r_out->len = idx_out;
  94.    
  95.    return r_out;
  96. }
  97.  
  98.  
  99. /* sprow_smlt -- sets r_out <- alpha*r1 
  100.    -- can be in situ
  101.    -- only for columns j0, j0+1, ...
  102.    -- returns r_out */
  103. SPROW    *sprow_smlt(r1,alpha,j0,r_out,type)
  104. SPROW    *r1, *r_out;
  105. double    alpha;
  106. int    j0, type;
  107. {
  108.    int    idx1, idx_out, len1;
  109.    row_elt    *elt1, *elt_out;
  110.    
  111.    if ( ! r1 )
  112.      error(E_NULL,"sprow_smlt");
  113.    if ( j0 < 0 )
  114.      error(E_BOUNDS,"sprow_smlt");
  115.    if ( ! r_out )
  116.      r_out = sprow_get(MINROWLEN);
  117.    
  118.    /* Initialise */
  119.    len1 = r1->len;
  120.    idx1    = sprow_idx(r1,j0);
  121.    idx_out = sprow_idx(r_out,j0);
  122.    idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
  123.    idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
  124.    elt1    = &(r1->elt[idx1]);
  125.  
  126.    r_out = sprow_resize(r_out,idx_out+len1-idx1,type);  
  127.    elt_out = &(r_out->elt[idx_out]);
  128.  
  129.    for ( ; idx1 < len1; elt1++,elt_out++,idx1++,idx_out++ )
  130.    {
  131.       elt_out->col = elt1->col;
  132.       elt_out->val = alpha*elt1->val;
  133.    }
  134.  
  135.    r_out->len = idx_out;
  136.  
  137.    return r_out;
  138. }
  139.  
  140.   
  141. /* sprow_foutput -- print a representation of r on stream fp */
  142. void    sprow_foutput(fp,r)
  143. FILE    *fp;
  144. SPROW    *r;
  145. {
  146.    int    i, len;
  147.    row_elt    *e;
  148.    
  149.    if ( ! r )
  150.    {
  151.       fprintf(fp,"SparseRow: **** NULL ****\n");
  152.       return;
  153.    }
  154.    len = r->len;
  155.    fprintf(fp,"SparseRow: length: %d\n",len);
  156.    for ( i = 0, e = r->elt; i < len; i++, e++ )
  157.      fprintf(fp,"Column %d: %g, next row: %d, next index %d\n",
  158.          e->col, e->val, e->nxt_row, e->nxt_idx);
  159. }
  160.  
  161.  
  162. /* sprow_set_val -- sets the j-th column entry of the sparse row r
  163.    -- Note: destroys the usual column & row access paths */
  164. double  sprow_set_val(r,j,val)
  165. SPROW   *r;
  166. int     j;
  167. double  val;
  168. {
  169.    int  idx, idx2, new_len;
  170.    
  171.    if ( ! r )
  172.      error(E_NULL,"sprow_set_val");
  173.    
  174.    idx = sprow_idx(r,j);
  175.    if ( idx >= 0 )
  176.    {    r->elt[idx].val = val;  return val;     }
  177.    /* else */ if ( idx < -1 )
  178.    {
  179.       /* shift & insert new value */
  180.       idx = -(idx+2);   /* this is the intended insertion index */
  181.       if ( r->len >= r->maxlen )
  182.       {
  183.          r->len = r->maxlen;
  184.          new_len = max(2*r->maxlen+1,5);
  185.          if (mem_info_is_on()) {
  186.             mem_bytes(TYPE_SPROW,r->maxlen*sizeof(row_elt),
  187.                         new_len*sizeof(row_elt)); 
  188.          }
  189.          
  190.          r->elt = RENEW(r->elt,new_len,row_elt);
  191.          if ( ! r->elt )        /* can't allocate  = 0  i < m; i++ )
  192.    {
  193.       r = A->row+i;
  194.       max_idx = r->len;
  195.       elts    = r->elt;
  196.       tmp = x_ve[i];
  197.       for ( j_idx = 0; j_idx < max_idx; j_idx++, elts++ )
  198.     out_ve[elts->col] += elts->val*tmp;
  199.    }
  200.    
  201.    return out;
  202. of
  203.  
  204. /* sp_ge;
  205.       /*******matr**
  206.    -- len is number of elements available for each row without
  207.    allocating further memory */
  208. SPMAT    *sp_get(m,n,maxlen)
  209. int    m, n, maxlen;
  210. {
  211.    SPMAT    *A;
  212.    SPROW    *rows;
  213.    int    i;
  214.    
  215.    if ( m < 0 || n < 0 )
  216.      error(E_NEG,"sp_get");
  217.  
  218.    maxlen = max(maxlen,1);
  219.    
  220.    A = NEW(SPMAT);
  221.    if ( ! A )        /* can't allocate */
  222.      error(E_MEM,"sp_get");
  223.    else if (mem_info_is_on()) {
  224.       mem_bytes(TYPE_SPMAT,0,sizeof(SPMAT));
  225.       mem_numvar(TYPE_SPMAT,1);
  226.    }
  227.    /* fprintf(stderr,"Have SPMAT structure\n"); */
  228.    
  229.    A->row = rows = NEW_A(m,SPROW);
  230.    if ( ! A->row )        /* can't allocate */
  231.      error(E_MEM,"sp_get");
  232.    else if (mem_info_is_on()) {
  233.       mem_bytes(TYPE_SPMAT,0,m*sizeof(SPROW));
  234.    }
  235.    /* fprintf(stderr,"Have row structure array\n"); */
  236.    
  237.    A->start_row = NEW_A(n,int);
  238.    A->start_idx = NEW_A(n,int);
  239.    if ( ! A->start_row || ! A->start_idx )    /* can't allocate */
  240.      error(E_MEM,"sp_get");
  241.    else if (mem_info_is_on()) {
  242.       mem_bytes(TYPE_SPMAT,0,2*n*sizeof(int));
  243.    }
  244.    for ( i = 0; i < n; i++ )
  245.      A->start_row[i] = A->start_idx[i] = -1;
  246.    /* fprintf(stderr,"Have start_row array\n"); */
  247.    
  248.    A->m = A->max_m = m;
  249.    A->n = A->max_n = n;
  250.    
  251.    for ( i = 0; i < m; i++, rows++ )
  252.    {
  253.       rows->elt = NEW_A(maxlen,row_elt);
  254.       if ( ! rows->elt )
  255.     error(E_MEM,"sp_get");
  256.       else if (mem_info_is_on()) {
  257.      mem_bytes(TYPE_SPMAT,0,maxlen*sizeof(row_elt));
  258.       }
  259.       /* fprintf(stderr,"Have row %d element array\n",i); */
  260.       rows->len = 0;
  261.      s a
  262. **      distribution fee considered a charge.
  263. **
  264. ***************************************************************************/
  265.  
  266.  
  267. /*
  268.     Sparse matrix swap and permutation routines
  269.     Modified Mon 09th Nov 1992, 08:50:54 PM
  270.     to use Karen George's suggestion to use unordered rows
  271. */
  272.  
  273. static    char    rcsid[] = "$Id: spswap.c,v 1.3 1994/01/13 05:44:43 des Exp $";
  274.  
  275. #include    <stdio.h>
  276. #include    <math.h>
  277. #include    "matrix.h"
  278. #include    "sparse.h"
  279. #include        "sparse2.h"
  280.  
  281.  
  282. #define    btos(x)    ((x) ? "TRUE" : "FALSE")
  283.  
  284. /* scan_to -- updates scan (int) vectors to point to the last row in each
  285.     column with row # <= max_row, if any */
  286. void    scan_to(A, scan_row, scan_idx, col_list, max_row)
  287. SPMAT    *A;
  288. IVEC    *scan_row, *scan_idx, *col_list;
  289. int    max_row;
  290. {
  291.     int        col, idx, j_idx, row_num;
  292.     SPROW    *r;
  293.     row_elt    *e;
  294.  
  295.     if ( ! A || ! scan_row || ! scan_idx || ! col_list )
  296.     error(E_NULL,"scan_to");
  297.     if ( scan_row->dim != scan_idx->dim || scan_idx->dim != col_list->dim )
  298.     error(E_SIZES,"scan_to");
  299.  
  300.     if ( max_row < 0 )
  301.     return;
  302.  
  303.     if ( ! A->flag_col )
  304.     sp_col_access(A);
  305.  
  306.     for ( j_idx = 0; j_idx < scan_row->dim; j_idx++ )
  307.     {
  308.     row_num = scan_row->ive[j_idx];
  309.     idx = scan_idx->ive[j_idx];
  310.     col = col_list->ive[j_idx];
  311.  
  312.     if ( col < 0 || col >= A->n )
  313.         error(E_BOUNDS,"scan_to");
  314.     if ( row_num < 0 )
  315.     {
  316.         idx = col;
  317.         continue;
  318.     }
  319.     r = &(A->row[row_num]);
  320.     if ( idx < 0 )
  321.         error(E_INTERN,"scan_to");
  322.     e = &(r->elt[idx]);
  323.     if ( e->col != col )
  324.         error(E_INTERN,"scan_to");
  325.     if ( idx < 0 )
  326.     {
  327.         printf("scan_to: row_num = %d, idx = %d, col = %d\n",
  328.            row_num, idx, col);
  329.         error(E_INTERN,"scan_to");
  330.     }
  331.     /* if ( e->nxt_row <= max_row )
  332.         chase_col(A, col, &row_num, &idx, max_row); */
  333.     while ( e->nxt_row >= 0 && e->nxt_row <= max_row )
  334.     {
  335.         row_num = e->nxt_row;
  336.         idx = e->nxt_idx;
  337.         e = &(A->row[row_num].elt[idx]);
  338.     }
  339.         
  340.     /* printf("scan_to: computed j_idx = %d, row_num = %d, idx = %d\n",
  341.            j_idx, row_num, idx); */
  342.     scan_row->ive[j_idx] = row_num;
  343.     scan_idx->ive[j_idx] = idx;
  344.     }
  345. }
  346.  
  347. /* patch_col -- patches column access paths for fill-in */
  348. void patch_col(A, col, old_row, old_idx, row_num, idx)
  349. SPMAT    *A;
  350. int    col, old_row, old_idx, row_num, idx;
  351. {
  352.     SPROW    *r;
  353.     row_elt    *e;
  354.     
  355.     if ( old_row >= 0 )
  356.     {
  357.     r = &(A->row[old_row]);
  358.     old_idx = sprow_idx2(r,col,old_idx);
  359.     e = &(r->elt[old_idx]);
  360.     e->nxt_row = row_num;
  361.     e->nxt_idx = idx;
  362.     }
  363.     else
  364.     {
  365.     A->start_row[col] = row_num;
  366.     A->start_idx[col] = idx;
  367.     }
  368. }
  369.  
  370. /* chase_col -- chases column access path in column col, starting with
  371.    row_num and idx, to find last row # in this column <= max_row
  372.    -- row_num is returned; idx is also set by this routine
  373.    -- assumes that the column access paths (possibly without the
  374.    nxt_idx fields) are set up */
  375. row_elt *chase_col(A, col, row_num, idx, max_row)
  376. SPMAT    *A;
  377. int    col, *row_num, *idx, max_row;
  378. {
  379.     int        old_idx, old_row, tmp_idx, tmp_row;
  380.     SPROW    *r;
  381.     row_elt    *e;
  382.     
  383.     if ( col < 0 || col >= A->n )
  384.     error(E_BOUNDS,"chase_col");
  385.     tmp_row = *row_num;
  386.     if ( tmp_row < 0 )
  387.     {
  388.     if ( A->start_row[col] > max_row )
  389.     {
  390.         tmp_row = -1;
  391.         tmp_idx = col;
  392.         return (row_elt *)NULL;
  393.     }
  394.     else
  395.     {
  396.         tmp_row = A->start_row[col];
  397.         tmp_idx = A->start_idx[col];
  398.     }
  399.     }
  400.     else
  401.     tmp_idx = *idx;
  402.     
  403.     old_row = tmp_row;
  404.     old_idx = tmp_idx;
  405.     while ( tmp_row >= 0 && tmp_row < max_row )
  406.     {
  407.     r = &(A->row[tmp_row]);
  408.     /* tmp_idx = sprow_idx2(r,col,tmp_idx); */
  409.     if ( tmp_idx < 0 || tmp_idx >= r->len ||
  410.          r->elt[tmp_idx].col != col )
  411.     {
  412. #ifdef DEBUG
  413.         printf("chase_col:error: col = %d, row # = %d, idx = %d\n",
  414.            col, tmp_row, tmp_idx);
  415.         printf("chase_col:error: old_row = %d, old_idx = %d\n",
  416.            old_row, old_idx);
  417.         printf("chase_col:error: A =\n");
  418.         sp_dump(stdout,A);
  419. #endif
  420.         error(E_INTERN,"chase_col");
  421.     }
  422.     e = &(r->elt[tmp_idx]);
  423.     old_row = tmp_row;
  424.     old_idx = tmp_idx;
  425.     tmp_row = e->nxt_row;
  426.     tmp_idx = e->nxt_idx;
  427.     }
  428. A == SMNULL )
  429.      error(E_NULL,"sp_diag_access");
  430.    
  431.    m = A->m;
  432.    
  433.    row = A->row;
  434.    for ( i = 0; i < m; i++, row++ )
  435.      row->diag = sprow_idx(row,i);
  436.    
  437.    A->flag_diag = TRUE;
  438.    
  439.    return A;
  440. }
  441.  
  442. /* sp_m2dense -- convert a sparse matrix to a dense one */
  443. MAT    *sp_m2dense(A,out)
  444. SPMAT    *A;
  445. MAT    *out;
  446. {
  447.    int    i, j_idx;
  448.    SPROW    *row;
  449.    row_elt    *elt;
  450.    
  451.    if ( ! A )
  452.      error(E_NULL,"sp_m2dense");
  453.    if ( ! out || out->m < A->m || out->n < A->n )
  454.      out = m_get(A->m,A->n);
  455.    
  456.    m_zero(out);
  457.    for ( i = 0; i < A->m; i++ )
  458.    {
  459.       row = &(A->row[i]);
  460.       elt = row->elt;
  461.       for ( j_idx = 0; j_idx < row->len; j_idx++, elt++ )
  462.     out->me[i][elt->col] = elt->val;
  463.    }
  464.    
  465.    return out;
  466. }
  467.  
  468.  
  469. /*  C = A+B, can be in situ */
  470. SPMAT *sp_add(A,B,C)
  471. SPMAT *A, *B, *C;
  472. {
  473.    int i, in_situ;
  474.    SPROW *rc;
  475.    static SPROW *tmp;
  476.  
  477.    if ( ! A || ! B )
  478.      error(E_NULL,"sp_add");
  479.    if ( A->m != B->m || A->n != B->n )
  480.      error(E_SIZES,"sp_add");
  481.    if (C == A || C == B)
  482.      in_situ = TRUE;
  483.    else in_situ = FALSE;
  484.  
  485.    if ( ! C )
  486.      C = sp_get(A->m,A->n,5);
  487.    else {
  488.       if ( C->m != A->m || C->n != A->n  )
  489.     error(E_SIZES,"sp_add");
  490.       if (!in_situ) sp_zero(C);
  491.    }
  492.  
  493.    if (tmp == (SPROW *)NULL && in_situ) {
  494.       tmp = sprow_get(MINROWLEN);
  495.       MEM_STAT_REG(tmp,TYPE_SPROW);
  496.    }
  497.  
  498.    if (in_situ)
  499.      for (i=0; i < A->m; i++) {
  500.     rc = &(C->row[i]);
  501.     sprow_add(&(A->row[i]),&(B->row[i]),0,tmp,TYPE_SPROW);
  502.     sprow_resize(rc,tmp->len,TYPE_SPMAT);
  503.     MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
  504.     rc->len = tmp->len;
  505.      }
  506.    else
  507.      for (i=0; i < A->m; i++) {
  508.     sprow_add(&(A->row[i]),&(B->row[i]),0,&(C->row[i]),TYPE_SPMAT);
  509.      }
  510.  
  511.    C->flag_col = C->flag_diag = FALSE;
  512.  
  513.    return C;
  514. }
  515.  
  516. /*  C = A-B, cannot be in situ */
  517. SPMAT *sp_sub(A,B,C)
  518. SPMAT *A, *B, *C;
  519. {
  520.    int i, in_situ;
  521.    SPROW *rc;
  522.    static SPROW *tmp;
  523.    
  524.    if ( ! A || ! B )
  525.      error(E_NULL,"sp_sub");
  526.    if ( A->m != B->m || A->n != B->n )
  527.      error(E_SIZES,"sp_sub");
  528.    if (C == A || C == B)
  529.      in_situ = TRUE;
  530.    else in_situ = FALSE;
  531.  
  532.    if ( ! C )
  533.      C = sp_get(A->m,A->n,5);
  534.    else {
  535.       if ( C->m != A->m || C->n != A->n  )
  536.     error(E_SIZES,"sp_sub");
  537.       if (!in_situ) sp_zero(C);
  538.    }
  539.  
  540.    if (tmp == (SPROW *)NULL && in_situ) {
  541.       tmp = sprow_get(MINROWLEN);
  542.       MEM_STAT_REG(tmp,TYPE_SPROW);
  543.    }
  544.  
  545.    if (in_situ)
  546.      for (i=0; i < A->m; i++) {
  547.     rc = &(C->row[i]);
  548.     sprow_sub(&(A->row[i]),&(B->row[i]),0,tmp,TYPE_SPROW);
  549.     sprow_resize(rc,tmp->len,TYPE_SPMAT);
  550.     MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
  551.     rc->len = tmp->len;
  552.      }
  553.    else
  554.      for (i=0; i < A->m; i++) {
  555.     sprow_sub(&(A->row[i]),&(B->row[i]),0,&(C->row[i]),TYPE_SPMAT);
  556.      }
  557.  
  558.    C->flag_col = C->flag_diag = FALSE;
  559.    
  560.    return C;
  561. }
  562.  
  563. /*  C = A+alpha*B, cannot be in situ */
  564. SPMAT *sp_mltadd(A,B,alpha,C)
  565. SPMAT *A, *B, *C;
  566. double alpha;
  567. {
  568.    int i, in_situ;
  569.    SPROW *rc;
  570.    static SPROW *tmp;
  571.  
  572.    if ( ! A || ! B )
  573.      error(E_NULL,"sp_mltadd");
  574.    if ( A->m != B->m || A->n != B->n )
  575.      error(E_SIZES,"sp_mltadd");
  576.    if (C == A || C == B)
  577.      in_situ = TRUE;
  578.    else in_situ = FALSE;
  579.  
  580.    if ( ! C )
  581.      C = sp_get(A->m,A->n,5);
  582.    else {
  583.       if ( C->m != A->m || C->n != A->n  )
  584.     error(E_SIZES,"sp_mltadd");
  585.       if (!in_situ) sp_zero(C);
  586.    }
  587.  
  588.    if (tmp == (SPROW *)NULL && in_situ) {
  589.       tmp = sprow_get(MINROWLEN);
  590.       MEM_STAT_REG(tmp,TYPE_SPROW);
  591.    }
  592.  
  593.    if (in_situ)
  594.      for (i=0; i < A->m; i++) {
  595.     rc = &(C->row[i]);
  596.     sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,tmp,TYPE_SPROW);
  597.     sprow_resize(rc,tmp->len,TYPE_SPMAT);
  598.     MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
  599.     rc->len = tmp->len;
  600.      }
  601.    else
  602.      for (i=0; i < A->m; i++) {
  603.     sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,
  604.              &(C->row[i]),TYPE_SPMAT);
  605.      }
  606.    
  607.    C->flag_col = C->flag_diag = FALSE;
  608.    
  609.    return C;
  610. }
  611.  
  612.  
  613.  
  614. /*  B = alpha*A, can be in situ */
  615. SPMAT *sp_smlt(A,alpha,B)
  616. SPMAT *A, *B;
  617. double alpha;
  618. {
  619.    int i;
  620.  
  621.    if ( ! A )
  622.      error(E_NULL,"sp_smlt");
  623.    if ( ! B )
  624.      B = sp_get(A->m,A->n,5);
  625.    else
  626.      if ( A->m != B->m || A->n != B->n )
  627.        error(E_SIZES,"sp_smlt");
  628.  
  629.    for (i=0; i < A->m; i++) {
  630.       sprow_smlt(&(A->row[i]),alpha,0,&(B->row[i]),TYPE_SPMAT);
  631.    }
  632.    return B;
  633. }
  634.  
  635.  
  636.  
  637. /* sp_zero -- zero all the (represented) elements of a sparse matrix */
  638. SPMAT    *sp_zero(A)
  639. SPMAT    *A;
  640. {
  641.    int    i, idx, len;
  642.    row_elt    *elt;
  643.    
  644.    if ( ! A )
  645.      error(E_NULL,"sp_zero");
  646.    
  647.    for ( i = 0; i < A->m; i++ )
  648.    {
  649.       elt = A->row[i].elt;
  650.       len = A->row[i].len;
  651.       for ( idx = 0; idx < len; idx++ )
  652.     (*elt++).val = 0.0;
  653.    }
  654.    
  655.    return A;
  656. }
  657.  
  658. /* sp_copy2 -- copy sparse matrix (type 2) 
  659.    -- keeps structure of the OUT matrix */
  660. SPMAT    *sp_copy2(A,OUT)
  661. SPMAT    *A, *OUT;
  662. {
  663.    int    i /* , idx, len1, len2 */;
  664.    SPROW    *r1, *r2;
  665.    static SPROW    *scratch = (SPROW *)NULL;
  666.    /* row_elt    *e1, *e2; */
  667.    
  668.    if ( ! A )
  669.      error(E_NULL,"sp_copy2");
  670.    if ( ! OUT )
  671.      OUT = sp_get(A->m,A->n,10);
  672.    if ( ! scratch ) {
  673.       scratch = sprow_xpd(scratch,MINROWLEN,TYPE_SPROW);
  674.       MEM_STAT_REG(scratch,TYPE_SPROW);
  675.    }
  676.  
  677.    if ( OUT->m < A->m )
  678.    {
  679.       if (mem_info_is_on()) {
  680.      mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
  681.               A->m*sizeof(SPROW));
  682.       }
  683.  
  684.       OUT->row = RENEW(OUT->row,A->m,SPROW);
  685.       if ( ! OUT->row )
  686.     error(E_MEM,"sp_copy2");
  687.       
  688.       for ( i = OUT->m; i < A->m; i++ )
  689.       {
  690.      OUT->row[i].elt = NEW_A(MINROWLEN,row_elt);
  691.      if ( ! OUT->row[i].elt )
  692.        error(E_MEM,"sp_copy2");
  693.      else if (mem_info_is_on()) {
  694.         mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt));
  695.      }
  696.      
  697.      OUT->row[i].maxlen = MINROWLEN;
  698.      OUT->row[i].len = 0;
  699.       }
  700.       OUT->m = A->m;
  701.    }
  702.    
  703.    OUT->flag_col = OUT->flag_diag = FALSE;
  704.    /* sp_zero(OUT); */
  705.  
  706.    for ( i = 0; i < A->m; i++ )
  707.    {
  708.       r1 = &(A->row[i]);    r2 = &(OUT->row[i]);
  709.       sprow_copy(r1,r2,scratch,TYPE_SPROW);
  710.       if ( r2->maxlen < scratch->len )
  711.     sprow_xpd(r2,scratch->len,TYPE_SPMAT);
  712.       MEM_COPY((char *)(scratch->elt),(char *)(r2->elt),
  713.            scratch->len*sizeof(row_elt));
  714.       r2->len = scratch->len;
  715.       /*******************************************************
  716.     e1 = r1->elt;        e2 = r2->elt;
  717.     len1 = r1->len;        len2 = r2->len;
  718.     for ( idx = 0; idx < len2; idx++, e2++ )
  719.     e2->val = 0.0;
  720.     for ( idx = 0; idx < len1; idx++, e1++ )
  721.     sprow_set_val(r2,e1->col,e1->val);
  722.     *******************************************************/
  723.    }
  724.  
  725.    sp_col_access(OUT);
  726.    return OUT;
  727. }
  728.  
  729. /* sp_resize -- resize a sparse matrix
  730.    -- don't destroying any contents if possible
  731.    -- returns resized matrix */
  732. SPMAT    *sp_resize(A,m,n)
  733. SPMAT    *A;
  734. int    m, n;
  735. {
  736.    int    i, len;
  737.    SPROW    *r;
  738.    
  739.    if (m < 0 || n < 0)
  740.      error(E_NEG,"sp_resize");
  741.  
  742.    if ( ! A )
  743.      return sp_get(m,n,10);
  744.  
  745.    if (m == A->m && n == A->n)
  746.      return A;
  747.  
  748.    if ( m <= A->max_m )
  749.    {
  750.       for ( i = A->m; i < m; i++ )
  751.     A->row[i].len = 0;
  752.       A->m = m;
  753.    }
  754.    else
  755.    {
  756.       if (mem_info_is_on()) {
  757.      mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
  758.              m*sizeof(SPROW));
  759.       }
  760.  
  761.       A->row = RENEW(A->row,(unsigned)m,SPROW);
  762.       if ( ! A->row )
  763.     error(E_MEM,"sp_resize");
  764.       for ( i = A->m; i < m; i++ )
  765.       {
  766.      if ( ! (A->row[i].elt = NEW_A(MINROWLEN,row_elt)) )
  767.        error(E_MEM,"sp_resize");
  768.      else if (mem_info_is_on()) {
  769.         mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt)