amesos_cholmod_complex.c

Go to the documentation of this file.

/* ========================================================================== */
/* === Core/cholmod_complex ================================================= */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/Core Module.  Copyright (C) 2005-2006,
 * Univ. of Florida.  Author: Timothy A. Davis
 * The CHOLMOD/Core Module is licensed under Version 2.1 of the GNU
 * Lesser General Public License.  See lesser.txt for a text of the license.
 * CHOLMOD is also available under other licenses; contact authors for details.
 * http://www.cise.ufl.edu/research/sparse
 * -------------------------------------------------------------------------- */

/* If you convert a matrix that contains uninitialized data, valgrind will
 * complain.  This can occur in a factor L which has gaps (a partial
 * factorization, or after updates that change the nonzero pattern), an
 * unpacked sparse matrix, a dense matrix with leading dimension d > # of rows,
 * or any matrix (dense, sparse, triplet, or factor) with more space allocated
 * than is used.  You can safely ignore any of these complaints by valgrind. */

#include "amesos_cholmod_internal.h"
#include "amesos_cholmod_core.h"

/* ========================================================================== */
/* === cholmod_hypot ======================================================== */
/* ========================================================================== */

/* There is an equivalent routine called hypot in <math.h>, which conforms
 * to ANSI C99.  However, CHOLMOD does not assume that ANSI C99 is available.
 * You can use the ANSI C99 hypot routine with:
 *
 *  #include <math.h>
 *  Common->hypotenuse = hypot ;
 *
 * Default value of the Common->hypotenuse pointer is cholmod_hypot.
 *
 * s = hypot (x,y) computes s = sqrt (x*x + y*y) but does so more accurately.
 * The NaN cases for the double relops x >= y and x+y == x are safely ignored.
 */

double CHOLMOD(hypot) (double x, double y)
{
    double s, r ;
    x = fabs (x) ;
    y = fabs (y) ;
    if (x >= y)
    {
  if (x + y == x)
  {
      s = x ;
  }
  else
  {
      r = y / x ;
      s = x * sqrt (1.0 + r*r) ;
  }
    }
    else
    {
  if (y + x == y)
  {
      s = y ;
  }
  else
  {
      r = x / y ;
      s = y * sqrt (1.0 + r*r) ;
  }
    } 
    return (s) ;
}


/* ========================================================================== */
/* === cholmod_divcomplex =================================================== */
/* ========================================================================== */

/* c = a/b where c, a, and b are complex.  The real and imaginary parts are
 * passed as separate arguments to this routine.  The NaN case is ignored
 * for the double relop br >= bi.  Returns 1 if the denominator is zero,
 * 0 otherwise.  Note that this return value is the single exception to the
 * rule that all CHOLMOD routines that return int return TRUE if successful
 * or FALSE otherise.
 *
 * This uses ACM Algo 116, by R. L. Smith, 1962, which tries to avoid
 * underflow and overflow.
 *
 * c can be the same variable as a or b.
 *
 * Default value of the Common->complex_divide pointer is cholmod_divcomplex.
 */

int CHOLMOD(divcomplex)
(
    double ar, double ai, /* real and imaginary parts of a */
    double br, double bi, /* real and imaginary parts of b */
    double *cr, double *ci  /* real and imaginary parts of c */
)
{
    double tr, ti, r, den ;
    if (fabs (br) >= fabs (bi))
    {
  r = bi / br ;
  den = br + r * bi ;
  tr = (ar + ai * r) / den ;
  ti = (ai - ar * r) / den ;
    }
    else
    {
  r = br / bi ;
  den = r * br + bi ;
  tr = (ar * r + ai) / den ;
  ti = (ai * r - ar) / den ;
    }
    *cr = tr ;
    *ci = ti ;
    return (IS_ZERO (den)) ;
}


/* ========================================================================== */
/* === change_complexity ==================================================== */
/* ========================================================================== */

