1 /* ========================================================================== */
2 /* === UMF_ltsolve ========================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* Copyright (c) 2005-2012 by Timothy A. Davis, http://www.suitesparse.com. */
7 /* All Rights Reserved. See ../Doc/License.txt for License. */
8 /* -------------------------------------------------------------------------- */
9
10 /* Solves L'x = b or L.'x=b, where L is the lower triangular factor of a */
11 /* matrix. B is overwritten with the solution X. */
12 /* Returns the floating point operation count */
13
14 #include "umf_internal.h"
15 #include "umf_ltsolve.h"
16
17 GLOBAL double
18 #ifdef CONJUGATE_SOLVE
UMF_lhsolve(NumericType * Numeric,Entry X[],Int Pattern[])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 }
226