Ifpack_IKLU_Utils.cpp

00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "Ifpack_ConfigDefs.h"
00031 #include "Ifpack_IKLU_Utils.h"
00032 
00033 //-----------------------------------------------------------------
00034 // Allocation and Destruction
00035 //-----------------------------------------------------------------
00036 
00037 /* allocate a sparse matrix (triplet form or compressed-column form) */
00038 csr *csr_spalloc (int m, int n, int nzmax, int values, int triplet)
00039 {
00040   csr *A = (csr*)calloc (1, sizeof (csr)) ;    /* allocate the csr struct */
00041   if (!A) return (NULL) ;                 /* out of memory */
00042   A->m = m ;                              /* define dimensions and nzmax */
00043   A->n = n ;
00044   A->nzmax = nzmax = CS_MAX (nzmax, 1) ;
00045   A->nz = triplet ? 0 : -1 ;              /* allocate triplet or comp.row */
00046   A->p = (int*)malloc (triplet ? CS_MAX(nzmax,1) * sizeof (int) : CS_MAX(m+1,1) * sizeof (int)) ;
00047   A->j = (int*)malloc (CS_MAX(nzmax,1) * sizeof (int)) ;
00048   A->x = values ? (double*)malloc (CS_MAX(nzmax,1) * sizeof (double)) : NULL ;
00049   return ((!A->p || !A->j || (values && !A->x)) ? csr_spfree (A) : A) ;
00050 }
00051 
00052 /* change the max # of entries sparse matrix */
00053 int csr_sprealloc (csr *A, int nzmax)
00054 {
00055     int ok, oki, okj = 1, okx = 1 ;
00056     if (!A) return (0) ;
00057     nzmax = (nzmax <= 0) ? (A->p [A->m]) : nzmax ;
00058     A->j = (int*)csr_realloc (A->j, nzmax, sizeof (int), &oki) ;
00059     if (CS_TRIPLET (A)) A->p = (int*)csr_realloc (A->p, nzmax, sizeof (int), &okj) ;
00060     if (A->x) A->x = (double*)csr_realloc (A->x, nzmax, sizeof (double), &okx) ;
00061     ok = (oki && okj && okx) ;
00062     if (ok) A->nzmax = nzmax ;
00063     return (ok) ;
00064 }
00065 
00066 /* wrapper for realloc */
00067 void *csr_realloc (void *p, int n, size_t size, int *ok)
00068 {
00069     void *pnew ;
00070     pnew = realloc (p, CS_MAX (n,1) * size) ; /* realloc the block */
00071     *ok = (pnew != NULL) ;                  /* realloc fails if pnew is NULL */
00072     return ((*ok) ? pnew : p) ;             /* return original p if failure */
00073 }
00074 
00075 /* free a sparse matrix */
00076 csr *csr_spfree (csr *A)
00077 {
00078   if (!A) return (NULL);     /* do nothing if A already NULL */
00079   if (A->p) free(A->p);
00080   if (A->j) free(A->j);
00081   if (A->x) free(A->x);
00082   if (A) free(A);
00083   return (NULL) ;      /* free the csr struct and return NULL */
00084 }
00085 
00086 /* free a symbolic factorization */
00087 css *csr_sfree (css *S)
00088 {
00089   if (!S) return (NULL) ;     /* do nothing if S already NULL */
00090   if (S->pinv) free(S->pinv);
00091   if (S->q) free(S->q);
00092   if (S->parent) free(S->parent);
00093   if (S->cp) free(S->cp);
00094   if (S->leftmost) free(S->leftmost);
00095   if (S) free(S);
00096   return (NULL) ;      /* free the css struct and return NULL */
00097 }
00098 
00099 csrn *csr_nfree (csrn *N)
00100 {
00101   if (!N) return (NULL) ;     /* do nothing if N already NULL */
00102   csr_spfree (N->L) ;
00103   csr_spfree (N->U) ;
00104   if (N->pinv) free(N->pinv);
00105   if (N->perm) free(N->perm);
00106   if (N->B) free(N->B);
00107   if (N) free(N);
00108   return (NULL) ;      /* free the csn struct and return NULL */
00109 }
00110 
00111 /* free workspace and return a sparse matrix result */
00112 csr *csr_done (csr *C, void *w, void *x, int ok)
00113 {
00114   if (w) free(w);                       /* free workspace */
00115   if (x) free(x);
00116   return (ok ? C : csr_spfree (C)) ;  /* return result if OK, else free it */
00117 }
00118 
00119 /* free workspace and return a numeric factorization (Cholesky, LU, or QR) */
00120 csrn *csr_ndone (csrn *N, csr *C, void *w, void *x, int ok)
00121 {
00122   csr_spfree (C) ;                    /* free temporary matrix */
00123   if (w) free(w);                     /* free workspace */
00124   if (x) free(x);
00125   return (ok ? N : csr_nfree (N)) ;   /* return result if OK, else free it */
00126 }
00127 
00128 /* free workspace and return int array result */
00129 int *csr_idone (int *p, csr *C, void *w, int ok)
00130 {
00131   csr_spfree (C) ;                    /* free temporary matrix */
00132   if (w) free(w);                     /* free workspace */
00133   /* return result if OK, else free it */
00134   if (ok)
00135     return p;
00136   else {
00137     if (p) free(p);
00138     return NULL;
00139   }
00140 }
00141 
00142 //-----------------------------------------------------------------
00143 // Basic CSR math routines
00144 //-----------------------------------------------------------------
00145 
00146 /* p [0..n] = cumulative sum of c [0..n-1], and then copy p [0..n-1] into c */
00147 double csr_cumsum (int *p, int *c, int n)
00148 {
00149   int i, nz = 0 ;
00150   double nz2 = 0 ;
00151   if (!p || !c) return (-1) ;     /* check inputs */
00152   for (i = 0 ; i < n ; i++)
00153     {
00154       p [i] = nz ;
00155       nz += c [i] ;
00156       nz2 += c [i] ;              /* also in double to avoid int overflow */
00157       c [i] = p [i] ;             /* also copy p[0..n-1] back into c[0..n-1]*/
00158     }
00159   p [n] = nz ;
00160   return (nz2) ;                  /* return sum (c [0..n-1]) */
00161 }
00162 
00163 /* x = x + alpha * B(i,:), where x is a dense vector and B(i,:) is sparse */
00164 int csr_scatter (const csr *B, int i, double alpha, int *w, double *x, int mark,
00165          csr *C, int nz)
00166 {
00167   int j, p, *Bp, *Bj, *Cj ;
00168   double *Bx ;
00169   if (!CS_CSC (B) || !w || !CS_CSC (C)) return (-1) ;     /* check inputs */
00170   Bp = B->p ; Bj = B->j ; Bx = B->x ; Cj = C->j ;
00171   for (p = Bp [i] ; p < Bp [i+1] ; p++)
00172     {
00173       j = Bj [p] ;                            /* B(i,j) is nonzero */
00174       if (w [j] < mark)
00175         {
00176       w [j] = mark ;                      /* j is new entry in row i */
00177       Cj [nz++] = j ;                     /* add j to pattern of C(i,:) */
00178       if (x) x [j] = alpha * Bx [p] ;     /* x(j) = alpha*B(i,j) */
00179         }
00180       else if (x) x [j] += alpha * Bx [p] ;   /* j exists in C(i,:) already */
00181     }
00182   return (nz) ;
00183 }
00184 
00185 /* C = alpha*A + beta*B */
00186 csr *csr_add (const csr *A, const csr *B, double alpha, double beta)
00187 {
00188   int p, j, nz = 0, anz, *Cp, *Cj, *Bp, m, n, bnz, *w, values ;
00189   double *x, *Bx, *Cx ;
00190   csr *C ;
00191   if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ;         /* check inputs */
00192   if ( (A->m != B->m) || (A->n != B->n) ) return (NULL);
00193   m = A->m ; anz = A->p [m] ;
00194   n = B->n ; Bp = B->p ; Bx = B->x ; bnz = Bp [m] ;
00195   w = (int*)calloc (CS_MAX(n,1), sizeof (int)) ;                       /* get workspace */
00196   values = (A->x != NULL) && (Bx != NULL) ;
00197   x = values ? (double*)malloc (n * sizeof (double)) : NULL ;    /* get workspace */
00198   C = csr_spalloc (m, n, anz + bnz, values, 0) ;          /* allocate result*/
00199   if (!C || !w || (values && !x)) return (csr_done (C, w, x, 0)) ;
00200   Cp = C->p ; Cj = C->j ; Cx = C->x ;
00201   for (j = 0 ; j < n ; j++)
00202     {
00203       Cp [j] = nz ;                   /* row j of C starts here */
00204       nz = csr_scatter (A, j, alpha, w, x, j+1, C, nz) ;   /* alpha*A(j,:)*/
00205       nz = csr_scatter (B, j, beta, w, x, j+1, C, nz) ;    /* beta*B(j,:) */
00206       if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Cj [p]] ;
00207     }
00208   Cp [m] = nz ;                       /* finalize the last row of C */
00209   csr_sprealloc (C, 0) ;              /* remove extra space from C */
00210   return (csr_done (C, w, x, 1)) ;    /* success; free workspace, return C */
00211 }
00212 
00213 /* C = A' */
00214 csr *csr_transpose (const csr *A, int values)
00215 {
00216   int p, q, i, *Cp, *Cj, n, m, *Ap, *Aj, *w ;
00217   double *Cx, *Ax ;
00218   csr *C ;
00219   if (!CS_CSC (A)) return (NULL) ;    /* check inputs */
00220   m = A->m ; n = A->n ; Ap = A->p ; Aj = A->j ; Ax = A->x ;
00221   C = csr_spalloc (n, m, Ap [m], values && Ax, 0) ;      /* allocate result */
00222   w = (int*)calloc (CS_MAX(n,1), sizeof (int)) ;         /* get workspace */
00223   if (!C || !w) return (csr_done (C, w, NULL, 0)) ;      /* out of memory */
00224   Cp = C->p ; Cj = C->j ; Cx = C->x ;
00225   for (p = 0 ; p < Ap [m] ; p++) w [Aj [p]]++ ;          /* col counts */
00226   csr_cumsum (Cp, w, n) ;                                 /* col pointers */
00227   for (i = 0 ; i < m ; i++)
00228     {
00229       for (p = Ap [i] ; p < Ap [i+1] ; p++)
00230         {
00231       Cj [q = w [Aj [p]]++] = i ; /* place A(i,j) as entry C(j,i) */
00232       if (Cx) Cx [q] = Ax [p] ;
00233         }
00234     }
00235   return (csr_done (C, w, NULL, 1)) ; /* success; free w and return C */
00236 }
00237 
00238 /* C = A*B */
00239 csr *csr_multiply (const csr *A, const csr *B)
00240 {
00241   int p, j, nz = 0, anz, *Cp, *Cj, *Ap, m, k, n, bnz, *w, values, *Aj;
00242   double *x, *Ax, *Cx ;
00243   csr *C ;
00244   if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ;      /* check inputs */
00245   k = A->n;
00246   if (k != B->m ) return(NULL);
00247   m = A->m ; anz = A->p [A->m] ;
00248   n = B->n ; Ap = A->p ; Aj = A->j ; Ax = A->x ; bnz = B->p [k] ;
00249   w = (int*)calloc (CS_MAX(n,1), sizeof (int)) ;            /* get workspace */
00250   values = (Ax != NULL) && (B->x != NULL) ;
00251   x = values ? (double*)malloc (n * sizeof (double)) : NULL ; /* get workspace */
00252   C = csr_spalloc (m, n, anz + bnz, values, 0) ;       /* allocate result */
00253   if (!C || !w || (values && !x)) return (csr_done (C, w, x, 0)) ;
00254   Cp = C->p ;
00255   for (j = 0 ; j < m ; j++)
00256     {
00257       if (nz + n > C->nzmax && !csr_sprealloc (C, 2*(C->nzmax)+m))
00258         {
00259       return (csr_done (C, w, x, 0)) ;            /* out of memory */
00260         }
00261       Cj = C->j ; Cx = C->x ;         /* C->j and C->x may be reallocated */
00262       Cp [j] = nz ;                   /* row j of C starts here */
00263       for (p = Ap [j] ; p < Ap [j+1] ; p++)
00264         {
00265       nz = csr_scatter (B, Aj [p], Ax ? Ax [p] : 1, w, x, j+1, C, nz) ;
00266         }
00267       if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Cj [p]] ;
00268     }
00269   Cp [m] = nz ;                       /* finalize the last column of C */
00270   csr_sprealloc (C, 0) ;              /* remove extra space from C */
00271   return (csr_done (C, w, x, 1)) ;    /* success; free workspace, return C */
00272 }
00273 
00274 
00275 //-----------------------------------------------------------------
00276 // Symbolic analysis
00277 //-----------------------------------------------------------------
00278 
00279 /* symbolic ordering and analysis for LU */
00280 css *csr_sqr (int order, const csr *A )
00281 {
00282     int n, ok = 1;
00283     css *S ;
00284     if (!CS_CSC (A)) return (NULL) ;        /* check inputs */
00285     n = A->n ;
00286     S = (css*)calloc(1, sizeof (css)) ;       /* allocate result S */
00287     if (!S) return (NULL) ;                 /* out of memory */
00288     S->q = csr_amd (order, A) ;             /* fill-reducing ordering */
00289     if (!S->q)
00290     {
00291         printf(" csr_sqr error no permutation\n");
00292     }
00293     if (order && !S->q) return (csr_sfree (S)) ;
00294 
00295     /* LU factorization */
00296     S->unz = (double) CS_MIN(4*(A->p [n]) + n, n * n );
00297     S->lnz = S->unz ;                       /* guess nnz(L) and nnz(U) */
00298     return (ok ? S : csr_sfree (S)) ;        /* return result S */
00299 }
00300 
00301 
00302 /* xi [top...n-1] = nodes reachable from graph of G*P' via nodes in B(:,k).
00303  * xi [n...2n-1] used as workspace */
00304 int csr_reach (csr *G, const csr *B, int k, int *xi, const int *pinv)
00305 {
00306     int p, m, top, *Bp, *Bj, *Gp ;
00307     if (!CS_CSC (G) || !CS_CSC (B) || !xi) return (-1) ;    /* check inputs */
00308     m = G->m ; Bp = B->p ; Bj = B->j ; Gp = G->p ;
00309     top = m ;
00310     for (p = Bp [k] ; p < Bp [k+1] ; p++)
00311     {
00312         if (!CS_MARKED (Gp, Bj [p]))    /* start a dfs at unmarked node i */
00313         {
00314             top = csr_dfs (Bj [p], G, top, xi, xi+m, pinv) ;
00315         }
00316     }
00317     for (p = top ; p < m ; p++) CS_MARK (Gp, xi [p]) ;  /* restore G */
00318     return (top) ;
00319 }
00320 
00321 /* depth-first-search of the graph of a csr matrix, starting at node j */
00322 int csr_dfs (int j, csr *G, int top, int *xi, int *pstack, const int *pinv)
00323 {
00324     int i, p, p2, done, jnew, head = 0, *Gp, *Gj ;
00325     if (!CS_CSC (G) || !xi || !pstack) return (-1) ;    /* check inputs */
00326     Gp = G->p ; Gj = G->j ;
00327     xi [0] = j ;                /* initialize the recursion stack */
00328     while (head >= 0)
00329     {
00330         j = xi [head] ;         /* get j from the top of the recursion stack */
00331         jnew = pinv ? (pinv [j]) : j ;
00332         if (!CS_MARKED (Gp, j))
00333         {
00334             CS_MARK (Gp, j) ;       /* mark node j as visited */
00335             pstack [head] = (jnew < 0) ? 0 : CS_UNFLIP (Gp [jnew]) ;
00336         }
00337         done = 1 ;                  /* node j done if no unvisited neighbors */
00338         p2 = (jnew < 0) ? 0 : CS_UNFLIP (Gp [jnew+1]) ;
00339         for (p = pstack [head] ; p < p2 ; p++)  /* examine all neighbors of j */
00340         {
00341             i = Gj [p] ;            /* consider neighbor node i */
00342             if (CS_MARKED (Gp, i)) continue ;   /* skip visited node i */
00343             pstack [head] = p ;     /* pause depth-first search of node j */
00344             xi [++head] = i ;       /* start dfs at node i */
00345             done = 0 ;              /* node j is not done */
00346             break ;                 /* break, to start dfs (i) */
00347         }
00348         if (done)               /* depth-first search at node j is done */
00349         {
00350             head-- ;            /* remove j from the recursion stack */
00351             xi [--top] = j ;    /* and place in the output stack */
00352         }
00353     }
00354     return (top) ;
00355 }
00356 
00357 /* depth-first search and postorder of a tree rooted at node j */
00358 int csr_tdfs (int j, int k, int *head, const int *next, int *post, int *stack)
00359 {
00360   int i, p, top = 0 ;
00361   if (!head || !next || !post || !stack) return (-1) ;    /* check inputs */
00362   stack [0] = j ;                 /* place j on the stack */
00363   while (top >= 0)                /* while (stack is not empty) */
00364     {
00365       p = stack [top] ;           /* p = top of stack */
00366       i = head [p] ;              /* i = youngest child of p */
00367       if (i == -1)
00368         {
00369       top-- ;                 /* p has no unordered children left */
00370       post [k++] = p ;        /* node p is the kth postordered node */
00371         }
00372       else
00373         {
00374       head [p] = next [i] ;   /* remove i from children of p */
00375       stack [++top] = i ;     /* start dfs on child node i */
00376         }
00377     }
00378   return (k) ;
00379 }
00380 
00381 //-----------------------------------------------------------------
00382 // LU factorization
00383 //-----------------------------------------------------------------
00384 
00385 /*
00386  * Given sparse A,  
00387  * [L,U,pinv]=lu(A, [q lnz unz]). lnz and unz can be guesses 
00388  * Hypotheses: m=n, Lj[Lp[i+1]-1]=i and Uj[Lp[i]]=i
00389  */
00390 csrn *csr_lu (const csr *A, const css *S, double tol)
00391 {
00392     csr *L, *U ;
00393     csrn *N ;
00394     double pivot, *Lx, *Ux, *x,  a, t;
00395     int *Lp, *Lj, *Up, *Uj, *pinv, *xi, *q, n, ipiv, k, top, p, i, row, lnz, unz;
00396     int debug = 0;
00397 
00398     if (!CS_CSC (A) ) printf(" error csrlu: A not csc\n");
00399     if (!CS_CSC (A) || !S) return (NULL) ;      /* check inputs */
00400     n = A->n ;
00401     if (n != A->m) return (NULL) ;      /* check inputs */
00402     q = S->q ; lnz = (int)S->lnz ; unz = (int)S->unz ;
00403     x = (double*)malloc(n * sizeof (double)) ;      /* get double workspace */
00404     xi = (int*)malloc (2*n * sizeof (int)) ;        /* get int workspace */
00405     N = (csrn*)calloc (1, sizeof (csrn)) ;          /* allocate result */
00406     if (!(A) ) printf(" error csrlu: allocation of N failed\n");
00407     if (!x || !xi || !N) return (csr_ndone (N, NULL, xi, x, 0)) ;
00408 
00409     N->L = L = csr_spalloc (n, n, lnz, 1, 0) ;      /* allocate result L */
00410     N->U = U = csr_spalloc (n, n, unz, 1, 0) ;      /* allocate result U */
00411     N->pinv = pinv = (int*)malloc (n * sizeof (int)) ;  /* allocate result pinv */
00412     N->perm = (int*)malloc (n * sizeof (int)) ;         /* allocate result perm */
00413     if (!L || !U || !pinv) return (csr_ndone (N, NULL, xi, x, 0)) ;
00414     Lp = L->p ; Up = U->p ;
00415     for (i = 0 ; i < n ; i++) x [i] = 0 ;       /* clear workspace */
00416     for (i = 0 ; i < n ; i++) pinv [i] = -1 ;       /* no rows pivotal yet */
00417     for (k = 1 ; k <= n ; k++) Up [k] = 0 ;     /* no rows of U yet */
00418     for (k = 1 ; k <= n ; k++) Lp [k] = 0 ;     /* no rows of L yet either */
00419     lnz = unz = 0 ;
00420     if( debug )
00421     {
00422        printf ("A:\n") ; csr_print (A, 0) ;  
00423     }
00424     for (k = 0 ; k < n ; k++)       /* compute L(:,k) and U(:,k) */
00425     {
00426     /* --- Triangular solve --------------------------------------------- */
00427     Lp [k] = lnz ;          /* L(:,k) starts here */
00428     Up [k] = unz ;          /* U(:,k) starts here */
00429     if ((lnz + n > L->nzmax && !csr_sprealloc (L, 2*L->nzmax + n)) ||
00430         (unz + n > U->nzmax && !csr_sprealloc (U, 2*U->nzmax + n)))
00431     {
00432         return (csr_ndone (N, NULL, xi, x, 0)) ;
00433     }
00434         Lj = L->j ; Lx = L->x ; Uj = U->j ; Ux = U->x ;
00435         row = q ? (q [k]) : k ;
00436         if( debug > 1 )
00437         {
00438             printf("--------------------------------\n");
00439             printf(" %d spsolve row=%d \n",k, row);
00440             printf(" pinv = %d %d %d %d \n", pinv[0], pinv[1], pinv[2], pinv[3]);
00441         }
00442         top = csr_spsolve (U, A, row, xi, x, pinv, 1) ; /* x = U\A(row,:) */
00443         if( debug > 1 ) printf("top=%d  x = %g %g %g %g \n", top,x[0],x[1],x[2],x[3]);
00444         /* --- Find pivot --------------------------------------------------- */
00445         ipiv = -1 ;
00446         a = -1 ;
00447         for (p = top ; p < n ; p++)
00448         {
00449             i = xi [p] ;            /* x(i) is nonzero */
00450             if (pinv [i] < 0)       /* row i is not yet pivotal */
00451             {
00452                 if ((t = fabs (x [i])) > a)
00453                 {
00454                     a = t ;         /* largest pivot candidate so far */
00455                     ipiv = i ;
00456                 }
00457             }
00458             else                    /* x(i) is the entry L(pinv[i],k) */
00459             {
00460                 Lj [lnz] = pinv [i] ;
00461                 Lx [lnz++] = x [i] ;
00462             }
00463         }
00464         if (ipiv == -1 || a <= 0) return (csr_ndone (N, NULL, xi, x, 0)) ;
00465         if (pinv [row] < 0 && fabs (x [row]) >= a*tol) ipiv = row;
00466         pivot = x [ipiv] ;          /* the chosen pivot */
00467         Lj [lnz] = k ;              /* last entry in L(:,k) is L(k,k) */
00468         Lx [lnz++] = pivot ;
00469         if( debug > 1 ) { printf ("L:") ; csr_print (L, 0) ;  }
00470 
00471         /* --- Divide by pivot ---------------------------------------------- */
00472         pinv [ipiv] = k ;           /* ipiv is the kth pivot row */
00473         Uj [unz] = ipiv ;           /* first entry in U(:,k) is U(k,k) = 1 */
00474         Ux [unz++] = 1 ;
00475         for (p = top ; p < n ; p++) /* U(k+1:n,k) = x / pivot */
00476         {
00477             i = xi [p] ;
00478             if (pinv [i] < 0)       /* x(i) is an entry in U(:,k) */
00479             {
00480                 Uj [unz] = i ;      /* save unpermuted row in U */
00481                 Ux [unz++] = x [i] / pivot ;    /* scale pivot row */
00482             }
00483             x [i] = 0 ;             /* x [0..n-1] = 0 for next k */
00484         }
00485         if( debug > 1 )
00486         {
00487             printf ("U:") ; csr_print (U, 0) ;  
00488             printf("------------------------------------\n");
00489         }
00490     }
00491     /* --- Finalize U and L ------------------------------------------------- */
00492     Lp [n] = lnz ;
00493     if( debug ) { printf ("L:") ; csr_print (L, 0) ;  }
00494     Up [n] = unz ;
00495     Uj = U->j ;                     /* fix column indices of U for final pinv */
00496     for (p = 0 ; p < unz ; p++) Uj [p] = pinv [Uj [p]] ;
00497 
00498     csr_sprealloc (L, 0) ;      /* remove extra space from L and U */
00499     csr_sprealloc (U, 0) ;
00500     if( debug ) { printf ("U:") ; csr_print (U, 0) ;  }
00501     return (csr_ndone (N, NULL, xi, x, 1)) ;    /* success */
00502 }
00503 
00504 //-----------------------------------------------------------------
00505 // Triangular solves
00506 //-----------------------------------------------------------------
00507 
00508 /* solve xG=b(k,:), where G is either upper (up=1) or lower (up=0) triangular */
00509 int csr_spsolve (csr *G, const csr *B, int k, int *xi, double *x, const int *pinv, int up)
00510 {
00511   int i, I, p, q, px, top, n, *Gp, *Gj, *Bp, *Bj ;
00512   int debug = 0;
00513   double *Gx, *Bx ;
00514   if (!CS_CSC (G) || !CS_CSC (B) || !xi || !x) return (-1) ;
00515   Gp = G->p ; Gj = G->j ; Gx = G->x ; n = G->n ;
00516   Bp = B->p ; Bj = B->j ; Bx = B->x ;
00517   top = csr_reach (G, B, k, xi, pinv) ;       /* xi[top..n-1]=Reach(B(:,k)) */
00518   for (p = top ; p < n ; p++) x [xi [p]] = 0 ;    /* clear x */
00519   for (p = Bp [k] ; p < Bp [k+1] ; p++) x [Bj [p]] = Bx [p] ; /* scatter B */
00520   if( debug )
00521     printf("solve k=%d   x= %g %g %g %g \n", k, x[0], x[1], x[2], x[3]);
00522   if( debug )
00523     printf("top=%d   xi= %d %d %d %d \n", top , xi[0], xi[1], xi[2], xi[3]);
00524   for (px = top ; px < n ; px++)
00525     {
00526       i = xi [px] ;                               /* x(i) is nonzero */
00527       I = pinv ? (pinv [i]) : i ;                 /* i maps to col I of G */
00528       if (I < 0) continue ;                       /* row I is empty */
00529       /* dmd */
00530       x [i] /= Gx [up ? (Gp [I]) : (Gp [I+1]-1)] ;/* x(i) /= G(i,i) */
00531       p = up ? (Gp [I]+1) : (Gp [I]) ;            /* up: L(i,i) 1st entry */
00532       q = up ? (Gp [I+1]) : (Gp [I+1]-1) ;        /* up: U(i,i) last entry */
00533       for ( ; p < q ; p++)
00534         {
00535       if( debug )
00536         printf("%d %d solve %d %g %g \n", px, i ,p,  Gx [p] , x [Gj [p]]  );
00537       
00538       x [Gj[p]] -= Gx [p] * x [i] ;           /* x(i) -= G(i,j) * x(j) */
00539         }
00540       if( debug )
00541     printf(" x= %g %g %g %g \n", x[0], x[1], x[2], x[3]);
00542     }
00543   return (top) ;                                  /* return top of stack */
00544 }
00545 
00546 //-----------------------------------------------------------------
00547 // AMD routine
00548 //-----------------------------------------------------------------
00549 
00550 /*
00551  * amd for csr matrices , and  order =1
00552  * Hypothesis:  m = n
00553  */
00554 /* clear w */
00555 static int csr_wclear (int mark, int lemax, int *w, int n)
00556 {
00557     int k ;
00558     if (mark < 2 || (mark + lemax < 0))
00559     {
00560         for (k = 0 ; k < n ; k++) if (w [k] != 0) w [k] = 1 ;
00561         mark = 2 ;
00562     }
00563     return (mark) ;     /* at this point, w [0..n-1] < mark holds */
00564 }
00565 
00566 /* keep off-diagonal entries; drop diagonal entries */
00567 static int csr_diag (int i, int j, double aij, void *other) { return (i != j) ; }
00568 
00569 /* p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */
00570 int *csr_amd (int order, const csr *A)  /* order 0:natural, 1:Chol, 2:LU, 3:QR */
00571 {
00572   csr *C, *A2, *AT ;
00573   int *Cp, *Cj, *last, *W, *len, *nv, *next, *P, *head, *elen, *degree, *w,
00574     *hhead, *ATp, *ATj, d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
00575     k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
00576     ok, cnz, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, n, m, t ;
00577   unsigned int h ;
00578   /* --- Construct matrix C ----------------------------------------------- */
00579   if (!CS_CSC (A) || order <= 0 || order > 3) return (NULL) ; /* check */
00580   AT = csr_transpose (A, 0) ;             /* compute A' */
00581   if (!AT) return (NULL) ;
00582   m = A->m ; n = A->n ;
00583   if ( n != m) return(NULL); /* For rectangular matrices, use csr_amd */
00584   dense = (int)CS_MAX (16, 10 * sqrt ((double) n)) ;   /* find dense threshold */
00585   dense = CS_MIN (n-2, dense) ;
00586   if (order == 1 && n == m)
00587     {
00588       C = csr_add (A, AT, 0, 0) ;         /* C = A+A' */
00589     }
00590   else if (order == 2)
00591     {
00592       ATp = AT->p ;                       /* drop dense columns from AT */
00593       ATj = AT->j ;
00594       for (p2 = 0, j = 0 ; j < m ; j++)
00595         {
00596       p = ATp [j] ;                   /* column j of AT starts here */
00597       ATp [j] = p2 ;                  /* new column j starts here */
00598       if (ATp [j+1] - p > dense) continue ;   /* skip dense col j */
00599       for ( ; p < ATp [j+1] ; p++) ATj [p2++] = ATj [p] ;
00600         }
00601       ATp [m] = p2 ;                      /* finalize AT */
00602       A2 = csr_transpose (AT, 0) ;        /* A2 = AT' */
00603       C = A2 ? csr_multiply (AT, A2) : NULL ; /* C=A'*A with no dense rows */
00604       csr_spfree (A2) ;
00605     }
00606   else
00607     {
00608       C = csr_multiply (AT, A) ;          /* C=A'*A */
00609     }
00610   csr_spfree (AT) ;
00611   if (!C) return (NULL) ;
00612   csr_fkeep (C, &csr_diag, NULL) ;         /* drop diagonal entries */
00613   Cp = C->p ;
00614   cnz = Cp [n] ;
00615   P = (int*)malloc (CS_MAX(n+1,1) * sizeof (int)) ;     /* allocate result */
00616   W = (int*)malloc (CS_MAX(8*(n+1),1) * sizeof (int)) ; /* get workspace */
00617   t = cnz + cnz/5 + 2*n ;                 /* add elbow room to C */
00618   if (!P || !W || !csr_sprealloc (C, t)) return (csr_idone (P, C, W, 0)) ;
00619   len  = W           ; nv     = W +   (n+1) ; next   = W + 2*(n+1) ;
00620   head = W + 3*(n+1) ; elen   = W + 4*(n+1) ; degree = W + 5*(n+1) ;
00621   w    = W + 6*(n+1) ; hhead  = W + 7*(n+1) ;
00622   last = P ;                              /* use P as workspace for last */
00623   /* --- Initialize quotient graph ---------------------------------------- */
00624   for (k = 0 ; k < n ; k++) len [k] = Cp [k+1] - Cp [k] ;
00625   len [n] = 0 ;
00626   nzmax = C->nzmax ;
00627   Cj = C->j ;
00628   for (i = 0 ; i <= n ; i++)
00629     {
00630       head [i] = -1 ;                     /* degree list i is empty */
00631       last [i] = -1 ;
00632       next [i] = -1 ;
00633       hhead [i] = -1 ;                    /* hash list i is empty */
00634       nv [i] = 1 ;                        /* node i is just one node */
00635       w [i] = 1 ;                         /* node i is alive */
00636       elen [i] = 0 ;                      /* Ek of node i is empty */
00637       degree [i] = len [i] ;              /* degree of node i */
00638     }
00639   mark = csr_wclear (0, 0, w, n) ;         /* clear w */
00640   elen [n] = -2 ;                         /* n is a dead element */
00641   Cp [n] = -1 ;                           /* n is a root of assembly tree */
00642   w [n] = 0 ;                             /* n is a dead element */
00643   /* --- Initialize degree lists ------------------------------------------ */
00644   for (i = 0 ; i < n ; i++)
00645     {
00646       d = degree [i] ;
00647       if (d == 0)                         /* node i is empty */
00648         {
00649       elen [i] = -2 ;                 /* element i is dead */
00650       nel++ ;
00651       Cp [i] = -1 ;                   /* i is a root of assemby tree */
00652       w [i] = 0 ;
00653         }
00654       else if (d > dense)                 /* node i is dense */
00655         {
00656       nv [i] = 0 ;                    /* absorb i into element n */
00657       elen [i] = -1 ;                 /* node i is dead */
00658       nel++ ;
00659       Cp [i] = CS_FLIP (n) ;
00660       nv [n]++ ;
00661         }
00662       else
00663         {
00664       if (head [d] != -1) last [head [d]] = i ;
00665       next [i] = head [d] ;           /* put node i in degree list d */
00666       head [d] = i ;
00667         }
00668     }
00669   while (nel < n)                         /* while (selecting pivots) do */
00670     {
00671       /* --- Select node of minimum approximate degree -------------------- */
00672       for (k = -1 ; mindeg < n && (k = head [mindeg]) == -1 ; mindeg++) ;
00673       if (next [k] != -1) last [next [k]] = -1 ;
00674       head [mindeg] = next [k] ;          /* remove k from degree list */
00675       elenk = elen [k] ;                  /* elenk = |Ek| */
00676       nvk = nv [k] ;                      /* # of nodes k represents */
00677       nel += nvk ;                        /* nv[k] nodes of A eliminated */
00678       /* --- Garbage collection ------------------------------------------- */
00679       if (elenk > 0 && cnz + mindeg >= nzmax)
00680         {
00681       for (j = 0 ; j < n ; j++)
00682             {
00683           if ((p = Cp [j]) >= 0)      /* j is a live node or element */
00684                 {
00685           Cp [j] = Cj [p] ;       /* save first entry of object */
00686           Cj [p] = CS_FLIP (j) ;  /* first entry is now CS_FLIP(j) */
00687                 }
00688             }
00689       for (q = 0, p = 0 ; p < cnz ; ) /* scan all of memory */
00690             {
00691           if ((j = CS_FLIP (Cj [p++])) >= 0)  /* found object j */
00692                 {
00693           Cj [q] = Cp [j] ;       /* restore first entry of object */
00694           Cp [j] = q++ ;          /* new pointer to object j */
00695           for (k3 = 0 ; k3 < len [j]-1 ; k3++) Cj [q++] = Cj [p++] ;
00696                 }
00697             }
00698       cnz = q ;                       /* Cj [cnz...nzmax-1] now free */
00699         }
00700       /* --- Construct new element ---------------------------------------- */
00701       dk = 0 ;
00702       nv [k] = -nvk ;                     /* flag k as in Lk */
00703       p = Cp [k] ;
00704       pk1 = (elenk == 0) ? p : cnz ;      /* do in place if elen[k] == 0 */
00705       pk2 = pk1 ;
00706       for (k1 = 1 ; k1 <= elenk + 1 ; k1++)
00707         {
00708       if (k1 > elenk)
00709             {
00710           e = k ;                     /* search the nodes in k */
00711           pj = p ;                    /* list of nodes starts at Cj[pj]*/
00712           ln = len [k] - elenk ;      /* length of list of nodes in k */
00713             }
00714       else
00715             {
00716           e = Cj [p++] ;              /* search the nodes in e */
00717           pj = Cp [e] ;
00718           ln = len [e] ;              /* length of list of nodes in e */
00719             }
00720       for (k2 = 1 ; k2 <= ln ; k2++)
00721             {
00722           i = Cj [pj++] ;
00723           if ((nvi = nv [i]) <= 0) continue ; /* node i dead, or seen */
00724           dk += nvi ;                 /* degree[Lk] += size of node i */
00725           nv [i] = -nvi ;             /* negate nv[i] to denote i in Lk*/
00726           Cj [pk2++] = i ;            /* place i in Lk */
00727           if (next [i] != -1) last [next [i]] = last [i] ;
00728           if (last [i] != -1)         /* remove i from degree list */
00729                 {
00730           next [last [i]] = next [i] ;
00731                 }
00732           else
00733                 {
00734           head [degree [i]] = next [i] ;
00735                 }
00736             }
00737       if (e != k)
00738             {
00739           Cp [e] = CS_FLIP (k) ;      /* absorb e into k */
00740           w [e] = 0 ;                 /* e is now a dead element */
00741             }
00742         }
00743       if (elenk != 0) cnz = pk2 ;         /* Cj [cnz...nzmax] is free */
00744       degree [k] = dk ;                   /* external degree of k - |Lk\i| */
00745       Cp [k] = pk1 ;                      /* element k is in Cj[pk1..pk2-1] */
00746       len [k] = pk2 - pk1 ;
00747       elen [k] = -2 ;                     /* k is now an element */
00748       /* --- Find set differences ----------------------------------------- */
00749       mark = csr_wclear (mark, lemax, w, n) ;  /* clear w if necessary */
00750       for (pk = pk1 ; pk < pk2 ; pk++)    /* scan 1: find |Le\Lk| */
00751         {
00752       i = Cj [pk] ;
00753       if ((eln = elen [i]) <= 0) continue ;/* skip if elen[i] empty */
00754       nvi = -nv [i] ;                      /* nv [i] was negated */
00755       wnvi = mark - nvi ;
00756       for (p = Cp [i] ; p <= Cp [i] + eln - 1 ; p++)  /* scan Ei */
00757             {
00758           e = Cj [p] ;
00759           if (w [e] >= mark)
00760                 {
00761           w [e] -= nvi ;          /* decrement |Le\Lk| */
00762                 }
00763           else if (w [e] != 0)        /* ensure e is a live element */
00764                 {
00765           w [e] = degree [e] + wnvi ; /* 1st time e seen in scan 1 */
00766                 }
00767             }
00768         }
00769       /* --- Degree update ------------------------------------------------ */
00770       for (pk = pk1 ; pk < pk2 ; pk++)    /* scan2: degree update */
00771         {
00772       i = Cj [pk] ;                   /* consider node i in Lk */
00773       p1 = Cp [i] ;
00774       p2 = p1 + elen [i] - 1 ;
00775       pn = p1 ;
00776       for (h = 0, d = 0, p = p1 ; p <= p2 ; p++)    /* scan Ei */
00777             {
00778           e = Cj [p] ;
00779           if (w [e] != 0)             /* e is an unabsorbed element */
00780                 {
00781           dext = w [e] - mark ;   /* dext = |Le\Lk| */
00782           if (dext > 0)
00783                     {
00784               d += dext ;         /* sum up the set differences */
00785               Cj [pn++] = e ;     /* keep e in Ei */
00786               h += e ;            /* compute the hash of node i */
00787                     }
00788           else
00789                     {
00790               Cp [e] = CS_FLIP (k) ;  /* aggressive absorb. e->k */
00791               w [e] = 0 ;             /* e is a dead element */
00792                     }
00793                 }
00794             }
00795       elen [i] = pn - p1 + 1 ;        /* elen[i] = |Ei| */
00796       p3 = pn ;
00797       p4 = p1 + len [i] ;
00798       for (p = p2 + 1 ; p < p4 ; p++) /* prune edges in Ai */
00799             {
00800           j = Cj [p] ;
00801           if ((nvj = nv [j]) <= 0) continue ; /* node j dead or in Lk */
00802           d += nvj ;                  /* degree(i) += |j| */
00803           Cj [pn++] = j ;             /* place j in node list of i */
00804           h += j ;                    /* compute hash for node i */
00805             }
00806       if (d == 0)                     /* check for mass elimination */
00807             {
00808           Cp [i] = CS_FLIP (k) ;      /* absorb i into k */
00809           nvi = -nv [i] ;
00810           dk -= nvi ;                 /* |Lk| -= |i| */
00811           nvk += nvi ;                /* |k| += nv[i] */
00812           nel += nvi ;
00813           nv [i] = 0 ;
00814           elen [i] = -1 ;             /* node i is dead */
00815             }
00816       else
00817             {
00818           degree [i] = CS_MIN (degree [i], d) ;   /* update degree(i) */
00819           Cj [pn] = Cj [p3] ;         /* move first node to end */
00820           Cj [p3] = Cj [p1] ;         /* move 1st el. to end of Ei */
00821           Cj [p1] = k ;               /* add k as 1st element in of Ei */
00822           
00823           len [i] = pn - p1 + 1 ;     /* new len of adj. list of node i */
00824           h %= n ;                    /* finalize hash of i */
00825           next [i] = hhead [h] ;      /* place i in hash bucket */
00826           hhead [h] = i ;
00827           last [i] = h ;              /* save hash of i in last[i] */
00828             }
00829         }                                   /* scan2 is done */
00830       degree [k] = dk ;                   /* finalize |Lk| */
00831       lemax = CS_MAX (lemax, dk) ;
00832       mark = csr_wclear (mark+lemax, lemax, w, n) ;    /* clear w */
00833       /* --- Supernode detection ------------------------------------------ */
00834       for (pk = pk1 ; pk < pk2 ; pk++)
00835         {
00836       i = Cj [pk] ;
00837       if (nv [i] >= 0) continue ;         /* skip if i is dead */
00838       h = last [i] ;                      /* scan hash bucket of node i */
00839       i = hhead [h] ;
00840       hhead [h] = -1 ;                    /* hash bucket will be empty */
00841       for ( ; i != -1 && next [i] != -1 ; i = next [i], mark++)
00842             {
00843           ln = len [i] ;
00844           eln = elen [i] ;
00845           for (p = Cp [i]+1 ; p <= Cp [i] + ln-1 ; p++) w [Cj [p]] = mark;
00846           jlast = i ;
00847           for (j = next [i] ; j != -1 ; ) /* compare i with all j */
00848                 {
00849           ok = (len [j] == ln) && (elen [j] == eln) ;
00850           for (p = Cp [j] + 1 ; ok && p <= Cp [j] + ln - 1 ; p++)
00851             {
00852               if (w [Cj [p]] != mark) ok = 0 ;    /* compare i and j*/
00853                     }
00854           if (ok)                     /* i and j are identical */
00855                     {
00856               Cp [j] = CS_FLIP (i) ;  /* absorb j into i */
00857               nv [i] += nv [j] ;
00858               nv [j] = 0 ;
00859               elen [j] = -1 ;         /* node j is dead */
00860               j = next [j] ;          /* delete j from hash bucket */
00861               next [jlast] = j ;
00862                     }
00863           else
00864                     {
00865               jlast = j ;             /* j and i are different */
00866               j = next [j] ;
00867                     }
00868                 }
00869             }
00870         }
00871       /* --- Finalize new element------------------------------------------ */
00872       for (p = pk1, pk = pk1 ; pk < pk2 ; pk++)   /* finalize Lk */
00873         {
00874       i = Cj [pk] ;
00875       if ((nvi = -nv [i]) <= 0) continue ;/* skip if i is dead */
00876       nv [i] = nvi ;                      /* restore nv[i] */
00877       d = degree [i] + dk - nvi ;         /* compute external degree(i) */
00878       d = CS_MIN (d, n - nel - nvi) ;
00879       if (head [d] != -1) last [head [d]] = i ;
00880       next [i] = head [d] ;               /* put i back in degree list */
00881       last [i] = -1 ;
00882       head [d] = i ;
00883       mindeg = CS_MIN (mindeg, d) ;       /* find new minimum degree */
00884       degree [i] = d ;
00885       Cj [p++] = i ;                      /* place i in Lk */
00886         }
00887       nv [k] = nvk ;                      /* # nodes absorbed into k */
00888       if ((len [k] = p-pk1) == 0)         /* length of adj list of element k*/
00889         {
00890       Cp [k] = -1 ;                   /* k is a root of the tree */
00891       w [k] = 0 ;                     /* k is now a dead element */
00892         }
00893       if (elenk != 0) cnz = p ;           /* free unused space in Lk */
00894     }
00895   /* --- Postordering ----------------------------------------------------- */
00896   for (i = 0 ; i < n ; i++) Cp [i] = CS_FLIP (Cp [i]) ;/* fix assembly tree */
00897   for (j = 0 ; j <= n ; j++) head [j] = -1 ;
00898   for (j = n ; j >= 0 ; j--)              /* place unordered nodes in lists */
00899     {
00900       if (nv [j] > 0) continue ;          /* skip if j is an element */
00901       next [j] = head [Cp [j]] ;          /* place j in list of its parent */
00902       head [Cp [j]] = j ;
00903     }
00904   for (e = n ; e >= 0 ; e--)              /* place elements in lists */
00905     {
00906       if (nv [e] <= 0) continue ;         /* skip unless e is an element */
00907       if (Cp [e] != -1)
00908         {
00909       next [e] = head [Cp [e]] ;      /* place e in list of its parent */
00910       head [Cp [e]] = e ;
00911         }
00912     }
00913   for (k = 0, i = 0 ; i <= n ; i++)       /* postorder the assembly tree */
00914     {
00915       if (Cp [i] == -1) k = csr_tdfs (i, k, head, next, P, w) ;
00916     }
00917   return (csr_idone (P, C, W, 1)) ;
00918 }
00919 
00920 
00921 //-----------------------------------------------------------------
00922 // Misc utilities
00923 //-----------------------------------------------------------------
00924 
00925 /* print a sparse matrix */
00926 int csr_print (const csr *A, int brief)
00927 {
00928   int p, j, m, n, nzmax, nz, nnz, *Ap, *Aj ;
00929   double *Ax ;
00930   if (!A) { printf ("(null)\n") ; return (0) ; }
00931   m = A->m ; n = A->n ; Ap = A->p ; Aj = A->j ; Ax = A->x ;
00932   nzmax = A->nzmax ; nz = A->nz ; nnz = Ap [m];
00933   if (nz < 0)
00934     {
00935       if( nnz == 0)
00936         {  /* print intermeidate matrices from csr_lu */
00937       while ((Ap[m] == 0) && (m > 0))
00938         {
00939               --m;
00940         }
00941       nnz = Ap [m];
00942         }
00943       if( nnz > 0)
00944         {
00945       printf ("%d-by-%d, nzmax: %d nnz: %d, mxnorm: %g\n", m, n, nzmax,
00946           Ap [m], csr_norm (A)) ;
00947       for (j = 0 ; j < m ; j++)
00948             {
00949           printf ("    row %d : locations %d to %d\n", j, Ap [j], Ap [j+1]-1);
00950           for (p = Ap [j] ; p < Ap [j+1] ; p++)
00951                 {
00952           printf ("      %d : %g\n", Aj [p], Ax ? Ax [p] : 1) ;
00953           if (brief && p > 20) { printf ("  ...\n") ; return (1) ; }
00954                 }
00955             }
00956         }else{
00957     printf ("%d-by-%d, zero matrix with nzmax: %d\n", m, n, nzmax);
00958       }
00959     }
00960   else
00961     {
00962       printf ("triplet: %d-by-%d, nzmax: %d nnz: %d\n", m, n, nzmax, nz) ;
00963       for (p = 0 ; p < nz ; p++)
00964         {
00965       printf ("    %d %d : %g\n", Aj [p], Ap [p], Ax ? Ax [p] : 1) ;
00966       if (brief && p > 20) { printf ("  ...\n") ; return (1) ; }
00967         }
00968     }
00969   return (1) ;
00970 }
00971 
00972 /* infinity-norm of a sparse matrix = norm(A,inf), max row sum */
00973 double csr_norm (const csr *A)
00974 {
00975     int p, j, m, *Ap ;
00976     double *Ax,  norm = 0, s ;
00977     if (!CS_CSC (A) || !A->x) return (-1) ;        /* check inputs */
00978     m = A->m ; Ap = A->p ; Ax = A->x ;
00979     for (j = 0 ; j < m ; j++)
00980     {
00981         for (s = 0, p = Ap [j] ; p < Ap [j+1] ; p++) s += fabs (Ax [p]);
00982         norm = CS_MAX (norm, s) ;
00983     }
00984     return (norm) ;
00985 }
00986 
00987 /* drop entries for which fkeep(A(i,j)) is false; return nz if OK, else -1 */
00988 int csr_fkeep (csr *A, int (*fkeep) (int, int, double, void *), void *other)
00989 {
00990     int j, p, nz = 0, m, *Ap, *Aj ;
00991     double *Ax ;
00992     if (!CS_CSC (A) || !fkeep) return (-1) ;    /* check inputs */
00993     m = A->m ; Ap = A->p ; Aj = A->j ; Ax = A->x ;
00994     for (j = 0 ; j < m ; j++)
00995     {
00996         p = Ap [j] ;                        /* get current location of col j */
00997         Ap [j] = nz ;                       /* record new location of col j */
00998         for ( ; p < Ap [j+1] ; p++)
00999         {
01000             if (fkeep (Aj [p], j, Ax ? Ax [p] : 1, other))
01001             {
01002                 if (Ax) Ax [nz] = Ax [p] ;  /* keep A(i,j) */
01003                 Aj [nz++] = Aj [p] ;
01004             }
01005         }
01006     }
01007     Ap [m] = nz ;                           /* finalize A */
01008     csr_sprealloc (A, 0) ;                  /* remove extra space from A */
01009     return (nz) ;
01010 }

Generated on Wed May 12 21:30:18 2010 for IFPACK by  doxygen 1.4.7