Amesos Package Browser (Single Doxygen Collection) Development
amesos_camd_aat.c
Go to the documentation of this file.
```00001 /* ========================================================================= */
00002 /* === CAMD_aat ============================================================ */
00003 /* ========================================================================= */
00004
00005 /* ------------------------------------------------------------------------- */
00006 /* CAMD, Copyright (c) Timothy A. Davis, Yanqing Chen,           */
00007 /* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
00008 /* email: davis at cise.ufl.edu    CISE Department, Univ. of Florida.        */
00009 /* web: http://www.cise.ufl.edu/research/sparse/camd                         */
00010 /* ------------------------------------------------------------------------- */
00011
00012 /* CAMD_aat:  compute the symmetry of the pattern of A, and count the number of
00013  * nonzeros each column of A+A' (excluding the diagonal).  Assumes the input
00014  * matrix has no errors, with sorted columns and no duplicates
00015  * (CAMD_valid (n, n, Ap, Ai) must be CAMD_OK, but this condition is not
00016  * checked).
00017  */
00018
00019 #include "amesos_camd_internal.h"
00020
00021 GLOBAL size_t CAMD_aat  /* returns nz in A+A' */
00022 (
00023     Int n,
00024     const Int Ap [ ],
00025     const Int Ai [ ],
00026     Int Len [ ],  /* Len [j]: length of column j of A+A', excl diagonal*/
00027     Int Tp [ ],   /* workspace of size n */
00028     double Info [ ]
00029 )
00030 {
00031     Int p1, p2, p, i, j, pj, pj2, k, nzdiag, nzboth, nz ;
00032     double sym ;
00033     size_t nzaat ;
00034
00035 #ifndef NDEBUG
00036     CAMD_debug_init ("CAMD AAT") ;
00037     for (k = 0 ; k < n ; k++) Tp [k] = EMPTY ;
00038     ASSERT (CAMD_valid (n, n, Ap, Ai) == CAMD_OK) ;
00039 #endif
00040
00041     if (Info != (double *) NULL)
00042     {
00043   /* clear the Info array, if it exists */
00044   for (i = 0 ; i < CAMD_INFO ; i++)
00045   {
00046       Info [i] = EMPTY ;
00047   }
00048   Info [CAMD_STATUS] = CAMD_OK ;
00049     }
00050
00051     for (k = 0 ; k < n ; k++)
00052     {
00053   Len [k] = 0 ;
00054     }
00055
00056     nzdiag = 0 ;
00057     nzboth = 0 ;
00058     nz = Ap [n] ;
00059
00060     for (k = 0 ; k < n ; k++)
00061     {
00062   p1 = Ap [k] ;
00063   p2 = Ap [k+1] ;
00064   CAMD_DEBUG2 (("\nAAT Column: "ID" p1: "ID" p2: "ID"\n", k, p1, p2)) ;
00065
00066   /* construct A+A' */
00067   for (p = p1 ; p < p2 ; )
00068   {
00069       /* scan the upper triangular part of A */
00070       j = Ai [p] ;
00071       if (j < k)
00072       {
00073     /* entry A (j,k) is in the strictly upper triangular part,
00074      * add both A (j,k) and A (k,j) to the matrix A+A' */
00075     Len [j]++ ;
00076     Len [k]++ ;
00077     CAMD_DEBUG3 (("    upper ("ID","ID") ("ID","ID")\n", j,k, k,j));
00078     p++ ;
00079       }
00080       else if (j == k)
00081       {
00082     /* skip the diagonal */
00083     p++ ;
00084     nzdiag++ ;
00085     break ;
00086       }
00087       else /* j > k */
00088       {
00089     /* first entry below the diagonal */
00090     break ;
00091       }
00092       /* scan lower triangular part of A, in column j until reaching
00093        * row k.  Start where last scan left off. */
00094       ASSERT (Tp [j] != EMPTY) ;
00095       ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
00096       pj2 = Ap [j+1] ;
00097       for (pj = Tp [j] ; pj < pj2 ; )
00098       {
00099     i = Ai [pj] ;
00100     if (i < k)
00101     {
00102         /* A (i,j) is only in the lower part, not in upper.
00103          * add both A (i,j) and A (j,i) to the matrix A+A' */
00104         Len [i]++ ;
00105         Len [j]++ ;
00106         CAMD_DEBUG3 (("    lower ("ID","ID") ("ID","ID")\n",
00107       i,j, j,i)) ;
00108         pj++ ;
00109     }
00110     else if (i == k)
00111     {
00112         /* entry A (k,j) in lower part and A (j,k) in upper */
00113         pj++ ;
00114         nzboth++ ;
00115         break ;
00116     }
00117     else /* i > k */
00118     {
00119         /* consider this entry later, when k advances to i */
00120         break ;
00121     }
00122       }
00123       Tp [j] = pj ;
00124   }
00125   /* Tp [k] points to the entry just below the diagonal in column k */
00126   Tp [k] = p ;
00127     }
00128
00129     /* clean up, for remaining mismatched entries */
00130     for (j = 0 ; j < n ; j++)
00131     {
00132   for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
00133   {
00134       i = Ai [pj] ;
00135       /* A (i,j) is only in the lower part, not in upper.
00136        * add both A (i,j) and A (j,i) to the matrix A+A' */
00137       Len [i]++ ;
00138       Len [j]++ ;
00139       CAMD_DEBUG3 (("    lower cleanup ("ID","ID") ("ID","ID")\n",
00140     i,j, j,i)) ;
00141   }
00142     }
00143
00144     /* --------------------------------------------------------------------- */
00145     /* compute the symmetry of the nonzero pattern of A */
00146     /* --------------------------------------------------------------------- */
00147
00148     /* Given a matrix A, the symmetry of A is:
00149      *  B = tril (spones (A), -1) + triu (spones (A), 1) ;
00150      *  sym = nnz (B & B') / nnz (B) ;
00151      *  or 1 if nnz (B) is zero.
00152      */
00153
00154     if (nz == nzdiag)
00155     {
00156   sym = 1 ;
00157     }
00158     else
00159     {
00160   sym = (2 * (double) nzboth) / ((double) (nz - nzdiag)) ;
00161     }
00162
00163     nzaat = 0 ;
00164     for (k = 0 ; k < n ; k++)
00165     {
00166   nzaat += Len [k] ;
00167     }
00168     CAMD_DEBUG1 (("CAMD nz in A+A', excluding diagonal (nzaat) = %g\n",
00169   (double) nzaat)) ;
00170     CAMD_DEBUG1 (("   nzboth: "ID" nz: "ID" nzdiag: "ID" symmetry: %g\n",
00171     nzboth, nz, nzdiag, sym)) ;
00172
00173     if (Info != (double *) NULL)
00174     {
00175   Info [CAMD_STATUS] = CAMD_OK ;
00176   Info [CAMD_N] = n ;
00177   Info [CAMD_NZ] = nz ;
00178   Info [CAMD_SYMMETRY] = sym ;      /* symmetry of pattern of A */
00179   Info [CAMD_NZDIAG] = nzdiag ;     /* nonzeros on diagonal of A */
00180   Info [CAMD_NZ_A_PLUS_AT] = nzaat ;   /* nonzeros in A+A' */
00181     }
00182
00183     return (nzaat) ;
00184 }
```