Amesos Package Browser (Single Doxygen Collection) Development
amesos_amd_l_aat.c
Go to the documentation of this file.
```00001 /* ========================================================================= */
00002 /* === AMD_aat ============================================================= */
00003 /* ========================================================================= */
00004
00005 /* ------------------------------------------------------------------------- */
00006 /* AMD, Copyright (c) Timothy A. Davis,              */
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/amd                          */
00010 /* ------------------------------------------------------------------------- */
00011
00012 /* AMD_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  * (AMD_valid (n, n, Ap, Ai) must be AMD_OK, but this condition is not
00016  * checked).
00017  */
00018
00019 /* This file should make the long int version of AMD */
00020 #define DLONG 1
00021
00022 #include "amesos_amd_internal.h"
00023
00024 GLOBAL size_t AMD_aat /* returns nz in A+A' */
00025 (
00026     Int n,
00027     const Int Ap [ ],
00028     const Int Ai [ ],
00029     Int Len [ ],  /* Len [j]: length of column j of A+A', excl diagonal*/
00030     Int Tp [ ],   /* workspace of size n */
00031     double Info [ ]
00032 )
00033 {
00034     Int p1, p2, p, i, j, pj, pj2, k, nzdiag, nzboth, nz ;
00035     double sym ;
00036     size_t nzaat ;
00037
00038 #ifndef NDEBUG
00039     AMD_debug_init ("AMD AAT") ;
00040     for (k = 0 ; k < n ; k++) Tp [k] = EMPTY ;
00041     ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
00042 #endif
00043
00044     if (Info != (double *) NULL)
00045     {
00046   /* clear the Info array, if it exists */
00047   for (i = 0 ; i < AMD_INFO ; i++)
00048   {
00049       Info [i] = EMPTY ;
00050   }
00051   Info [AMD_STATUS] = AMD_OK ;
00052     }
00053
00054     for (k = 0 ; k < n ; k++)
00055     {
00056   Len [k] = 0 ;
00057     }
00058
00059     nzdiag = 0 ;
00060     nzboth = 0 ;
00061     nz = Ap [n] ;
00062
00063     for (k = 0 ; k < n ; k++)
00064     {
00065   p1 = Ap [k] ;
00066   p2 = Ap [k+1] ;
00067   AMD_DEBUG2 (("\nAAT Column: "ID" p1: "ID" p2: "ID"\n", k, p1, p2)) ;
00068
00069   /* construct A+A' */
00070   for (p = p1 ; p < p2 ; )
00071   {
00072       /* scan the upper triangular part of A */
00073       j = Ai [p] ;
00074       if (j < k)
00075       {
00076     /* entry A (j,k) is in the strictly upper triangular part,
00077      * add both A (j,k) and A (k,j) to the matrix A+A' */
00078     Len [j]++ ;
00079     Len [k]++ ;
00080     AMD_DEBUG3 (("    upper ("ID","ID") ("ID","ID")\n", j,k, k,j));
00081     p++ ;
00082       }
00083       else if (j == k)
00084       {
00085     /* skip the diagonal */
00086     p++ ;
00087     nzdiag++ ;
00088     break ;
00089       }
00090       else /* j > k */
00091       {
00092     /* first entry below the diagonal */
00093     break ;
00094       }
00095       /* scan lower triangular part of A, in column j until reaching
00096        * row k.  Start where last scan left off. */
00097       ASSERT (Tp [j] != EMPTY) ;
00098       ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
00099       pj2 = Ap [j+1] ;
00100       for (pj = Tp [j] ; pj < pj2 ; )
00101       {
00102     i = Ai [pj] ;
00103     if (i < k)
00104     {
00105         /* A (i,j) is only in the lower part, not in upper.
00106          * add both A (i,j) and A (j,i) to the matrix A+A' */
00107         Len [i]++ ;
00108         Len [j]++ ;
00109         AMD_DEBUG3 (("    lower ("ID","ID") ("ID","ID")\n",
00110       i,j, j,i)) ;
00111         pj++ ;
00112     }
00113     else if (i == k)
00114     {
00115         /* entry A (k,j) in lower part and A (j,k) in upper */
00116         pj++ ;
00117         nzboth++ ;
00118         break ;
00119     }
00120     else /* i > k */
00121     {
00122         /* consider this entry later, when k advances to i */
00123         break ;
00124     }
00125       }
00126       Tp [j] = pj ;
00127   }
00128   /* Tp [k] points to the entry just below the diagonal in column k */
00129   Tp [k] = p ;
00130     }
00131
00132     /* clean up, for remaining mismatched entries */
00133     for (j = 0 ; j < n ; j++)
00134     {
00135   for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
00136   {
00137       i = Ai [pj] ;
00138       /* A (i,j) is only in the lower part, not in upper.
00139        * add both A (i,j) and A (j,i) to the matrix A+A' */
00140       Len [i]++ ;
00141       Len [j]++ ;
00142       AMD_DEBUG3 (("    lower cleanup ("ID","ID") ("ID","ID")\n",
00143     i,j, j,i)) ;
00144   }
00145     }
00146
00147     /* --------------------------------------------------------------------- */
00148     /* compute the symmetry of the nonzero pattern of A */
00149     /* --------------------------------------------------------------------- */
00150
00151     /* Given a matrix A, the symmetry of A is:
00152      *  B = tril (spones (A), -1) + triu (spones (A), 1) ;
00153      *  sym = nnz (B & B') / nnz (B) ;
00154      *  or 1 if nnz (B) is zero.
00155      */
00156
00157     if (nz == nzdiag)
00158     {
00159   sym = 1 ;
00160     }
00161     else
00162     {
00163   sym = (2 * (double) nzboth) / ((double) (nz - nzdiag)) ;
00164     }
00165
00166     nzaat = 0 ;
00167     for (k = 0 ; k < n ; k++)
00168     {
00169   nzaat += Len [k] ;
00170     }
00171
00172     AMD_DEBUG1 (("AMD nz in A+A', excluding diagonal (nzaat) = %g\n",
00173   (double) nzaat)) ;
00174     AMD_DEBUG1 (("   nzboth: "ID" nz: "ID" nzdiag: "ID" symmetry: %g\n",
00175     nzboth, nz, nzdiag, sym)) ;
00176
00177     if (Info != (double *) NULL)
00178     {
00179   Info [AMD_STATUS] = AMD_OK ;
00180   Info [AMD_N] = n ;
00181   Info [AMD_NZ] = nz ;
00182   Info [AMD_SYMMETRY] = sym ;     /* symmetry of pattern of A */
00183   Info [AMD_NZDIAG] = nzdiag ;      /* nonzeros on diagonal of A */
00184   Info [AMD_NZ_A_PLUS_AT] = nzaat ;   /* nonzeros in A+A' */
00185     }
00186
00187     return (nzaat) ;
00188 }
```