← Back to branch summary

~ubuntu-branches/ubuntu/natty/suitesparse/natty

« back to all changes in this revision

Viewing changes to CXSparse/Demo/cs_di_demo.c

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2006-12-22 10:16:15 UTC
  • Revision ID: james.westby@ubuntu.com-20061222101615-2ohaj8902oix2rnk
Tags: upstream-2.3.1
Import upstream version 2.3.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
CXSparse/Demo/cs_di_demo.c
 
1
#include "cs_di_demo.h"
 
2
#include <time.h>
 
3
/* 1 if A is square & upper tri., -1 if square & lower tri., 0 otherwise */
 
4
static int is_sym (cs_di *A)
 
5
{
 
6
    int is_upper, is_lower, j, p, n = A->n, m = A->m, *Ap = A->p, *Ai = A->i ;
 
7
    if (m != n) return (0) ;
 
8
    is_upper = 1 ;
 
9
    is_lower = 1 ;
 
10
    for (j = 0 ; j < n ; j++)
 
11
    {
 
12
        for (p = Ap [j] ; p < Ap [j+1] ; p++)
 
13
        {
 
14
            if (Ai [p] > j) is_upper = 0 ;
 
15
            if (Ai [p] < j) is_lower = 0 ;
 
16
        }
 
17
    }
 
18
    return (is_upper ? 1 : (is_lower ? -1 : 0)) ;
 
19
}
 
20
 
 
21
/* true for off-diagonal entries */
 
22
static int dropdiag (int i, int j, double aij, void *other) { return (i != j) ;}
 
23
 
 
24
/* C = A + triu(A,1)' */
 
25
static cs_di *make_sym (cs_di *A)
 
26
{
 
27
    cs_di *AT, *C ;
 
28
    AT = cs_di_transpose (A, 1) ;               /* AT = A' */
 
29
    cs_di_fkeep (AT, &dropdiag, NULL) ; /* drop diagonal entries from AT */
 
30
    C = cs_di_add (A, AT, 1, 1) ;               /* C = A+AT */
 
31
    cs_di_spfree (AT) ;
 
32
    return (C) ;
 
33
}
 
34
 
 
35
/* create a right-hand side */
 
36
static void rhs (double *x, double *b, int m)
 
37
{
 
38
    int i ;
 
39
    for (i = 0 ; i < m ; i++) b [i] = 1 + ((double) i) / m ;
 
40
    for (i = 0 ; i < m ; i++) x [i] = b [i] ;
 
41
}
 
42
 
 
43
/* infinity-norm of x */
 
44
static double norm (double *x, int n)
 
45
{
 
46
    int i ;
 
47
    double normx = 0 ;
 
48
    for (i = 0 ; i < n ; i++) normx = CS_MAX (normx, fabs (x [i])) ;
 
49
    return (normx) ;
 
50
}
 
51
 
 
52
/* compute residual, norm(A*x-b,inf) / (norm(A,1)*norm(x,inf) + norm(b,inf)) */
 
53
static void print_resid (int ok, cs_di *A, double *x, double *b, double *resid)
 
54
{
 
55
    int i, m, n ;
 
56
    if (!ok) { printf ("    (failed)\n") ; return ; }
 
57
    m = A->m ; n = A->n ;
 
58
    for (i = 0 ; i < m ; i++) resid [i] = -b [i] ;  /* resid = -b */
 
59
    cs_di_gaxpy (A, x, resid) ;                     /* resid = resid + A*x  */
 
60
    printf ("resid: %8.2e\n", norm (resid,m) / ((n == 0) ? 1 :
 
61
        (cs_di_norm (A) * norm (x,n) + norm (b,m)))) ;
 
62
}
 
63
 
 
64
static double tic (void) { return (clock () / (double) CLOCKS_PER_SEC) ; }
 
65
static double toc (double t) { double s = tic () ; return (CS_MAX (0, s-t)) ; }
 
66
 
 
67
static void print_order (int order)
 
68
{
 
69
    switch (order)
 
70
    {
 
71
        case 0: printf ("natural    ") ; break ;
 
72
        case 1: printf ("amd(A+A')  ") ; break ;
 
73
        case 2: printf ("amd(S'*S)  ") ; break ;
 
74
        case 3: printf ("amd(A'*A)  ") ; break ;
 
75
    }
 
76
}
 
77
 
 
78
/* read a problem from a file */
 
79
problem *get_problem (FILE *f, double tol)
 
