1 /* ========================================================================== */
2 /* === Modify/t_cholmod_updown ============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Modify Module. Copyright (C) 2005-2006,
7 * Timothy A. Davis and William W. Hager.
8 * The CHOLMOD/Modify Module is licensed under Version 2.0 of the GNU
9 * General Public License. See gpl.txt for a text of the license.
10 * CHOLMOD is also available under other licenses; contact authors for details.
11 * http://www.suitesparse.com
12 * -------------------------------------------------------------------------- */
13
14 /* Updates/downdates the LDL' factorization, by computing a new factorization of
15 *
16 * Lnew * Dnew * Lnew' = Lold * Dold * Lold' +/- C*C'
17 *
18 * This file is not compiled separately. It is included into
19 * cholmod_updown.c. There are no user-callable routines in this file.
20 *
21 * The next include statements, below, create the numerical update/downdate
22 * kernels from t_cholmod_updown_numkr.c. There are 4 compiled versions of this
23 * file, one for each value of WDIM in the set 1, 2, 4, and 8. Each calls
24 * multiple versions of t_cholmod_updown_numkr; the number of versions of each
25 * is equal to WDIM. Each t_cholmod_updown_numkr version is included as a
26 * static function within its t_cholmod_updown.c caller routine. Thus:
27 *
28 * t*_updown.c creates these versions of t_cholmod_updown_numkr.c:
29 * --------- ---------------------------------------------------
30 *
31 * updown_1_r updown_1_1
32 *
33 * updown_2_r updown_2_1 updown_2_2
34 *
35 * updown_4_r updown_4_1 updown_4_2 updown_4_3 updown_4_4
36 *
37 * updown_8_r updown_8_1 updown_8_2 updown_8_3 updown_8_4
38 * updown_8_5 updown_8_6 updown_8_7 updown_8_8
39 *
40 * workspace: Xwork (nrow*wdim)
41 */
42
43 /* ========================================================================== */
44 /* === routines for numeric update/downdate along one path ================== */
45 /* ========================================================================== */
46
47 #undef FORM_NAME
48 #undef NUMERIC
49
50 #define FORM_NAME(k,rank) updown_ ## k ## _ ## rank
51 #define NUMERIC(k,rank) FORM_NAME(k,rank)
52
53 #define RANK 1
54 #include "t_cholmod_updown_numkr.c"
55
56 #if WDIM >= 2
57 #define RANK 2
58 #include "t_cholmod_updown_numkr.c"
59 #endif
60
61 #if WDIM >= 4
62 #define RANK 3
63 #include "t_cholmod_updown_numkr.c"
64 #define RANK 4
65 #include "t_cholmod_updown_numkr.c"
66 #endif
67
68 #if WDIM == 8
69 #define RANK 5
70 #include "t_cholmod_updown_numkr.c"
71 #define RANK 6
72 #include "t_cholmod_updown_numkr.c"
73 #define RANK 7
74 #include "t_cholmod_updown_numkr.c"
75 #define RANK 8
76 #include "t_cholmod_updown_numkr.c"
77 #endif
78
79
80 /* ========================================================================== */
81 /* === numeric update/downdate for all paths ================================ */
82 /* ========================================================================== */
83
NUMERIC(WDIM,r)84 static void NUMERIC (WDIM, r)
85 (
86 int update, /* TRUE for update, FALSE for downdate */
87 cholmod_sparse *C, /* in packed or unpacked, and sorted form */
88 /* no empty columns */
89 Int rank, /* rank of the update/downdate */
90 cholmod_factor *L, /* with unit diagonal (diagonal not stored) */
91 /* temporary workspaces: */
92 double W [ ], /* n-by-WDIM dense matrix, initially zero */
93 Path_type Path [ ],
94 Int npaths,
95 Int mask [ ], /* size n */
96 Int maskmark,
97 cholmod_common *Common
98 )
99 {
100 double Alpha [8] ;
101 double *Cx, *Wpath, *W1, *a ;
102 Int i, j, p, ccol, pend, wfirst, e, path, packed ;
103 Int *Ci, *Cp, *Cnz ;
104
105 /* ---------------------------------------------------------------------- */
106 /* get inputs */
107 /* ---------------------------------------------------------------------- */
108
109 Ci = C->i ;
110 Cx = C->x ;
111 Cp = C->p ;
112 Cnz = C->nz ;
113 packed = C->packed ;
114 ASSERT (IMPLIES (!packed, Cnz != NULL)) ;
115 ASSERT (L->n == C->nrow) ;
116 DEBUG (CHOLMOD(dump_real) ("num_d: in W:", W, WDIM, L->n, FALSE, 1,Common));
117
118 /* ---------------------------------------------------------------------- */
119 /* scatter C into W */
120 /* ---------------------------------------------------------------------- */
121
122 for (path = 0 ; path < rank ; path++)
123 {
124 /* W (:, path) = C (:, Path [path].col) */
125 ccol = Path [path].ccol ;
126 Wpath = W + path ;
127 PRINT1 (("Ordered Columns [path = "ID"] = "ID"\n", path, ccol)) ;
128 p = Cp [ccol] ;
129 pend = (packed) ? (Cp [ccol+1]) : (p + Cnz [ccol]) ;
130 /* column C can be empty */
131 for ( ; p < pend ; p++)
132 {
133 i = Ci [p] ;
134 ASSERT (i >= 0 && i < (Int) (C->nrow)) ;
135 if (mask == NULL || mask [i] < maskmark)
136 {
137 Wpath [WDIM * i] = Cx [p] ;
138 }
139 PRINT1 ((" row "ID" : %g mask "ID"\n", i, Cx [p],
140 (mask) ? mask [i] : 0)) ;
141 }
142 Alpha [path] = 1.0 ;
143 }
144 DEBUG (CHOLMOD(dump_real) ("num_d: W:", W, WDIM, L->n, FALSE, 1,Common)) ;
145
146 /* ---------------------------------------------------------------------- */
147 /* numeric update/downdate of the paths */
148 /* ---------------------------------------------------------------------- */
149
150 /* for each disjoint subpath in Tbar in DFS order do */
151 for (path = rank ; path < npaths ; path++)
152 {
153
154 /* determine which columns of W to use */
155 wfirst = Path [path].wfirst ;
156 e = Path [path].end ;
157 j = Path [path].start ;
158 ASSERT (e >= 0 && e < (Int) (L->n)) ;
159 ASSERT (j >= 0 && j < (Int) (L->n)) ;
160
161 W1 = W + wfirst ; /* pointer to row 0, column wfirst of W */
162 a = Alpha + wfirst ; /* pointer to Alpha [wfirst] */
163
164 PRINT1 (("Numerical update/downdate of path "ID"\n", path)) ;
165 PRINT1 (("start "ID" end "ID" wfirst "ID" rank "ID" ccol "ID"\n", j, e,
166 wfirst, Path [path].rank, Path [path].ccol)) ;
167
168 #if WDIM == 1
169 NUMERIC (WDIM,1) (update, j, e, a, W1, L, Common) ;
170 #else
171
172 switch (Path [path].rank)
173 {
174 case 1:
175 NUMERIC (WDIM,1) (update, j, e, a, W1, L, Common) ;
176 break ;
177
178 #if WDIM >= 2
179 case 2:
180 NUMERIC (WDIM,2) (update, j, e, a, W1, L, Common) ;
181 break ;
182 #endif
183
184 #if WDIM >= 4
185 case 3:
186 NUMERIC (WDIM,3) (update, j, e, a, W1, L, Common) ;
187 break ;
188 case 4:
189 NUMERIC (WDIM,4) (update, j, e, a, W1, L, Common) ;
190 break ;
191 #endif
192
193 #if WDIM == 8
194 case 5:
195 NUMERIC (WDIM,5) (update, j, e, a, W1, L, Common) ;
196 break ;
197 case 6:
198 NUMERIC (WDIM,6) (update, j, e, a, W1, L, Common) ;
199 break ;
200 case 7:
201 NUMERIC (WDIM,7) (update, j, e, a, W1, L, Common) ;
202 break ;
203 case 8:
204 NUMERIC (WDIM,8) (update, j, e, a, W1, L, Common) ;
205 break ;
206 #endif
207
208 }
209 #endif
210
211 }
212 }
213
214 /* prepare for the next inclusion of this file in cholmod_updown.c */
215 #undef WDIM
216