/* X and Z represent an array of size nz, with numeric xtype given by xtype_in.
 *
 * If xtype_in is:
 * CHOLMOD_PATTERN: X and Z must be NULL.
 * CHOLMOD_REAL:    X is of size nz, Z must be NULL.
 * CHOLMOD_COMPLEX: X is of size 2*nz, Z must be NULL.
 * CHOLMOD_ZOMPLEX: X is of size nz, Z is of size nz.
 *
 * The array is changed into the numeric xtype given by xtype_out, with the
 * same definitions of X and Z above.  Note that the input conditions, above,
 * are not checked.  These are checked in the caller routine.
 *
 * Returns TRUE if successful, FALSE otherwise.  X and Z are not modified if
 * not successful.
 */

static int change_complexity
(
    /* ---- input ---- */
    Int nz,   /* size of X and/or Z */
    int xtype_in, /* xtype of X and Z on input */
    int xtype_out,  /* requested xtype of X and Z on output */
    int xtype1,   /* xtype_out must be in the range [xtype1 .. xtype2] */
    int xtype2,
    /* ---- in/out --- */
    void **XX,    /* old X on input, new X on output */
    void **ZZ,    /* old Z on input, new Z on output */
    /* --------------- */
    cholmod_common *Common
)
{
    double *Xold, *Zold, *Xnew, *Znew ;
    Int k ;
    size_t nz2 ;

    if (xtype_out < xtype1 || xtype_out > xtype2)
    {
  ERROR (CHOLMOD_INVALID, "invalid xtype") ;
  return (FALSE) ;
    }

    Common->status = CHOLMOD_OK ;
    Xold = *XX ;
    Zold = *ZZ ;

    switch (xtype_in)
    {

  /* ------------------------------------------------------------------ */
  /* converting from pattern */
  /* ------------------------------------------------------------------ */

  case CHOLMOD_PATTERN:

      switch (xtype_out)
      {

    /* ---------------------------------------------------------- */
    /* pattern -> real */
    /* ---------------------------------------------------------- */

    case CHOLMOD_REAL:
        /* allocate X and set to all ones */
        Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
        if (Common->status < CHOLMOD_OK)
        {
      return (FALSE) ;
        }
        for (k = 0 ; k < nz ; k++)
        {
      Xnew [k] = 1 ;
        }
        *XX = Xnew ;
        break ;

    /* ---------------------------------------------------------- */
    /* pattern -> complex */
    /* ---------------------------------------------------------- */

    case CHOLMOD_COMPLEX:
        /* allocate X and set to all ones */
        Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
        if (Common->status < CHOLMOD_OK)
        {
      return (FALSE) ;
        }
        for (k = 0 ; k < nz ; k++)
        {
      Xnew [2*k  ] = 1 ;
      Xnew [2*k+1] = 0 ;
        }
        *XX = Xnew ;
        break ;

    /* ---------------------------------------------------------- */
    /* pattern -> zomplex */
    /* ---------------------------------------------------------- */

    case CHOLMOD_ZOMPLEX:
        /* allocate X and Z and set to all ones */
        Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
        Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
        if (Common->status < CHOLMOD_OK)
        {
      CHOLMOD(free) (nz, sizeof (double), Xnew, Common) ;
      CHOLMOD(free) (nz, sizeof (double), Znew, Common) ;
      return (FALSE) ;
        }
        for (k = 0 ; k < nz ; k++)
        {
      Xnew [k] = 1 ;
      Znew [k] = 0 ;
        }
        *XX = Xnew ;
        *ZZ = Znew ;
        break ;
      }
      break ;

  /* ------------------------------------------------------------------ */
  /* converting from real */
  /* ------------------------------------------------------------------ */

  case CHOLMOD_REAL:

      switch (xtype_out)
      {

    /* ---------------------------------------------------------- */
    /* real -> pattern */
    /* ---------------------------------------------------------- */

    case CHOLMOD_PATTERN:
        /* free X */
        *XX = CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
        break ;

    /* ---------------------------------------------------------- */
    /* real -> complex */
    /* ---------------------------------------------------------- */

    case CHOLMOD_COMPLEX:
        /* allocate a new X and copy the old X */
        Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
        if (Common->status < CHOLMOD_OK)
        {
      return (FALSE) ;
        }
        for (k = 0 ; k < nz ; k++)
        {
      Xnew [2*k  ] = Xold [k] ;
      Xnew [2*k+1] = 0 ;
        }
        CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
        *XX = Xnew ;
        break ;

    /* ---------------------------------------------------------- */
    /* real -> zomplex */
    /* ---------------------------------------------------------- */

    case CHOLMOD_ZOMPLEX:
        /* allocate a new Z and set it to zero */
        Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
        if (Common->status < CHOLMOD_OK)
        {
      return (FALSE) ;
        }
        for (k = 0 ; k < nz ; k++)
        {
      Znew [k] = 0 ;
        }
        *ZZ = Znew ;
        break ;
      }
      break ;

  /* ------------------------------------------------------------------ */
  /* converting from complex */
  /* ------------------------------------------------------------------ */

  case CHOLMOD_COMPLEX:

      switch (xtype_out)
      {

    /* ---------------------------------------------------------- */
    /* complex -> pattern */
    /* ---------------------------------------------------------- */

    case CHOLMOD_PATTERN:
        /* free X */
        *XX = CHOLMOD(free) (nz, 2*sizeof (double), *XX, Common) ;
        break ;

    /* ---------------------------------------------------------- */
    /* complex -> real */
    /* ---------------------------------------------------------- */

    case CHOLMOD_REAL:
        /* pack the real part of X, discarding the imaginary part */
        for (k = 0 ; k < nz ; k++)
        {
      Xold [k] = Xold [2*k] ;
        }
        /* shrink X in half (this cannot fail) */
        nz2 = 2*nz ;
        *XX = CHOLMOD(realloc) (nz, sizeof (double), *XX, &nz2,
          Common) ;
        break ;

    /* ---------------------------------------------------------- */
    /* complex -> zomplex */
    /* ---------------------------------------------------------- */

    case CHOLMOD_ZOMPLEX:
        /* allocate X and Z and copy the old X into them */
        Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
        Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
        if (Common->status < CHOLMOD_OK)
        {
      CHOLMOD(free) (nz, sizeof (double), Xnew, Common) ;
      CHOLMOD(free) (nz, sizeof (double), Znew, Common) ;
      return (FALSE) ;
        }
        for (k = 0 ; k < nz ; k++)
        {
      Xnew [k] = Xold [2*k  ] ;
      Znew [k] = Xold [2*k+1] ;
        }
        CHOLMOD(free) (nz, 2*sizeof (double), *XX, Common) ;
        *XX = Xnew ;
        *ZZ = Znew ;
        break ;
      }
      break ;

  /* ------------------------------------------------------------------ */
  /* converting from zomplex */
  /* ------------------------------------------------------------------ */

  case CHOLMOD_ZOMPLEX:

      switch (xtype_out)
      {

    /* ---------------------------------------------------------- */
    /* zomplex -> pattern */
    /* ---------------------------------------------------------- */

    case CHOLMOD_PATTERN:
        /* free X and Z */
        *XX = CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
        *ZZ = CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
        break ;

    /* ---------------------------------------------------------- */
    /* zomplex -> real */
    /* ---------------------------------------------------------- */

    case CHOLMOD_REAL:
        /* free the imaginary part */
        *ZZ = CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
        break ;

    /* ---------------------------------------------------------- */
    /* zomplex -> complex */
    /* ---------------------------------------------------------- */

    case CHOLMOD_COMPLEX:
        Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
        if (Common->status < CHOLMOD_OK)
        {
      return (FALSE) ;
        }
        for (k = 0 ; k < nz ; k++)
        {
      Xnew [2*k  ] = Xold [k] ;
      Xnew [2*k+1] = Zold [k] ;
        }
        CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
        CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
        *XX = Xnew ;
        *ZZ = NULL ;
        break ;

      }
      break ;
    }

    return (TRUE) ;
}