80
{
 
81
    cs_di *T, *A, *C ;
 
82
    int sym, m, n, mn, nz1, nz2 ;
 
83
    problem *Prob ;
 
84
    Prob = cs_di_calloc (1, sizeof (problem)) ;
 
85
    if (!Prob) return (NULL) ;
 
86
    T = cs_di_load (f) ;                        /* load triplet matrix T from a file */
 
87
    Prob->A = A = cs_di_compress (T) ;  /* A = compressed-column form of T */
 
88
    cs_di_spfree (T) ;                  /* clear T */
 
89
    if (!cs_di_dupl (A)) return (free_problem (Prob)) ; /* sum up duplicates */
 
90
    Prob->sym = sym = is_sym (A) ;      /* determine if A is symmetric */
 
91
    m = A->m ; n = A->n ;
 
92
    mn = CS_MAX (m,n) ;
 
93
    nz1 = A->p [n] ;
 
94
    cs_di_dropzeros (A) ;                       /* drop zero entries */
 
95
    nz2 = A->p [n] ;
 
96
    if (tol > 0) cs_di_droptol (A, tol) ;       /* drop tiny entries (just to test) */
 
97
    Prob->C = C = sym ? make_sym (A) : A ;  /* C = A + triu(A,1)', or C=A */
 
98
    if (!C) return (free_problem (Prob)) ;
 
99
    printf ("\n--- Matrix: %d-by-%d, nnz: %d (sym: %d: nnz %d), norm: %8.2e\n",
 
100
            m, n, A->p [n], sym, sym ? C->p [n] : 0, cs_di_norm (C)) ;
 
101
    if (nz1 != nz2) printf ("zero entries dropped: %d\n", nz1 - nz2) ;
 
102
    if (nz2 != A->p [n]) printf ("tiny entries dropped: %d\n", nz2 - A->p [n]) ;
 
103
    Prob->b = cs_di_malloc (mn, sizeof (double)) ;
 
104
    Prob->x = cs_di_malloc (mn, sizeof (double)) ;
 
105
    Prob->resid = cs_di_malloc (mn, sizeof (double)) ;
 
106
    return ((!Prob->b || !Prob->x || !Prob->resid) ? free_problem (Prob) : Prob) ;
 
107
}
 
108
 
 
109
/* free a problem */
 
110
problem *free_problem (problem *Prob)
 
111
{
 
112
    if (!Prob) return (NULL) ;
 
113
    cs_di_spfree (Prob->A) ;
 
114
    if (Prob->sym) cs_di_spfree (Prob->C) ;
 
115
    cs_di_free (Prob->b) ;
 
116
    cs_di_free (Prob->x) ;
 
117
    cs_di_free (Prob->resid) ;
 
118
    return (cs_di_free (Prob)) ;
 
119
}
 
120
 
 
121
/* solve a linear system using Cholesky, LU, and QR, with various orderings */
 
122
int demo2 (problem *Prob)
 
