Mercurial > octave-nkf

comparison liboctave/UMFPACK/UMFPACK/Source/umf_lsolve.c @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[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_lsolve =========================================================== */
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 Lx = b, where L is the lower 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_lsolve
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 Entry xk ;
25 Entry *xp, *Lval ;
26 Int k, deg, *ip, j, row, *Lpos, *Lilen, *Lip, llen, lp, newLchain,
27 pos, npiv, n1, *Li ;
28
29 /* ---------------------------------------------------------------------- */
30
31 if (Numeric->n_row != Numeric->n_col) return (0.) ;
32 npiv = Numeric->npiv ;
33 Lpos = Numeric->Lpos ;
34 Lilen = Numeric->Lilen ;
35 Lip = Numeric->Lip ;
36 n1 = Numeric->n1 ;
37
38 #ifndef NDEBUG
39 DEBUG4 (("Lsolve start:\n")) ;
40 for (j = 0 ; j < Numeric->n_row ; j++)
41 {
42 DEBUG4 (("Lsolve start "ID": ", j)) ;
43 EDEBUG4 (X [j]) ;
44 DEBUG4 (("\n")) ;
45 }
46 #endif
47
48 /* ---------------------------------------------------------------------- */
49 /* singletons */
50 /* ---------------------------------------------------------------------- */
51
52 for (k = 0 ; k < n1 ; k++)
53 {
54 DEBUG4 (("Singleton k "ID"\n", k)) ;
55 xk = X [k] ;
56 deg = Lilen [k] ;
57 if (deg > 0 && IS_NONZERO (xk))
58 {
59 lp = Lip [k] ;
60 Li = (Int *) (Numeric->Memory + lp) ;
61 lp += UNITS (Int, deg) ;
62 Lval = (Entry *) (Numeric->Memory + lp) ;
63 for (j = 0 ; j < deg ; j++)
64 {
65 DEBUG4 ((" row "ID" k "ID" value", Li [j], k)) ;
66 EDEBUG4 (Lval [j]) ;
67 DEBUG4 (("\n")) ;
68 /* X [Li [j]] -= xk * Lval [j] ; */
69 MULT_SUB (X [Li [j]], xk, Lval [j]) ;
70 }
71 }
72 }
73
74 /* ---------------------------------------------------------------------- */
75 /* rest of L */
76 /* ---------------------------------------------------------------------- */
77
78 deg = 0 ;
79
80 for (k = n1 ; k < npiv ; k++)
81 {
82
83 /* ------------------------------------------------------------------ */
84 /* make column of L in Pattern [0..deg-1] */
85 /* ------------------------------------------------------------------ */
86
87 lp = Lip [k] ;
88 newLchain = (lp < 0) ;
89 if (newLchain)
90 {
91 lp = -lp ;
92 deg = 0 ;
93 DEBUG4 (("start of chain for column of L\n")) ;
94 }
95
96 /* remove pivot row */
97 pos = Lpos [k] ;
98 if (pos != EMPTY)
99 {
100 DEBUG4 ((" k "ID" removing row "ID" at position "ID"\n",
101 k, Pattern [pos], pos)) ;
102 ASSERT (!newLchain) ;
103 ASSERT (deg > 0) ;
104 ASSERT (pos >= 0 && pos < deg) ;
105 ASSERT (Pattern [pos] == k) ;
106 Pattern [pos] = Pattern [--deg] ;
107 }
108
109 /* concatenate the pattern */
110 ip = (Int *) (Numeric->Memory + lp) ;
111 llen = Lilen [k] ;
112 for (j = 0 ; j < llen ; j++)
113 {
114 row = *ip++ ;
115 DEBUG4 ((" row "ID" k "ID"\n", row, k)) ;
116 ASSERT (row > k) ;
117 Pattern [deg++] = row ;
118 }
119
120 /* ------------------------------------------------------------------ */
121 /* use column k of L */
122 /* ------------------------------------------------------------------ */
123
124 xk = X [k] ;
125 if (IS_NONZERO (xk))
126 {
127 xp = (Entry *) (Numeric->Memory + lp + UNITS (Int, llen)) ;
128 for (j = 0 ; j < deg ; j++)
129 {
130 DEBUG4 ((" row "ID" k "ID" value", Pattern [j], k)) ;
131 EDEBUG4 (*xp) ;
132 DEBUG4 (("\n")) ;
133 /* X [Pattern [j]] -= xk * (*xp) ; */
134 MULT_SUB (X [Pattern [j]], xk, *xp) ;
135 xp++ ;
136 }
137 }
138 }
139
140 #ifndef NDEBUG
141 for (j = 0 ; j < Numeric->n_row ; j++)
142 {
143 DEBUG4 (("Lsolve done "ID": ", j)) ;
144 EDEBUG4 (X [j]) ;
145 DEBUG4 (("\n")) ;
146 }
147 DEBUG4 (("Lsolve done.\n")) ;
148 #endif
149
150 return (MULTSUB_FLOPS * ((double) Numeric->lnz)) ;
151 }