/* ========================================================================== */
/* === cholmod_sparse_xtype ================================================= */
/* ========================================================================== */

/* Change the numeric xtype of a sparse matrix.  Supports any type on input
 * and output (pattern, real, complex, or zomplex). */

int CHOLMOD(sparse_xtype)
(
    /* ---- input ---- */
    int to_xtype, /* requested xtype */
    /* ---- in/out --- */
    cholmod_sparse *A,  /* sparse matrix to change */
    /* --------------- */
    cholmod_common *Common
)
{
    Int ok ;
    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;

    ok = change_complexity (A->nzmax, A->xtype, to_xtype,
      CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, &(A->x), &(A->z), Common) ;
    if (ok)
    {
  A->xtype = to_xtype ;
    }
    return (ok) ;
}


/* ========================================================================== */
/* === cholmod_triplet_xtype ================================================ */
/* ========================================================================== */

/* Change the numeric xtype of a triplet matrix.  Supports any type on input
 * and output (pattern, real, complex, or zomplex). */

int CHOLMOD(triplet_xtype)
(
    /* ---- input ---- */
    int to_xtype, /* requested xtype */
    /* ---- in/out --- */
    cholmod_triplet *T, /* triplet matrix to change */
    /* --------------- */
    cholmod_common *Common
)
{
    Int ok ;
    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (T, FALSE) ;
    RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    ok = change_complexity (T->nzmax, T->xtype, to_xtype,
      CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, &(T->x), &(T->z), Common) ;
    if (ok)
    {
  T->xtype = to_xtype ;
    }
    return (ok) ;
}


