Mercurial > octave-nkf
comparison liboctave/UMFPACK/UMFPACK/Source/umf_usolve.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_usolve =========================================================== */ | |
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 Ux = b, where U is the upper triangular factor of a matrix. */ | |
12 /* B is overwritten with the solution X. */ | |
13 /* Returns the floating point operation count */ | |
14 | |
15 #include "umf_internal.h" | |
16 | |
17 GLOBAL double UMF_usolve | |
18 ( | |
19 NumericType *Numeric, | |
20 Entry X [ ], /* b on input, solution x on output */ | |
21 Int Pattern [ ] /* a work array of size n */ | |
22 ) | |
23 { | |
24 /* ---------------------------------------------------------------------- */ | |
25 /* local variables */ | |
26 /* ---------------------------------------------------------------------- */ | |
27 | |
28 Entry xk ; | |
29 Entry *xp, *D, *Uval ; | |
30 Int k, deg, j, *ip, col, *Upos, *Uilen, pos, | |
31 *Uip, n, ulen, up, newUchain, npiv, n1, *Ui ; | |
32 | |
33 /* ---------------------------------------------------------------------- */ | |
34 /* get parameters */ | |
35 /* ---------------------------------------------------------------------- */ | |
36 | |
37 if (Numeric->n_row != Numeric->n_col) return (0.) ; | |
38 n = Numeric->n_row ; | |
39 npiv = Numeric->npiv ; | |
40 Upos = Numeric->Upos ; | |
41 Uilen = Numeric->Uilen ; | |
42 Uip = Numeric->Uip ; | |
43 D = Numeric->D ; | |
44 n1 = Numeric->n1 ; | |
45 | |
46 #ifndef NDEBUG | |
47 DEBUG4 (("Usolve start: npiv = "ID" n = "ID"\n", npiv, n)) ; | |
48 for (j = 0 ; j < n ; j++) | |
49 { | |
50 DEBUG4 (("Usolve start "ID": ", j)) ; | |
51 EDEBUG4 (X [j]) ; | |
52 DEBUG4 (("\n")) ; | |
53 } | |
54 #endif | |
55 | |
56 /* ---------------------------------------------------------------------- */ | |
57 /* singular case */ | |
58 /* ---------------------------------------------------------------------- */ | |
59 | |
60 #ifndef NO_DIVIDE_BY_ZERO | |
61 /* handle the singular part of D, up to just before the last pivot */ | |
62 for (k = n-1 ; k >= npiv ; k--) | |
63 { | |
64 /* This is an *** intentional *** divide-by-zero, to get Inf or Nan, | |
65 * as appropriate. It is not a bug. */ | |
66 ASSERT (IS_ZERO (D [k])) ; | |
67 xk = X [k] ; | |
68 /* X [k] = xk / D [k] ; */ | |
69 DIV (X [k], xk, D [k]) ; | |
70 } | |
71 #else | |
72 /* Do not divide by zero */ | |
73 #endif | |
74 | |
75 deg = Numeric->ulen ; | |
76 if (deg > 0) | |
77 { | |
78 /* :: make last pivot row of U (singular matrices only) :: */ | |
79 for (j = 0 ; j < deg ; j++) | |
80 { | |
81 DEBUG1 (("Last row of U: j="ID"\n", j)) ; | |
82 DEBUG1 (("Last row of U: Upattern[j]="ID"\n", | |
83 Numeric->Upattern [j]) ); | |
84 Pattern [j] = Numeric->Upattern [j] ; | |
85 } | |
86 } | |
87 | |
88 /* ---------------------------------------------------------------------- */ | |
89 /* nonsingletons */ | |
90 /* ---------------------------------------------------------------------- */ | |
91 | |
92 for (k = npiv-1 ; k >= n1 ; k--) | |
93 { | |
94 | |
95 /* ------------------------------------------------------------------ */ | |
96 /* use row k of U */ | |
97 /* ------------------------------------------------------------------ */ | |
98 | |
99 up = Uip [k] ; | |
100 ulen = Uilen [k] ; | |
101 newUchain = (up < 0) ; | |
102 if (newUchain) | |
103 { | |
104 up = -up ; | |
105 xp = (Entry *) (Numeric->Memory + up + UNITS (Int, ulen)) ; | |
106 } | |
107 else | |
108 { | |
109 xp = (Entry *) (Numeric->Memory + up) ; | |
110 } | |
111 | |
112 xk = X [k] ; | |
113 for (j = 0 ; j < deg ; j++) | |
114 { | |
115 DEBUG4 ((" k "ID" col "ID" value", k, Pattern [j])) ; | |
116 EDEBUG4 (*xp) ; | |
117 DEBUG4 (("\n")) ; | |
118 /* xk -= X [Pattern [j]] * (*xp) ; */ | |
119 MULT_SUB (xk, X [Pattern [j]], *xp) ; | |
120 xp++ ; | |
121 } | |
122 | |
123 #ifndef NO_DIVIDE_BY_ZERO | |
124 /* Go ahead and divide by zero if D [k] is zero */ | |
125 /* X [k] = xk / D [k] ; */ | |
126 DIV (X [k], xk, D [k]) ; | |
127 #else | |
128 /* Do not divide by zero */ | |
129 if (IS_NONZERO (D [k])) | |
130 { | |
131 /* X [k] = xk / D [k] ; */ | |
132 DIV (X [k], xk, D [k]) ; | |
133 } | |
134 #endif | |
135 | |
136 /* ------------------------------------------------------------------ */ | |
137 /* make row k-1 of U in Pattern [0..deg-1] */ | |
138 /* ------------------------------------------------------------------ */ | |
139 | |
140 if (k == n1) break ; | |
141 | |
142 if (newUchain) | |
143 { | |
144 /* next row is a new Uchain */ | |
145 deg = ulen ; | |
146 ASSERT (IMPLIES (k == 0, deg == 0)) ; | |
147 DEBUG4 (("end of chain for row of U "ID" deg "ID"\n", k-1, deg)) ; | |
148 ip = (Int *) (Numeric->Memory + up) ; | |
149 for (j = 0 ; j < deg ; j++) | |
150 { | |
151 col = *ip++ ; | |
152 DEBUG4 ((" k "ID" col "ID"\n", k-1, col)) ; | |
153 ASSERT (k <= col) ; | |
154 Pattern [j] = col ; | |
155 } | |
156 } | |
157 else | |
158 { | |
159 deg -= ulen ; | |
160 DEBUG4 (("middle of chain for row of U "ID" deg "ID"\n", k, deg)) ; | |
161 ASSERT (deg >= 0) ; | |
162 pos = Upos [k] ; | |
163 if (pos != EMPTY) | |
164 { | |
165 /* add the pivot column */ | |
166 DEBUG4 (("k "ID" add pivot entry at pos "ID"\n", k, pos)) ; | |
167 ASSERT (pos >= 0 && pos <= deg) ; | |
168 Pattern [deg++] = Pattern [pos] ; | |
169 Pattern [pos] = k ; | |
170 } | |
171 } | |
172 } | |
173 | |
174 /* ---------------------------------------------------------------------- */ | |
175 /* singletons */ | |
176 /* ---------------------------------------------------------------------- */ | |
177 | |
178 for (k = n1 - 1 ; k >= 0 ; k--) | |
179 { | |
180 deg = Uilen [k] ; | |
181 xk = X [k] ; | |
182 DEBUG4 (("Singleton k "ID"\n", k)) ; | |
183 if (deg > 0) | |
184 { | |
185 up = Uip [k] ; | |
186 Ui = (Int *) (Numeric->Memory + up) ; | |
187 up += UNITS (Int, deg) ; | |
188 Uval = (Entry *) (Numeric->Memory + up) ; | |
189 for (j = 0 ; j < deg ; j++) | |
190 { | |
191 DEBUG4 ((" k "ID" col "ID" value", k, Ui [j])) ; | |
192 EDEBUG4 (Uval [j]) ; | |
193 DEBUG4 (("\n")) ; | |
194 /* xk -= X [Ui [j]] * Uval [j] ; */ | |
195 ASSERT (Ui [j] >= 0 && Ui [j] < n) ; | |
196 MULT_SUB (xk, X [Ui [j]], Uval [j]) ; | |
197 } | |
198 } | |
199 | |
200 #ifndef NO_DIVIDE_BY_ZERO | |
201 /* Go ahead and divide by zero if D [k] is zero */ | |
202 /* X [k] = xk / D [k] ; */ | |
203 DIV (X [k], xk, D [k]) ; | |
204 #else | |
205 /* Do not divide by zero */ | |
206 if (IS_NONZERO (D [k])) | |
207 { | |
208 /* X [k] = xk / D [k] ; */ | |
209 DIV (X [k], xk, D [k]) ; | |
210 } | |
211 #endif | |
212 | |
213 } | |
214 | |
215 #ifndef NDEBUG | |
216 for (j = 0 ; j < n ; j++) | |
217 { | |
218 DEBUG4 (("Usolve done "ID": ", j)) ; | |
219 EDEBUG4 (X [j]) ; | |
220 DEBUG4 (("\n")) ; | |
221 } | |
222 DEBUG4 (("Usolve done.\n")) ; | |
223 #endif | |
224 | |
225 return (DIV_FLOPS * ((double) n) + MULTSUB_FLOPS * ((double) Numeric->unz)); | |
226 } |