← Back to branch summary

~logan/ubuntu/trusty/suitesparse/4.2.1-3ubuntu1

« back to all changes in this revision

Viewing changes to KLU/Demo/kludemo.c

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2007-05-29 09:36:29 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070529093629-zowquo0b7slkk6nc
Tags: 3.0.0-2
http://bugs.debian.org/426349
http://bugs.debian.org/426354
* suitesparse builds properly twice in a row
* Bug fix: "suitesparse - FTBFS: Broken build depens: libgfortran1-dev",
  thanks to Bastian Blank (Closes: #426349).
* Bug fix: "suitesparse_3.0.0-1: FTBFS: build-depends on
  libgfortran1-dev", thanks to Steve Langasek (Closes: #426354).

Show diffs side-by-side

added added

removed removed

Lines of Context:
KLU/Demo/kludemo.c
1
1
/* ========================================================================== */
2
 
/* === kludemo ============================================================== */
 
2
/* === KLU DEMO ============================================================= */
3
3
/* ========================================================================== */
4
4
 
5
 
/* Demo program for KLU.  Reads in a column-form 0-based matrix in tmp/Ap,
6
 
 * tmp/Ai, and tmp/Ax, whose size and # of nonzeros are in the file tmp/Asize.
7
 
 * Then calls KLU to analyze, factor, and solve the system.
8
 
 *
9
 
 * Example:
10
 
 *
11
 
 *      kludemo
12
 
 *
13
 
 * The right-hand-side can be provided in the optional tmp/b file.  The solution
14
 
 * is written to tmp/x.
15
 
 */
 
5
/* Read in a Matrix Market matrix (using CHOLMOD) and solve a linear system. */
16
6
 
 
7
#include <math.h>
17
8
#include <stdio.h>
18
 
#include <stdlib.h>
19
 
#include <limits.h>
20
 
#include <math.h>
21
9
#include "klu.h"
22
 
#include "klu_cholmod.h"
23
 
 
24
 
#define ABS(x) ((x) >= 0 ? (x) : -(x))
 
10
 
 
11
/* for handling complex matrices */
 
12
#define REAL(X,i) (X [2*(i)])
 
13
#define IMAG(X,i) (X [2*(i)+1])
 
14
#define CABS(X,i) (sqrt (REAL (X,i) * REAL (X,i) + IMAG (X,i) * IMAG (X,i)))
 
15
 
25
16
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
26
 
#define XTRUE(i,n) (1.0 + ((double) i) / ((double) n))
27
 
#define ISNAN(x) ((x) != (x))
28
 
 
29
 
#ifndef FALSE
30
 
#define FALSE 0
31
 
#endif
32
 
 
33
 
#ifndef TRUE
34
 
#define TRUE 1
35
 
#endif
36
 
 
37
 
#define NTRIAL 1
38
 
#define NRHS 24
39
 
 
40
 
#include "dsecnd.h"
41
 
 
42
 
#include <time.h>
43
 
#define CPUTIME (((double) clock ( )) / CLOCKS_PER_SEC)
44
 
 
45
 
/* -------------------------------------------------------------------------- */
46
 
/* err: compute the relative error, ||x-xtrue||/||xtrue|| */
47
 
/* -------------------------------------------------------------------------- */
48
 
 
49
 
static double err
50
 
(
51
 
    int n,
52
 
    double x [ ]
53
 
)
54
 
{
55
 
    int i  ;
56
 
    double enorm, e, abse, absxtrue, xnorm ;
57
 
    enorm = 0 ;
58
 
    xnorm = 0 ;
59
 
 
60
 
    for (i = 0 ; i < n ; i++)
61
 
    {
62
 
        if (ISNAN (x [i]))
63
 
        {
64
 
            enorm = x [i] ;
65
 
            break ;
66
 
        }
67
 
        e = x [i] - XTRUE (i,n) ;
68
 
        abse = ABS (e) ;
69
 
        enorm = MAX (enorm, abse) ;
70
 
    }
71
 
 
72
 
    for (i = 0 ; i < n ; i++)
73
 
    {
74
 
        /* XTRUE is positive, but do this in case XTRUE is redefined */
75
 
        absxtrue = ABS (XTRUE (i,n)) ;
76
 
        xnorm = MAX (xnorm, absxtrue) ;
77
 
    }
78
 
 
79
 
    if (xnorm == 0)
80
 
    {
81
 
        xnorm = 1 ;
82
 
    }
83
 
    return (enorm / xnorm) ;
84
 
}
85
 
 
86
 
 
87
 
/* -------------------------------------------------------------------------- */
88
 
/* resid: compute the relative residual, ||Ax-b||/||b|| or ||A'x-b||/||b|| */
89
 
/* -------------------------------------------------------------------------- */
90
 
 
91
 
static double resid
92
 
(
93
 
    int n,
94
 
    int Ap [ ],
95
 
    int Ai [ ],
96
 
    double Ax [ ],
97
 
    double x [ ],
98
 
    double r [ ],
99
 
    double b [ ],
100
 
    int transpose
101
 
)
102
 
{
103
 
    int i, j, p ;
104
 
    double rnorm, absr, absb, bnorm ;
105
 
    for (i = 0 ; i < n ; i++)
106
 
    {
107
 
        r [i] = 0 ;
108
 
    }
109
 
 
110
 
    if (transpose)
111
 
    {
112
 
        for (j = 0 ; j < n ; j++)
113
 
        {
114
 
            for (p = Ap [j] ; p < Ap [j+1] ; p++)
115
 
            {
116
 
                i = Ai [p] ;
117
 
                r [j] += Ax [p] * x [i] ;
118
 
            }
119
 
        }
120
 
    }
121
 
    else
122
 
    {
123
 
        for (j = 0 ; j < n ; j++)
124
 
        {
125
 
            for (p = Ap [j] ; p < Ap [j+1] ; p++)
126
 
            {
127
 
                i = Ai [p] ;
128
 
                r [i] += Ax [p] * x [j] ;
129
 
            }
130
 
        }
131
 
    }
132
 
 
133
 
    for (i = 0 ; i < n ; i++)
134
 
    {
135
 
        r [i] -= b [i] ;
136
 
    }
137
 
    rnorm = 0. ;
138
 
    bnorm = 0. ;
139
 
    for (i = 0 ; i < n ; i++)
140
 
    {
141
 
        if (ISNAN (r [i]))
142
 
        {
143
 
            rnorm = r [i] ;
144
 
            break ;
145
 
        }
146
 
        absr = ABS (r [i]) ;
147
 
        rnorm = MAX (rnorm, absr) ;
148
 
    }
149
 
    for (i = 0 ; i < n ; i++)
150
 
    {
151
 
        if (ISNAN (b [i]))
152
 
        {
153
 
            bnorm = b [i] ;
154
 
            break ;
155
 
        }
156
 
        absb = ABS (b [i]) ;
157
 
        bnorm = MAX (bnorm, absb) ;
158
 
    }
159
 
    if (bnorm == 0)
160
 
    {
161
 
        bnorm = 1 ;
162
 
    }
163
 
    return (rnorm / bnorm) ;
164
 
}
165
 
 
166
 
 
167
 
/* -------------------------------------------------------------------------- */
168
 
/* Atimesx: compute y = A*x  or A'*x, where x (i) = 1 + i/n */
169
 
/* -------------------------------------------------------------------------- */
170
 
 
171
 
static void Atimesx
172
 
(
173
 
    int n,
174
 
    int Ap [ ],
175
 
    int Ai [ ],
176
 
    double Ax [ ],
177
 
    double y [ ],
178
 
    int transpose
179
 
)
180
 
{
181
 
    int i, j, p ;
182
 
    for (i = 0 ; i < n ; i++)
183
 
    {
184
 
        y [i] = 0 ;
185
 
    }
186
 
    if (transpose)
187
 
    {
188
 
        for (j = 0 ; j < n ; j++)
189
 
        {
190
 
            for (p = Ap [j] ; p < Ap [j+1] ; p++)
191
 
            {
192
 
                i = Ai [p] ;
193
 
                y [j] += Ax [p] * XTRUE (i,n) ;
194
 
            }
195
 
        }
196
 
    }
197
 
    else
198
 
    {
199
 
        for (j = 0 ; j < n ; j++)
200
 
        {
201
 
            for (p = Ap [j] ; p < Ap [j+1] ; p++)
202
 
            {
203
 
                i = Ai [p] ;
204
 
                y [i] += Ax [p] * XTRUE (j,n) ;
205
 
            }
206
 
        }
207
 
    }
208
 
}
209
 
 
210
 
/* -------------------------------------------------------------------------- */
211
 
/* main program */
212
 
/* -------------------------------------------------------------------------- */
 
17
 
 
18
/* ========================================================================== */
 
19
/* === klu_backslash ======================================================== */
 
20
/* ========================================================================== */
 
21
 
 
22
static int klu_backslash    /* return 1 if successful, 0 otherwise */
 
23
(
 
24
    /* --- input ---- */
 
25
    int n,              /* A is n-by-n */
 
26
    int *Ap,            /* size n+1, column pointers */
 
27
    int *Ai,            /* size nz = Ap [n], row indices */
 
28
    double *Ax,         /* size nz, numerical values */
 
29
    int isreal,         /* nonzero if A is real, 0 otherwise */
 
30
    double *B,          /* size n, right-hand-side */
 
31
 
 
32
    /* --- output ---- */
 
33
    double *X,          /* size n, solution to Ax=b */
 
34
    double *R,          /* size n, residual r = b-A*x */
 
35
 
 
36
    /* --- scalar output --- */
 
37
    int *lunz,          /* nnz (L+U+F) */
 
38
    double *rnorm,      /* norm (b-A*x,1) / norm (A,1) */
 
39
 
 
40
    /* --- workspace - */
 
41
 
 
42
    klu_common *Common  /* default parameters and statistics */
 
43
)
 
44
{
 
45
    double anorm = 0, asum ;
 
46
    klu_symbolic *Symbolic ;
 
47
    klu_numeric *Numeric ;
 
48
    int i, j, p ;
 
49
 
 
50
    if (!Ap || !Ai || !Ax || !B || !X || !B) return (0) ;
 
51
 
 
52
    /* ---------------------------------------------------------------------- */
 
53
    /* symbolic ordering and analysis */
 
54
    /* ---------------------------------------------------------------------- */
 
55
 
 
56
    Symbolic = klu_analyze (n, Ap, Ai, Common) ;
 
57
    if (!Symbolic) return (0) ;
 
58
 
 
59
    if (isreal)
 
60
    {
 
61
 
 
62
        /* ------------------------------------------------------------------ */
 
63
        /* factorization */
 
64
        /* ------------------------------------------------------------------ */
 
65
 
 
66
        Numeric = klu_factor (Ap, Ai, Ax, Symbolic, Common) ;
 
67
        if (!Numeric)
 
68
        {
 
69
            klu_free_symbolic (&Symbolic, Common) ;
 
70
            return (0) ;
 
71
        }
 
72
 
 
73
        /* ------------------------------------------------------------------ */
 
74
        /* statistics (not required to solve Ax=b) */
 
75
        /* ------------------------------------------------------------------ */
 
76
 
 
77
        klu_rgrowth (Ap, Ai, Ax, Symbolic, Numeric, Common) ;
 
78
        klu_condest (Ap, Ax, Symbolic, Numeric, Common) ;
 
79
        klu_rcond (Symbolic, Numeric, Common) ;
 
80
        klu_flops (Symbolic, Numeric, Common) ;
 
81
        *lunz = Numeric->lnz + Numeric->unz - n + 
 
82
            ((Numeric->Offp) ? (Numeric->Offp [n]) : 0) ;
 
83
 
 
84
        /* ------------------------------------------------------------------ */
 
85
        /* solve Ax=b */
 
86
        /* ------------------------------------------------------------------ */
 
87
 
 
88
        for (i = 0 ; i < n ; i++)
 
89
        {
 
90
            X [i] = B [i] ;
 
91
        }
 
92
        klu_solve (Symbolic, Numeric, n, 1, X, Common) ;
 
93
 
 
94
        /* ------------------------------------------------------------------ */
 
95
        /* compute residual, rnorm = norm(b-Ax,1) / norm(A,1) */
 
96
        /* ------------------------------------------------------------------ */
 
97
 
 
98
        for (i = 0 ; i < n ; i++)
 
99
        {
 
100
            R [i] = B [i] ;
 
101
        }
 
102
        for (j = 0 ; j < n ; j++)
 
103
        {
 
104
            asum = 0 ;
 
105
            for (p = Ap [j] ; p < Ap [j+1] ; p++)
 
106
            {
 
107
                /* R (i) -= A (i,j) * X (j) */
 
108
                R [Ai [p]] -= Ax [p] * X [j] ;
 
109
                asum += fabs (Ax [p]) ;
 
110
            }
 
111
            anorm = MAX (anorm, asum) ;
 
112
        }
 
113
        *rnorm = 0 ;
 
114
        for (i = 0 ; i < n ; i++)
 
115
        {
 
116
            *rnorm = MAX (*rnorm, fabs (R [i])) ;
 
117
        }
 
118
 
 
119
        /* ------------------------------------------------------------------ */
 
120
        /* free numeric factorization */
 
121
        /* ------------------------------------------------------------------ */
 
122
 
 
123
        klu_free_numeric (&Numeric, Common) ;
 
124
 
 
125
    }
 
126
    else
 
127
    {
 
128
 
 
129
        /* ------------------------------------------------------------------ */
 
130
        /* statistics (not required to solve Ax=b) */
 
131
        /* ------------------------------------------------------------------ */
 
132
 
 
133
        Numeric = klu_z_factor (Ap, Ai, Ax, Symbolic, Common) ;
 
134
        if (!Numeric)
 
135
        {
 
136
            klu_free_symbolic (&Symbolic, Common) ;
 
137
            return (0) ;
 
138
        }
 
139
 
 
140
        /* ------------------------------------------------------------------ */
 
141
        /* statistics */
 
142
        /* ------------------------------------------------------------------ */
 
143
 
 
144
        klu_z_rgrowth (Ap, Ai, Ax, Symbolic, Numeric, Common) ;
 
145
        klu_z_condest (Ap, Ax, Symbolic, Numeric, Common) ;
 
146
        klu_z_rcond (Symbolic, Numeric, Common) ;
 
147
        klu_z_flops (Symbolic, Numeric, Common) ;
 
148
        *lunz = Numeric->lnz + Numeric->unz - n + 
 
149
            ((Numeric->Offp) ? (Numeric->Offp [n]) : 0) ;
 
150
 
 
151
        /* ------------------------------------------------------------------ */
 
152
        /* solve Ax=b */
 
153
        /* ------------------------------------------------------------------ */
 
154
 
 
155
        for (i = 0 ; i < 2*n ; i++)
 
156
        {
 
157
            X [i] = B [i] ;
 
158
        }
 
159
        klu_z_solve (Symbolic, Numeric, n, 1, X, Common) ;
 
160
 
 
161
        /* ------------------------------------------------------------------ */
 
162
        /* compute residual, rnorm = norm(b-Ax,1) / norm(A,1) */
 
163
        /* ------------------------------------------------------------------ */
 
164
 
 
165
        for (i = 0 ; i < 2*n ; i++)
 
166
        {
 
167
            R [i] = B [i] ;
 
168
        }
 
169
        for (j = 0 ; j < n ; j++)
 
170
        {
 
171
            asum = 0 ;
 
172
            for (p = Ap [j] ; p < Ap [j+1] ; p++)
 
173
            {
 
174
                /* R (i) -= A (i,j) * X (j) */
 
175
                i = Ai [p] ;
 
176
                REAL (R,i) -= REAL(Ax,p) * REAL(X,j) - IMAG(Ax,p) * IMAG(X,j) ;
 
177
                IMAG (R,i) -= IMAG(Ax,p) * REAL(X,j) + REAL(Ax,p) * IMAG(X,j) ;
 
178
                asum += CABS (Ax, p) ;
 
179
            }
 
180
            anorm = MAX (anorm, asum) ;
 
181
        }
 
182
        *rnorm = 0 ;
 
183
        for (i = 0 ; i < n ; i++)
 
184
        {
 
185
            *rnorm = MAX (*rnorm, CABS (R, i)) ;
 
186
        }
 
187
 
 
188
        /* ------------------------------------------------------------------ */
 
189
        /* free numeric factorization */
 
190
        /* ------------------------------------------------------------------ */
 
191
 
 
192
        klu_z_free_numeric (&Numeric, Common) ;
 
193
    }
 
194
 
 
195
    /* ---------------------------------------------------------------------- */
 
196
    /* free symbolic analysis, and residual */
 
197
    /* ---------------------------------------------------------------------- */
 
198
 
 
199
    klu_free_symbolic (&Symbolic, Common) ;
 
200
    return (1) ;
 
201
}
 
202
 
 
203
 
 
204
/* ========================================================================== */
 
205
/* === klu_demo ============================================================= */
 
206
/* ========================================================================== */
 
207
 
 
208
/* Given a sparse matrix A, set up a right-hand-side and solve X = A\b */
 
209
 
 
210
static void klu_demo (int n, int *Ap, int *Ai, double *Ax, int isreal)
 
211
{
 
212
    double rnorm ;
 
213
    klu_common Common ;
 
214
    double *B, *X, *R ;
 
215
    int i, lunz ;
 
216
 
 
217
    printf ("KLU: %s, version: %d.%d.%d\n", KLU_DATE, KLU_MAIN_VERSION,
 
218
        KLU_SUB_VERSION, KLU_SUBSUB_VERSION) ;
 
219
 
 
220
    /* ---------------------------------------------------------------------- */
 
221
    /* set defaults */
 
222
    /* ---------------------------------------------------------------------- */
 
223
 
 
224
    klu_defaults (&Common) ;
 
225
 
 
226
    /* ---------------------------------------------------------------------- */
 
227
    /* create a right-hand-side */
 
228
    /* ---------------------------------------------------------------------- */
 
229
 
 
230
    if (isreal)
 
231
    {
 
232
        /* B = 1 + (1:n)/n */
 
233
        B = klu_malloc (n, sizeof (double), &Common) ;
 
234
        X = klu_malloc (n, sizeof (double), &Common) ;
 
235
        R = klu_malloc (n, sizeof (double), &Common) ;
 
236
        if (B)
 
237
        {
 
238
            for (i = 0 ; i < n ; i++)
 
239
            {
 
240
                B [i] = 1 + ((double) i+1) / ((double) n) ;
 
241
            }
 
242
        }
 
243
    }
 
244
    else
 
245
    {
 
246
        /* real (B) = 1 + (1:n)/n, imag(B) = (n:-1:1)/n */
 
247
        B = klu_malloc (n, 2 * sizeof (double), &Common) ;
 
248
        X = klu_malloc (n, 2 * sizeof (double), &Common) ;
 
249
        R = klu_malloc (n, 2 * sizeof (double), &Common) ;
 
250
        if (B)
 
251
        {
 
252
            for (i = 0 ; i < n ; i++)
 
253
            {
 
254
                REAL (B, i) = 1 + ((double) i+1) / ((double) n) ;
 
255
                IMAG (B, i) = ((double) n-i) / ((double) n) ;
 
256
            }
 
257
        }
 
258
    }
 
259
 
 
260
    /* ---------------------------------------------------------------------- */
 
261
    /* X = A\b using KLU and print statistics */
 
262
    /* ---------------------------------------------------------------------- */
 
263
 
 
264
    if (!klu_backslash (n, Ap, Ai, Ax, isreal, B, X, R, &lunz, &rnorm, &Common))
 
265
    {
 
266
        printf ("KLU failed\n") ;
 
267
    }
 
268
    else
 
269
    {
 
270
        printf ("n %d nnz(A) %d nnz(L+U+F) %d resid %g\n"
 
271
            "recip growth %g condest %g rcond %g flops %g\n",
 
272
            n, Ap [n], lunz, rnorm, Common.rgrowth, Common.condest,
 
273
            Common.rcond, Common.flops) ;
 
274
    }
 
275
 
 
276
    /* ---------------------------------------------------------------------- */
 
277
    /* free the problem */
 
278
    /* ---------------------------------------------------------------------- */
 
279
 
 
280
    if (isreal)
 
281
    {
 
282
        klu_free (B, n, sizeof (double), &Common) ;
 
283
        klu_free (X, n, sizeof (double), &Common) ;
 
284
        klu_free (R, n, sizeof (double), &Common) ;
 
285
    }
 
286
    else
 
287
    {
 
288
        klu_free (B, 2*n, sizeof (double), &Common) ;
 
289
        klu_free (X, 2*n, sizeof (double), &Common) ;
 
290
        klu_free (R, 2*n, sizeof (double), &Common) ;
 
291
    }
 
292
    printf ("peak memory usage: %g bytes\n\n", (double) (Common.mempeak)) ;
 
293
}
 
294
 
 
295
 
 
296
/* ========================================================================== */
 
297
/* === main ================================================================= */
 
298
/* ========================================================================== */
 
299
 
 
300
/* Read in a sparse matrix in Matrix Market format using CHOLMOD, and then
 
301
 * solve Ax=b with KLU.  Note that CHOLMOD is only used to read the matrix. */
 
302
 
 
303
#include "cholmod.h"
213
304
 
214
305
int main (void)
215
306
{
216
 
    double *Ax, *b, *x, *r, *Y ;
217
 
    double stats [2], condest, growth, t, c, rcond ;
218
 
    /* atime [2], ftime [2], rtime [2], stime [2], stime2 [2], */
219
 
    int p, i, j, n, nz, *Ap, *Ai, nrow, ncol, rhs, trial, s, k, nrhs, do_sort,
220
 
        ldim, ordering, do_btf, scale, rescale, lnz, unz, halt_if_singular ;
221
 
    klu_symbolic *Symbolic ;
222
 
    klu_numeric  *Numeric ;
223
 
    klu_common Kommon, *km ;
224
 
    FILE *f ;
225
 
 
226
 
    printf ("\nKLU stand-alone demo:\n") ;
227
 
    km = &Kommon ;
228
 
    klu_defaults (km) ;
229
 
 
230
 
    /* ---------------------------------------------------------------------- */
231
 
    /* get n and nz */
232
 
    /* ---------------------------------------------------------------------- */
233
 
 
234
 
    printf ("File: tmp/Asize\n") ;
235
 
    f = fopen ("tmp/Asize", "r") ;
236
 
    if (!f)
237
 
    {
238
 
        printf ("Unable to open file\n") ;
239
 
        exit (1) ;
240
 
    }
241
 
    fscanf (f, "%d %d %d\n", &nrow, &ncol, &nz) ;
242
 
    fclose (f) ;
243
 
    n  = MAX (1, MAX (nrow, ncol)) ;
244
 
    nz = MAX (1, nz) ;
245
 
    ldim = n + 13 ;
246
 
    nrhs = NRHS ;
247
 
 
248
 
    /* ---------------------------------------------------------------------- */
249
 
    /* allocate space for the matrix A, and for the vectors b, r, and x */
250
 
    /* ---------------------------------------------------------------------- */
251
 
 
252
 
    Ap = (int *) malloc ((n+1) * sizeof (int)) ;
253
 
    Ai = (int *) malloc (nz * sizeof (int)) ;
254
 
    Ax = (double *) malloc (nz * sizeof (double)) ;
255
 
    b = (double *) malloc (ldim * NRHS * sizeof (double)) ;
256
 
    r = (double *) malloc (n * sizeof (double)) ;
257
 
    x = (double *) malloc (ldim * NRHS * sizeof (double)) ;
258
 
 
259
 
    if (!Ap || !Ai || !Ax || !b || !r || !x)
260
 
    {
261
 
        printf ("out of memory for input matrix\n") ;
262
 
        exit (1) ;
263
 
    }
264
 
 
265
 
    /* ---------------------------------------------------------------------- */
266
 
    /* get the matrix (tmp/Ap, tmp/Ai, and tmp/Ax) */
267
 
    /* ---------------------------------------------------------------------- */
268
 
 
269
 
    printf ("File: tmp/Ap\n") ;
270
 
    f = fopen ("tmp/Ap", "r") ;
271
 
    if (!f)
272
 
    {
273
 
        printf ("Unable to open file\n") ;
274
 
        exit (1) ;
275
 
    }
276
 
    for (j = 0 ; j <= n ; j++)
277
 
    {
278
 
        fscanf (f, "%d\n", &Ap [j]) ;
279
 
    }
280
 
    fclose (f) ;
281
 
 
282
 
    printf ("File: tmp/Ai\n") ;
283
 
    f = fopen ("tmp/Ai", "r") ;
284
 
    if (!f)
285
 
    {
286
 
        printf ("Unable to open file\n") ;
287
 
        exit (1) ;
288
 
    }
289
 
    for (p = 0 ; p < nz ; p++)
290
 
    {
291
 
        fscanf (f, "%d\n", &Ai [p]) ;
292
 
    }
293
 
    fclose (f) ;
294
 
 
295
 
    printf ("File: tmp/Ax\n") ;
296
 
    f = fopen ("tmp/Ax", "r") ;
297
 
    if (!f)
298
 
    {
299
 
        printf ("Unable to open file\n") ;
300
 
        exit (1) ;
301
 
    }
302
 
    for (p = 0 ; p < nz ; p++)
303
 
    {
304
 
        fscanf (f, "%lg\n", &Ax [p]) ;
305
 
    }
306
 
    fclose (f) ;
307
 
 
308
 
    printf ("n %d nrow %d ncol %d nz %d\n", n, nrow, ncol, nz) ;
309
 
 
310
 
    /* ---------------------------------------------------------------------- */
311
 
    /* b = A * xtrue, or read b from a file */
312
 
    /* ---------------------------------------------------------------------- */
313
 
 
314
 
    for (i = 0 ; i < ldim * NRHS ; i++) b [i] = 0 ;
315
 
 
316
 
    rhs = FALSE ;
317
 
    f = fopen ("tmp/b", "r") ;
318
 
    if (f != (FILE *) NULL)
319
 
    {
320
 
        printf ("Reading tmp/b\n") ;
321
 
        rhs = TRUE ;
322
 
        for (i = 0 ; i < n ; i++)
323
 
        {
324
 
            fscanf (f, "%lg\n", &b [i]) ;
325
 
        }
326
 
        fclose (f) ;
327
 
    }
328
 
    else
329
 
    {
330
 
        Atimesx (n, Ap, Ai, Ax, b, FALSE) ;
331
 
    }
332
 
 
333
 
    /* create b (:,2:NRHS) */
334
 
    Y = b + ldim ;
335
 
    for (s = 1 ; s < NRHS ; s++)
336
 
    {
337
 
        for (k = 0 ; k < n ; k++)
338
 
        {
339
 
            Y [k] = s + ((double) n) / 1000 ;
340
 
        }
341
 
        Y += ldim ;
342
 
    }
343
 
 
344
 
    /*
345
 
    atime [0] = 0 ;
346
 
    atime [1] = 0 ;
347
 
    ftime [0] = 0 ;
348
 
    ftime [1] = 0 ;
349
 
    rtime [0] = 0 ;
350
 
    rtime [1] = 0 ;
351
 
    stime [0] = 0 ;
352
 
    stime [1] = 0 ;
353
 
    stime2 [0] = 0 ;
354
 
    stime2 [1] = 0 ;
355
 
    */
356
 
    stats [0] = 0 ;
357
 
    stats [1] = 0 ;
358
 
 
359
 
    do_btf = 1 ;
360
 
    ordering = 0 ;
361
 
    scale = -1 ;
362
 
    halt_if_singular = TRUE ;
363
 
    do_sort = FALSE ;
364
 
 
365
 
    /*
366
 
for (do_btf = 0 ; do_btf <= 1 ; do_btf++)
367
 
{
368
 
for (scale = -1 ; scale <= 2 ; scale++)
369
 
{
370
 
for (halt_if_singular = TRUE ; halt_if_singular >= FALSE ; halt_if_singular--)
371
 
{
372
 
*/
373
 
 
374
 
for (ordering = 0 ; ordering <= 3 ; ordering++)
375
 
{
376
 
if (ordering == 2 || ordering == 1) continue ;
377
 
 
378
 
/*
379
 
for (do_sort = FALSE ; do_sort <= TRUE ; do_sort++)
380
 
{
381
 
*/
382
 
 
383
 
klu_defaults (km) ;
384
 
km->ordering = ordering ;           /* 0: amd, 1: colamd, 2: given, 3: user */
385
 
km->btf = do_btf ;                  /* 1: BTF , 0: no BTF */
386
 
km->scale = scale ;                 /* -1: none (and no error check)
387
 
                                    0: none, 1: row, 2: max */
388
 
km->user_order = klu_cholmod ;
389
 
km->halt_if_singular = halt_if_singular ;
390
 
 
391
 
printf ("\nKLUDEMO: ordering: %d  BTF: %d scale: %d halt_if_singular: %d\n",
392
 
        ordering, do_btf, scale, halt_if_singular) ;
393
 
 
394
 
    /* ---------------------------------------------------------------------- */
395
 
    /* symbolic factorization */
396
 
    /* ---------------------------------------------------------------------- */
397
 
 
398
 
for (trial = 0 ; trial < NTRIAL ; trial++)
399
 
{
400
 
printf ("---------------- trial %d\n", trial) ;
401
 
 
402
 
    c = CPUTIME ;
403
 
    t = dsecnd_ ( ) ;
404
 
    Symbolic = klu_analyze (n, Ap, Ai, km) ;
405
 
    stats [0] = dsecnd_ ( ) - t ;
406
 
    stats [1] = CPUTIME - c ;
407
 
 
408
 
    if (Symbolic == NULL)
409
 
    {
410
 
        printf ("klu_analyze failed\n") ;
411
 
        exit (1) ;
412
 
    }
413
 
 
414
 
    printf ("\nKLU Symbolic analysis:\n"
415
 
            "dimension of matrix:          %d\n"
416
 
            "nz in off diagonal blocks:    %d\n"
417
 
            "# of off diagonal blocks:     %d\n"
418
 
            "largest block dimension:      %d\n"
419
 
            "estimated nz in L:            %g\n"
420
 
            "estimated nz in U:            %g\n"
421
 
            "symbolic analysis time:       %10.6f (wall) %10.6f (cpu)\n",
422
 
            n, Symbolic->nzoff, Symbolic->nblocks,
423
 
            Symbolic->maxblock, Symbolic->lnz, Symbolic->unz,
424
 
            stats [0], stats [1]) ;
425
 
 
426
 
    /* ---------------------------------------------------------------------- */
427
 
    /* numeric factorization */
428
 
    /* ---------------------------------------------------------------------- */
429
 
 
430
 
    c = CPUTIME ;
431
 
    t = dsecnd_ ( ) ;
432
 
    Numeric = klu_factor (Ap, Ai, Ax, Symbolic, km) ;
433
 
    stats [0] = dsecnd_ ( ) - t ;
434
 
    stats [1] = CPUTIME - c ;
435
 
 
436
 
    if (Numeric == NULL)
437
 
    {
438
 
        printf ("klu_factor failed\n") ;
439
 
    }
440
 
 
441
 
    lnz = (Numeric != NULL) ? Numeric->lnz : 0 ;
442
 
    unz = (Numeric != NULL) ? Numeric->unz : 0 ;
443
 
 
444
 
    printf ("\nKLU Numeric factorization:\n"
445
 
            "actual nz in L:               %d\n"
446
 
            "actual nz in U:               %d\n"
447
 
            "actual nz in L+U:             %d\n"
448
 
            "# of offdiagonal pivots:      %d\n"
449
 
            "numeric factorization time:   %10.6f (wall) %10.6f (cpu)\n",
450
 
            lnz, unz, lnz + unz - n, km->noffdiag, stats [0], stats [1]) ;
451
 
 
452
 
    c = CPUTIME ;
453
 
    t = dsecnd_ ( ) ;
454
 
    klu_flops (Symbolic, Numeric, km) ;
455
 
    t = dsecnd_ ( ) - t ;
456
 
    c = CPUTIME - c ;
457
 
    printf ("flop count %g time %10.6f %10.6f\n", km->flops, t,c) ;
458
 
 
459
 
    c = CPUTIME ;
460
 
    t = dsecnd_ ( ) ;
461
 
    klu_condest (Ap, Ax, Symbolic, Numeric, &condest, km) ;
462
 
    t = dsecnd_ ( ) - t ;
463
 
    c = CPUTIME - c ;
464
 
 
465
 
    printf ("condition number estimate %g time %10.6f %10.6f\n", condest,t,c) ;
466
 
 
467
 
    c = CPUTIME ;
468
 
    t = dsecnd_ ( ) ;
469
 
    klu_growth (Ap, Ai, Ax, Symbolic, Numeric, &growth, km) ;
470
 
    t = dsecnd_ ( ) - t ;
471
 
    c = CPUTIME - c ;
472
 
 
473
 
    printf ("reciprocal pivot growth:  %g time %10.6f %10.6f\n", growth,t,c) ;
474
 
 
475
 
    c = CPUTIME ;
476
 
    t = dsecnd_ ( ) ;
477
 
    klu_rcond (Symbolic, Numeric, &rcond, km) ;
478
 
    t = dsecnd_ ( ) - t ;
479
 
    c = CPUTIME - c ;
480
 
 
481
 
    printf ("cheap reciprocal cond:    %g time %10.6f %10.6f\n", rcond,t,c) ;
482
 
 
483
 
    if (do_sort)
484
 
    {
485
 
        c = CPUTIME ;
486
 
        t = dsecnd_ ( ) ;
487
 
        klu_sort (Symbolic, Numeric, km) ;
488
 
        t = dsecnd_ ( ) - t ;
489
 
        c = CPUTIME - c ;
490
 
        printf (" -------------------- sort time: %10.6f (wall) %10.6f (cpu)\n",
491
 
            t, c) ;
492
 
    }
493
 
 
494
 
    /* ---------------------------------------------------------------------- */
495
 
    /* solve Ax=b */
496
 
    /* ---------------------------------------------------------------------- */
497
 
 
498
 
    for (i = 0 ; i < n ; i++) x [i] = b [i] ;
499
 
 
500
 
    c = CPUTIME ;
501
 
    t = dsecnd_ ( ) ;
502
 
    klu_solve (Symbolic, Numeric, n, 1, x, km) ;  /* does x = A\x */
503
 
    stats [0] = dsecnd_ ( ) - t ;
504
 
    stats [1] = CPUTIME - c ;
505
 
 
506
 
    printf ("\nKLU solve:\n"
507
 
            "solve time:                   %10.6f (wall) %10.6f (cpu)\n",
508
 
            stats [0], stats [1]) ;
509
 
 
510
 
    printf ("\nrelative maxnorm of residual, ||Ax-b||/||b||:  %8.2e\n",
511
 
        resid (n, Ap, Ai, Ax, x, r, b, FALSE)) ;
512
 
    if (!rhs)
513
 
    {
514
 
        printf ("relative maxnorm of error, ||x-xtrue||/||xtrue||: %8.2e\n\n",
515
 
            err (n, x)) ;
516
 
    }
517
 
 
518
 
    if (trial == 0)
519
 
    {
520
 
        printf ("Writing solution to file: tmp/x\n") ;
521
 
        f = fopen ("tmp/x", "w") ;
522
 
        if (!f)
523
 
        {
524
 
            printf ("Unable to open file\n") ;
525
 
            exit (1) ;
526
 
        }
527
 
        for (i = 0 ; i < n ; i++)
528
 
        {
529
 
            fprintf (f, "%24.16e\n", x [i]) ;
530
 
        }
531
 
        fclose (f) ;
532
 
    }
533
 
 
534
 
    /* ---------------------------------------------------------------------- */
535
 
    /* solve A'x=b */
536
 
    /* ---------------------------------------------------------------------- */
537
 
 
538
 
    for (i = 0 ; i < n ; i++) x [i] = b [i] ;
539
 
 
540
 
    c = CPUTIME ;
541
 
    t = dsecnd_ ( ) ;
542
 
    klu_tsolve (Symbolic, Numeric, n, 1, x, km) ;
543
 
    stats [0] = dsecnd_ ( ) - t ;
544
 
    stats [1] = CPUTIME - c ;
545
 
 
546
 
    printf ("\nKLU tsolve:\n"
547
 
            "solve time:                   %10.6f (wall) %10.6f (cpu)\n",
548
 
            stats [0], stats [1]) ;
549
 
 
550
 
    printf ("\nrelative maxnorm of residual, ||A'x-b||/||b||: %8.2e\n",
551
 
        resid (n, Ap, Ai, Ax, x, r, b, TRUE)) ;
552
 
 
553
 
    /* ---------------------------------------------------------------------- */
554
 
    /* solve AX=B */
555
 
    /* ---------------------------------------------------------------------- */
556
 
 
557
 
#if 0
558
 
 
559
 
    for (i = 0 ; i < ldim*NRHS ; i++) x [i] = b [i] ;
560
 
 
561
 
    c = CPUTIME ;
562
 
    t = dsecnd_ ( ) ;
563
 
    klu_solve (Symbolic, Numeric, ldim, NRHS, x, km) ;  /* does x = A\x */
564
 
    stats [0] = dsecnd_ ( ) - t ;
565
 
    stats [1] = CPUTIME - c ;
566
 
 
567
 
    printf ("\nKLU solve (nrhs = %d:\n"
568
 
            "solve time:                   %10.6f (wall) %10.6f (cpu)\n",
569
 
            NRHS, stats [0], stats [1]) ;
570
 
 
571
 
    for (s = 0 ; s < NRHS ; s++)
572
 
    {
573
 
        printf ("relative maxnorm of residual, ||Ax-b||/||b||:  %8.2e\n",
574
 
            resid (n, Ap, Ai, Ax, x + s*ldim, r, b + s*ldim, FALSE)) ;
575
 
    }
576
 
 
577
 
    for (i = 0 ; i < ldim*NRHS ; i++) x [i] = b [i] ;
578
 
 
579
 
    for (s = 0 ; s < NRHS ; s++)
580
 
    {
581
 
        klu_solve (Symbolic, Numeric, n, 1, x+s*ldim, km) ;  /* does x = A\x */
582
 
    }
583
 
    for (s = 0 ; s < NRHS ; s++)
584
 
    {
585
 
        printf ("relative maxnorm of residual, ||Ax-b||/||b||:  %8.2e\n",
586
 
            resid (n, Ap, Ai, Ax, x + s*ldim, r, b + s*ldim, FALSE)) ;
587
 
    }
588
 
 
589
 
    for (nrhs = 1 ; nrhs <= NRHS ; nrhs++)
590
 
    {
591
 
        for (i = 0 ; i < ldim*NRHS ; i++) x [i] = b [i] ;
592
 
 
593
 
        klu_solve (Symbolic, Numeric, ldim, nrhs, x, km) ;  /* does x = A\x */
594
 
 
595
 
        printf ("\nKLU  solve (nrhs = %2d : ", nrhs) ;
596
 
 
597
 
        for (s = 0 ; s < nrhs ; s++)
598
 
        {
599
 
            printf ("relative maxnorm of residual, ||Ax-b||/||b||:  %8.2e\n",
600
 
                resid (n, Ap, Ai, Ax, x + s*ldim, r, b + s*ldim, FALSE)) ;
601
 
        }
602
 
    }
603
 
 
604
 
    /* ---------------------------------------------------------------------- */
605
 
    /* solve A'X=B */
606
 
    /* ---------------------------------------------------------------------- */
607
 
 
608
 
    for (i = 0 ; i < ldim*NRHS ; i++) x [i] = b [i] ;
609
 
 
610
 
    c = CPUTIME ;
611
 
    t = dsecnd_ ( ) ;
612
 
    klu_tsolve (Symbolic, Numeric, ldim, NRHS, x, km) ;  /* does x = A'\x */
613
 
    stats [0] = dsecnd_ ( ) - t ;
614
 
    stats [1] = CPUTIME - c ;
615
 
 
616
 
    printf ("\nKLU tsolve (nrhs = %d:\n"
617
 
            "solve time:                   %10.6f (wall) %10.6f (cpu)\n",
618
 
            NRHS, stats [0], stats [1]) ;
619
 
 
620
 
    for (s = 0 ; s < NRHS ; s++)
621
 
    {
622
 
        printf ("relative maxnorm of residual, ||A'x-b||/||b||: %8.2e\n",
623
 
            resid (n, Ap, Ai, Ax, x + s*ldim, r, b + s*ldim, TRUE)) ;
624
 
    }
625
 
 
626
 
    for (i = 0 ; i < ldim*NRHS ; i++) x [i] = b [i] ;
627
 
 
628
 
    for (s = 0 ; s < NRHS ; s++)
629
 
    {
630
 
        klu_tsolve (Symbolic, Numeric, n, 1, x+s*ldim, 0, km) ;
631
 
    }
632
 
    for (s = 0 ; s < NRHS ; s++)
633
 
    {
634
 
        printf ("relative maxnorm of residual, ||A'x-b||/||b||: %8.2e\n",
635
 
            resid (n, Ap, Ai, Ax, x + s*ldim, r, b + s*ldim, TRUE)) ;
636
 
    }
637
 
 
638
 
    for (nrhs = 1 ; nrhs <= NRHS ; nrhs++)
639
 
    {
640
 
        for (i = 0 ; i < ldim*NRHS ; i++) x [i] = b [i] ;
641
 
 
642
 
        klu_tsolve (Symbolic, Numeric, ldim, nrhs, x, 0, km) ;
643
 
 
644
 
        printf ("\nKLU tsolve (nrhs = %2d : ", nrhs) ;
645
 
 
646
 
        for (s = 0 ; s < nrhs ; s++)
647
 
        {
648
 
            printf ("relative maxnorm of residual, ||A'x-b||/||b||: %8.2e\n",
649
 
                resid (n, Ap, Ai, Ax, x + s*ldim, r, b + s*ldim, TRUE)) ;
650
 
        }
651
 
    }
652
 
 
653
 
#endif
654
 
 
655
 
    /* ---------------------------------------------------------------------- */
656
 
    /* numeric refactorization, no pivot (just to test the code) */
657
 
    /* ---------------------------------------------------------------------- */
658
 
 
659
 
for (rescale = -1 ; rescale <= 2 ; rescale++)
660
 
{
661
 
    km->scale = rescale ;
662
 
 
663
 
    c = CPUTIME ;
664
 
    t = dsecnd_ ( ) ;
665
 
    if (!klu_refactor (Ap, Ai, Ax, Symbolic, Numeric, km))
666
 
    {
667
 
        printf ("klu_refactor failed\n") ;
668
 
    }
669
 
    stats [0] = dsecnd_ ( ) - t ;
670
 
    stats [1] = CPUTIME - c ;
671
 
 
672
 
    printf ("\nKLU Numeric refactorization: scale %d\n"
673
 
            "numeric refactorization time: %10.6f (wall) %10.6f (cpu)\n",
674
 
            km->scale, stats [0], stats [1]) ;
675
 
 
676
 
    klu_condest (Ap, Ax, Symbolic, Numeric, &condest, km) ;
677
 
    printf ("condition number estimate %g\n", condest) ;
678
 
    klu_growth (Ap, Ai, Ax, Symbolic, Numeric, &growth, km) ;
679
 
    printf ("reciprocal pivot growth:  %g\n", growth) ;
680
 
    klu_rcond (Symbolic, Numeric, &rcond, km) ;
681
 
    printf ("cheap reciprocal cond:    %g\n", rcond) ;
682
 
 
683
 
    /* ---------------------------------------------------------------------- */
684
 
    /* solve Ax=b again */
685
 
    /* ---------------------------------------------------------------------- */
686
 
 
687
 
    for (i = 0 ; i < n ; i++) x [i] = b [i] ;
688
 
 
689
 
    c = CPUTIME ;
690
 
    t = dsecnd_ ( ) ;
691
 
    klu_solve (Symbolic, Numeric, n, 1, x, km) ;
692
 
    stats [0] = dsecnd_ ( ) - t ;
693
 
    stats [1] = CPUTIME - c ;
694
 
 
695
 
    printf ("\nKLU solve:\n"
696
 
            "solve time:                   %10.6f (wall) %10.6f (cpu)\n",
697
 
            stats [0], stats [1]) ;
698
 
 
699
 
    printf ("\nrelative maxnorm of residual, ||Ax-b||/||b||:  %8.2e\n",
700
 
        resid (n, Ap, Ai, Ax, x, r, b, FALSE)) ;
701
 
    if (!rhs)
702
 
    {
703
 
        printf ("relative maxnorm of error, ||x-xtrue||/||xtrue||: %8.2e\n\n",
704
 
            err (n, x)) ;
705
 
    }
706
 
 
707
 
}
708
 
    /* restore scale parameter for next trial */
709
 
    km->scale = scale ;
710
 
 
711
 
    /* ---------------------------------------------------------------------- */
712
 
    /* free the Symbolic and Numeric factorization, and all workspace */
713
 
    /* ---------------------------------------------------------------------- */
714
 
 
715
 
    klu_free_symbolic (&Symbolic, km) ;
716
 
    klu_free_numeric (&Numeric, km) ;
717
 
}
718
 
 
719
 
/*
720
 
    printf ("analysis: %12.4f  %12.4f\n", atime[0] / NTRIAL, atime [1] / NTRIAL) ;
721
 
    printf ("factor:   %12.4f  %12.4f\n", ftime[0] / NTRIAL, ftime [1] / NTRIAL) ;
722
 
    printf ("refactor: %12.4f  %12.4f\n", rtime[0] / NTRIAL, rtime [1] / NTRIAL) ;
723
 
    printf ("solve:    %12.4f  %12.4f\n", stime[0] / (4*NTRIAL), stime [1] / (4*NTRIAL)) ;
724
 
    printf ("solve2:   %12.4f  %12.4f (per rhs)\n", stime2[0] / (NTRIAL*NRHS), stime2 [1] / (NTRIAL*NRHS)) ;
725
 
    */
726
 
 
727
 
 
728
 
fflush (stdout) ;
729
 
 
730
 
/* } */ }
731
 
 
732
 
/*
733
 
} } }
734
 
*/
735
 
 
736
 
    free (Ap) ;
737
 
    free (Ai) ;
738
 
    free (Ax) ;
739
 
    free (b) ;
740
 
    free (r) ;
741
 
    free (x) ;
742
 
 
743
 
    printf ("kludemo done\n") ;
 
307
    cholmod_sparse *A ;
 
308
    cholmod_common ch ;
 
309
    cholmod_start (&ch) ;
 
310
    A = cholmod_read_sparse (stdin, &ch) ;
 
311
    if (A)
 
312
    {
 
313
        if (A->nrow != A->ncol || A->stype != 0
 
314
            || (!(A->xtype == CHOLMOD_REAL || A->xtype == CHOLMOD_COMPLEX)))
 
315
        {
 
316
            printf ("invalid matrix\n") ;
 
317
        }
 
318
        else
 
319
        {
 
320
            klu_demo (A->nrow, A->p, A->i, A->x, A->xtype == CHOLMOD_REAL) ;
 
321
        }
 
322
        cholmod_free_sparse (&A, &ch) ;
 
323
    }
 
324
    cholmod_finish (&ch) ;
744
325
    return (0) ;
745
326
}