Mercurial > octave-nkf
comparison liboctave/UMFPACK/UMFPACK/Source/umf_utsolve.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 |
comparison
equal
deleted
inserted
replaced
5163:9f3299378193 | 5164:57077d0ddc8e |
---|---|
1 /* ========================================================================== */ | |
2 /* === UMF_utsolve ========================================================== */ | |
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 /* solves U'x = b or U.'x=b, where U is the upper triangular factor of a */ | |
12 /* matrix. B is overwritten with the solution X. */ | |
13 /* Returns the floating point operation count */ | |
14 | |
15 #include "umf_internal.h" | |
16 | |
17 GLOBAL double | |
18 #ifdef CONJUGATE_SOLVE | |
19 UMF_uhsolve /* solve U'x=b (complex conjugate transpose) */ | |
20 #else | |
21 UMF_utsolve /* solve U.'x=b (array transpose) */ | |
22 #endif | |
23 ( | |
24 NumericType *Numeric, | |
25 Entry X [ ], /* b on input, solution x on output */ | |
26 Int Pattern [ ] /* a work array of size n */ | |
27 ) | |
28 { | |
29 /* ---------------------------------------------------------------------- */ | |
30 /* local variables */ | |
31 /* ---------------------------------------------------------------------- */ | |
32 | |
33 Entry xk ; | |
34 Entry *xp, *D, *Uval ; | |
35 Int k, deg, j, *ip, col, *Upos, *Uilen, kstart, kend, up, | |
36 *Uip, n, uhead, ulen, pos, npiv, n1, *Ui ; | |
37 | |
38 /* ---------------------------------------------------------------------- */ | |
39 /* get parameters */ | |
40 /* ---------------------------------------------------------------------- */ | |
41 | |
42 if (Numeric->n_row != Numeric->n_col) return (0.) ; | |
43 n = Numeric->n_row ; | |
44 npiv = Numeric->npiv ; | |
45 Upos = Numeric->Upos ; | |
46 Uilen = Numeric->Uilen ; | |
47 Uip = Numeric->Uip ; | |
48 D = Numeric->D ; | |
49 kend = 0 ; | |
50 n1 = Numeric->n1 ; | |
51 | |
52 #ifndef NDEBUG | |
53 DEBUG4 (("Utsolve start: npiv "ID" n "ID"\n", npiv, n)) ; | |
54 for (j = 0 ; j < n ; j++) | |
55 { | |
56 DEBUG4 (("Utsolve start "ID": ", j)) ; | |
57 EDEBUG4 (X [j]) ; | |
58 DEBUG4 (("\n")) ; | |
59 } | |
60 #endif | |
61 | |
62 /* ---------------------------------------------------------------------- */ | |
63 /* singletons */ | |
64 /* ---------------------------------------------------------------------- */ | |
65 | |
66 for (k = 0 ; k < n1 ; k++) | |
67 { | |
68 DEBUG4 (("Singleton k "ID"\n", k)) ; | |
69 | |
70 #ifndef NO_DIVIDE_BY_ZERO | |
71 /* Go ahead and divide by zero if D [k] is zero. */ | |
72 #ifdef CONJUGATE_SOLVE | |
73 /* xk = X [k] / conjugate (D [k]) ; */ | |
74 DIV_CONJ (xk, X [k], D [k]) ; | |
75 #else | |
76 /* xk = X [k] / D [k] ; */ | |
77 DIV (xk, X [k], D [k]) ; | |
78 #endif | |
79 #else | |
80 /* Do not divide by zero */ | |
81 if (IS_NONZERO (D [k])) | |
82 { | |
83 #ifdef CONJUGATE_SOLVE | |
84 /* xk = X [k] / conjugate (D [k]) ; */ | |
85 DIV_CONJ (xk, X [k], D [k]) ; | |
86 #else | |
87 /* xk = X [k] / D [k] ; */ | |
88 DIV (xk, X [k], D [k]) ; | |
89 #endif | |
90 } | |
91 #endif | |
92 | |
93 X [k] = xk ; | |
94 deg = Uilen [k] ; | |
95 if (deg > 0 && IS_NONZERO (xk)) | |
96 { | |
97 up = Uip [k] ; | |
98 Ui = (Int *) (Numeric->Memory + up) ; | |
99 up += UNITS (Int, deg) ; | |
100 Uval = (Entry *) (Numeric->Memory + up) ; | |
101 for (j = 0 ; j < deg ; j++) | |
102 { | |
103 DEBUG4 ((" k "ID" col "ID" value", k, Ui [j])) ; | |
104 EDEBUG4 (Uval [j]) ; | |
105 DEBUG4 (("\n")) ; | |
106 #ifdef CONJUGATE_SOLVE | |
107 /* X [Ui [j]] -= xk * conjugate (Uval [j]) ; */ | |
108 MULT_SUB_CONJ (X [Ui [j]], xk, Uval [j]) ; | |
109 #else | |
110 /* X [Ui [j]] -= xk * Uval [j] ; */ | |
111 MULT_SUB (X [Ui [j]], xk, Uval [j]) ; | |
112 #endif | |
113 } | |
114 } | |
115 } | |
116 | |
117 /* ---------------------------------------------------------------------- */ | |
118 /* nonsingletons */ | |
119 /* ---------------------------------------------------------------------- */ | |
120 | |
121 for (kstart = n1 ; kstart < npiv ; kstart = kend + 1) | |
122 { | |
123 | |
124 /* ------------------------------------------------------------------ */ | |
125 /* find the end of this Uchain */ | |
126 /* ------------------------------------------------------------------ */ | |
127 | |
128 DEBUG4 (("kstart "ID" kend "ID"\n", kstart, kend)) ; | |
129 /* for (kend = kstart ; kend < npiv && Uip [kend+1] > 0 ; kend++) ; */ | |
130 kend = kstart ; | |
131 while (kend < npiv && Uip [kend+1] > 0) | |
132 { | |
133 kend++ ; | |
134 } | |
135 | |
136 /* ------------------------------------------------------------------ */ | |
137 /* scan the whole Uchain to find the pattern of the first row of U */ | |
138 /* ------------------------------------------------------------------ */ | |
139 | |
140 k = kend+1 ; | |
141 DEBUG4 (("\nKend "ID" K "ID"\n", kend, k)) ; | |
142 | |
143 /* ------------------------------------------------------------------ */ | |
144 /* start with last row in Uchain of U in Pattern [0..deg-1] */ | |
145 /* ------------------------------------------------------------------ */ | |
146 | |
147 if (k == npiv) | |
148 { | |
149 deg = Numeric->ulen ; | |
150 if (deg > 0) | |
151 { | |
152 /* :: make last pivot row of U (singular matrices only) :: */ | |
153 for (j = 0 ; j < deg ; j++) | |
154 { | |
155 Pattern [j] = Numeric->Upattern [j] ; | |
156 } | |
157 } | |
158 } | |
159 else | |
160 { | |
161 ASSERT (k >= 0 && k < npiv) ; | |
162 up = -Uip [k] ; | |
163 ASSERT (up > 0) ; | |
164 deg = Uilen [k] ; | |
165 DEBUG4 (("end of chain for row of U "ID" deg "ID"\n", k-1, deg)) ; | |
166 ip = (Int *) (Numeric->Memory + up) ; | |
167 for (j = 0 ; j < deg ; j++) | |
168 { | |
169 col = *ip++ ; | |
170 DEBUG4 ((" k "ID" col "ID"\n", k-1, col)) ; | |
171 ASSERT (k <= col) ; | |
172 Pattern [j] = col ; | |
173 } | |
174 } | |
175 | |
176 /* empty the stack at the bottom of Pattern */ | |
177 uhead = n ; | |
178 | |
179 for (k = kend ; k > kstart ; k--) | |
180 { | |
181 /* Pattern [0..deg-1] is the pattern of row k of U */ | |
182 | |
183 /* -------------------------------------------------------------- */ | |
184 /* make row k-1 of U in Pattern [0..deg-1] */ | |
185 /* -------------------------------------------------------------- */ | |
186 | |
187 ASSERT (k >= 0 && k < npiv) ; | |
188 ulen = Uilen [k] ; | |
189 /* delete, and push on the stack */ | |
190 for (j = 0 ; j < ulen ; j++) | |
191 { | |
192 ASSERT (uhead >= deg) ; | |
193 Pattern [--uhead] = Pattern [--deg] ; | |
194 } | |
195 DEBUG4 (("middle of chain for row of U "ID" deg "ID"\n", k, deg)) ; | |
196 ASSERT (deg >= 0) ; | |
197 | |
198 pos = Upos [k] ; | |
199 if (pos != EMPTY) | |
200 { | |
201 /* add the pivot column */ | |
202 DEBUG4 (("k "ID" add pivot entry at position "ID"\n", k, pos)) ; | |
203 ASSERT (pos >= 0 && pos <= deg) ; | |
204 Pattern [deg++] = Pattern [pos] ; | |
205 Pattern [pos] = k ; | |
206 } | |
207 } | |
208 | |
209 /* Pattern [0..deg-1] is now the pattern of the first row in Uchain */ | |
210 | |
211 /* ------------------------------------------------------------------ */ | |
212 /* solve using this Uchain, in reverse order */ | |
213 /* ------------------------------------------------------------------ */ | |
214 | |
215 DEBUG4 (("Unwinding Uchain\n")) ; | |
216 for (k = kstart ; k <= kend ; k++) | |
217 { | |
218 | |
219 /* -------------------------------------------------------------- */ | |
220 /* construct row k */ | |
221 /* -------------------------------------------------------------- */ | |
222 | |
223 ASSERT (k >= 0 && k < npiv) ; | |
224 pos = Upos [k] ; | |
225 if (pos != EMPTY) | |
226 { | |
227 /* remove the pivot column */ | |
228 DEBUG4 (("k "ID" add pivot entry at position "ID"\n", k, pos)) ; | |
229 ASSERT (k > kstart) ; | |
230 ASSERT (pos >= 0 && pos < deg) ; | |
231 ASSERT (Pattern [pos] == k) ; | |
232 Pattern [pos] = Pattern [--deg] ; | |
233 } | |
234 | |
235 up = Uip [k] ; | |
236 ulen = Uilen [k] ; | |
237 if (k > kstart) | |
238 { | |
239 /* concatenate the deleted pattern; pop from the stack */ | |
240 for (j = 0 ; j < ulen ; j++) | |
241 { | |
242 ASSERT (deg <= uhead && uhead < n) ; | |
243 Pattern [deg++] = Pattern [uhead++] ; | |
244 } | |
245 DEBUG4 (("middle of chain, row of U "ID" deg "ID"\n", k, deg)) ; | |
246 ASSERT (deg >= 0) ; | |
247 } | |
248 | |
249 /* -------------------------------------------------------------- */ | |
250 /* use row k of U */ | |
251 /* -------------------------------------------------------------- */ | |
252 | |
253 #ifndef NO_DIVIDE_BY_ZERO | |
254 /* Go ahead and divide by zero if D [k] is zero. */ | |
255 #ifdef CONJUGATE_SOLVE | |
256 /* xk = X [k] / conjugate (D [k]) ; */ | |
257 DIV_CONJ (xk, X [k], D [k]) ; | |
258 #else | |
259 /* xk = X [k] / D [k] ; */ | |
260 DIV (xk, X [k], D [k]) ; | |
261 #endif | |
262 #else | |
263 /* Do not divide by zero */ | |
264 if (IS_NONZERO (D [k])) | |
265 { | |
266 #ifdef CONJUGATE_SOLVE | |
267 /* xk = X [k] / conjugate (D [k]) ; */ | |
268 DIV_CONJ (xk, X [k], D [k]) ; | |
269 #else | |
270 /* xk = X [k] / D [k] ; */ | |
271 DIV (xk, X [k], D [k]) ; | |
272 #endif | |
273 } | |
274 #endif | |
275 | |
276 X [k] = xk ; | |
277 if (IS_NONZERO (xk)) | |
278 { | |
279 if (k == kstart) | |
280 { | |
281 up = -up ; | |
282 xp = (Entry *) (Numeric->Memory + up + UNITS (Int, ulen)) ; | |
283 } | |
284 else | |
285 { | |
286 xp = (Entry *) (Numeric->Memory + up) ; | |
287 } | |
288 for (j = 0 ; j < deg ; j++) | |
289 { | |
290 DEBUG4 ((" k "ID" col "ID" value", k, Pattern [j])) ; | |
291 EDEBUG4 (*xp) ; | |
292 DEBUG4 (("\n")) ; | |
293 #ifdef CONJUGATE_SOLVE | |
294 /* X [Pattern [j]] -= xk * conjugate (*xp) ; */ | |
295 MULT_SUB_CONJ (X [Pattern [j]], xk, *xp) ; | |
296 #else | |
297 /* X [Pattern [j]] -= xk * (*xp) ; */ | |
298 MULT_SUB (X [Pattern [j]], xk, *xp) ; | |
299 #endif | |
300 xp++ ; | |
301 } | |
302 } | |
303 } | |
304 ASSERT (uhead == n) ; | |
305 } | |
306 | |
307 #ifndef NO_DIVIDE_BY_ZERO | |
308 for (k = npiv ; k < n ; k++) | |
309 { | |
310 /* This is an *** intentional *** divide-by-zero, to get Inf or Nan, | |
311 * as appropriate. It is not a bug. */ | |
312 ASSERT (IS_ZERO (D [k])) ; | |
313 /* For conjugate solve, D [k] == conjugate (D [k]), in this case */ | |
314 /* xk = X [k] / D [k] ; */ | |
315 DIV (xk, X [k], D [k]) ; | |
316 X [k] = xk ; | |
317 } | |
318 #endif | |
319 | |
320 #ifndef NDEBUG | |
321 for (j = 0 ; j < n ; j++) | |
322 { | |
323 DEBUG4 (("Utsolve done "ID": ", j)) ; | |
324 EDEBUG4 (X [j]) ; | |
325 DEBUG4 (("\n")) ; | |
326 } | |
327 DEBUG4 (("Utsolve done.\n")) ; | |
328 #endif | |
329 | |
330 return (DIV_FLOPS * ((double) n) + MULTSUB_FLOPS * ((double) Numeric->unz)); | |
331 } |