1 /* ========================================================================== */
2 /* === Tcov/aug ============================================================= */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Tcov Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * -------------------------------------------------------------------------- */
9
10 /* Create the augmented system S = [-I A' ; A alpha*I], solve Sx=b, and return
11 * the residual. The system solved is (alpha*I + A*A')*x=b.
12 *
13 * r1 = norm (Sx-b)
14 * r2 = norm ((alpha*I+AA')x-b)
15 * alpha = norm(A)
16 */
17
18 #include "cm.h"
19
20
21 /* ========================================================================== */
22 /* === aug ================================================================== */
23 /* ========================================================================== */
24
aug(cholmod_sparse * A)25 double aug (cholmod_sparse *A)
26 {
27 double r, maxerr = 0, bnorm, anorm ;
28 cholmod_sparse *S, *Im, *In, *At, *A1, *A2, *Sup ;
29 cholmod_dense *Alpha, *B, *Baug, *X, *W1, *W2, *R, *X2, X2mat ;
30 cholmod_factor *L ;
31 double *b, *baug, *rx, *w, *x ;
32 Int nrow, ncol, nrhs, i, j, d, d2, save, save2, save3 ;
33
34 if (A == NULL)
35 {
36 ERROR (CHOLMOD_INVALID, "cm: no A for aug") ;
37 return (1) ;
38 }
39
40 if (A->xtype != CHOLMOD_REAL)
41 {
42 return (0) ;
43 }
44
45 /* ---------------------------------------------------------------------- */
46 /* A is m-by-n, B must be m-by-nrhs */
47 /* ---------------------------------------------------------------------- */
48
49 nrow = A->nrow ;
50 ncol = A->ncol ;
51 B = rhs (A, 5, A->nrow + 7) ;
52
53 /* ---------------------------------------------------------------------- */
54 /* create scalars */
55 /* ---------------------------------------------------------------------- */
56
57 bnorm = CHOLMOD(norm_dense) (B, 0, cm) ;
58 anorm = CHOLMOD(norm_sparse) (A, 1, cm) ;
59
60 Alpha = CHOLMOD(eye) (1, 1, CHOLMOD_REAL, cm) ;
61 if (Alpha != NULL)
62 {
63 ((double *) (Alpha->x)) [0] = anorm ;
64 }
65
66 CHOLMOD(print_dense) (M1, "MinusOne", cm) ;
67 CHOLMOD(print_dense) (Alpha, "Alpha = norm(A)", cm) ;
68
69 /* ---------------------------------------------------------------------- */
70 /* create augmented system, S = [-I A' ; A anorm*I] */
71 /* ---------------------------------------------------------------------- */
72
73 Im = CHOLMOD(speye) (nrow, nrow, CHOLMOD_REAL, cm) ;
74 In = CHOLMOD(speye) (ncol, ncol, CHOLMOD_REAL, cm) ;
75 CHOLMOD(scale) (Alpha, CHOLMOD_SCALAR, Im, cm) ;
76 CHOLMOD(scale) (M1, CHOLMOD_SCALAR, In, cm) ;
77 At = CHOLMOD(transpose) (A, 2, cm) ;
78
79 /* use one of two equivalent methods */
80 if (nrow % 2)
81 {
82 /* S = [[-In A'] ; [A alpha*Im]] */
83 A1 = CHOLMOD(horzcat) (In, At, TRUE, cm) ;
84 A2 = CHOLMOD(horzcat) (A, Im, TRUE, cm) ;
85 S = CHOLMOD(vertcat) (A1, A2, TRUE, cm) ;
86 }
87 else
88 {
89 /* S = [[-In ; A] [A' ; alpha*Im]] */
90 A1 = CHOLMOD(vertcat) (In, A, TRUE, cm) ;
91 A2 = CHOLMOD(vertcat) (At, Im, TRUE, cm) ;
92 S = CHOLMOD(horzcat) (A1, A2, TRUE, cm) ;
93 }
94
95 CHOLMOD(free_sparse) (&Im, cm) ;
96 CHOLMOD(free_sparse) (&In, cm) ;
97
98 CHOLMOD(print_sparse) (S, "S, augmented system", cm) ;
99
100 /* make a symmetric (upper) copy of S */
101 Sup = CHOLMOD(copy) (S, 1, 1, cm) ;
102
103 CHOLMOD(print_sparse) (S, "S, augmented system (upper)", cm) ;
104 CHOLMOD(print_sparse) (Sup, "Sup", cm) ;
105
106 /* ---------------------------------------------------------------------- */
107 /* create augmented right-hand-side, Baug = [ zeros(ncol,nrhs) ; B ] */
108 /* ---------------------------------------------------------------------- */
109
110 b = NULL ;
111 d = 0 ;
112 nrhs = 0 ;
113 d2 = 0 ;
114 if (B != NULL)
115 {
116 nrhs = B->ncol ;
117 d = B->d ;
118 b = B->x ;
119 Baug = CHOLMOD(zeros) (nrow+ncol, nrhs, CHOLMOD_REAL, cm) ;
120 if (Baug != NULL)
121 {
122 d2 = Baug->d ;
123 baug = Baug->x ;
124 for (j = 0 ; j < nrhs ; j++)
125 {
126 for (i = 0 ; i < nrow ; i++)
127 {
128 baug [(i+ncol)+j*d2] = b [i+j*d] ;
129 }
130 }
131 }
132 }
133 else
134 {
135 Baug = NULL ;
136 }
137
138 /* ---------------------------------------------------------------------- */
139 /* solve Sx=baug */
140 /* ---------------------------------------------------------------------- */
141
142 /* S is symmetric indefinite, so do not use a supernodal LL' */
143 save = cm->supernodal ;
144 save2 = cm->final_asis ;
145 cm->supernodal = CHOLMOD_SIMPLICIAL ;
146 cm->final_asis = TRUE ;
147 save3 = cm->metis_memory ;
148 cm->metis_memory = 2.0 ;
149 L = CHOLMOD(analyze) (Sup, cm) ;
150 CHOLMOD(factorize) (Sup, L, cm) ;
151 X = CHOLMOD(solve) (CHOLMOD_A, L, Baug, cm) ;
152 cm->supernodal = save ;
153 cm->final_asis = save2 ;
154 cm->metis_memory = save3 ;
155
156 /* ---------------------------------------------------------------------- */
157 /* compute the residual */
158 /* ---------------------------------------------------------------------- */
159
160 r = resid (Sup, X, Baug) ;
161 MAXERR (maxerr, r, 1) ;
162
163 /* ---------------------------------------------------------------------- */
164 /* create a shallow submatrix of X, X2 = X (ncol:end, :) */
165 /* ---------------------------------------------------------------------- */
166
167 if (X == NULL)
168 {
169 X2 = NULL ;
170 }
171 else
172 {
173 X2 = &X2mat ;
174 X2->nrow = nrow ;
175 X2->ncol = nrhs ;
176 X2->nzmax = X->nzmax ;
177 X2->d = X->d ;
178 X2->x = ((double *) X->x) + ncol ;
179 X2->z = NULL ;
180 X2->xtype = X->xtype ;
181 X2->dtype = X->dtype ;
182 }
183
184 CHOLMOD(print_dense) (X, "X", cm) ;
185 CHOLMOD(print_dense) (X2, "X2 = X (ncol:end,:)", cm) ;
186
187 /* ---------------------------------------------------------------------- */
188 /* compute norm ((alpha*I + A*A')*x-b) */
189 /* ---------------------------------------------------------------------- */
190
191 /* W1 = A'*X2 */
192 W1 = CHOLMOD(zeros) (ncol, nrhs, CHOLMOD_REAL, cm) ;
193 CHOLMOD(sdmult) (A, TRUE, one, zero, X2, W1, cm) ;
194
195 /* W2 = A*W1 */
196 W2 = CHOLMOD(zeros) (nrow, nrhs, CHOLMOD_REAL, cm) ;
197 CHOLMOD(sdmult) (A, FALSE, one, zero, W1, W2, cm) ;
198
199 /* R = alpha*x + w2 - b */
200 R = CHOLMOD(zeros) (nrow, nrhs, CHOLMOD_REAL, cm) ;
201
202 if (R != NULL && W2 != NULL && X != NULL)
203 {
204 w = W2->x ;
205 rx = R->x ;
206 x = X2->x ;
207 for (j = 0 ; j < nrhs ; j++)
208 {
209 for (i = 0 ; i < nrow ; i++)
210 {
211 rx [i+j*nrow] = anorm * x [i+j*d2] + w [i+j*nrow] - b [i+j*d] ;
212 }
213 }
214 }
215
216 r = CHOLMOD(norm_dense) (R, 1, cm) ;
217 MAXERR (maxerr, r, bnorm) ;
218
219 /* ---------------------------------------------------------------------- */
220 /* free everything */
221 /* ---------------------------------------------------------------------- */
222
223 CHOLMOD(free_sparse) (&At, cm) ;
224 CHOLMOD(free_sparse) (&A1, cm) ;
225 CHOLMOD(free_sparse) (&A2, cm) ;
226 CHOLMOD(free_sparse) (&S, cm) ;
227 CHOLMOD(free_sparse) (&Sup, cm) ;
228 CHOLMOD(free_factor) (&L, cm) ;
229 CHOLMOD(free_dense) (&R, cm) ;
230 CHOLMOD(free_dense) (&W1, cm) ;
231 CHOLMOD(free_dense) (&W2, cm) ;
232 CHOLMOD(free_dense) (&B, cm) ;
233 CHOLMOD(free_dense) (&Baug, cm) ;
234 CHOLMOD(free_dense) (&X, cm) ;
235 CHOLMOD(free_dense) (&Alpha, cm) ;
236
237 progress (0, '.') ;
238 return (maxerr) ;
239 }
240