Amesos Package Browser (Single Doxygen Collection) Development
amesos_klu_extract.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === KLU_extract ========================================================== */
00003 /* ========================================================================== */
00004 
00005 /* Extract KLU factorization into conventional compressed-column matrices.
00006  * If any output array is NULL, that part of the LU factorization is not
00007  * extracted (this is not an error condition).
00008  *
00009  * nnz(L) = Numeric->lnz, nnz(U) = Numeric->unz, and nnz(F) = Numeric->Offp [n]
00010  */
00011 
00012 #include "amesos_klu_internal.h"
00013 
00014 Int KLU_extract     /* returns TRUE if successful, FALSE otherwise */
00015 (
00016     /* inputs: */
00017     KLU_numeric *Numeric,
00018     KLU_symbolic *Symbolic,
00019 
00020     /* outputs, all of which must be allocated on input */
00021 
00022     /* L */
00023     Int *Lp,      /* size n+1 */
00024     Int *Li,      /* size nnz(L) */
00025     double *Lx,     /* size nnz(L) */
00026 #ifdef COMPLEX
00027     double *Lz,     /* size nnz(L) for the complex case, ignored if real */
00028 #endif
00029 
00030     /* U */
00031     Int *Up,      /* size n+1 */
00032     Int *Ui,      /* size nnz(U) */
00033     double *Ux,     /* size nnz(U) */
00034 #ifdef COMPLEX
00035     double *Uz,     /* size nnz(U) for the complex case, ignored if real */
00036 #endif
00037 
00038     /* F */
00039     Int *Fp,      /* size n+1 */
00040     Int *Fi,      /* size nnz(F) */
00041     double *Fx,     /* size nnz(F) */
00042 #ifdef COMPLEX
00043     double *Fz,     /* size nnz(F) for the complex case, ignored if real */
00044 #endif
00045 
00046     /* P, row permutation */
00047     Int *P,     /* size n */
00048 
00049     /* Q, column permutation */
00050     Int *Q,     /* size n */
00051 
00052     /* Rs, scale factors */
00053     double *Rs,     /* size n */
00054 
00055     /* R, block boundaries */
00056     Int *R,     /* size nblocks+1 */
00057 
00058     KLU_common *Common
00059 )
00060 {
00061     Int *Lip, *Llen, *Uip, *Ulen, *Li2, *Ui2 ;
00062     Unit *LU ;
00063     Entry *Lx2, *Ux2, *Ukk ;
00064     Int i, k, block, nblocks, n, nz, k1, k2, nk, len, kk, p ;
00065 
00066     if (Common == NULL)
00067     {
00068   return (FALSE) ;
00069     }
00070 
00071     if (Symbolic == NULL || Numeric == NULL)
00072     {
00073   Common->status = KLU_INVALID ;
00074   return (FALSE) ;
00075     }
00076 
00077     Common->status = KLU_OK ;
00078     n = Symbolic->n ;
00079     nblocks = Symbolic->nblocks ;
00080 
00081     /* ---------------------------------------------------------------------- */
00082     /* extract scale factors */
00083     /* ---------------------------------------------------------------------- */
00084 
00085     if (Rs != NULL)
00086     {
00087   if (Numeric->Rs != NULL)
00088   {
00089       for (i = 0 ; i < n ; i++)
00090       {
00091     Rs [i] = Numeric->Rs [i] ;
00092       }
00093   }
00094   else
00095   {
00096       /* no scaling */
00097       for (i = 0 ; i < n ; i++)
00098       {
00099     Rs [i] = 1 ;
00100       }
00101   }
00102     }
00103 
00104     /* ---------------------------------------------------------------------- */
00105     /* extract block boundaries */
00106     /* ---------------------------------------------------------------------- */
00107 
00108     if (R != NULL)
00109     {
00110   for (block = 0 ; block <= nblocks ; block++)
00111   {
00112       R [block] = Symbolic->R [block] ;
00113   }
00114     }
00115 
00116     /* ---------------------------------------------------------------------- */
00117     /* extract final row permutation */
00118     /* ---------------------------------------------------------------------- */
00119 
00120     if (P != NULL)
00121     {
00122   for (k = 0 ; k < n ; k++)
00123   {
00124       P [k] = Numeric->Pnum [k] ;
00125   }
00126     }
00127 
00128     /* ---------------------------------------------------------------------- */
00129     /* extract column permutation */
00130     /* ---------------------------------------------------------------------- */
00131 
00132     if (Q != NULL)
00133     {
00134   for (k = 0 ; k < n ; k++)
00135   {
00136       Q [k] = Symbolic->Q [k] ;
00137   }
00138     }
00139 
00140     /* ---------------------------------------------------------------------- */
00141     /* extract each block of L */
00142     /* ---------------------------------------------------------------------- */
00143 
00144     if (Lp != NULL && Li != NULL && Lx != NULL
00145 #ifdef COMPLEX
00146   && Lz != NULL
00147 #endif
00148     )
00149     {
00150   nz = 0 ;
00151   for (block = 0 ; block < nblocks ; block++)
00152   {
00153       k1 = Symbolic->R [block] ;
00154       k2 = Symbolic->R [block+1] ;
00155       nk = k2 - k1 ;
00156       if (nk == 1)
00157       {
00158     /* singleton block */
00159     Lp [k1] = nz ;
00160     Li [nz] = k1 ;
00161     Lx [nz] = 1 ;
00162 #ifdef COMPLEX
00163     Lz [nz] = 0 ;
00164 #endif
00165     nz++ ;
00166       }
00167       else
00168       {
00169     /* non-singleton block */
00170     LU = Numeric->LUbx [block] ;
00171     Lip = Numeric->Lip + k1 ;
00172     Llen = Numeric->Llen + k1 ;
00173     for (kk = 0 ; kk < nk ; kk++)
00174     {
00175         Lp [k1+kk] = nz ;
00176         /* add the unit diagonal entry */
00177         Li [nz] = k1 + kk ;
00178         Lx [nz] = 1 ;
00179 #ifdef COMPLEX
00180         Lz [nz] = 0 ;
00181 #endif
00182         nz++ ;
00183         GET_POINTER (LU, Lip, Llen, Li2, Lx2, kk, len) ;
00184         for (p = 0 ; p < len ; p++)
00185         {
00186       Li [nz] = k1 + Li2 [p] ;
00187       Lx [nz] = REAL (Lx2 [p]) ;
00188 #ifdef COMPLEX
00189       Lz [nz] = IMAG (Lx2 [p]) ;
00190 #endif
00191       nz++ ;
00192         }
00193     }
00194       }
00195   }
00196   Lp [n] = nz ;
00197   ASSERT (nz == Numeric->lnz) ;
00198     }
00199 
00200     /* ---------------------------------------------------------------------- */
00201     /* extract each block of U */
00202     /* ---------------------------------------------------------------------- */
00203 
00204     if (Up != NULL && Ui != NULL && Ux != NULL
00205 #ifdef COMPLEX
00206   && Uz != NULL
00207 #endif
00208     )
00209     {
00210   nz = 0 ;
00211   for (block = 0 ; block < nblocks ; block++)
00212   {
00213       k1 = Symbolic->R [block] ;
00214       k2 = Symbolic->R [block+1] ;
00215       nk = k2 - k1 ;
00216       Ukk = ((Entry *) Numeric->Udiag) + k1 ;
00217       if (nk == 1)
00218       {
00219     /* singleton block */
00220     Up [k1] = nz ;
00221     Ui [nz] = k1 ;
00222     Ux [nz] = REAL (Ukk [0]) ;
00223 #ifdef COMPLEX
00224     Uz [nz] = IMAG (Ukk [0]) ;
00225 #endif
00226     nz++ ;
00227       }
00228       else
00229       {
00230     /* non-singleton block */
00231     LU = Numeric->LUbx [block] ;
00232     Uip = Numeric->Uip + k1 ;
00233     Ulen = Numeric->Ulen + k1 ;
00234     for (kk = 0 ; kk < nk ; kk++)
00235     {
00236         Up [k1+kk] = nz ;
00237         GET_POINTER (LU, Uip, Ulen, Ui2, Ux2, kk, len) ;
00238         for (p = 0 ; p < len ; p++)
00239         {
00240       Ui [nz] = k1 + Ui2 [p] ;
00241       Ux [nz] = REAL (Ux2 [p]) ;
00242 #ifdef COMPLEX
00243       Uz [nz] = IMAG (Ux2 [p]) ;
00244 #endif
00245       nz++ ;
00246         }
00247         /* add the diagonal entry */
00248         Ui [nz] = k1 + kk ;
00249         Ux [nz] = REAL (Ukk [kk]) ;
00250 #ifdef COMPLEX
00251         Uz [nz] = IMAG (Ukk [kk]) ;
00252 #endif
00253         nz++ ;
00254     }
00255       }
00256   }
00257   Up [n] = nz ;
00258   ASSERT (nz == Numeric->unz) ;
00259     }
00260 
00261     /* ---------------------------------------------------------------------- */
00262     /* extract the off-diagonal blocks, F */
00263     /* ---------------------------------------------------------------------- */
00264 
00265     if (Fp != NULL && Fi != NULL && Fx != NULL
00266 #ifdef COMPLEX
00267   && Fz != NULL
00268 #endif
00269     )
00270     {
00271   for (k = 0 ; k <= n ; k++)
00272   {
00273       Fp [k] = Numeric->Offp [k] ;
00274   }
00275   nz = Fp [n] ;
00276   for (k = 0 ; k < nz ; k++)
00277   {
00278       Fi [k] = Numeric->Offi [k] ;
00279   }
00280   for (k = 0 ; k < nz ; k++)
00281   {
00282       Fx [k] = REAL (((Entry *) Numeric->Offx) [k]) ;
00283 #ifdef COMPLEX
00284       Fz [k] = IMAG (((Entry *) Numeric->Offx) [k]) ;
00285 #endif
00286   }
00287     }
00288 
00289     return (TRUE) ;
00290 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines