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