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