home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / spchfctr.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  17KB  |  630 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*
  28.     Sparse Cholesky factorisation code
  29.     To be used with sparse.h, sparse.c etc
  30.  
  31. */
  32.  
  33. static char    rcsid[] = "$Id: spchfctr.c,v 1.4 1994/01/13 05:31:32 des Exp $";
  34.  
  35. #include    <stdio.h>
  36. #include    <math.h>
  37. #include    "matrix.h"
  38. #include    "sparse.h"
  39. #include        "sparse2.h"
  40.  
  41.  
  42. #ifndef MALLOCDECL
  43. #ifndef ANSI_C
  44. extern    char    *calloc(), *realloc();
  45. #endif
  46. #endif
  47.  
  48.  
  49.  
  50. /* sprow_ip -- finds the (partial) inner product of a pair of sparse rows
  51.     -- uses a "merging" approach & assumes column ordered rows
  52.     -- row indices for inner product are all < lim */
  53. double    sprow_ip(row1, row2, lim)
  54. SPROW    *row1, *row2;
  55. int    lim;
  56. {
  57.     int            idx1, idx2, len1, len2, tmp;
  58.     int            sprow_idx();
  59.     register row_elt    *elts1, *elts2;
  60.     register Real        sum;
  61.  
  62.     elts1 = row1->elt;    elts2 = row2->elt;
  63.     len1 = row1->len;    len2 = row2->len;
  64.  
  65.     sum = 0.0;
  66.  
  67.     if ( len1 <= 0 || len2 <= 0 )
  68.         return 0.0;
  69.     if ( elts1->col >= lim || elts2->col >= lim )
  70.         return 0.0;
  71.  
  72.     /* use sprow_idx() to speed up inner product where one row is
  73.         much longer than the other */
  74.     idx1 = idx2 = 0;
  75.     if ( len1 > 2*len2 )
  76.     {
  77.         idx1 = sprow_idx(row1,elts2->col);
  78.         idx1 = (idx1 < 0) ? -(idx1+2) : idx1;
  79.         if ( idx1 < 0 )
  80.             error(E_UNKNOWN,"sprow_ip");
  81.         len1 -= idx1;
  82.     }
  83.     else if ( len2 > 2*len1 )
  84.     {
  85.         idx2 = sprow_idx(row2,elts1->col);
  86.         idx2 = (idx2 < 0) ? -(idx2+2) : idx2;
  87.         if ( idx2 < 0 )
  88.             error(E_UNKNOWN,"sprow_ip");
  89.         len2 -= idx2;
  90.     }
  91.     if ( len1 <= 0 || len2 <= 0 )
  92.         return 0.0;
  93.  
  94.     elts1 = &(elts1[idx1]);        elts2 = &(elts2[idx2]);
  95.  
  96.  
  97.     for ( ; ; )    /* forever do... */
  98.     {
  99.         if ( (tmp=elts1->col-elts2->col) < 0 )
  100.         {
  101.             len1--;        elts1++;
  102.             if ( ! len1 || elts1->col >= lim )
  103.             break;
  104.         }
  105.         else if ( tmp > 0 )
  106.         {
  107.             len2--;        elts2++;
  108.             if ( ! len2 || elts2->col >= lim )
  109.             break;
  110.         }
  111.         else
  112.         {
  113.             sum += elts1->val * elts2->val;
  114.             len1--;        elts1++;
  115.             len2--;        elts2++;
  116.             if ( ! len1 || ! len2 ||
  117.                 elts1->col >= lim || elts2->col >= lim )
  118.             break;
  119.         }
  120.     }
  121.  
  122.     return sum;
  123. }
  124.  
  125. /* sprow_sqr -- returns same as sprow_ip(row, row, lim) */
  126. double    sprow_sqr(row, lim)
  127. SPROW    *row;
  128. int    lim;
  129. {
  130.     register    row_elt    *elts;
  131.     int        idx, len;
  132.     register    Real    sum, tmp;
  133.  
  134.     sum = 0.0;
  135.     elts = row->elt;    len = row->len;
  136.     for ( idx = 0; idx < len; idx++, elts++ )
  137.     {
  138.         if ( elts->col >= lim )
  139.             break;
  140.         tmp = elts->val;
  141.         sum += tmp*tmp;
  142.     }
  143.  
  144.     return sum;
  145. }
  146.  
  147. static    int    *scan_row = (int *)NULL, *scan_idx = (int *)NULL,
  148.             *col_list = (int *)NULL;
  149. static    int    scan_len = 0;
  150.  
  151. /* set_scan -- expand scan_row and scan_idx arrays
  152.     -- return new length */
  153. int    set_scan(new_len)
  154. int    new_len;
  155. {
  156.     if ( new_len <= scan_len )
  157.         return scan_len;
  158.     if ( new_len <= scan_len+5 )
  159.         new_len += 5;
  160.  
  161.     if ( ! scan_row || ! scan_idx || ! col_list )
  162.     {
  163.         scan_row = (int *)calloc(new_len,sizeof(int));
  164.         scan_idx = (int *)calloc(new_len,sizeof(int));
  165.         col_list = (int *)calloc(new_len,sizeof(int));
  166.     }
  167.     else
  168.     {
  169.         scan_row = (int *)realloc((char *)scan_row,new_len*sizeof(int));
  170.         scan_idx = (int *)realloc((char *)scan_idx,new_len*sizeof(int));
  171.         col_list = (int *)realloc((char *)col_list,new_len*sizeof(int));
  172.     }
  173.  
  174.     if ( ! scan_row || ! scan_idx || ! col_list )
  175.         error(E_MEM,"set_scan");
  176.     return new_len;
  177. }
  178.  
  179. /* spCHfactor -- sparse Cholesky factorisation
  180.     -- only the lower triangular part of A (incl. diagonal) is used */
  181. SPMAT    *spCHfactor(A)
  182. SPMAT    *A;
  183. {
  184.     register     int    i;
  185.     int    idx, k, m, minim, n, num_scan, diag_idx, tmp1;
  186.     Real    pivot, tmp2;
  187.     SPROW    *r_piv, *r_op;
  188.     row_elt    *elt_piv, *elt_op, *old_elt;
  189.  
  190.     if ( A == SMNULL )
  191.         error(E_NULL,"spCHfactor");
  192.     if ( A->m != A->n )
  193.         error(E_SQUARE,"spCHfactor");
  194.  
  195.     /* set up access paths if not already done so */
  196.     sp_col_access(A);
  197.     sp_diag_access(A);
  198.  
  199.     /* printf("spCHfactor() -- checkpoint 1\n"); */
  200.     m = A->m;    n = A->n;
  201.     for ( k = 0; k < m; k++ )
  202.     {
  203.         r_piv = &(A->row[k]);
  204.         if ( r_piv->len > scan_len )
  205.             set_scan(r_piv->len);
  206.         elt_piv = r_piv->elt;
  207.         diag_idx = sprow_idx2(r_piv,k,r_piv->diag);
  208.         if ( diag_idx < 0 )
  209.             error(E_POSDEF,"spCHfactor");
  210.         old_elt = &(elt_piv[diag_idx]);
  211.         for ( i = 0; i < r_piv->len; i++ )
  212.         {
  213.             if ( elt_piv[i].col > k )
  214.                 break;
  215.             col_list[i] = elt_piv[i].col;
  216.             scan_row[i] = elt_piv[i].nxt_row;
  217.             scan_idx[i] = elt_piv[i].nxt_idx;
  218.         }
  219.         /* printf("spCHfactor() -- checkpoint 2\n"); */
  220.         num_scan = i;    /* number of actual entries in scan_row etc. */
  221.         /* printf("num_scan = %d\n",num_scan); */
  222.  
  223.         /* set diagonal entry of Cholesky factor */
  224.         tmp2 = elt_piv[diag_idx].val - sprow_sqr(r_piv,k);
  225.         if ( tmp2 <= 0.0 )
  226.             error(E_POSDEF,"spCHfactor");
  227.         elt_piv[diag_idx].val = pivot = sqrt(tmp2);
  228.  
  229.         /* now set the k-th column of the Cholesky factors */
  230.         /* printf("k = %d\n",k); */
  231.         for ( ; ; )    /* forever do... */
  232.         {
  233.             /* printf("spCHfactor() -- checkpoint 3\n"); */
  234.             /* find next row where something (non-trivial) happens
  235.             i.e. find min(scan_row) */
  236.             /* printf("scan_row: "); */
  237.             minim = n;
  238.             for ( i = 0; i < num_scan; i++ )
  239.             {
  240.             tmp1 = scan_row[i];
  241.             /* printf("%d ",tmp1); */
  242.             minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
  243.             }
  244.             /* printf("minim = %d\n",minim); */
  245.             /* printf("col_list: "); */
  246. /**********************************************************************
  247.             for ( i = 0; i < num_scan; i++ )
  248.             printf("%d ",col_list[i]);
  249.             printf("\n");
  250. **********************************************************************/
  251.  
  252.             if ( minim >= n )
  253.             break;    /* nothing more to do for this column */
  254.             r_op = &(A->row[minim]);
  255.             elt_op = r_op->elt;
  256.  
  257.             /* set next entry in column k of Cholesky factors */
  258.             idx = sprow_idx2(r_op,k,scan_idx[num_scan-1]);
  259.             if ( idx < 0 )
  260.             {    /* fill-in */
  261.             sp_set_val(A,minim,k,
  262.                     -sprow_ip(r_piv,r_op,k)/pivot);
  263.             /* in case a realloc() has occurred... */
  264.             elt_op = r_op->elt;
  265.             /* now set up column access path again */
  266.             idx = sprow_idx2(r_op,k,-(idx+2));
  267.             tmp1 = old_elt->nxt_row;
  268.             old_elt->nxt_row = minim;
  269.             r_op->elt[idx].nxt_row = tmp1;
  270.             tmp1 = old_elt->nxt_idx;
  271.             old_elt->nxt_idx = idx;
  272.             r_op->elt[idx].nxt_idx = tmp1;
  273.             }
  274.             else
  275.                 elt_op[idx].val = (elt_op[idx].val -
  276.                 sprow_ip(r_piv,r_op,k))/pivot;
  277.  
  278.             /* printf("spCHfactor() -- checkpoint 4\n"); */
  279.  
  280.             /* remember current element in column k for column chain */
  281.             idx = sprow_idx2(r_op,k,idx);
  282.             old_elt = &(r_op->elt[idx]);
  283.  
  284.             /* update scan_row */
  285.             /* printf("spCHfactor() -- checkpoint 5\n"); */
  286.             /* printf("minim = %d\n",minim); */
  287.             for ( i = 0; i < num_scan; i++ )
  288.             {
  289.             if ( scan_row[i] != minim )
  290.                 continue;
  291.             idx = sprow_idx2(r_op,col_list[i],scan_idx[i]);
  292.             if ( idx < 0 )
  293.             {    scan_row[i] = -1;    continue;    }
  294.             scan_row[i] = elt_op[idx].nxt_row;
  295.             scan_idx[i] = elt_op[idx].nxt_idx;
  296.             /* printf("scan_row[%d] = %d\n",i,scan_row[i]); */
  297.             /* printf("scan_idx[%d] = %d\n",i,scan_idx[i]); */
  298.             }
  299.             
  300.         }
  301.         /* printf("spCHfactor() -- checkpoint 6\n"); */
  302.         /* sp_dump(stdout,A); */
  303.         /* printf("\n\n\n"); */
  304.     }
  305.  
  306.     return A;
  307. }
  308.  
  309. /* spCHsolve -- solve L.L^T.out=b where L is a sparse matrix,
  310.     -- out, b dense vectors
  311.     -- returns out; operation may be in-situ */
  312. VEC    *spCHsolve(L,b,out)
  313. SPMAT    *L;
  314. VEC    *b, *out;
  315. {
  316.     int    i, j_idx, n, scan_idx, scan_row;
  317.     SPROW    *row;
  318.     row_elt    *elt;
  319.     Real    diag_val, sum, *out_ve;
  320.  
  321.     if ( L == SMNULL || b == VNULL )
  322.         error(E_NULL,"spCHsolve");
  323.     if ( L->m != L->n )
  324.         error(E_SQUARE,"spCHsolve");
  325.     if ( b->dim != L->m )
  326.         error(E_SIZES,"spCHsolve");
  327.  
  328.     if ( ! L->flag_col )
  329.         sp_col_access(L);
  330.     if ( ! L->flag_diag )
  331.         sp_diag_access(L);
  332.  
  333.     out = v_copy(b,out);
  334.     out_ve = out->ve;
  335.  
  336.     /* forward substitution: solve L.x=b for x */
  337.     n = L->n;
  338.     for ( i = 0; i < n; i++ )
  339.     {
  340.         sum = out_ve[i];
  341.         row = &(L->row[i]);
  342.         elt = row->elt;
  343.         for ( j_idx = 0; j_idx < row->len; j_idx++, elt++ )
  344.         {
  345.             if ( elt->col >= i )
  346.             break;
  347.             sum -= elt->val*out_ve[elt->col];
  348.         }
  349.         if ( row->diag >= 0 )
  350.             out_ve[i] = sum/(row->elt[row->diag].val);
  351.         else
  352.             error(E_SING,"spCHsolve");
  353.     }
  354.  
  355.     /* backward substitution: solve L^T.out = x for out */
  356.     for ( i = n-1; i >= 0; i-- )
  357.     {
  358.         sum = out_ve[i];
  359.         row = &(L->row[i]);
  360.         /* Note that row->diag >= 0 by above loop */
  361.         elt = &(row->elt[row->diag]);
  362.         diag_val = elt->val;
  363.  
  364.         /* scan down column */
  365.         scan_idx = elt->nxt_idx;
  366.         scan_row = elt->nxt_row;
  367.         while ( scan_row >= 0 /* && scan_idx >= 0 */ )
  368.         {
  369.             row = &(L->row[scan_row]);
  370.             elt = &(row->elt[scan_idx]);
  371.             sum -= elt->val*out_ve[scan_row];
  372.             scan_idx = elt->nxt_idx;
  373.             scan_row = elt->nxt_row;
  374.         }
  375.         out_ve[i] = sum/diag_val;
  376.     }
  377.  
  378.     return out;
  379. }
  380.  
  381. /* spICHfactor -- sparse Incomplete Cholesky factorisation
  382.     -- does a Cholesky factorisation assuming NO FILL-IN
  383.     -- as for spCHfactor(), only the lower triangular part of A is used */
  384. SPMAT    *spICHfactor(A)
  385. SPMAT    *A;
  386. {
  387.     int    k, m, n, nxt_row, nxt_idx, diag_idx;
  388.     Real    pivot, tmp2;
  389.     SPROW    *r_piv, *r_op;
  390.     row_elt    *elt_piv, *elt_op;
  391.  
  392.     if ( A == SMNULL )
  393.         error(E_NULL,"spICHfactor");
  394.     if ( A->m != A->n )
  395.         error(E_SQUARE,"spICHfactor");
  396.  
  397.     /* set up access paths if not already done so */
  398.     if ( ! A->flag_col )
  399.         sp_col_access(A);
  400.     if ( ! A->flag_diag )
  401.         sp_diag_access(A);
  402.  
  403.     m = A->m;    n = A->n;
  404.     for ( k = 0; k < m; k++ )
  405.     {
  406.         r_piv = &(A->row[k]);
  407.  
  408.         diag_idx = r_piv->diag;
  409.         if ( diag_idx < 0 )
  410.             error(E_POSDEF,"spICHfactor");
  411.  
  412.         elt_piv = r_piv->elt;
  413.  
  414.         /* set diagonal entry of Cholesky factor */
  415.         tmp2 = elt_piv[diag_idx].val - sprow_sqr(r_piv,k);
  416.         if ( tmp2 <= 0.0 )
  417.             error(E_POSDEF,"spICHfactor");
  418.         elt_piv[diag_idx].val = pivot = sqrt(tmp2);
  419.  
  420.         /* find next row where something (non-trivial) happens */
  421.         nxt_row = elt_piv[diag_idx].nxt_row;
  422.         nxt_idx = elt_piv[diag_idx].nxt_idx;
  423.  
  424.         /* now set the k-th column of the Cholesky factors */
  425.         while ( nxt_row >= 0 && nxt_idx >= 0 )
  426.         {
  427.             /* nxt_row and nxt_idx give next next row (& index)
  428.             of the entry to be modified */
  429.             r_op = &(A->row[nxt_row]);
  430.             elt_op = r_op->elt;
  431.             elt_op[nxt_idx].val = (elt_op[nxt_idx].val -
  432.                 sprow_ip(r_piv,r_op,k))/pivot;
  433.  
  434.             nxt_row = elt_op[nxt_idx].nxt_row;
  435.             nxt_idx = elt_op[nxt_idx].nxt_idx;
  436.         }
  437.     }
  438.  
  439.     return A;
  440. }
  441.  
  442.  
  443. /* spCHsymb -- symbolic sparse Cholesky factorisation
  444.     -- does NOT do any floating point arithmetic; just sets up the structure
  445.     -- only the lower triangular part of A (incl. diagonal) is used */
  446. SPMAT    *spCHsymb(A)
  447. SPMAT    *A;
  448. {
  449.     register     int    i;
  450.     int    idx, k, m, minim, n, num_scan, diag_idx, tmp1;
  451.     SPROW    *r_piv, *r_op;
  452.     row_elt    *elt_piv, *elt_op, *old_elt;
  453.  
  454.     if ( A == SMNULL )
  455.         error(E_NULL,"spCHsymb");
  456.     if ( A->m != A->n )
  457.         error(E_SQUARE,"spCHsymb");
  458.  
  459.     /* set up access paths if not already done so */
  460.     if ( ! A->flag_col )
  461.         sp_col_access(A);
  462.     if ( ! A->flag_diag )
  463.         sp_diag_access(A);
  464.  
  465.     /* printf("spCHsymb() -- checkpoint 1\n"); */
  466.     m = A->m;    n = A->n;
  467.     for ( k = 0; k < m; k++ )
  468.     {
  469.         r_piv = &(A->row[k]);
  470.         if ( r_piv->len > scan_len )
  471.             set_scan(r_piv->len);
  472.         elt_piv = r_piv->elt;
  473.         diag_idx = sprow_idx2(r_piv,k,r_piv->diag);
  474.         if ( diag_idx < 0 )
  475.             error(E_POSDEF,"spCHsymb");
  476.         old_elt = &(elt_piv[diag_idx]);
  477.         for ( i = 0; i < r_piv->len; i++ )
  478.         {
  479.             if ( elt_piv[i].col > k )
  480.                 break;
  481.             col_list[i] = elt_piv[i].col;
  482.             scan_row[i] = elt_piv[i].nxt_row;
  483.             scan_idx[i] = elt_piv[i].nxt_idx;
  484.         }
  485.         /* printf("spCHsymb() -- checkpoint 2\n"); */
  486.         num_scan = i;    /* number of actual entries in scan_row etc. */
  487.         /* printf("num_scan = %d\n",num_scan); */
  488.  
  489.         /* now set the k-th column of the Cholesky factors */
  490.         /* printf("k = %d\n",k); */
  491.         for ( ; ; )    /* forever do... */
  492.         {
  493.             /* printf("spCHsymb() -- checkpoint 3\n"); */
  494.             /* find next row where something (non-trivial) happens
  495.             i.e. find min(scan_row) */
  496.             minim = n;
  497.             for ( i = 0; i < num_scan; i++ )
  498.             {
  499.             tmp1 = scan_row[i];
  500.             /* printf("%d ",tmp1); */
  501.             minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
  502.             }
  503.  
  504.             if ( minim >= n )
  505.             break;    /* nothing more to do for this column */
  506.             r_op = &(A->row[minim]);
  507.             elt_op = r_op->elt;
  508.  
  509.             /* set next entry in column k of Cholesky factors */
  510.             idx = sprow_idx2(r_op,k,scan_idx[num_scan-1]);
  511.             if ( idx < 0 )
  512.             {    /* fill-in */
  513.             sp_set_val(A,minim,k,0.0);
  514.             /* in case a realloc() has occurred... */
  515.             elt_op = r_op->elt;
  516.             /* now set up column access path again */
  517.             idx = sprow_idx2(r_op,k,-(idx+2));
  518.             tmp1 = old_elt->nxt_row;
  519.             old_elt->nxt_row = minim;
  520.             r_op->elt[idx].nxt_row = tmp1;
  521.             tmp1 = old_elt->nxt_idx;
  522.             old_elt->nxt_idx = idx;
  523.             r_op->elt[idx].nxt_idx = tmp1;
  524.             }
  525.  
  526.             /* printf("spCHsymb() -- checkpoint 4\n"); */
  527.  
  528.             /* remember current element in column k for column chain */
  529.             idx = sprow_idx2(r_op,k,idx);
  530.             old_elt = &(r_op->elt[idx]);
  531.  
  532.             /* update scan_row */
  533.             /* printf("spCHsymb() -- checkpoint 5\n"); */
  534.             /* printf("minim = %d\n",minim); */
  535.             for ( i = 0; i < num_scan; i++ )
  536.             {
  537.             if ( scan_row[i] != minim )
  538.                 continue;
  539.             idx = sprow_idx2(r_op,col_list[i],scan_idx[i]);
  540.             if ( idx < 0 )
  541.             {    scan_row[i] = -1;    continue;    }
  542.             scan_row[i] = elt_op[idx].nxt_row;
  543.             scan_idx[i] = elt_op[idx].nxt_idx;
  544.             /* printf("scan_row[%d] = %d\n",i,scan_row[i]); */
  545.             /* printf("scan_idx[%d] = %d\n",i,scan_idx[i]); */
  546.             }
  547.             
  548.         }
  549.         /* printf("spCHsymb() -- checkpoint 6\n"); */
  550.     }
  551.  
  552.     return A;
  553. }
  554.  
  555. /* comp_AAT -- compute A.A^T where A is a given sparse matrix */
  556. SPMAT    *comp_AAT(A)
  557. SPMAT    *A;
  558. {
  559.     SPMAT    *AAT;
  560.     SPROW    *r, *r2;
  561.     row_elt    *elts, *elts2;
  562.     int    i, idx, idx2, j, m, minim, n, num_scan, tmp1;
  563.     Real    ip;
  564.  
  565.     if ( ! A )
  566.         error(E_NULL,"comp_AAT");
  567.     m = A->m;    n = A->n;
  568.  
  569.     /* set up column access paths */
  570.     if ( ! A->flag_col )
  571.         sp_col_access(A);
  572.  
  573.     AAT = sp_get(m,m,10);
  574.  
  575.     for ( i = 0; i < m; i++ )
  576.     {
  577.         /* initialisation */
  578.         r = &(A->row[i]);
  579.         elts = r->elt;
  580.  
  581.         /* set up scan lists for this row */
  582.         if ( r->len > scan_len )
  583.             set_scan(r->len);
  584.         for ( j = 0; j < r->len; j++ )
  585.         {
  586.             col_list[j] = elts[j].col;
  587.             scan_row[j] = elts[j].nxt_row;
  588.             scan_idx[j] = elts[j].nxt_idx;
  589.         }
  590.         num_scan = r->len;
  591.  
  592.         /* scan down the rows for next non-zero not
  593.             associated with a diagonal entry */
  594.         for ( ; ; )
  595.         {
  596.             minim = m;
  597.             for ( idx = 0; idx < num_scan; idx++ )
  598.             {
  599.             tmp1 = scan_row[idx];
  600.             minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
  601.             }
  602.             if ( minim >= m )
  603.              break;
  604.             r2 = &(A->row[minim]);
  605.             if ( minim > i )
  606.             {
  607.             ip = sprow_ip(r,r2,n);
  608.                 sp_set_val(AAT,minim,i,ip);
  609.                 sp_set_val(AAT,i,minim,ip);
  610.             }
  611.             /* update scan entries */
  612.             elts2 = r2->elt;
  613.             for ( idx = 0; idx < num_scan; idx++ )
  614.             {
  615.             if ( scan_row[idx] != minim || scan_idx[idx] < 0 )
  616.                 continue;
  617.             idx2 = scan_idx[idx];
  618.             scan_row[idx] = elts2[idx2].nxt_row;
  619.             scan_idx[idx] = elts2[idx2].nxt_idx;
  620.             }
  621.         }
  622.  
  623.         /* set the diagonal entry */
  624.         sp_set_val(AAT,i,i,sprow_sqr(r,n));
  625.     }
  626.  
  627.     return AAT;
  628. }
  629.  
  630.