Mercurial > octave-nkf
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 + + } +}