|
5164
|
1 /* ========================================================================== */ |
|
|
2 /* === UMF_ltsolve ========================================================== */ |
|
|
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 L'x = b or L.'x=b, where L is the lower 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_lhsolve /* solve L'x=b (complex conjugate transpose) */ |
|
|
20 #else |
|
|
21 UMF_ltsolve /* solve L.'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 Entry xk ; |
|
|
30 Entry *xp, *Lval ; |
|
|
31 Int k, deg, *ip, j, row, *Lpos, *Lilen, kstart, kend, *Lip, llen, |
|
|
32 lp, pos, npiv, n1, *Li ; |
|
|
33 |
|
|
34 /* ---------------------------------------------------------------------- */ |
|
|
35 |
|
|
36 if (Numeric->n_row != Numeric->n_col) return (0.) ; |
|
|
37 npiv = Numeric->npiv ; |
|
|
38 Lpos = Numeric->Lpos ; |
|
|
39 Lilen = Numeric->Lilen ; |
|
|
40 Lip = Numeric->Lip ; |
|
|
41 kstart = npiv ; |
|
|
42 n1 = Numeric->n1 ; |
|
|
43 |
|
|
44 #ifndef NDEBUG |
|
|
45 DEBUG4 (("Ltsolve start:\n")) ; |
|
|
46 for (j = 0 ; j < Numeric->n_row ; j++) |
|
|
47 { |
|
|
48 DEBUG4 (("Ltsolve start "ID": ", j)) ; |
|
|
49 EDEBUG4 (X [j]) ; |
|
|
50 DEBUG4 (("\n")) ; |
|
|
51 } |
|
|
52 #endif |
|
|
53 |
|
|
54 /* ---------------------------------------------------------------------- */ |
|
|
55 /* non-singletons */ |
|
|
56 /* ---------------------------------------------------------------------- */ |
|
|
57 |
|
|
58 for (kend = npiv-1 ; kend >= n1 ; kend = kstart-1) |
|
|
59 { |
|
|
60 |
|
|
61 /* ------------------------------------------------------------------ */ |
|
|
62 /* find the start of this Lchain */ |
|
|
63 /* ------------------------------------------------------------------ */ |
|
|
64 |
|
|
65 /* for (kstart = kend ; kstart >= 0 && Lip [kstart] > 0 ; kstart--) ; */ |
|
|
66 kstart = kend ; |
|
|
67 while (kstart >= 0 && Lip [kstart] > 0) |
|
|
68 { |
|
|
69 kstart-- ; |
|
|
70 } |
|
|
71 |
|
|
72 /* the Lchain goes from kstart to kend */ |
|
|
73 |
|
|
74 /* ------------------------------------------------------------------ */ |
|
|
75 /* scan the whole chain to find the pattern of the last column of L */ |
|
|
76 /* ------------------------------------------------------------------ */ |
|
|
77 |
|
|
78 deg = 0 ; |
|
|
79 DEBUG4 (("start of chain for column of L\n")) ; |
|
|
80 for (k = kstart ; k <= kend ; k++) |
|
|
81 { |
|
|
82 ASSERT (k >= 0 && k < npiv) ; |
|
|
83 |
|
|
84 /* -------------------------------------------------------------- */ |
|
|
85 /* make column k of L in Pattern [0..deg-1] */ |
|
|
86 /* -------------------------------------------------------------- */ |
|
|
87 |
|
|
88 /* remove pivot row */ |
|
|
89 pos = Lpos [k] ; |
|
|
90 if (pos != EMPTY) |
|
|
91 { |
|
|
92 DEBUG4 ((" k "ID" removing row "ID" at position "ID"\n", |
|
|
93 k, Pattern [pos], pos)) ; |
|
|
94 ASSERT (k != kstart) ; |
|
|
95 ASSERT (deg > 0) ; |
|
|
96 ASSERT (pos >= 0 && pos < deg) ; |
|
|
97 ASSERT (Pattern [pos] == k) ; |
|
|
98 Pattern [pos] = Pattern [--deg] ; |
|
|
99 } |
|
|
100 |
|
|
101 /* concatenate the pattern */ |
|
|
102 lp = Lip [k] ; |
|
|
103 if (k == kstart) |
|
|
104 { |
|
|
105 lp = -lp ; |
|
|
106 } |
|
|
107 ASSERT (lp > 0) ; |
|
|
108 ip = (Int *) (Numeric->Memory + lp) ; |
|
|
109 llen = Lilen [k] ; |
|
|
110 for (j = 0 ; j < llen ; j++) |
|
|
111 { |
|
|
112 row = *ip++ ; |
|
|
113 DEBUG4 ((" row "ID" k "ID"\n", row, k)) ; |
|
|
114 ASSERT (row > k) ; |
|
|
115 Pattern [deg++] = row ; |
|
|
116 } |
|
|
117 |
|
|
118 } |
|
|
119 /* Pattern [0..deg-1] is now the pattern of column kend */ |
|
|
120 |
|
|
121 /* ------------------------------------------------------------------ */ |
|
|
122 /* solve using this chain, in reverse order */ |
|
|
123 /* ------------------------------------------------------------------ */ |
|
|
124 |
|
|
125 DEBUG4 (("Unwinding Lchain\n")) ; |
|
|
126 for (k = kend ; k >= kstart ; k--) |
|
|
127 { |
|
|
128 |
|
|
129 /* -------------------------------------------------------------- */ |
|
|
130 /* use column k of L */ |
|
|
131 /* -------------------------------------------------------------- */ |
|
|
132 |
|
|
133 ASSERT (k >= 0 && k < npiv) ; |
|
|
134 lp = Lip [k] ; |
|
|
135 if (k == kstart) |
|
|
136 { |
|
|
137 lp = -lp ; |
|
|
138 } |
|
|
139 ASSERT (lp > 0) ; |
|
|
140 llen = Lilen [k] ; |
|
|
141 xp = (Entry *) (Numeric->Memory + lp + UNITS (Int, llen)) ; |
|
|
142 xk = X [k] ; |
|
|
143 for (j = 0 ; j < deg ; j++) |
|
|
144 { |
|
|
145 DEBUG4 ((" row "ID" k "ID" value", Pattern [j], k)) ; |
|
|
146 EDEBUG4 (*xp) ; |
|
|
147 DEBUG4 (("\n")) ; |
|
|
148 |
|
|
149 #ifdef CONJUGATE_SOLVE |
|
|
150 /* xk -= X [Pattern [j]] * conjugate (*xp) ; */ |
|
|
151 MULT_SUB_CONJ (xk, X [Pattern [j]], *xp) ; |
|
|
152 #else |
|
|
153 /* xk -= X [Pattern [j]] * (*xp) ; */ |
|
|
154 MULT_SUB (xk, X [Pattern [j]], *xp) ; |
|
|
155 #endif |
|
|
156 |
|
|
157 xp++ ; |
|
|
158 } |
|
|
159 X [k] = xk ; |
|
|
160 |
|
|
161 /* -------------------------------------------------------------- */ |
|
|
162 /* construct column k-1 of L */ |
|
|
163 /* -------------------------------------------------------------- */ |
|
|
164 |
|
|
165 /* un-concatenate the pattern */ |
|
|
166 deg -= llen ; |
|
|
167 |
|
|
168 /* add pivot row */ |
|
|
169 pos = Lpos [k] ; |
|
|
170 if (pos != EMPTY) |
|
|
171 { |
|
|
172 DEBUG4 ((" k "ID" adding row "ID" at position "ID"\n", |
|
|
173 k, k, pos)) ; |
|
|
174 ASSERT (k != kstart) ; |
|
|
175 ASSERT (pos >= 0 && pos <= deg) ; |
|
|
176 Pattern [deg++] = Pattern [pos] ; |
|
|
177 Pattern [pos] = k ; |
|
|
178 } |
|
|
179 } |
|
|
180 } |
|
|
181 |
|
|
182 /* ---------------------------------------------------------------------- */ |
|
|
183 /* singletons */ |
|
|
184 /* ---------------------------------------------------------------------- */ |
|
|
185 |
|
|
186 for (k = n1 - 1 ; k >= 0 ; k--) |
|
|
187 { |
|
|
188 DEBUG4 (("Singleton k "ID"\n", k)) ; |
|
|
189 deg = Lilen [k] ; |
|
|
190 if (deg > 0) |
|
|
191 { |
|
|
192 xk = X [k] ; |
|
|
193 lp = Lip [k] ; |
|
|
194 Li = (Int *) (Numeric->Memory + lp) ; |
|
|
195 lp += UNITS (Int, deg) ; |
|
|
196 Lval = (Entry *) (Numeric->Memory + lp) ; |
|
|
197 for (j = 0 ; j < deg ; j++) |
|
|
198 { |
|
|
199 DEBUG4 ((" row "ID" k "ID" value", Li [j], k)) ; |
|
|
200 EDEBUG4 (Lval [j]) ; |
|
|
201 DEBUG4 (("\n")) ; |
|
|
202 #ifdef CONJUGATE_SOLVE |
|
|
203 /* xk -= X [Li [j]] * conjugate (Lval [j]) ; */ |
|
|
204 MULT_SUB_CONJ (xk, X [Li [j]], Lval [j]) ; |
|
|
205 #else |
|
|
206 /* xk -= X [Li [j]] * Lval [j] ; */ |
|
|
207 MULT_SUB (xk, X [Li [j]], Lval [j]) ; |
|
|
208 #endif |
|
|
209 } |
|
|
210 X [k] = xk ; |
|
|
211 } |
|
|
212 } |
|
|
213 |
|
|
214 #ifndef NDEBUG |
|
|
215 for (j = 0 ; j < Numeric->n_row ; j++) |
|
|
216 { |
|
|
217 DEBUG4 (("Ltsolve done "ID": ", j)) ; |
|
|
218 EDEBUG4 (X [j]) ; |
|
|
219 DEBUG4 (("\n")) ; |
|
|
220 } |
|
|
221 DEBUG4 (("Ltsolve done.\n")) ; |
|
|
222 #endif |
|
|
223 |
|
|
224 return (MULTSUB_FLOPS * ((double) Numeric->lnz)) ; |
|
|
225 } |
