Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_l_analyze_given.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === klu_analyze_given ==================================================== */
00003 /* ========================================================================== */
00004 
00005 /* Given an input permutation P and Q, create the Symbolic object.  BTF can
00006  * be done to modify the user's P and Q (does not perform the max transversal;
00007  * just finds the strongly-connected components). */
00008 
00009 /* This file should make the long int version of KLU */
00010 #define DLONG 1
00011 
00012 #include "amesos_klu_internal.h"
00013 
00014 /* ========================================================================== */
00015 /* === klu_alloc_symbolic =================================================== */
00016 /* ========================================================================== */
00017 
00018 /* Allocate Symbolic object, and check input matrix.  Not user callable. */
00019 
00020 KLU_symbolic *KLU_alloc_symbolic
00021 (
00022     Int n,
00023     Int *Ap,
00024     Int *Ai,
00025     KLU_common *Common
00026 )
00027 {
00028     KLU_symbolic *Symbolic ;
00029     Int *P, *Q, *R ;
00030     double *Lnz ;
00031     Int nz, i, j, p, pend ;
00032 
00033     if (Common == NULL)
00034     {
00035   return (NULL) ;
00036     }
00037     Common->status = KLU_OK ;
00038 
00039     /* A is n-by-n, with n > 0.  Ap [0] = 0 and nz = Ap [n] >= 0 required.
00040      * Ap [j] <= Ap [j+1] must hold for all j = 0 to n-1.  Row indices in Ai
00041      * must be in the range 0 to n-1, and no duplicate entries can be present.
00042      * The list of row indices in each column of A need not be sorted.
00043      */
00044 
00045     if (n <= 0 || Ap == NULL || Ai == NULL)
00046     {
00047   /* Ap and Ai must be present, and n must be > 0 */
00048   Common->status = KLU_INVALID ;
00049   return (NULL) ;
00050     }
00051 
00052     nz = Ap [n] ;
00053     if (Ap [0] != 0 || nz < 0)
00054     {
00055   /* nz must be >= 0 and Ap [0] must equal zero */
00056   Common->status = KLU_INVALID ;
00057   return (NULL) ;
00058     }
00059 
00060     for (j = 0 ; j < n ; j++)
00061     {
00062   if (Ap [j] > Ap [j+1])
00063   {
00064       /* column pointers must be non-decreasing */
00065       Common->status = KLU_INVALID ;
00066       return (NULL) ;
00067   }
00068     }
00069     P = KLU_malloc (n, sizeof (Int), Common) ;
00070     if (Common->status < KLU_OK)
00071     {
00072   /* out of memory */
00073   Common->status = KLU_OUT_OF_MEMORY ;
00074   return (NULL) ;
00075     }
00076     for (i = 0 ; i < n ; i++)
00077     {
00078   P [i] = EMPTY ;
00079     }
00080     for (j = 0 ; j < n ; j++)
00081     {
00082   pend = Ap [j+1] ;
00083   for (p = Ap [j] ; p < pend ; p++)
00084   {
00085       i = Ai [p] ;
00086       if (i < 0 || i >= n || P [i] == j)
00087       {
00088     /* row index out of range, or duplicate entry */
00089     KLU_free (P, n, sizeof (Int), Common) ;
00090     Common->status = KLU_INVALID ;
00091     return (NULL) ;
00092       }
00093       /* flag row i as appearing in column j */
00094       P [i] = j ;
00095   }
00096     }
00097 
00098     /* ---------------------------------------------------------------------- */
00099     /* allocate the Symbolic object */
00100     /* ---------------------------------------------------------------------- */
00101 
00102     Symbolic = KLU_malloc (sizeof (KLU_symbolic), 1, Common) ;
00103     if (Common->status < KLU_OK)
00104     {
00105   /* out of memory */
00106   KLU_free (P, n, sizeof (Int), Common) ;
00107   Common->status = KLU_OUT_OF_MEMORY ;
00108   return (NULL) ;
00109     }
00110 
00111     Q = KLU_malloc (n, sizeof (Int), Common) ;
00112     R = KLU_malloc (n+1, sizeof (Int), Common) ;
00113     Lnz = KLU_malloc (n, sizeof (double), Common) ;
00114 
00115     Symbolic->n = n ;
00116     Symbolic->nz = nz ;
00117     Symbolic->P = P ;
00118     Symbolic->Q = Q ;
00119     Symbolic->R = R ;
00120     Symbolic->Lnz = Lnz ;
00121 
00122     if (Common->status < KLU_OK)
00123     {
00124   /* out of memory */
00125   KLU_free_symbolic (&Symbolic, Common) ;
00126   Common->status = KLU_OUT_OF_MEMORY ;
00127   return (NULL) ;
00128     }
00129 
00130     return (Symbolic) ;
00131 }
00132 
00133 
00134 /* ========================================================================== */
00135 /* === KLU_analyze_given ==================================================== */
00136 /* ========================================================================== */
00137 
00138 KLU_symbolic *KLU_analyze_given     /* returns NULL if error, or a valid
00139                KLU_symbolic object if successful */
00140 (
00141     /* inputs, not modified */
00142     Int n,    /* A is n-by-n */
00143     Int Ap [ ],   /* size n+1, column pointers */
00144     Int Ai [ ],   /* size nz, row indices */
00145     Int Puser [ ],  /* size n, user's row permutation (may be NULL) */
00146     Int Quser [ ],  /* size n, user's column permutation (may be NULL) */
00147     /* -------------------- */
00148     KLU_common *Common
00149 )
00150 {
00151     KLU_symbolic *Symbolic ;
00152     double *Lnz ;
00153     Int nblocks, nz, block, maxblock, *P, *Q, *R, nzoff, p, pend, do_btf, k ;
00154 
00155     /* ---------------------------------------------------------------------- */
00156     /* determine if input matrix is valid, and get # of nonzeros */
00157     /* ---------------------------------------------------------------------- */
00158 
00159     Symbolic = KLU_alloc_symbolic (n, Ap, Ai, Common) ;
00160     if (Symbolic == NULL)
00161     {
00162   return (NULL) ;
00163     }
00164     P = Symbolic->P ;
00165     Q = Symbolic->Q ;
00166     R = Symbolic->R ;
00167     Lnz = Symbolic->Lnz ;
00168     nz = Symbolic->nz ;
00169 
00170     /* ---------------------------------------------------------------------- */
00171     /* Q = Quser, or identity if Quser is NULL */
00172     /* ---------------------------------------------------------------------- */
00173 
00174     if (Quser == (Int *) NULL)
00175     {
00176   for (k = 0 ; k < n ; k++)
00177   {
00178       Q [k] = k ;
00179   }
00180     }
00181     else
00182     {
00183   for (k = 0 ; k < n ; k++)
00184   {
00185       Q [k] = Quser [k] ;
00186   }
00187     }
00188 
00189     /* ---------------------------------------------------------------------- */
00190     /* get the control parameters for BTF and ordering method */
00191     /* ---------------------------------------------------------------------- */
00192 
00193     do_btf = Common->btf ;
00194     do_btf = (do_btf) ? TRUE : FALSE ;
00195     Symbolic->ordering = 2 ;
00196     Symbolic->do_btf = do_btf ;
00197 
00198     /* ---------------------------------------------------------------------- */
00199     /* find the block triangular form, if requested */
00200     /* ---------------------------------------------------------------------- */
00201 
00202     if (do_btf)
00203     {
00204 
00205   /* ------------------------------------------------------------------ */
00206   /* get workspace for BTF_strongcomp */
00207   /* ------------------------------------------------------------------ */
00208 
00209   Int *Pinv, *Work, *Bi, k1, k2, nk, oldcol ;
00210 
00211   Work = KLU_malloc (4*n, sizeof (Int), Common) ;
00212   Pinv = KLU_malloc (n, sizeof (Int), Common) ;
00213   if (Puser != (Int *) NULL)
00214   {
00215       Bi = KLU_malloc (nz+1, sizeof (Int), Common) ;
00216   }
00217   else
00218   {
00219       Bi = Ai ;
00220   }
00221 
00222   if (Common->status < KLU_OK)
00223   {
00224       /* out of memory */
00225       KLU_free (Work, 4*n, sizeof (Int), Common) ;
00226       KLU_free (Pinv, n, sizeof (Int), Common) ;
00227       if (Puser != (Int *) NULL)
00228       {
00229     KLU_free (Bi, nz+1, sizeof (Int), Common) ;
00230       }
00231       KLU_free_symbolic (&Symbolic, Common) ;
00232       Common->status = KLU_OUT_OF_MEMORY ;
00233       return (NULL) ;
00234   }
00235 
00236   /* ------------------------------------------------------------------ */
00237   /* B = Puser * A */
00238   /* ------------------------------------------------------------------ */
00239 
00240   if (Puser != (Int *) NULL)
00241   {
00242       for (k = 0 ; k < n ; k++)
00243       {
00244     Pinv [Puser [k]] = k ;
00245       }
00246       for (p = 0 ; p < nz ; p++)
00247       {
00248     Bi [p] = Pinv [Ai [p]] ;
00249       }
00250   }
00251 
00252   /* ------------------------------------------------------------------ */
00253   /* find the strongly-connected components */
00254   /* ------------------------------------------------------------------ */
00255 
00256   /* modifies Q, and determines P and R */
00257   nblocks = BTF_strongcomp (n, Ap, Bi, Q, P, R, Work) ;
00258 
00259   /* ------------------------------------------------------------------ */
00260   /* P = P * Puser */
00261   /* ------------------------------------------------------------------ */
00262 
00263   if (Puser != (Int *) NULL)
00264   {
00265       for (k = 0 ; k < n ; k++)
00266       {
00267     Work [k] = Puser [P [k]] ;
00268       }
00269       for (k = 0 ; k < n ; k++)
00270       {
00271     P [k] = Work [k] ;
00272       }
00273   }
00274 
00275   /* ------------------------------------------------------------------ */
00276   /* Pinv = inverse of P */
00277   /* ------------------------------------------------------------------ */
00278 
00279   for (k = 0 ; k < n ; k++)
00280   {
00281       Pinv [P [k]] = k ;
00282   }
00283 
00284   /* ------------------------------------------------------------------ */
00285   /* analyze each block */
00286   /* ------------------------------------------------------------------ */
00287 
00288   nzoff = 0 ;     /* nz in off-diagonal part */
00289   maxblock = 1 ;      /* size of the largest block */
00290 
00291   for (block = 0 ; block < nblocks ; block++)
00292   {
00293 
00294       /* -------------------------------------------------------------- */
00295       /* the block is from rows/columns k1 to k2-1 */
00296       /* -------------------------------------------------------------- */
00297 
00298       k1 = R [block] ;
00299       k2 = R [block+1] ;
00300       nk = k2 - k1 ;
00301       PRINTF (("BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1, k2-1, nk)) ;
00302       maxblock = MAX (maxblock, nk) ;
00303 
00304       /* -------------------------------------------------------------- */
00305       /* scan the kth block, C */
00306       /* -------------------------------------------------------------- */
00307 
00308       for (k = k1 ; k < k2 ; k++)
00309       {
00310     oldcol = Q [k] ;
00311     pend = Ap [oldcol+1] ;
00312     for (p = Ap [oldcol] ; p < pend ; p++)
00313     {
00314         if (Pinv [Ai [p]] < k1)
00315         {
00316       nzoff++ ;
00317         }
00318     }
00319       }
00320 
00321       /* fill-in not estimated */
00322       Lnz [block] = EMPTY ;
00323   }
00324 
00325   /* ------------------------------------------------------------------ */
00326   /* free all workspace */
00327   /* ------------------------------------------------------------------ */
00328 
00329   KLU_free (Work, 4*n, sizeof (Int), Common) ;
00330   KLU_free (Pinv, n, sizeof (Int), Common) ;
00331   if (Puser != (Int *) NULL)
00332   {
00333       KLU_free (Bi, nz+1, sizeof (Int), Common) ;
00334   }
00335 
00336     }
00337     else
00338     {
00339 
00340   /* ------------------------------------------------------------------ */
00341   /* BTF not requested */
00342   /* ------------------------------------------------------------------ */
00343 
00344   nzoff = 0 ;
00345   nblocks = 1 ;
00346   maxblock = n ;
00347   R [0] = 0 ;
00348   R [1] = n ;
00349   Lnz [0] = EMPTY ;
00350 
00351   /* ------------------------------------------------------------------ */
00352   /* P = Puser, or identity if Puser is NULL */
00353   /* ------------------------------------------------------------------ */
00354 
00355   for (k = 0 ; k < n ; k++)
00356   {
00357       P [k] = (Puser == NULL) ? k : Puser [k] ;
00358   }
00359     }
00360 
00361     /* ---------------------------------------------------------------------- */
00362     /* return the symbolic object */
00363     /* ---------------------------------------------------------------------- */
00364 
00365     Symbolic->nblocks = nblocks ;
00366     Symbolic->maxblock = maxblock ;
00367     Symbolic->lnz = EMPTY ;
00368     Symbolic->unz = EMPTY ;
00369     Symbolic->nzoff = nzoff ;
00370 
00371     return (Symbolic) ;
00372 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines