Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_sort.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === KLU_sort ============================================================= */
00003 /* ========================================================================== */
00004 
00005 /* sorts the columns of L and U so that the row indices appear in strictly
00006  * increasing order.
00007  */
00008 
00009 #include "amesos_klu_internal.h"
00010 
00011 /* ========================================================================== */
00012 /* === sort ================================================================= */
00013 /* ========================================================================== */
00014 
00015 /* Sort L or U using a double-transpose */
00016 
00017 static void sort (Int n, Int *Xip, Int *Xlen, Unit *LU, Int *Tp, Int *Tj,
00018     Entry *Tx, Int *W)
00019 {
00020     Int *Xi ;
00021     Entry *Xx ;
00022     Int p, i, j, len, nz, tp, xlen, pend ;
00023 
00024     ASSERT (KLU_valid_LU (n, FALSE, Xip, Xlen, LU)) ;
00025 
00026     /* count the number of entries in each row of L or U */ 
00027     for (i = 0 ; i < n ; i++)
00028     {
00029   W [i] = 0 ;
00030     }
00031     for (j = 0 ; j < n ; j++)
00032     {
00033   GET_POINTER (LU, Xip, Xlen, Xi, Xx, j, len) ;
00034   for (p = 0 ; p < len ; p++)
00035   {
00036       W [Xi [p]]++ ;
00037   }
00038     }
00039 
00040     /* construct the row pointers for T */
00041     nz = 0 ;
00042     for (i = 0 ; i < n ; i++)
00043     {
00044   Tp [i] = nz ;
00045   nz += W [i] ;
00046     }
00047     Tp [n] = nz ;
00048     for (i = 0 ; i < n ; i++)
00049     {
00050   W [i] = Tp [i] ;
00051     }
00052 
00053     /* transpose the matrix into Tp, Ti, Tx */
00054     for (j = 0 ; j < n ; j++)
00055     {
00056   GET_POINTER (LU, Xip, Xlen, Xi, Xx, j, len) ;
00057   for (p = 0 ; p < len ; p++)
00058   {
00059       tp = W [Xi [p]]++ ;
00060       Tj [tp] = j ;
00061       Tx [tp] = Xx [p] ;
00062   }
00063     }
00064 
00065     /* transpose the matrix back into Xip, Xlen, Xi, Xx */
00066     for (j = 0 ; j < n ; j++)
00067     {
00068   W [j] = 0 ;
00069     }
00070     for (i = 0 ; i < n ; i++)
00071     {
00072   pend = Tp [i+1] ;
00073   for (p = Tp [i] ; p < pend ; p++)
00074   {
00075       j = Tj [p] ;
00076       GET_POINTER (LU, Xip, Xlen, Xi, Xx, j, len) ;
00077       xlen = W [j]++ ;
00078       Xi [xlen] = i ;
00079       Xx [xlen] = Tx [p] ;
00080   }
00081     }
00082 
00083     ASSERT (KLU_valid_LU (n, FALSE, Xip, Xlen, LU)) ;
00084 }
00085 
00086 
00087 /* ========================================================================== */
00088 /* === KLU_sort ============================================================= */
00089 /* ========================================================================== */
00090 
00091 Int KLU_sort
00092 (
00093     KLU_symbolic *Symbolic,
00094     KLU_numeric *Numeric,
00095     KLU_common *Common
00096 )
00097 {
00098     Int *R, *W, *Tp, *Ti, *Lip, *Uip, *Llen, *Ulen ;
00099     Entry *Tx ;
00100     Unit **LUbx ;
00101     Int n, nk, nz, block, nblocks, maxblock, k1 ;
00102     size_t m1 ;
00103 
00104     if (Common == NULL)
00105     {
00106   return (FALSE) ;
00107     }
00108     Common->status = KLU_OK ;
00109 
00110     n = Symbolic->n ;
00111     R = Symbolic->R ;
00112     nblocks = Symbolic->nblocks ;
00113     maxblock = Symbolic->maxblock ;
00114 
00115     Lip  = Numeric->Lip ;
00116     Llen = Numeric->Llen ;
00117     Uip  = Numeric->Uip ;
00118     Ulen = Numeric->Ulen ;
00119     LUbx = (Unit **) Numeric->LUbx ;
00120 
00121     m1 = ((size_t) maxblock) + 1 ;
00122 
00123     /* allocate workspace */
00124     nz = MAX (Numeric->max_lnz_block, Numeric->max_unz_block) ;
00125     W  = KLU_malloc (maxblock, sizeof (Int), Common) ;
00126     Tp = KLU_malloc (m1, sizeof (Int), Common) ;
00127     Ti = KLU_malloc (nz, sizeof (Int), Common) ;
00128     Tx = KLU_malloc (nz, sizeof (Entry), Common) ;
00129 
00130     PRINTF (("\n======================= Start sort:\n")) ;
00131 
00132     if (Common->status == KLU_OK)
00133     {
00134   /* sort each block of L and U */
00135   for (block = 0 ; block < nblocks ; block++)
00136   {
00137       k1 = R [block] ;
00138       nk = R [block+1] - k1 ;
00139       if (nk > 1)
00140       {
00141     PRINTF (("\n-------------------block: %d nk %d\n", block, nk)) ;
00142     sort (nk, Lip + k1, Llen + k1, LUbx [block], Tp, Ti, Tx, W) ;
00143     sort (nk, Uip + k1, Ulen + k1, LUbx [block], Tp, Ti, Tx, W) ;
00144       }
00145   }
00146     }
00147 
00148     PRINTF (("\n======================= sort done.\n")) ;
00149 
00150     /* free workspace */
00151     KLU_free (W, maxblock, sizeof (Int), Common) ;
00152     KLU_free (Tp, m1, sizeof (Int), Common) ;
00153     KLU_free (Ti, nz, sizeof (Int), Common) ;
00154     KLU_free (Tx, nz, sizeof (Entry), Common) ;
00155     return (Common->status == KLU_OK) ;
00156 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines