Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_l_complex.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Core/cholmod_complex ================================================= */
00003 /* ========================================================================== */
00004 
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Core Module.  Copyright (C) 2005-2006,
00007  * Univ. of Florida.  Author: Timothy A. Davis
00008  * The CHOLMOD/Core Module is licensed under Version 2.1 of the GNU
00009  * Lesser General Public License.  See lesser.txt for a text of the license.
00010  * CHOLMOD is also available under other licenses; contact authors for details.
00011  * http://www.cise.ufl.edu/research/sparse
00012  * -------------------------------------------------------------------------- */
00013 
00014 /* If you convert a matrix that contains uninitialized data, valgrind will
00015  * complain.  This can occur in a factor L which has gaps (a partial
00016  * factorization, or after updates that change the nonzero pattern), an
00017  * unpacked sparse matrix, a dense matrix with leading dimension d > # of rows,
00018  * or any matrix (dense, sparse, triplet, or factor) with more space allocated
00019  * than is used.  You can safely ignore any of these complaints by valgrind. */
00020 
00021 /* This file should make the long int version of CHOLMOD */
00022 #define DLONG 1
00023 
00024 #include "amesos_cholmod_internal.h"
00025 #include "amesos_cholmod_core.h"
00026 
00027 /* ========================================================================== */
00028 /* === cholmod_hypot ======================================================== */
00029 /* ========================================================================== */
00030 
00031 /* There is an equivalent routine called hypot in <math.h>, which conforms
00032  * to ANSI C99.  However, CHOLMOD does not assume that ANSI C99 is available.
00033  * You can use the ANSI C99 hypot routine with:
00034  *
00035  *  #include <math.h>
00036  *  Common->hypotenuse = hypot ;
00037  *
00038  * Default value of the Common->hypotenuse pointer is cholmod_hypot.
00039  *
00040  * s = hypot (x,y) computes s = sqrt (x*x + y*y) but does so more accurately.
00041  * The NaN cases for the double relops x >= y and x+y == x are safely ignored.
00042  */
00043 
00044 double CHOLMOD(hypot) (double x, double y)
00045 {
00046     double s, r ;
00047     x = fabs (x) ;
00048     y = fabs (y) ;
00049     if (x >= y)
00050     {
00051   if (x + y == x)
00052   {
00053       s = x ;
00054   }
00055   else
00056   {
00057       r = y / x ;
00058       s = x * sqrt (1.0 + r*r) ;
00059   }
00060     }
00061     else
00062     {
00063   if (y + x == y)
00064   {
00065       s = y ;
00066   }
00067   else
00068   {
00069       r = x / y ;
00070       s = y * sqrt (1.0 + r*r) ;
00071   }
00072     } 
00073     return (s) ;
00074 }
00075 
00076 
00077 /* ========================================================================== */
00078 /* === cholmod_divcomplex =================================================== */
00079 /* ========================================================================== */
00080 
00081 /* c = a/b where c, a, and b are complex.  The real and imaginary parts are
00082  * passed as separate arguments to this routine.  The NaN case is ignored
00083  * for the double relop br >= bi.  Returns 1 if the denominator is zero,
00084  * 0 otherwise.  Note that this return value is the single exception to the
00085  * rule that all CHOLMOD routines that return int return TRUE if successful
00086  * or FALSE otherise.
00087  *
00088  * This uses ACM Algo 116, by R. L. Smith, 1962, which tries to avoid
00089  * underflow and overflow.
00090  *
00091  * c can be the same variable as a or b.
00092  *
00093  * Default value of the Common->complex_divide pointer is cholmod_divcomplex.
00094  */
00095 
00096 int CHOLMOD(divcomplex)
00097 (
00098     double ar, double ai, /* real and imaginary parts of a */
00099     double br, double bi, /* real and imaginary parts of b */
00100     double *cr, double *ci  /* real and imaginary parts of c */
00101 )
00102 {
00103     double tr, ti, r, den ;
00104     if (fabs (br) >= fabs (bi))
00105     {
00106   r = bi / br ;
00107   den = br + r * bi ;
00108   tr = (ar + ai * r) / den ;
00109   ti = (ai - ar * r) / den ;
00110     }
00111     else
00112     {
00113   r = br / bi ;
00114   den = r * br + bi ;
00115   tr = (ar * r + ai) / den ;
00116   ti = (ai * r - ar) / den ;
00117     }
00118     *cr = tr ;
00119     *ci = ti ;
00120     return (IS_ZERO (den)) ;
00121 }
00122 
00123 
00124 /* ========================================================================== */
00125 /* === change_complexity ==================================================== */
00126 /* ========================================================================== */
00127 
00128 /* X and Z represent an array of size nz, with numeric xtype given by xtype_in.
00129  *
00130  * If xtype_in is:
00131  * CHOLMOD_PATTERN: X and Z must be NULL.
00132  * CHOLMOD_REAL:    X is of size nz, Z must be NULL.
00133  * CHOLMOD_COMPLEX: X is of size 2*nz, Z must be NULL.
00134  * CHOLMOD_ZOMPLEX: X is of size nz, Z is of size nz.
00135  *
00136  * The array is changed into the numeric xtype given by xtype_out, with the
00137  * same definitions of X and Z above.  Note that the input conditions, above,
00138  * are not checked.  These are checked in the caller routine.
00139  *
00140  * Returns TRUE if successful, FALSE otherwise.  X and Z are not modified if
00141  * not successful.
00142  */
00143 
00144 static int change_complexity
00145 (
00146     /* ---- input ---- */
00147     Int nz,   /* size of X and/or Z */
00148     int xtype_in, /* xtype of X and Z on input */
00149     int xtype_out,  /* requested xtype of X and Z on output */
00150     int xtype1,   /* xtype_out must be in the range [xtype1 .. xtype2] */
00151     int xtype2,
00152     /* ---- in/out --- */
00153     void **XX,    /* old X on input, new X on output */
00154     void **ZZ,    /* old Z on input, new Z on output */
00155     /* --------------- */
00156     cholmod_common *Common
00157 )
00158 {
00159     double *Xold, *Zold, *Xnew, *Znew ;
00160     Int k ;
00161     size_t nz2 ;
00162 
00163     if (xtype_out < xtype1 || xtype_out > xtype2)
00164     {
00165   ERROR (CHOLMOD_INVALID, "invalid xtype") ;
00166   return (FALSE) ;
00167     }
00168 
00169     Common->status = CHOLMOD_OK ;
00170     Xold = *XX ;
00171     Zold = *ZZ ;
00172 
00173     switch (xtype_in)
00174     {
00175 
00176   /* ------------------------------------------------------------------ */
00177   /* converting from pattern */
00178   /* ------------------------------------------------------------------ */
00179 
00180   case CHOLMOD_PATTERN:
00181 
00182       switch (xtype_out)
00183       {
00184 
00185     /* ---------------------------------------------------------- */
00186     /* pattern -> real */
00187     /* ---------------------------------------------------------- */
00188 
00189     case CHOLMOD_REAL:
00190         /* allocate X and set to all ones */
00191         Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
00192         if (Common->status < CHOLMOD_OK)
00193         {
00194       return (FALSE) ;
00195         }
00196         for (k = 0 ; k < nz ; k++)
00197         {
00198       Xnew [k] = 1 ;
00199         }
00200         *XX = Xnew ;
00201         break ;
00202 
00203     /* ---------------------------------------------------------- */
00204     /* pattern -> complex */
00205     /* ---------------------------------------------------------- */
00206 
00207     case CHOLMOD_COMPLEX:
00208         /* allocate X and set to all ones */
00209         Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
00210         if (Common->status < CHOLMOD_OK)
00211         {
00212       return (FALSE) ;
00213         }
00214         for (k = 0 ; k < nz ; k++)
00215         {
00216       Xnew [2*k  ] = 1 ;
00217       Xnew [2*k+1] = 0 ;
00218         }
00219         *XX = Xnew ;
00220         break ;
00221 
00222     /* ---------------------------------------------------------- */
00223     /* pattern -> zomplex */
00224     /* ---------------------------------------------------------- */
00225 
00226     case CHOLMOD_ZOMPLEX:
00227         /* allocate X and Z and set to all ones */
00228         Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
00229         Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
00230         if (Common->status < CHOLMOD_OK)
00231         {
00232       CHOLMOD(free) (nz, sizeof (double), Xnew, Common) ;
00233       CHOLMOD(free) (nz, sizeof (double), Znew, Common) ;
00234       return (FALSE) ;
00235         }
00236         for (k = 0 ; k < nz ; k++)
00237         {
00238       Xnew [k] = 1 ;
00239       Znew [k] = 0 ;
00240         }
00241         *XX = Xnew ;
00242         *ZZ = Znew ;
00243         break ;
00244       }
00245       break ;
00246 
00247   /* ------------------------------------------------------------------ */
00248   /* converting from real */
00249   /* ------------------------------------------------------------------ */
00250 
00251   case CHOLMOD_REAL:
00252 
00253       switch (xtype_out)
00254       {
00255 
00256     /* ---------------------------------------------------------- */
00257     /* real -> pattern */
00258     /* ---------------------------------------------------------- */
00259 
00260     case CHOLMOD_PATTERN:
00261         /* free X */
00262         *XX = CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
00263         break ;
00264 
00265     /* ---------------------------------------------------------- */
00266     /* real -> complex */
00267     /* ---------------------------------------------------------- */
00268 
00269     case CHOLMOD_COMPLEX:
00270         /* allocate a new X and copy the old X */
00271         Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
00272         if (Common->status < CHOLMOD_OK)
00273         {
00274       return (FALSE) ;
00275         }
00276         for (k = 0 ; k < nz ; k++)
00277         {
00278       Xnew [2*k  ] = Xold [k] ;
00279       Xnew [2*k+1] = 0 ;
00280         }
00281         CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
00282         *XX = Xnew ;
00283         break ;
00284 
00285     /* ---------------------------------------------------------- */
00286     /* real -> zomplex */
00287     /* ---------------------------------------------------------- */
00288 
00289     case CHOLMOD_ZOMPLEX:
00290         /* allocate a new Z and set it to zero */
00291         Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
00292         if (Common->status < CHOLMOD_OK)
00293         {
00294       return (FALSE) ;
00295         }
00296         for (k = 0 ; k < nz ; k++)
00297         {
00298       Znew [k] = 0 ;
00299         }
00300         *ZZ = Znew ;
00301         break ;
00302       }
00303       break ;
00304 
00305   /* ------------------------------------------------------------------ */
00306   /* converting from complex */
00307   /* ------------------------------------------------------------------ */
00308 
00309   case CHOLMOD_COMPLEX:
00310 
00311       switch (xtype_out)
00312       {
00313 
00314     /* ---------------------------------------------------------- */
00315     /* complex -> pattern */
00316     /* ---------------------------------------------------------- */
00317 
00318     case CHOLMOD_PATTERN:
00319         /* free X */
00320         *XX = CHOLMOD(free) (nz, 2*sizeof (double), *XX, Common) ;
00321         break ;
00322 
00323     /* ---------------------------------------------------------- */
00324     /* complex -> real */
00325     /* ---------------------------------------------------------- */
00326 
00327     case CHOLMOD_REAL:
00328         /* pack the real part of X, discarding the imaginary part */
00329         for (k = 0 ; k < nz ; k++)
00330         {
00331       Xold [k] = Xold [2*k] ;
00332         }
00333         /* shrink X in half (this cannot fail) */
00334         nz2 = 2*nz ;
00335         *XX = CHOLMOD(realloc) (nz, sizeof (double), *XX, &nz2,
00336           Common) ;
00337         break ;
00338 
00339     /* ---------------------------------------------------------- */
00340     /* complex -> zomplex */
00341     /* ---------------------------------------------------------- */
00342 
00343     case CHOLMOD_ZOMPLEX:
00344         /* allocate X and Z and copy the old X into them */
00345         Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
00346         Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
00347         if (Common->status < CHOLMOD_OK)
00348         {
00349       CHOLMOD(free) (nz, sizeof (double), Xnew, Common) ;
00350       CHOLMOD(free) (nz, sizeof (double), Znew, Common) ;
00351       return (FALSE) ;
00352         }
00353         for (k = 0 ; k < nz ; k++)
00354         {
00355       Xnew [k] = Xold [2*k  ] ;
00356       Znew [k] = Xold [2*k+1] ;
00357         }
00358         CHOLMOD(free) (nz, 2*sizeof (double), *XX, Common) ;
00359         *XX = Xnew ;
00360         *ZZ = Znew ;
00361         break ;
00362       }
00363       break ;
00364 
00365   /* ------------------------------------------------------------------ */
00366   /* converting from zomplex */
00367   /* ------------------------------------------------------------------ */
00368 
00369   case CHOLMOD_ZOMPLEX:
00370 
00371       switch (xtype_out)
00372       {
00373 
00374     /* ---------------------------------------------------------- */
00375     /* zomplex -> pattern */
00376     /* ---------------------------------------------------------- */
00377 
00378     case CHOLMOD_PATTERN:
00379         /* free X and Z */
00380         *XX = CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
00381         *ZZ = CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
00382         break ;
00383 
00384     /* ---------------------------------------------------------- */
00385     /* zomplex -> real */
00386     /* ---------------------------------------------------------- */
00387 
00388     case CHOLMOD_REAL:
00389         /* free the imaginary part */
00390         *ZZ = CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
00391         break ;
00392 
00393     /* ---------------------------------------------------------- */
00394     /* zomplex -> complex */
00395     /* ---------------------------------------------------------- */
00396 
00397     case CHOLMOD_COMPLEX:
00398         Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
00399         if (Common->status < CHOLMOD_OK)
00400         {
00401       return (FALSE) ;
00402         }
00403         for (k = 0 ; k < nz ; k++)
00404         {
00405       Xnew [2*k  ] = Xold [k] ;
00406       Xnew [2*k+1] = Zold [k] ;
00407         }
00408         CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
00409         CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
00410         *XX = Xnew ;
00411         *ZZ = NULL ;
00412         break ;
00413 
00414       }
00415       break ;
00416     }
00417 
00418     return (TRUE) ;
00419 }
00420 
00421 
00422 /* ========================================================================== */
00423 /* === cholmod_sparse_xtype ================================================= */
00424 /* ========================================================================== */
00425 
00426 /* Change the numeric xtype of a sparse matrix.  Supports any type on input
00427  * and output (pattern, real, complex, or zomplex). */
00428 
00429 int CHOLMOD(sparse_xtype)
00430 (
00431     /* ---- input ---- */
00432     int to_xtype, /* requested xtype */
00433     /* ---- in/out --- */
00434     cholmod_sparse *A,  /* sparse matrix to change */
00435     /* --------------- */
00436     cholmod_common *Common
00437 )
00438 {
00439     Int ok ;
00440     RETURN_IF_NULL_COMMON (FALSE) ;
00441     RETURN_IF_NULL (A, FALSE) ;
00442     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00443 
00444     ok = change_complexity (A->nzmax, A->xtype, to_xtype,
00445       CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, &(A->x), &(A->z), Common) ;
00446     if (ok)
00447     {
00448   A->xtype = to_xtype ;
00449     }
00450     return (ok) ;
00451 }
00452 
00453 
00454 /* ========================================================================== */
00455 /* === cholmod_triplet_xtype ================================================ */
00456 /* ========================================================================== */
00457 
00458 /* Change the numeric xtype of a triplet matrix.  Supports any type on input
00459  * and output (pattern, real, complex, or zomplex). */
00460 
00461 int CHOLMOD(triplet_xtype)
00462 (
00463     /* ---- input ---- */
00464     int to_xtype, /* requested xtype */
00465     /* ---- in/out --- */
00466     cholmod_triplet *T, /* triplet matrix to change */
00467     /* --------------- */
00468     cholmod_common *Common
00469 )
00470 {
00471     Int ok ;
00472     RETURN_IF_NULL_COMMON (FALSE) ;
00473     RETURN_IF_NULL (T, FALSE) ;
00474     RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00475     ok = change_complexity (T->nzmax, T->xtype, to_xtype,
00476       CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, &(T->x), &(T->z), Common) ;
00477     if (ok)
00478     {
00479   T->xtype = to_xtype ;
00480     }
00481     return (ok) ;
00482 }
00483 
00484 
00485 /* ========================================================================== */
00486 /* === cholmod_dense_xtype ================================================= */
00487 /* ========================================================================== */
00488 
00489 /* Change the numeric xtype of a dense matrix.  Supports real, complex or
00490  * zomplex on input and output */
00491 
00492 int CHOLMOD(dense_xtype)
00493 (
00494     /* ---- input ---- */
00495     int to_xtype, /* requested xtype */
00496     /* ---- in/out --- */
00497     cholmod_dense *X, /* dense matrix to change */
00498     /* --------------- */
00499     cholmod_common *Common
00500 )
00501 {
00502     Int ok ;
00503     RETURN_IF_NULL_COMMON (FALSE) ;
00504     RETURN_IF_NULL (X, FALSE) ;
00505     RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
00506     ok = change_complexity (X->nzmax, X->xtype, to_xtype,
00507       CHOLMOD_REAL, CHOLMOD_ZOMPLEX, &(X->x), &(X->z), Common) ;
00508     if (ok)
00509     {
00510   X->xtype = to_xtype ;
00511     }
00512     return (ok) ;
00513 }
00514 
00515 
00516 /* ========================================================================== */
00517 /* === cholmod_factor_xtype ================================================= */
00518 /* ========================================================================== */
00519 
00520 /* Change the numeric xtype of a factor.  Supports real, complex or zomplex on
00521  * input and output.   Supernodal zomplex factors are not supported. */
00522 
00523 int CHOLMOD(factor_xtype)
00524 (
00525     /* ---- input ---- */
00526     int to_xtype, /* requested xtype */
00527     /* ---- in/out --- */
00528     cholmod_factor *L,  /* factor to change */
00529     /* --------------- */
00530     cholmod_common *Common
00531 )
00532 {
00533     Int ok ;
00534     RETURN_IF_NULL_COMMON (FALSE) ;
00535     RETURN_IF_NULL (L, FALSE) ;
00536     RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
00537     if (L->is_super &&
00538       (L->xtype == CHOLMOD_ZOMPLEX || to_xtype == CHOLMOD_ZOMPLEX))
00539     {
00540   ERROR (CHOLMOD_INVALID, "invalid xtype for supernodal L") ;
00541   return (FALSE) ;
00542     }
00543     ok = change_complexity ((L->is_super ? L->xsize : L->nzmax), L->xtype,
00544       to_xtype, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, &(L->x), &(L->z), Common) ;
00545     if (ok)
00546     {
00547   L->xtype = to_xtype ;
00548     }
00549     return (ok) ;
00550 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines