home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / sparse < prev    next >
Encoding:
Text File  |  1994-03-08  |  23.2 KB  |  1,035 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 matrix 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: sparse.c,v 1.10 1994/03/08 05:46:07 des Exp $";
  38.  
  39. #define    MINROWLEN    10
  40.  
  41.  
  42.  
  43. /* sp_get_val -- returns the (i,j) entry of the sparse matrix A */
  44. double    sp_get_val(A,i,j)
  45. SPMAT    *A;
  46. int    i, j;
  47. {
  48.    SPROW    *r;
  49.    int    idx;
  50.    
  51.    if ( A == SMNULL )
  52.      error(E_NULL,"sp_get_val");
  53.    if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
  54.      error(E_SIZES,"sp_get_val");
  55.    
  56.    r = A->row+i;
  57.    idx = sprow_idx(r,j);
  58.    if ( idx < 0 )
  59.      return 0.0;
  60.    /* else */
  61.    return r->elt[idx].val;
  62. }
  63.  
  64. /* sp_set_val -- sets the (i,j) entry of the sparse matrix A */
  65. double    sp_set_val(A,i,j,val)
  66. SPMAT    *A;
  67. int    i, j;
  68. double    val;
  69. {
  70.    SPROW    *r;
  71.    int    idx, idx2, new_len;
  72.    
  73.    if ( A == SMNULL )
  74.      error(E_NULL,"sp_set_val");
  75.    if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
  76.      error(E_SIZES,"sp_set_val");
  77.    
  78.    r = A->row+i;
  79.    idx = sprow_idx(r,j);
  80.    /* printf("sp_set_val: idx = %d\n",idx); */
  81.    if ( idx >= 0 )
  82.    {    r->elt[idx].val = val;    return val;    }
  83.    /* else */ if ( idx < -1 )
  84.    {
  85.       /* Note: this destroys the column & diag access paths */
  86.       A->flag_col = A->flag_diag = FALSE;
  87.       /* shift & insert new value */
  88.       idx = -(idx+2);    /* this is the intended insertion index */
  89.       if ( r->len >= r->maxlen )
  90.       {
  91.      r->len = r->maxlen;
  92.      new_len = max(2*r->maxlen+1,5);
  93.      if (mem_info_is_on()) {
  94.         mem_bytes(TYPE_SPMAT,A->row[i].maxlen*sizeof(row_elt),
  95.                 new_len*sizeof(row_elt));
  96.      }
  97.  
  98.      r->elt = RENEW(r->elt,new_len,row_elt);
  99.      if ( ! r->elt )    /* can't allocate */
  100.        error(E_MEM,"sp_set_val");
  101.      r->maxlen = 2*r->maxlen+1;
  102.       }
  103.       for ( idx2 = r->len-1; idx2 >= idx; idx2-- )
  104.     MEM_COPY((char *)(&(r->elt[idx2])),
  105.          (char *)(&(r->elt[idx2+1])),sizeof(row_elt));
  106.       /************************************************************
  107.     if ( idx < r->len )
  108.     MEM_COPY((char *)(&(r->elt[idx])),(char *)(&(r->elt[idx+1])),
  109.     (r->len-idx)*sizeof(row_elt));
  110.     ************************************************************/
  111.       r->len++;
  112.       r->elt[idx].col = j;
  113.       return r->elt[idx].val = val;
  114.    }
  115.    /* else -- idx == -1, error in index/matrix! */
  116.    return 0.0;
  117. }
  118.  
  119. /* sp_mv_mlt -- sparse matrix/dense vector multiply
  120.    -- result is in out, which is returned unless out==NULL on entry
  121.    --  if out==NULL on entry then the result vector is created */
  122. VEC    *sp_mv_mlt(A,x,out)
  123. SPMAT    *A;
  124. VEC    *x, *out;
  125. {
  126.    int    i, j_idx, m, n, max_idx;
  127.    Real    sum, *x_ve;
  128.    SPROW    *r;
  129.    row_elt    *elts;
  130.    
  131.    if ( ! A || ! x )
  132.      error(E_NULL,"sp_mv_mlt");
  133.    if ( x->dim != A->n )
  134.      error(E_SIZES,"sp_mv_mlt");
  135.    if ( ! out || out->dim < A->m )
  136.      out = v_resize(out,A->m);
  137.    if ( out == x )
  138.      error(E_INSITU,"sp_mv_mlt");
  139.    m = A->m;    n = A->n;
  140.    x_ve = x->ve;
  141.    
  142.    for ( i = 0; i < m; i++ )
  143.    {
  144.       sum = 0.0;
  145.       r = &(A->row[i]);
  146.       max_idx = r->len;
  147.       elts    = r->elt;
  148.       for ( j_idx = 0; j_idx < max_idx; j_idx++, elts++ )
  149.     sum += elts->val*x_ve[elts->col];
  150.       out->ve[i] = sum;
  151.    }
  152.    return out;
  153. }
  154.  
  155. /* sp_vm_mlt -- sparse matrix/dense vector multiply from left
  156.    -- result is in out, which is returned unless out==NULL on entry
  157.    -- if out==NULL on entry then result vector is created & returned */
  158. VEC    *sp_vm_mlt(A,x,out)
  159. SPMAT    *A;
  160. VEC    *x, *out;
  161. {
  162.    int    i, j_idx, m, n, max_idx;
  163.    Real    tmp, *x_ve, *out_ve;
  164.    SPROW    *r;
  165.    row_elt    *elts;
  166.    
  167.    if ( ! A || ! x )
  168.      error(E_NULL,"sp_vm_mlt");
  169.    if ( x->dim != A->m )
  170.      error(E_SIZES,"sp_vm_mlt");
  171.    if ( ! out || out->dim < A->n )
  172.      out = v_resize(out,A->n);
  173.    if ( out == x )
  174.      error(E_INSITU,"sp_vm_mlt");
  175.    
  176.    m = A->m;    n = A->n;
  177.    v_zero(out);
  178.    x_ve = x->ve;    out_ve = out->ve;
  179.    
  180.    for ( i = 0; i < m; i++ )
  181.    {
  182.       r = A->row+i;
  183.       max_idx = r->len;
  184.       elts    = r->elt;
  185.       tmp = x_ve[i];
  186.       for ( j_idx = 0; j_idx < max_idx; j_idx++, elts++ )
  187.     out_ve[elts->col] += elts->val*tmp;
  188.    }
  189.    
  190.    return out;
  191. }
  192.  
  193.  
  194. /* sp_get -- get sparse matrix
  195.    -- len is number of elements available for each row without
  196.    allocating further memory */
  197. SPMAT    *sp_get(m,n,maxlen)
  198. int    m, n, maxlen;
  199. {
  200.    SPMAT    *A;
  201.    SPROW    *rows;
  202.    int    i;
  203.    
  204.    if ( m < 0 || n < 0 )
  205.      error(E_NEG,"sp_get");
  206.  
  207.    maxlen = max(maxlen,1);
  208.    
  209.    A = NEW(SPMAT);
  210.    if ( ! A )        /* can't allocate */
  211.      error(E_MEM,"sp_get");
  212.    else if (mem_info_is_on()) {
  213.       mem_bytes(TYPE_SPMAT,0,sizeof(SPMAT));
  214.       mem_numvar(TYPE_SPMAT,1);
  215.    }
  216.    /* fprintf(stderr,"Have SPMAT structure\n"); */
  217.    
  218.    A->row = rows = NEW_A(m,SPROW);
  219.    if ( ! A->row )        /* can't allocate */
  220.      error(E_MEM,"sp_get");
  221.    else if (mem_info_is_on()) {
  222.       mem_bytes(TYPE_SPMAT,0,m*sizeof(SPROW));
  223.    }
  224.    /* fprintf(stderr,"Have row structure array\n"); */
  225.    
  226.    A->start_row = NEW_A(n,int);
  227.    A->start_idx = NEW_A(n,int);
  228.    if ( ! A->start_row || ! A->start_idx )    /* can't allocate */
  229.      error(E_MEM,"sp_get");
  230.    else if (mem_info_is_on()) {
  231.       mem_bytes(TYPE_SPMAT,0,2*n*sizeof(int));
  232.    }
  233.    for ( i = 0; i < n; i++ )
  234.      A->start_row[i] = A->start_idx[i] = -1;
  235.    /* fprintf(stderr,"Have start_row array\n"); */
  236.    
  237.    A->m = A->max_m = m;
  238.    A->n = A->max_n = n;
  239.    
  240.    for ( i = 0; i < m; i++, rows++ )
  241.    {
  242.       rows->elt = NEW_A(maxlen,row_elt);
  243.       if ( ! rows->elt )
  244.     error(E_MEM,"sp_get");
  245.       else if (mem_info_is_on()) {
  246.      mem_bytes(TYPE_SPMAT,0,maxlen*sizeof(row_elt));
  247.       }
  248.       /* fprintf(stderr,"Have row %d element array\n",i); */
  249.       rows->len = 0;
  250.       rows->maxlen = maxlen;
  251.       rows->diag = -1;
  252.    }
  253.    
  254.    return A;
  255. }
  256.  
  257.  
  258. /* sp_free -- frees up the memory for a sparse matrix */
  259. int    sp_free(A)
  260. SPMAT    *A;
  261. {
  262.    SPROW    *r;
  263.    int    i;
  264.    
  265.    if ( ! A )
  266.      return -1;
  267.    if ( A->start_row != (int *)NULL ) {
  268.       if (mem_info_is_on()) {
  269.      mem_bytes(TYPE_SPMAT,A->max_n*sizeof(int),0);
  270.       }
  271.       free((char *)(A->start_row));
  272.    }
  273.    if ( A->start_idx != (int *)NULL ) {
  274.       if (mem_info_is_on()) {
  275.      mem_bytes(TYPE_SPMAT,A->max_n*sizeof(int),0);
  276.       }
  277.       
  278.       free((char *)(A->start_idx));
  279.    }
  280.    if ( ! A->row )
  281.    {
  282.       if (mem_info_is_on()) {
  283.      mem_bytes(TYPE_SPMAT,sizeof(SPMAT),0);
  284.      mem_numvar(TYPE_SPMAT,-1);
  285.       }
  286.       
  287.       free((char *)A);
  288.       return 0;
  289.    }
  290.    for ( i = 0; i < A->m; i++ )
  291.    {
  292.       r = &(A->row[i]);
  293.       if ( r->elt != (row_elt *)NULL ) {
  294.      if (mem_info_is_on()) {
  295.         mem_bytes(TYPE_SPMAT,A->row[i].maxlen*sizeof(row_elt),0);
  296.      }
  297.      free((char *)(r->elt));
  298.       }
  299.    }
  300.    
  301.    if (mem_info_is_on()) {
  302.       if (A->row) 
  303.     mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),0);
  304.       mem_bytes(TYPE_SPMAT,sizeof(SPMAT),0);
  305.       mem_numvar(TYPE_SPMAT,-1);
  306.    }
  307.    
  308.    free((char *)(A->row));
  309.    free((char *)A);
  310.  
  311.    return 0;
  312. }
  313.  
  314.  
  315. /* sp_copy -- constructs a copy of a given matrix
  316.    -- note that the max_len fields (etc) are no larger in the copy
  317.    than necessary
  318.    -- result is returned */
  319. SPMAT    *sp_copy(A)
  320. SPMAT    *A;
  321. {
  322.    SPMAT    *out;
  323.    SPROW    *row1, *row2;
  324.    int    i;
  325.    
  326.    if ( A == SMNULL )
  327.      error(E_NULL,"sp_copy");
  328.    if ( ! (out=NEW(SPMAT)) )
  329.      error(E_MEM,"sp_copy");
  330.    else if (mem_info_is_on()) {
  331.       mem_bytes(TYPE_SPMAT,0,sizeof(SPMAT));
  332.       mem_numvar(TYPE_SPMAT,1);
  333.    }
  334.    out->m = out->max_m = A->m;    out->n = out->max_n = A->n;
  335.    
  336.    /* set up rows */
  337.    if ( ! (out->row=NEW_A(A->m,SPROW)) )
  338.      error(E_MEM,"sp_copy");
  339.    else if (mem_info_is_on()) {
  340.       mem_bytes(TYPE_SPMAT,0,A->m*sizeof(SPROW));
  341.    }
  342.    for ( i = 0; i < A->m; i++ )
  343.    {
  344.       row1 = &(A->row[i]);
  345.       row2 = &(out->row[i]);
  346.       if ( ! (row2->elt=NEW_A(max(row1->len,3),row_elt)) )
  347.     error(E_MEM,"sp_copy");
  348.       else if (mem_info_is_on()) {
  349.      mem_bytes(TYPE_SPMAT,0,max(row1->len,3)*sizeof(row_elt));
  350.       }
  351.       row2->len = row1->len;
  352.       row2->maxlen = max(row1->len,3);
  353.       row2->diag = row1->diag;
  354.       MEM_COPY((char *)(row1->elt),(char *)(row2->elt),
  355.            row1->len*sizeof(row_elt));
  356.    }
  357.    
  358.    /* set up start arrays -- for column access */
  359.    if ( ! (out->start_idx=NEW_A(A->n,int)) ||
  360.        ! (out->start_row=NEW_A(A->n,int)) )
  361.      error(E_MEM,"sp_copy");
  362.    else if (mem_info_is_on()) {
  363.       mem_bytes(TYPE_SPMAT,0,2*A->n*sizeof(int));
  364.    }
  365.    MEM_COPY((char *)(A->start_idx),(char *)(out->start_idx),
  366.         A->n*sizeof(int));
  367.    MEM_COPY((char *)(A->start_row),(char *)(out->start_row),
  368.         A->n*sizeof(int));
  369.    
  370.    return out;
  371. }
  372.  
  373. /* sp_col_access -- set column access path; i.e. nxt_row, nxt_idx fields
  374.    -- returns A */
  375. SPMAT    *sp_col_access(A)
  376. SPMAT    *A;
  377. {
  378.    int    i, j, j_idx, len, m, n;
  379.    SPROW    *row;
  380.    row_elt    *r_elt;
  381.    int    *start_row, *start_idx;
  382.    
  383.    if ( A == SMNULL )
  384.      error(E_NULL,"sp_col_access");
  385.    
  386.    m = A->m;    n = A->n;
  387.    
  388.    /* initialise start_row and start_idx */
  389.    start_row = A->start_row;    start_idx = A->start_idx;
  390.    for ( j = 0; j < n; j++ )
  391.    {    *start_row++ = -1;    *start_idx++ = -1;    }
  392.    
  393.    start_row = A->start_row;    start_idx = A->start_idx;
  394.    
  395.    /* now work UP the rows, setting nxt_row, nxt_idx fields */
  396.    for ( i = m-1; i >= 0; i-- )
  397.    {
  398.       row = &(A->row[i]);
  399.       r_elt = row->elt;
  400.       len   = row->len;
  401.       for ( j_idx = 0; j_idx < len; j_idx++, r_elt++ )
  402.       {
  403.      j = r_elt->col;
  404.      r_elt->nxt_row = start_row[j];
  405.      r_elt->nxt_idx = start_idx[j];
  406.      start_row[j] = i;
  407.      start_idx[j] = j_idx;
  408.       }
  409.    }
  410.    
  411.    A->flag_col = TRUE;
  412.    return A;
  413. }
  414.  
  415. /* sp_diag_access -- set diagonal access path(s) */
  416. SPMAT    *sp_diag_access(A)
  417. SPMAT    *A;
  418. {
  419.    int    i, m;
  420.    SPROW    *row;
  421.    
  422.    if ( A == SMNULL )
  423.      error(E_NULL,"sp_diag_access");
  424.    
  425.    m = A->m;
  426.    
  427.    row = A->row;
  428.    for ( i = 0; i < m; i++, row++ )
  429.      row->diag = sprow_idx(row,i);
  430.    
  431.    A->flag_diag = TRUE;
  432.    
  433.    return A;
  434. }
  435.  
  436. /* sp_m2dense -- convert a sparse matrix to a dense one */
  437. MAT    *sp_m2dense(A,out)
  438. SPMAT    *A;
  439. MAT    *out;
  440. {
  441.    int    i, j_idx;
  442.    SPROW    *row;
  443.    row_elt    *elt;
  444.    
  445.    if ( ! A )
  446.      error(E_NULL,"sp_m2dense");
  447.    if ( ! out || out->m < A->m || out->n < A->n )
  448.      out = m_get(A->m,A->n);
  449.    
  450.    m_zero(out);
  451.    for ( i = 0; i < A->m; i++ )
  452.    {
  453.       row = &(A->row[i]);
  454.       elt = row->elt;
  455.       for ( j_idx = 0; j_idx < row->len; j_idx++, elt++ )
  456.     out->me[i][elt->col] = elt->val;
  457.    }
  458.    
  459.    return out;
  460. }
  461.  
  462.  
  463. /*  C = A+B, can be in situ */
  464. SPMAT *sp_add(A,B,C)
  465. SPMAT *A, *B, *C;
  466. {
  467.    int i, in_situ;
  468.    SPROW *rc;
  469.    static SPROW *tmp;
  470.  
  471.    if ( ! A || ! B )
  472.      error(E_NULL,"sp_add");
  473.    if ( A->m != B->m || A->n != B->n )
  474.      error(E_SIZES,"sp_add");
  475.    if (C == A || C == B)
  476.      in_situ = TRUE;
  477.    else in_situ = FALSE;
  478.  
  479.    if ( ! C )
  480.      C = sp_get(A->m,A->n,5);
  481.    else {
  482.       if ( C->m != A->m || C->n != A->n  )
  483.     error(E_SIZES,"sp_add");
  484.       if (!in_situ) sp_zero(C);
  485.    }
  486.  
  487.    if (tmp == (SPROW *)NULL && in_situ) {
  488.       tmp = sprow_get(MINROWLEN);
  489.       MEM_STAT_REG(tmp,TYPE_SPROW);
  490.    }
  491.  
  492.    if (in_situ)
  493.      for (i=0; i < A->m; i++) {
  494.     rc = &(C->row[i]);
  495.     sprow_add(&(A->row[i]),&(B->row[i]),0,tmp,TYPE_SPROW);
  496.     sprow_resize(rc,tmp->len,TYPE_SPMAT);
  497.     MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
  498.     rc->len = tmp->len;
  499.      }
  500.    else
  501.      for (i=0; i < A->m; i++) {
  502.     sprow_add(&(A->row[i]),&(B->row[i]),0,&(C->row[i]),TYPE_SPMAT);
  503.      }
  504.  
  505.    C->flag_col = C->flag_diag = FALSE;
  506.  
  507.    return C;
  508. }
  509.  
  510. /*  C = A-B, cannot be in situ */
  511. SPMAT *sp_sub(A,B,C)
  512. SPMAT *A, *B, *C;
  513. {
  514.    int i, in_situ;
  515.    SPROW *rc;
  516.    static SPROW *tmp;
  517.    
  518.    if ( ! A || ! B )
  519.      error(E_NULL,"sp_sub");
  520.    if ( A->m != B->m || A->n != B->n )
  521.      error(E_SIZES,"sp_sub");
  522.    if (C == A || C == B)
  523.      in_situ = TRUE;
  524.    else in_situ = FALSE;
  525.  
  526.    if ( ! C )
  527.      C = sp_get(A->m,A->n,5);
  528.    else {
  529.       if ( C->m != A->m || C->n != A->n  )
  530.     error(E_SIZES,"sp_sub");
  531.       if (!in_situ) sp_zero(C);
  532.    }
  533.  
  534.    if (tmp == (SPROW *)NULL && in_situ) {
  535.       tmp = sprow_get(MINROWLEN);
  536.       MEM_STAT_REG(tmp,TYPE_SPROW);
  537.    }
  538.  
  539.    if (in_situ)
  540.      for (i=0; i < A->m; i++) {
  541.     rc = &(C->row[i]);
  542.     sprow_sub(&(A->row[i]),&(B->row[i]),0,tmp,TYPE_SPROW);
  543.     sprow_resize(rc,tmp->len,TYPE_SPMAT);
  544.     MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
  545.     rc->len = tmp->len;
  546.      }
  547.    else
  548.      for (i=0; i < A->m; i++) {
  549.     sprow_sub(&(A->row[i]),&(B->row[i]),0,&(C->row[i]),TYPE_SPMAT);
  550.      }
  551.  
  552.    C->flag_col = C->flag_diag = FALSE;
  553.    
  554.    return C;
  555. }
  556.  
  557. /*  C = A+alpha*B, cannot be in situ */
  558. SPMAT *sp_mltadd(A,B,alpha,C)
  559. SPMAT *A, *B, *C;
  560. double alpha;
  561. {
  562.    int i, in_situ;
  563.    SPROW *rc;
  564.    static SPROW *tmp;
  565.  
  566.    if ( ! A || ! B )
  567.      error(E_NULL,"sp_mltadd");
  568.    if ( A->m != B->m || A->n != B->n )
  569.      error(E_SIZES,"sp_mltadd");
  570.    if (C == A || C == B)
  571.      in_situ = TRUE;
  572.    else in_situ = FALSE;
  573.  
  574.    if ( ! C )
  575.      C = sp_get(A->m,A->n,5);
  576.    else {
  577.       if ( C->m != A->m || C->n != A->n  )
  578.     error(E_SIZES,"sp_mltadd");
  579.       if (!in_situ) sp_zero(C);
  580.    }
  581.  
  582.    if (tmp == (SPROW *)NULL && in_situ) {
  583.       tmp = sprow_get(MINROWLEN);
  584.       MEM_STAT_REG(tmp,TYPE_SPROW);
  585.    }
  586.  
  587.    if (in_situ)
  588.      for (i=0; i < A->m; i++) {
  589.     rc = &(C->row[i]);
  590.     sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,tmp,TYPE_SPROW);
  591.     sprow_resize(rc,tmp->len,TYPE_SPMAT);
  592.     MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
  593.     rc->len = tmp->len;
  594.      }
  595.    else
  596.      for (i=0; i < A->m; i++) {
  597.     sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,
  598.              &(C->row[i]),TYPE_SPMAT);
  599.      }
  600.    
  601.    C->flag_col = C->flag_diag = FALSE;
  602.    
  603.    return C;
  604. }
  605.  
  606.  
  607.  
  608. /*  B = alpha*A, can be in situ */
  609. SPMAT *sp_smlt(A,alpha,B)
  610. SPMAT *A, *B;
  611. double alpha;
  612. {
  613.    int i;
  614.  
  615.    if ( ! A )
  616.      error(E_NULL,"sp_smlt");
  617.    if ( ! B )
  618.      B = sp_get(A->m,A->n,5);
  619.    else
  620.      if ( A->m != B->m || A->n != B->n )
  621.        error(E_SIZES,"sp_smlt");
  622.  
  623.    for (i=0; i < A->m; i++) {
  624.       sprow_smlt(&(A->row[i]),alpha,0,&(B->row[i]),TYPE_SPMAT);
  625.    }
  626.    return B;
  627. }
  628.  
  629.  
  630.  
  631. /* sp_zero -- zero all the (represented) elements of a sparse matrix */
  632. SPMAT    *sp_zero(A)
  633. SPMAT    *A;
  634. {
  635.    int    i, idx, len;
  636.    row_elt    *elt;
  637.    
  638.    if ( ! A )
  639.      error(E_NULL,"sp_zero");
  640.    
  641.    for ( i = 0; i < A->m; i++ )
  642.    {
  643.       elt = A->row[i].elt;
  644.       len = A->row[i].len;
  645.       for ( idx = 0; idx < len; idx++ )
  646.     (*elt++).val = 0.0;
  647.    }
  648.    
  649.    return A;
  650. }
  651.  
  652. /* sp_copy2 -- copy sparse matrix (type 2) 
  653.    -- keeps structure of the OUT matrix */
  654. SPMAT    *sp_copy2(A,OUT)
  655. SPMAT    *A, *OUT;
  656. {
  657.    int    i /* , idx, len1, len2 */;
  658.    SPROW    *r1, *r2;
  659.    static SPROW    *scratch = (SPROW *)NULL;
  660.    /* row_elt    *e1, *e2; */
  661.    
  662.    if ( ! A )
  663.      error(E_NULL,"sp_copy2");
  664.    if ( ! OUT )
  665.      OUT = sp_get(A->m,A->n,10);
  666.    if ( ! scratch ) {
  667.       scratch = sprow_xpd(scratch,MINROWLEN,TYPE_SPROW);
  668.       MEM_STAT_REG(scratch,TYPE_SPROW);
  669.    }
  670.  
  671.    if ( OUT->m < A->m )
  672.    {
  673.       if (mem_info_is_on()) {
  674.      mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
  675.               A->m*sizeof(SPROW));
  676.       }
  677.  
  678.       OUT->row = RENEW(OUT->row,A->m,SPROW);
  679.       if ( ! OUT->row )
  680.     error(E_MEM,"sp_copy2");
  681.       
  682.       for ( i = OUT->m; i < A->m; i++ )
  683.       {
  684.      OUT->row[i].elt = NEW_A(MINROWLEN,row_elt);
  685.      if ( ! OUT->row[i].elt )
  686.        error(E_MEM,"sp_copy2");
  687.      else if (mem_info_is_on()) {
  688.         mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt));
  689.      }
  690.      
  691.      OUT->row[i].maxlen = MINROWLEN;
  692.      OUT->row[i].len = 0;
  693.       }
  694.       OUT->m = A->m;
  695.    }
  696.    
  697.    OUT->flag_col = OUT->flag_diag = FALSE;
  698.    /* sp_zero(OUT); */
  699.  
  700.    for ( i = 0; i < A->m; i++ )
  701.    {
  702.       r1 = &(A->row[i]);    r2 = &(OUT->row[i]);
  703.       sprow_copy(r1,r2,scratch,TYPE_SPROW);
  704.       if ( r2->maxlen < scratch->len )
  705.     sprow_xpd(r2,scratch->len,TYPE_SPMAT);
  706.       MEM_COPY((char *)(scratch->elt),(char *)(r2->elt),
  707.            scratch->len*sizeof(row_elt));
  708.       r2->len = scratch->len;
  709.       /*******************************************************
  710.     e1 = r1->elt;        e2 = r2->elt;
  711.     len1 = r1->len;        len2 = r2->len;
  712.     for ( idx = 0; idx < len2; idx++, e2++ )
  713.     e2->val = 0.0;
  714.     for ( idx = 0; idx < len1; idx++, e1++ )
  715.     sprow_set_val(r2,e1->col,e1->val);
  716.     *******************************************************/
  717.    }
  718.  
  719.    sp_col_access(OUT);
  720.    return OUT;
  721. }
  722.  
  723. /* sp_resize -- resize a sparse matrix
  724.    -- don't destroying any contents if possible
  725.    -- returns resized matrix */
  726. SPMAT    *sp_resize(A,m,n)
  727. SPMAT    *A;
  728. int    m, n;
  729. {
  730.    int    i, len;
  731.    SPROW    *r;
  732.    
  733.    if (m < 0 || n < 0)
  734.      error(E_NEG,"sp_resize");
  735.  
  736.    if ( ! A )
  737.      return sp_get(m,n,10);
  738.  
  739.    if (m == A->m && n == A->n)
  740.      return A;
  741.  
  742.    if ( m <= A->max_m )
  743.    {
  744.       for ( i = A->m; i < m; i++ )
  745.     A->row[i].len = 0;
  746.       A->m = m;
  747.    }
  748.    else
  749.    {
  750.       if (mem_info_is_on()) {
  751.      mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
  752.              m*sizeof(SPROW));
  753.       }
  754.  
  755.       A->row = RENEW(A->row,(unsigned)m,SPROW);
  756.       if ( ! A->row )
  757.     error(E_MEM,"sp_resize");
  758.       for ( i = A->m; i < m; i++ )
  759.       {
  760.      if ( ! (A->row[i].elt = NEW_A(MINROWLEN,row_elt)) )
  761.        error(E_MEM,"sp_resize");
  762.      else if (mem_info_is_on()) {
  763.         mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt));
  764.      }
  765.      A->row[i].len = 0;    A->row[i].maxlen = MINROWLEN;
  766.       }
  767.       A->m = A->max_m = m;
  768.    }
  769.  
  770.    /* update number of rows */
  771.    A->n = n;
  772.  
  773.    /* do we need to increase the size of start_idx[] and start_row[] ? */
  774.    if ( n > A->max_n )
  775.    {    /* only have to update the start_idx & start_row arrays */
  776.       if (mem_info_is_on())
  777.       {
  778.       mem_bytes(TYPE_SPMAT,2*A->max_n*sizeof(int),
  779.             2*n*sizeof(int));
  780.       }
  781.  
  782.       A->start_row = RENEW(A->start_row,(unsigned)n,int);
  783.       A->start_idx = RENEW(A->start_idx,(unsigned)n,int);
  784.       if ( ! A->start_row || ! A->start_idx )
  785.     error(E_MEM,"sp_resize");
  786.       A->max_n = n;    /* ...and update max_n */
  787.  
  788.       return A;
  789.    }
  790.  
  791.    if ( n <= A->n )
  792.        /* make sure that all rows are truncated just before column n */
  793.        for ( i = 0; i < A->m; i++ )
  794.        {
  795.        r = &(A->row[i]);
  796.        len = sprow_idx(r,n);
  797.        if ( len < 0 )
  798.            len = -(len+2);
  799.        if ( len < 0 )
  800.            error(E_MEM,"sp_resize");
  801.        r->len = len;
  802.        }
  803.    
  804.    return A;
  805. }
  806.  
  807.  
  808. /* sp_compact -- removes zeros and near-zeros from a sparse matrix */
  809. SPMAT    *sp_compact(A,tol)
  810. SPMAT    *A;
  811. double    tol;
  812. {
  813.    int    i, idx1, idx2;
  814.    SPROW    *r;
  815.    row_elt    *elt1, *elt2;
  816.    
  817.    if (  ! A )
  818.      error(E_NULL,"sp_compact");
  819.    if ( tol < 0.0 )
  820.      error(E_RANGE,"sp_compact");
  821.    
  822.    A->flag_col = A->flag_diag = FALSE;
  823.    
  824.    for ( i = 0; i < A->m; i++ )
  825.    {
  826.       r = &(A->row[i]);
  827.       elt1 = elt2 = r->elt;
  828.       idx1 = idx2 = 0;
  829.       while ( idx1 < r->len )
  830.       {
  831.      /* printf("# sp_compact: idx1 = %d, idx2 = %d\n",idx1,idx2); */
  832.      if ( fabs(elt1->val) <= tol )
  833.      {    idx1++;    elt1++;    continue;    }
  834.      if ( elt1 != elt2 )
  835.        MEM_COPY(elt1,elt2,sizeof(row_elt));
  836.      idx1++;    elt1++;
  837.      idx2++;    elt2++;
  838.       }
  839.       r->len = idx2;
  840.    }
  841.    
  842.    return A;
  843. }
  844.  
  845. /* varying number of arguments */
  846.  
  847. #ifdef ANSI_C
  848.  
  849. /* To allocate memory to many arguments. 
  850.    The function should be called:
  851.    sp_get_vars(m,n,deg,&x,&y,&z,...,NULL);
  852.    where 
  853.      int m,n,deg;
  854.      SPMAT *x, *y, *z,...;
  855.      The last argument should be NULL ! 
  856.      m x n is the dimension of matrices x,y,z,...
  857.      returned value is equal to the number of allocated variables
  858. */
  859.  
  860. int sp_get_vars(int m,int n,int deg,...) 
  861. {
  862.    va_list ap;
  863.    int i=0;
  864.    SPMAT **par;
  865.    
  866.    va_start(ap, deg);
  867.    while (par = va_arg(ap,SPMAT **)) {   /* NULL ends the list*/
  868.       *par = sp_get(m,n,deg);
  869.       i++;
  870.    } 
  871.  
  872.    va_end(ap);
  873.    return i;
  874. }
  875.  
  876.  
  877. /* To resize memory for many arguments. 
  878.    The function should be called:
  879.    sp_resize_vars(m,n,&x,&y,&z,...,NULL);
  880.    where 
  881.      int m,n;
  882.      SPMAT *x, *y, *z,...;
  883.      The last argument should be NULL ! 
  884.      m X n is the resized dimension of matrices x,y,z,...
  885.      returned value is equal to the number of allocated variables.
  886.      If one of x,y,z,.. arguments is NULL then memory is allocated to this 
  887.      argument. 
  888. */
  889.   
  890. int sp_resize_vars(int m,int n,...) 
  891. {
  892.    va_list ap;
  893.    int i=0;
  894.    SPMAT **par;
  895.    
  896.    va_start(ap, n);
  897.    while (par = va_arg(ap,SPMAT **)) {   /* NULL ends the list*/
  898.       *par = sp_resize(*par,m,n);
  899.       i++;
  900.    } 
  901.  
  902.    va_end(ap);
  903.    return i;
  904. }
  905.  
  906. /* To deallocate memory for many arguments. 
  907.    The function should be called:
  908.    sp_free_vars(&x,&y,&z,...,NULL);
  909.    where 
  910.      SPMAT *x, *y, *z,...;
  911.      The last argument should be NULL ! 
  912.      There must be at least one not NULL argument.
  913.      returned value is equal to the number of allocated variables.
  914.      Returned value of x,y,z,.. is VNULL.
  915. */
  916.  
  917. int sp_free_vars(SPMAT **va,...)
  918. {
  919.    va_list ap;
  920.    int i=1;
  921.    SPMAT **par;
  922.    
  923.    sp_free(*va);
  924.    *va = (SPMAT *) NULL;
  925.    va_start(ap, va);
  926.    while (par = va_arg(ap,SPMAT **)) {   /* NULL ends the list*/
  927.       sp_free(*par); 
  928.       *par = (SPMAT *)NULL;
  929.       i++;
  930.    } 
  931.  
  932.    va_end(ap);
  933.    return i;
  934. }
  935.  
  936.  
  937. #elif VARARGS
  938.  
  939. /* To allocate memory to many arguments. 
  940.    The function should be called:
  941.    sp_get_vars(m,n,deg,&x,&y,&z,...,NULL);
  942.    where 
  943.      int m,n,deg;
  944.      SPMAT *x, *y, *z,...;
  945.      The last argument should be NULL ! 
  946.      m x n is the dimension of matrices x,y,z,...
  947.      returned value is equal to the number of allocated variables
  948. */
  949.  
  950. int sp_get_vars(va_alist) va_dcl
  951. {
  952.    va_list ap;
  953.    int i=0, m, n, deg;
  954.    SPMAT **par;
  955.    
  956.    va_start(ap);
  957.    m = va_arg(ap,int);
  958.    n = va_arg(ap,int);
  959.    deg = va_arg(ap,int);
  960.    while (par = va_arg(ap,SPMAT **)) {   /* NULL ends the list*/
  961.       *par = sp_get(m,n,deg);
  962.       i++;
  963.    } 
  964.  
  965.    va_end(ap);
  966.    return i;
  967. }
  968.  
  969.  
  970. /* To resize memory for many arguments. 
  971.    The function should be called:
  972.    sp_resize_vars(m,n,&x,&y,&z,...,NULL);
  973.    where 
  974.      int m,n;
  975.      SPMAT *x, *y, *z,...;
  976.      The last argument should be NULL ! 
  977.      m X n is the resized dimension of matrices x,y,z,...
  978.      returned value is equal to the number of allocated variables.
  979.      If one of x,y,z,.. arguments is NULL then memory is allocated to this 
  980.      argument. 
  981. */
  982.  
  983. int sp_resize_vars(va_alist) va_dcl
  984. {
  985.    va_list ap;
  986.    int i=0, m, n;
  987.    SPMAT **par;
  988.    
  989.    va_start(ap);
  990.    m = va_arg(ap,int);
  991.    n = va_arg(ap,int);
  992.    while (par = va_arg(ap,SPMAT **)) {   /* NULL ends the list*/
  993.       *par = sp_resize(*par,m,n);
  994.       i++;
  995.    } 
  996.  
  997.    va_end(ap);
  998.    return i;
  999. }
  1000.  
  1001.  
  1002.  
  1003. /* To deallocate memory for many arguments. 
  1004.    The function should be called:
  1005.    sp_free_vars(&x,&y,&z,...,NULL);
  1006.    where 
  1007.      SPMAT *x, *y, *z,...;
  1008.      The last argument should be NULL ! 
  1009.      There must be at least one not NULL argument.
  1010.      returned value is equal to the number of allocated variables.
  1011.      Returned value of x,y,z,.. is VNULL.
  1012. */
  1013.  
  1014. int sp_free_vars(va_alist) va_dcl
  1015. {
  1016.    va_list ap;
  1017.    int i=0;
  1018.    SPMAT **par;
  1019.    
  1020.    va_start(ap);
  1021.    while (par = va_arg(ap,SPMAT **)) {   /* NULL ends the list*/
  1022.       sp_free(*par); 
  1023.       *par = (SPMAT *)NULL;
  1024.       i++;
  1025.    } 
  1026.  
  1027.    va_end(ap);
  1028.    return i;
  1029. }
  1030.  
  1031.  
  1032.  
  1033. #endif
  1034.  
  1035.