123
{
 
124
    cs_di *A, *C ;
 
125
    double *b, *x, *resid,  t, tol ;
 
126
    int k, m, n, ok, order, nb, ns, *r, *s, *rr, sprank ;
 
127
    cs_did *D ;
 
128
    if (!Prob) return (0) ;
 
129
    A = Prob->A ; C = Prob->C ; b = Prob->b ; x = Prob->x ; resid = Prob->resid;
 
130
    m = A->m ; n = A->n ;
 
131
    tol = Prob->sym ? 0.001 : 1 ;               /* partial pivoting tolerance */
 
132
    D = cs_di_dmperm (C, 1) ;                   /* randomized dmperm analysis */
 
133
    if (!D) return (0) ;
 
134
    nb = D->nb ; r = D->r ; s = D->s ; rr = D->rr ;
 
135
    sprank = rr [3] ;
 
136
    for (ns = 0, k = 0 ; k < nb ; k++)
 
137
    {
 
138
        ns += ((r [k+1] == r [k]+1) && (s [k+1] == s [k]+1)) ;
 
139
    }
 
140
    printf ("blocks: %d singletons: %d structural rank: %d\n", nb, ns, sprank) ;
 
141
    cs_di_dfree (D) ;
 
142
    for (order = 0 ; order <= 3 ; order += 3)   /* natural and amd(A'*A) */
 
143
    {
 
144
        if (!order && m > 1000) continue ;
 
145
        printf ("QR   ") ;
 
146
        print_order (order) ;
 
147
        rhs (x, b, m) ;                         /* compute right-hand side */
 
148
        t = tic () ;
 
149
        ok = cs_di_qrsol (order, C, x) ;                /* min norm(Ax-b) with QR */
 
150
        printf ("time: %8.2f ", toc (t)) ;
 
151
        print_resid (ok, C, x, b, resid) ;      /* print residual */
 
152
    }
 
153
    if (m != n || sprank < n) return (1) ;      /* return if rect. or singular*/
 
154
    for (order = 0 ; order <= 3 ; order++)      /* try all orderings */
 
155
    {
 
156
        if (!order && m > 1000) continue ;
 
157
        printf ("LU   ") ;
 
158
        print_order (order) ;
 
159
        rhs (x, b, m) ;                         /* compute right-hand side */
 
160
        t = tic () ;
 
161
        ok = cs_di_lusol (order, C, x, tol) ;   /* solve Ax=b with LU */
 
162
        printf ("time: %8.2f ", toc (t)) ;
 
163
        print_resid (ok, C, x, b, resid) ;      /* print residual */
 
164
    }
 
165
    if (!Prob->sym) return (1) ;
 
166
    for (order = 0 ; order <= 1 ; order++)      /* natural and amd(A+A') */
 
167
    {
 
168
        if (!order && m > 1000) continue ;
 
169
        printf ("Chol ") ;
 
170
        print_order (order) ;
 
171
        rhs (x, b, m) ;                         /* compute right-hand side */
 
172
        t = tic () ;
 
173
        ok = cs_di_cholsol (order, C, x) ;              /* solve Ax=b with Cholesky */
 
174
        printf ("time: %8.2f ", toc (t)) ;
 
175
        print_resid (ok, C, x, b, resid) ;      /* print residual */
 
176
    }
 
177
    return (1) ;
 
178
} 
 
179
 
 
180
/* free workspace for demo3 */
 
181
static int done3 (int ok, cs_dis *S, cs_din *N, double *y, cs_di *W, cs_di *E, int *p)
 
182
{
 
183
    cs_di_sfree (S) ;
 
184
    cs_di_nfree (N) ;
 
185
    cs_di_free (y) ;
 
186
    cs_di_spfree (W) ;
 
187
    cs_di_spfree (E) ;
 
188
    cs_di_free (p) ;
 
189
    return (ok) ;
 
190
}
 
191
 
 
192
/* Cholesky update/downdate */
 
193
int demo3 (problem *Prob)
 
