diff liboctave/UMFPACK/UMFPACK/Source/umf_scale.c @ 5164:57077d0ddc8e

[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/UMFPACK/UMFPACK/Source/umf_scale.c	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,127 @@
+/* ========================================================================== */
+/* === UMF_scale ============================================================ */
+/* ========================================================================== */
+
+/* -------------------------------------------------------------------------- */
+/* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis.  CISE Dept,   */
+/* Univ. of Florida.  All Rights Reserved.  See ../Doc/License for License.   */
+/* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
+/* -------------------------------------------------------------------------- */
+
+/* Divide a vector of stride 1 by the pivot value. */
+
+#include "umf_internal.h"
+
+GLOBAL void UMF_scale
+(
+    Int n,
+    Entry pivot,
+    Entry X [ ]
+)
+{
+    Entry x ;
+    double s ;
+    Int i ;
+
+    /* ---------------------------------------------------------------------- */
+    /* compute the approximate absolute value of the pivot, and select method */
+    /* ---------------------------------------------------------------------- */
+
+    APPROX_ABS (s, pivot) ;
+
+    if (s < RECIPROCAL_TOLERANCE || IS_NAN (pivot))
+    {
+	/* ------------------------------------------------------------------ */
+	/* tiny, or zero, pivot case */
+	/* ------------------------------------------------------------------ */
+
+	/* The pivot is tiny, or NaN.  Do not divide zero by the pivot value,
+	 * and do not multiply by 1/pivot, either. */
+
+	for (i = 0 ; i < n ; i++)
+	{
+	    /* X [i] /= pivot ; */
+	    x = X [i] ;
+
+#ifndef NO_DIVIDE_BY_ZERO
+	    if (IS_NONZERO (x))
+	    {
+		DIV (X [i], x, pivot) ;
+	    }
+#else
+	    /* Do not divide by zero */
+	    if (IS_NONZERO (x) && IS_NONZERO (pivot))
+	    {
+		DIV (X [i], x, pivot) ;
+	    }
+#endif
+
+	}
+
+    }
+    else
+    {
+
+	/* ------------------------------------------------------------------ */
+	/* normal case.  select the x/pivot or x * (1/pivot) method */
+	/* ------------------------------------------------------------------ */
+
+	/* The pivot is not tiny, and is not NaN.   Don't bother to check for
+	 * zeros in the pivot column, X. */
+
+#if !defined (NRECIPROCAL) && !(defined (__GNUC__) && defined (COMPLEX))
+
+	    /* -------------------------------------------------------------- */
+	    /* multiply x by (1/pivot) */
+	    /* -------------------------------------------------------------- */
+
+	    /* Slightly less accurate, but faster.  It allows the use of
+	     * the level-1 BLAS dscal or zscal routine.  This not used when
+	     * UMFPACK is used in MATLAB (either as a built-in routine, or as
+	     * a mexFunction).
+	     *
+	     * Using gcc version 3.2 can cause the following code to fail for
+	     * some complex matrices (not all), with or without the BLAS.  This
+	     * was found in Red Hat Linux 7.3 on a Dell Latitude C840 with a
+	     * Pentium 4M.  Thus, this code is not used when gcc is used, for
+	     * the complex case.
+	     *
+	     * It works just fine with Intel's icc compiler, version 7.0.
+	     */
+
+	    /* pivot = 1 / pivot */
+	    RECIPROCAL (pivot) ;
+
+#if defined (USE_NO_BLAS)
+	    for (i = 0 ; i < n ; i++)
+	    {
+		/* X [i] *= pivot ; */
+		x = X [i] ;
+		MULT (X [i], x, pivot) ;
+	    }
+#else
+	    BLAS_SCAL (n, pivot, X) ;
+#endif
+
+#else
+
+	    /* -------------------------------------------------------------- */
+	    /* divide x by the pivot */
+	    /* -------------------------------------------------------------- */
+
+	    /* This is slightly more accurate, particularly if the pivot column
+	     * consists of only IEEE subnormals.  Always do this if UMFPACK is
+	     * being compiled as a built-in routine or mexFunction in MATLAB,
+	     * or if gcc is being used with complex matrices. */
+
+	    for (i = 0 ; i < n ; i++)
+	    {
+		/* X [i] /= pivot ; */
+		x = X [i] ;
+		DIV (X [i], x, pivot) ;
+	    }
+
+#endif
+
+    }
+}