1 /* ========================================================================== */
2 /* === UMFPACK_scale ======================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* Copyright (c) 2005-2012 by Timothy A. Davis, http://www.suitesparse.com. */
7 /* All Rights Reserved. See ../Doc/License.txt for License. */
8 /* -------------------------------------------------------------------------- */
9
10 /*
11 User-callable. Applies the scale factors computed during numerical
12 factorization to a vector. See umfpack_scale.h for more details.
13
14 The LU factorization is L*U = P*R*A*Q, where P and Q are permutation
15 matrices, and R is diagonal. This routine computes X = R * B using the
16 matrix R stored in the Numeric object.
17
18 Returns FALSE if any argument is invalid, TRUE otherwise.
19
20 If R not present in the Numeric object, then R = I and no floating-point
21 work is done. B is simply copied into X.
22 */
23
24 #include "umf_internal.h"
25 #include "umf_valid_numeric.h"
26
UMFPACK_scale(double Xx[],double Xz[],const double Bx[],const double Bz[],void * NumericHandle)27 GLOBAL Int UMFPACK_scale
28 (
29 double Xx [ ],
30 #ifdef COMPLEX
31 double Xz [ ],
32 #endif
33 const double Bx [ ],
34 #ifdef COMPLEX
35 const double Bz [ ],
36 #endif
37 void *NumericHandle
38 )
39 {
40 /* ---------------------------------------------------------------------- */
41 /* local variables */
42 /* ---------------------------------------------------------------------- */
43
44 NumericType *Numeric ;
45 Int n, i ;
46 double *Rs ;
47 #ifdef COMPLEX
48 Int split = SPLIT (Xz) && SPLIT (Bz) ;
49 #endif
50
51 Numeric = (NumericType *) NumericHandle ;
52 if (!UMF_valid_numeric (Numeric))
53 {
54 return (UMFPACK_ERROR_invalid_Numeric_object) ;
55 }
56
57 n = Numeric->n_row ;
58 Rs = Numeric->Rs ;
59
60 if (!Xx || !Bx)
61 {
62 return (UMFPACK_ERROR_argument_missing) ;
63 }
64
65 /* ---------------------------------------------------------------------- */
66 /* X = R*B or R\B */
67 /* ---------------------------------------------------------------------- */
68
69 if (Rs != (double *) NULL)
70 {
71 #ifndef NRECIPROCAL
72 if (Numeric->do_recip)
73 {
74 /* multiply by the scale factors */
75 #ifdef COMPLEX
76 if (split)
77 {
78 for (i = 0 ; i < n ; i++)
79 {
80 Xx [i] = Bx [i] * Rs [i] ;
81 Xz [i] = Bz [i] * Rs [i] ;
82 }
83 }
84 else
85 {
86 for (i = 0 ; i < n ; i++)
87 {
88 Xx [2*i ] = Bx [2*i ] * Rs [i] ;
89 Xx [2*i+1] = Bx [2*i+1] * Rs [i] ;
90 }
91 }
92 #else
93 for (i = 0 ; i < n ; i++)
94 {
95 Xx [i] = Bx [i] * Rs [i] ;
96 }
97 #endif
98 }
99 else
100 #endif
101 {
102 /* divide by the scale factors */
103 #ifdef COMPLEX
104 if (split)
105 {
106 for (i = 0 ; i < n ; i++)
107 {
108 Xx [i] = Bx [i] / Rs [i] ;
109 Xz [i] = Bz [i] / Rs [i] ;
110 }
111 }
112 else
113 {
114 for (i = 0 ; i < n ; i++)
115 {
116 Xx [2*i ] = Bx [2*i ] / Rs [i] ;
117 Xx [2*i+1] = Bx [2*i+1] / Rs [i] ;
118 }
119 }
120 #else
121 for (i = 0 ; i < n ; i++)
122 {
123 Xx [i] = Bx [i] / Rs [i] ;
124 }
125 #endif
126 }
127 }
128 else
129 {
130 /* no scale factors, just copy B into X */
131 #ifdef COMPLEX
132 if (split)
133 {
134 for (i = 0 ; i < n ; i++)
135 {
136 Xx [i] = Bx [i] ;
137 Xz [i] = Bz [i] ;
138 }
139 }
140 else
141 {
142 for (i = 0 ; i < n ; i++)
143 {
144 Xx [2*i ] = Bx [2*i ] ;
145 Xx [2*i+1] = Bx [2*i+1] ;
146 }
147 }
148 #else
149 for (i = 0 ; i < n ; i++)
150 {
151 Xx [i] = Bx [i] ;
152 }
153 #endif
154 }
155
156 return (UMFPACK_OK) ;
157 }
158