/* ========================================================================== */
/* === cholmod_dense_xtype ================================================= */
/* ========================================================================== */

/* Change the numeric xtype of a dense matrix.  Supports real, complex or
 * zomplex on input and output */

int CHOLMOD(dense_xtype)
(
    /* ---- input ---- */
    int to_xtype, /* requested xtype */
    /* ---- in/out --- */
    cholmod_dense *X, /* dense matrix to change */
    /* --------------- */
    cholmod_common *Common
)
{
    Int ok ;
    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (X, FALSE) ;
    RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
    ok = change_complexity (X->nzmax, X->xtype, to_xtype,
      CHOLMOD_REAL, CHOLMOD_ZOMPLEX, &(X->x), &(X->z), Common) ;
    if (ok)
    {
  X->xtype = to_xtype ;
    }
    return (ok) ;
}


/* ========================================================================== */
/* === cholmod_factor_xtype ================================================= */
/* ========================================================================== */

/* Change the numeric xtype of a factor.  Supports real, complex or zomplex on
 * input and output.   Supernodal zomplex factors are not supported. */

int CHOLMOD(factor_xtype)
(
    /* ---- input ---- */
    int to_xtype, /* requested xtype */
    /* ---- in/out --- */
    cholmod_factor *L,  /* factor to change */
    /* --------------- */
    cholmod_common *Common
)
{
    Int ok ;
    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (L, FALSE) ;
    RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
    if (L->is_super &&
      (L->xtype == CHOLMOD_ZOMPLEX || to_xtype == CHOLMOD_ZOMPLEX))
    {
  ERROR (CHOLMOD_INVALID, "invalid xtype for supernodal L") ;
  return (FALSE) ;
    }
    ok = change_complexity ((L->is_super ? L->xsize : L->nzmax), L->xtype,
      to_xtype, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, &(L->x), &(L->z), Common) ;
    if (ok)
    {
  L->xtype = to_xtype ;
    }
    return (ok) ;
}