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