194
{
 
195
    cs_di *A, *C, *W = NULL, *WW, *WT, *E = NULL, *W2 ;
 
196
    int n, k, *Li, *Lp, *Wi, *Wp, p1, p2, *p = NULL, ok ;
 
197
    double *b, *x, *resid, *y = NULL, *Lx, *Wx, s,  t, t1 ;
 
198
    cs_dis *S = NULL ;
 
199
    cs_din *N = NULL ;
 
200
    if (!Prob || !Prob->sym || Prob->A->n == 0) return (0) ;
 
201
    A = Prob->A ; C = Prob->C ; b = Prob->b ; x = Prob->x ; resid = Prob->resid;
 
202
    n = A->n ;
 
203
    if (!Prob->sym || n == 0) return (1) ;
 
204
    rhs (x, b, n) ;                             /* compute right-hand side */
 
205
    printf ("\nchol then update/downdate ") ;
 
206
    print_order (1) ;
 
207
    y = cs_di_malloc (n, sizeof (double)) ;
 
208
    t = tic () ;
 
209
    S = cs_di_schol (1, C) ;                    /* symbolic Chol, amd(A+A') */
 
210
    printf ("\nsymbolic chol time %8.2f\n", toc (t)) ;
 
211
    t = tic () ;
 
212
    N = cs_di_chol (C, S) ;                     /* numeric Cholesky */
 
213
    printf ("numeric  chol time %8.2f\n", toc (t)) ;
 
214
    if (!S || !N || !y) return (done3 (0, S, N, y, W, E, p)) ;
 
215
    t = tic () ;
 
216
    cs_di_ipvec (S->pinv, b, y, n) ;            /* y = P*b */
 
217
    cs_di_lsolve (N->L, y) ;                    /* y = L\y */
 
218
    cs_di_ltsolve (N->L, y) ;                   /* y = L'\y */
 
219
    cs_di_pvec (S->pinv, y, x, n) ;             /* x = P'*y */
 
220
    printf ("solve    chol time %8.2f\n", toc (t)) ;
 
221
    printf ("original: ") ;
 
222
    print_resid (1, C, x, b, resid) ;           /* print residual */
 
223
    k = n/2 ;                                   /* construct W  */
 
224
    W = cs_di_spalloc (n, 1, n, 1, 0) ;
 
225
    if (!W) return (done3 (0, S, N, y, W, E, p)) ;
 
226
    Lp = N->L->p ; Li = N->L->i ; Lx = N->L->x ;
 
227
    Wp = W->p ; Wi = W->i ; Wx = W->x ;
 
228
    Wp [0] = 0 ;
 
229
    p1 = Lp [k] ;
 
230
    Wp [1] = Lp [k+1] - p1 ;
 
231
    s = Lx [p1] ;
 
232
    srand (1) ;
 
233
    for ( ; p1 < Lp [k+1] ; p1++)
 
234
    {
 
235
        p2 = p1 - Lp [k] ;
 
236
        Wi [p2] = Li [p1] ;
 
237
        Wx [p2] = s * rand () / ((double) RAND_MAX) ;
 
238
    }
 
239
    t = tic () ;
 
240
    ok = cs_di_updown (N->L, +1, W, S->parent) ;        /* update: L*L'+W*W' */
 
241
    t1 = toc (t) ;
 
242
    printf ("update:   time: %8.2f\n", t1) ;
 
243
    if (!ok) return (done3 (0, S, N, y, W, E, p)) ;
 
244
    t = tic () ;
 
245
    cs_di_ipvec (S->pinv, b, y, n) ;            /* y = P*b */
 
246
    cs_di_lsolve (N->L, y) ;                    /* y = L\y */
 
247
    cs_di_ltsolve (N->L, y) ;                   /* y = L'\y */
 
248
    cs_di_pvec (S->pinv, y, x, n) ;             /* x = P'*y */
 
249
    t = toc (t) ;
 
250
    p = cs_di_pinv (S->pinv, n) ;
 
251
    W2 = cs_di_permute (W, p, NULL, 1) ;                /* E = C + (P'W)*(P'W)' */
 
252
    WT = cs_di_transpose (W2,1) ;
 
253
    WW = cs_di_multiply (W2, WT) ;
 
254
    cs_di_spfree (WT) ;
 
255
    cs_di_spfree (W2) ;
 
256
    E = cs_di_add (C, WW, 1, 1) ;
 
257
    cs_di_spfree (WW) ;
 
258
    if (!E || !p) return (done3 (0, S, N, y, W, E, p)) ;
 
259
    printf ("update:   time: %8.2f (incl solve) ", t1+t) ;
 
260
    print_resid (1, E, x, b, resid) ;           /* print residual */
 
261
    cs_di_nfree (N) ;                           /* clear N */
 
262
    t = tic () ;
 
263
    N = cs_di_chol (E, S) ;                     /* numeric Cholesky */
 
264
    if (!N) return (done3 (0, S, N, y, W, E, p)) ;
 
265
    cs_di_ipvec (S->pinv, b, y, n) ;            /* y = P*b */
 
266
    cs_di_lsolve (N->L, y) ;                    /* y = L\y */
 
267
    cs_di_ltsolve (N->L, y) ;                   /* y = L'\y */
 
268
    cs_di_pvec (S->pinv, y, x, n) ;             /* x = P'*y */
 
269
    t = toc (t) ;
 
270
    printf ("rechol:   time: %8.2f (incl solve) ", t) ;
 
271
    print_resid (1, E, x, b, resid) ;           /* print residual */
 
272
    t = tic () ;
 
273
    ok = cs_di_updown (N->L, -1, W, S->parent) ;        /* downdate: L*L'-W*W' */
 
274
    t1 = toc (t) ;
 
275
    if (!ok) return (done3 (0, S, N, y, W, E, p)) ;
 
276
    printf ("downdate: time: %8.2f\n", t1) ;
 
277
    t = tic () ;
 
278
    cs_di_ipvec (S->pinv, b, y, n) ;            /* y = P*b */
 
279
    cs_di_lsolve (N->L, y) ;                    /* y = L\y */
 
280
    cs_di_ltsolve (N->L, y) ;                   /* y = L'\y */
 
281
    cs_di_pvec (S->pinv, y, x, n) ;             /* x = P'*y */
 
282
    t = toc (t) ;
 
283
    printf ("downdate: time: %8.2f (incl solve) ", t1+t) ;
 
284
    print_resid (1, C, x, b, resid) ;           /* print residual */
 
285
    return (done3 (1, S, N, y, W, E, p)) ;
 
286
}