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
UMF_usolve(NumericType * Numeric,Entry X[],Int Pattern[])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 }
227