1 /* ========================================================================== */
2 /* === Cholesky/t_cholmod_solve ============================================= */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2013, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Template routine for cholmod_solve. Supports any numeric xtype (real,
10 * complex, or zomplex). The xtypes of all matrices (L and Y) must match.
11 */
12
13 #include "cholmod_template.h"
14
15 /* ========================================================================== */
16 /* === simplicial template solvers ========================================== */
17 /* ========================================================================== */
18
19 /* LL': solve Lx=b with non-unit diagonal */
20 #define LL
21 #include "t_cholmod_lsolve.c"
22
23 /* LDL': solve LDx=b */
24 #define LD
25 #include "t_cholmod_lsolve.c"
26
27 /* LDL': solve Lx=b with unit diagonal */
28 #include "t_cholmod_lsolve.c"
29
30 /* LL': solve L'x=b with non-unit diagonal */
31 #define LL
32 #include "t_cholmod_ltsolve.c"
33
34 /* LDL': solve DL'x=b */
35 #define LD
36 #include "t_cholmod_ltsolve.c"
37
38 /* LDL': solve L'x=b with unit diagonal */
39 #include "t_cholmod_ltsolve.c"
40
41
42 /* ========================================================================== */
43 /* === t_ldl_dsolve ========================================================= */
44 /* ========================================================================== */
45
46 /* Solve Dx=b for an LDL' factorization, where Y holds b' on input and x' on
47 * output.
48 *
49 * The number of right-hand-sides (nrhs) is not restricted, even if Yseti
50 * is present.
51 */
52
TEMPLATE(ldl_dsolve)53 static void TEMPLATE (ldl_dsolve)
54 (
55 cholmod_factor *L,
56 cholmod_dense *Y, /* nr-by-n with leading dimension nr */
57 Int *Yseti, Int ysetlen
58 )
59 {
60 double d [1] ;
61 double *Lx, *Yx, *Yz ;
62 Int *Lp ;
63 Int n, nrhs, k, p, k1, k2, kk, kkiters ;
64
65 ASSERT (L->xtype == Y->xtype) ; /* L and Y must have the same xtype */
66 ASSERT (L->n == Y->ncol) ; /* dimensions must match */
67 ASSERT (Y->nrow == Y->d) ; /* leading dimension of Y = # rows of Y */
68 ASSERT (L->xtype != CHOLMOD_PATTERN) ; /* L is not symbolic */
69 ASSERT (!(L->is_super) && !(L->is_ll)) ; /* L is simplicial LDL' */
70
71 nrhs = Y->nrow ;
72 n = L->n ;
73 Lp = L->p ;
74 Lx = L->x ;
75 Yx = Y->x ;
76 Yz = Y->z ;
77 kkiters = Yseti ? ysetlen : n ;
78 for (kk = 0 ; kk < kkiters ; kk++)
79 {
80 k = Yseti ? Yseti [kk] : kk ;
81 k1 = k*nrhs ;
82 k2 = (k+1)*nrhs ;
83 ASSIGN_REAL (d,0, Lx,Lp[k]) ;
84 for (p = k1 ; p < k2 ; p++)
85 {
86 DIV_REAL (Yx,Yz,p, Yx,Yz,p, d,0) ;
87 }
88 }
89 }
90
91
92 /* ========================================================================== */
93 /* === t_simplicial_solver ================================================== */
94 /* ========================================================================== */
95
96 /* Solve a linear system, where Y' contains the (array-transposed) right-hand
97 * side on input, and the solution on output. No permutations are applied;
98 * these must have already been applied to Y on input.
99 *
100 * Yseti [0..ysetlen-1] is an optional list of indices from
101 * cholmod_lsolve_pattern. The solve is performed only on the columns of L
102 * corresponding to entries in Yseti. Ignored if NULL. If present, most
103 * functions require that Y' consist of a single dense column.
104 */
105
TEMPLATE(simplicial_solver)106 static void TEMPLATE (simplicial_solver)
107 (
108 int sys, /* system to solve */
109 cholmod_factor *L, /* factor to use, a simplicial LL' or LDL' */
110 cholmod_dense *Y, /* right-hand-side on input, solution on output */
111 Int *Yseti, Int ysetlen
112 )
113 {
114 if (L->is_ll)
115 {
116 /* The factorization is LL' */
117 if (sys == CHOLMOD_A || sys == CHOLMOD_LDLt)
118 {
119 /* Solve Ax=b or LL'x=b */
120 TEMPLATE (ll_lsolve_k) (L, Y, Yseti, ysetlen) ;
121 TEMPLATE (ll_ltsolve_k) (L, Y, Yseti, ysetlen) ;
122 }
123 else if (sys == CHOLMOD_L || sys == CHOLMOD_LD)
124 {
125 /* Solve Lx=b */
126 TEMPLATE (ll_lsolve_k) (L, Y, Yseti, ysetlen) ;
127 }
128 else if (sys == CHOLMOD_Lt || sys == CHOLMOD_DLt)
129 {
130 /* Solve L'x=b */
131 TEMPLATE (ll_ltsolve_k) (L, Y, Yseti, ysetlen) ;
132 }
133 }
134 else
135 {
136 /* The factorization is LDL' */
137 if (sys == CHOLMOD_A || sys == CHOLMOD_LDLt)
138 {
139 /* Solve Ax=b or LDL'x=b */
140 TEMPLATE (ldl_lsolve_k) (L, Y, Yseti, ysetlen) ;
141 TEMPLATE (ldl_dltsolve_k) (L, Y, Yseti, ysetlen) ;
142 }
143 else if (sys == CHOLMOD_LD)
144 {
145 /* Solve LDx=b */
146 TEMPLATE (ldl_ldsolve_k) (L, Y, Yseti, ysetlen) ;
147 }
148 else if (sys == CHOLMOD_L)
149 {
150 /* Solve Lx=b */
151 TEMPLATE (ldl_lsolve_k) (L, Y, Yseti, ysetlen) ;
152 }
153 else if (sys == CHOLMOD_Lt)
154 {
155 /* Solve L'x=b */
156 TEMPLATE (ldl_ltsolve_k) (L, Y, Yseti, ysetlen) ;
157 }
158 else if (sys == CHOLMOD_DLt)
159 {
160 /* Solve DL'x=b */
161 TEMPLATE (ldl_dltsolve_k) (L, Y, Yseti, ysetlen) ;
162 }
163 else if (sys == CHOLMOD_D)
164 {
165 /* Solve Dx=b */
166 TEMPLATE (ldl_dsolve) (L, Y, Yseti, ysetlen) ;
167 }
168 }
169 }
170
171 #undef PATTERN
172 #undef REAL
173 #undef COMPLEX
174 #undef ZOMPLEX
175