annotate 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
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
1 /* ========================================================================== */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
2 /* === UMF_scale ============================================================ */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
3 /* ========================================================================== */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
4
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
5 /* -------------------------------------------------------------------------- */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
6 /* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE Dept, */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
9 /* -------------------------------------------------------------------------- */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
10
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
11 /* Divide a vector of stride 1 by the pivot value. */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
12
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
13 #include "umf_internal.h"
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
14
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
15 GLOBAL void UMF_scale
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
16 (
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
17 Int n,
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
18 Entry pivot,
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
19 Entry X [ ]
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
20 )
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
21 {
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
22 Entry x ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
23 double s ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
24 Int i ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
25
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
26 /* ---------------------------------------------------------------------- */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
27 /* compute the approximate absolute value of the pivot, and select method */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
28 /* ---------------------------------------------------------------------- */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
29
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
30 APPROX_ABS (s, pivot) ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
31
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
32 if (s < RECIPROCAL_TOLERANCE || IS_NAN (pivot))
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
33 {
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
34 /* ------------------------------------------------------------------ */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
35 /* tiny, or zero, pivot case */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
36 /* ------------------------------------------------------------------ */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
37
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
38 /* The pivot is tiny, or NaN. Do not divide zero by the pivot value,
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
39 * and do not multiply by 1/pivot, either. */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
40
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
41 for (i = 0 ; i < n ; i++)
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
42 {
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
43 /* X [i] /= pivot ; */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
44 x = X [i] ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
45
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
46 #ifndef NO_DIVIDE_BY_ZERO
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
47 if (IS_NONZERO (x))
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
48 {
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
49 DIV (X [i], x, pivot) ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
50 }
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
51 #else
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
52 /* Do not divide by zero */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
53 if (IS_NONZERO (x) && IS_NONZERO (pivot))
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
54 {
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
55 DIV (X [i], x, pivot) ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
56 }
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
57 #endif
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
58
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
59 }
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
60
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
61 }
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
62 else
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
63 {
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
64
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
65 /* ------------------------------------------------------------------ */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
66 /* normal case. select the x/pivot or x * (1/pivot) method */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
67 /* ------------------------------------------------------------------ */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
68
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
69 /* The pivot is not tiny, and is not NaN. Don't bother to check for
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
70 * zeros in the pivot column, X. */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
71
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
72 #if !defined (NRECIPROCAL) && !(defined (__GNUC__) && defined (COMPLEX))
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
73
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
74 /* -------------------------------------------------------------- */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
75 /* multiply x by (1/pivot) */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
76 /* -------------------------------------------------------------- */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
77
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
78 /* Slightly less accurate, but faster. It allows the use of
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
79 * the level-1 BLAS dscal or zscal routine. This not used when
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
80 * UMFPACK is used in MATLAB (either as a built-in routine, or as
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
81 * a mexFunction).
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
82 *
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
83 * Using gcc version 3.2 can cause the following code to fail for
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
84 * some complex matrices (not all), with or without the BLAS. This
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
85 * was found in Red Hat Linux 7.3 on a Dell Latitude C840 with a
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
86 * Pentium 4M. Thus, this code is not used when gcc is used, for
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
87 * the complex case.
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
88 *
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
89 * It works just fine with Intel's icc compiler, version 7.0.
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
90 */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
91
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
92 /* pivot = 1 / pivot */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
93 RECIPROCAL (pivot) ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
94
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
95 #if defined (USE_NO_BLAS)
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
96 for (i = 0 ; i < n ; i++)
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
97 {
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
98 /* X [i] *= pivot ; */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
99 x = X [i] ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
100 MULT (X [i], x, pivot) ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
101 }
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
102 #else
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
103 BLAS_SCAL (n, pivot, X) ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
104 #endif
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
105
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
106 #else
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
107
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
108 /* -------------------------------------------------------------- */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
109 /* divide x by the pivot */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
110 /* -------------------------------------------------------------- */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
111
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
112 /* This is slightly more accurate, particularly if the pivot column
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
113 * consists of only IEEE subnormals. Always do this if UMFPACK is
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
114 * being compiled as a built-in routine or mexFunction in MATLAB,
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
115 * or if gcc is being used with complex matrices. */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
116
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
117 for (i = 0 ; i < n ; i++)
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
118 {
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
119 /* X [i] /= pivot ; */
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
120 x = X [i] ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
121 DIV (X [i], x, pivot) ;
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
122 }
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
123
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
124 #endif
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
125
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
126 }
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
127 }