← Back to branch summary

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

« back to all changes in this revision

Viewing changes to KLU/MATLAB/Test/test4.m

  • Committer: Bazaar Package Importer
  • Author(s): Rafael Laboissiere, Rafael Laboissiere, Ondrej Certik
  • Date: 2008-02-21 14:46:50 UTC
  • mfrom: (1.1.2 upstream) (5.1.1 hardy)
  • Revision ID: james.westby@ubuntu.com-20080221144650-tgeppgj0t7s759i8
Tags: 3.1.0-3
[ Rafael Laboissiere ]
* Upload to unstable

[ Ondrej Certik ]
* XS-DM-Upload-Allowed: yes field added
* Ondrej Certik added to uploaders

Show diffs side-by-side

added added

removed removed

Lines of Context:
KLU/MATLAB/Test/test4.m
1
 
%test4: KLU test script
 
1
function test4 (nmat)
 
2
%test4: KLU test
2
3
% Example:
3
4
%   test4
4
5
% See also klu
6
7
% Copyright 2004-2007 Timothy A. Davis, Univ. of Florida
7
8
% http://www.cise.ufl.edu/research/sparse
8
9
 
9
 
% clear all
10
 
% clear functions
11
10
% rand ('state', 0) ;
12
11
% warning ('off', 'MATLAB:singularMatrix') ;
13
12
% warning ('off', 'MATLAB:nearlySingularMatrix') ;
20
19
f = f (i) ;
21
20
% f = f (1:100) ;
22
21
 
23
 
nmat = length (f) ;
 
22
if (nargin < 1)
 
23
    nmat = 500 ;
 
24
end
 
25
nmat = min (nmat, length (f)) ;
 
26
f = f (1:nmat) ;
24
27
 
25
28
if (1)
26
29
    Tlu = -ones (nmat,1) ;
38
41
 
39
42
% 589: Sieber
40
43
 
41
 
for kk = 1:nmat
42
 
%% for kk = 527:nmat
43
 
 
44
 
    Prob = UFget (f (kk)) ;
45
 
    disp (Prob) ;
46
 
    if (isfield (Prob, 'kind'))
47
 
        if (~isempty (strfind (Prob.kind, 'subsequent')))
48
 
            fprintf ('skip ...\n') ;
49
 
            continue
50
 
        end
51
 
        if (~isempty (strfind (Prob.kind, 'random')))
52
 
            fprintf ('skip ...\n') ;
53
 
            continue
54
 
        end
55
 
    end
56
 
    k = k + 1 ;
57
 
    A = Prob.A ;
58
 
    n = size (A,1) ;
59
 
    err1 = 0 ;
60
 
    err2 = 0 ;
61
 
    err4 = 0 ;
62
 
    terr1 = 0 ;
63
 
    terr2 = 0 ;
64
 
    terr4 = 0 ;
65
 
 
66
 
 
67
 
    for do_imag = 0:1
68
 
        if (do_imag)
69
 
            A = sprand (A) + 1i * sprand (A) ;
70
 
        end
71
 
 
72
 
        % compare with UMFPACK
73
 
        tic
74
 
        [L,U,p,q,R1] = lu (A, 'vector') ;
75
 
        t1 = toc ;
76
 
        if (Tlu (k) == -1)
77
 
            Tlu (k) = t1 ;
78
 
            LUnz (k) = nnz (L) + nnz (U) ;
79
 
        end
80
 
 
81
 
        % note that the scaling R1 and R2 are different with KLU and UMFPACK.
82
 
        % UMFPACK:  L*U-P*(R1\A)*Q
83
 
        % KLU:      L*U-R2\(P*A*Q)
84
 
        %
85
 
        % R1 and R2 are related, via P, where R2 = P*R*P', or equivalently
86
 
        % R2 = R1 (p,p).
87
 
 
88
 
        if (any (diag (U) == 0))
89
 
            fprintf ('skip...\n') ;
90
 
            break ;
91
 
        end
92
 
 
93
 
        F.L = L ;
94
 
        F.U = U ;
95
 
        if (is64)
96
 
            F.p = int64(p) ;
97
 
            F.q = int64(q) ;
98
 
        else
99
 
            F.p = int32(p) ;
100
 
            F.q = int32(q) ;
101
 
        end
102
 
 
103
 
        F.R = R1(p,p) ;
104
 
        b = rand (n,1) ;
105
 
        x = klu (F, '\', b) ;
106
 
        y = klu (b', '/', F) ;
107
 
% y = (b'/A);
108
 
        fprintf ('solve with klu %g\n', ...
109
 
            norm (A*x-b,1)/norm(A,1)) ;
110
 
        fprintf ('solve with klu %g transpose\n', ...
111
 
            norm (y*A-b',1)/norm(A,1)) ;
112
 
 
113
 
 
114
 
        for nrhs = 1:10
115
 
            for do_b_imag = 0:1
116
 
                b = rand (n, nrhs) ;
117
 
                if (do_b_imag)
118
 
                    b = b + 1i * rand (n, nrhs) ;
119
 
                end
120
 
 
121
 
                % KLU backslash
122
 
                tic ;
123
 
                x = klu (A,'\',b) ;
124
 
                t2 = toc ;
125
 
 
126
 
                % KLU slash
127
 
                xt = klu (b','/',A) ;
128
 
 
129
 
                % KLU backslash with precomputed LU
130
 
                tic
131
 
                LU = klu (A) ;
132
 
                z = klu (LU,'\',b) ;
133
 
                t4 = toc ;
134
 
 
135
 
                % KLU slash with precomputed LU
136
 
                zt = klu (b','/',LU) ;
137
 
%               fprintf ('zt err: %g Aimag %d nrhs %d bimag %d\n', ...
138
 
%                   norm (zt*A-b',1), do_imag, nrhs, do_b_imag) ;
139
 
 
140
 
                % UMFPACK
141
 
                tic
142
 
                rb = R1 \ b ;
143
 
                y = U \ (L \ rb (p,:)) ;
144
 
                y (q,:) = y ;
145
 
                t3 = toc ;
146
 
 
147
 
%               yt = b' / A ;
148
 
%               yt = (A' \ b)' ;
149
 
 
150
 
                yt = (L' \ (U' \ b (q,:))) ;
151
 
                yt (p,:) = yt ;
152
 
                yt = R1 \ yt ;
153
 
                yt = yt' ;
154
 
 
155
 
%               norm (yt*A-b',1)
156
 
%               pause
157
 
 
158
 
                if (Tklu (k) == -1)
159
 
                    Tlu (k) = Tlu (k) + t3 ;
160
 
                    Tklu (k) = t2 ;
161
 
                    Tklu2 (k) = t4 ;
162
 
                end
163
 
 
164
 
                err1 = max (err1, norm (A*x-b,1) / norm (A,1)) ;
165
 
                err2 = max (err2, norm (A*y-b,1) / norm (A,1)) ;
166
 
                err4 = max (err4, norm (A*z-b,1) / norm (A,1)) ;
167
 
 
168
 
                terr1 = max (terr1, norm (xt*A-b',1) / norm (A,1)) ;
169
 
                terr2 = max (terr2, norm (yt*A-b',1) / norm (A,1)) ;
170
 
                terr4 = max (terr4, norm (zt*A-b',1) / norm (A,1)) ;
171
 
            end
172
 
        end
173
 
    end
174
 
 
175
 
    fprintf ('err %g %g %g\n', err1, err2, err4) ;
176
 
    if (err1 > 1e4*err2 | err4 > 1e4*err2)                                  %#ok
177
 
        fprintf ('warning: KLU inaccurate!\n')
178
 
    end
179
 
 
180
 
    fprintf ('terr %g %g %g\n', terr1, terr2, terr4) ;
181
 
    if (terr1 > 1e4*terr2 | terr4 > 1e4*terr2)                              %#ok
182
 
        fprintf ('warning: KLU T inaccurate!\n')
183
 
    end
184
 
 
185
 
    lunzmax = max (LUnz (1:k)) ;
186
 
    loglog ( ...
187
 
        LUnz (1:k), Tklu (1:k) ./ Tlu (1:k), 'o', ...
188
 
        LUnz (1:k), Tklu2 (1:k) ./ Tlu (1:k), 'x', ...
189
 
        [10 lunzmax], [1 1], 'r-') ;
190
 
    drawnow
191
 
 
 
44
h = waitbar (0, 'KLU test 4 of 5') ;
 
45
 
 
46
figure (1)
 
47
clf
 
48
 
 
49
try
 
50
 
 
51
    for kk = 1:nmat
 
52
 
 
53
        Prob = UFget (f (kk), index) ;
 
54
 
 
55
        waitbar (kk/nmat, h) ;
 
56
 
 
57
        disp (Prob) ;
 
58
        if (isfield (Prob, 'kind'))
 
59
            if (~isempty (strfind (Prob.kind, 'subsequent')))
 
60
                fprintf ('skip ...\n') ;
 
61
                continue
 
62
            end
 
63
            if (~isempty (strfind (Prob.kind, 'random')))
 
64
                fprintf ('skip ...\n') ;
 
65
                continue
 
66
            end
 
67
        end
 
68
        k = k + 1 ;
 
69
        A = Prob.A ;
 
70
        n = size (A,1) ;
 
71
        err1 = 0 ;
 
72
        err2 = 0 ;
 
73
        err4 = 0 ;
 
74
        terr1 = 0 ;
 
75
        terr2 = 0 ;
 
76
        terr4 = 0 ;
 
77
 
 
78
 
 
79
        for do_imag = 0:1
 
80
            if (do_imag)
 
81
                A = sprand (A) + 1i * sprand (A) ;
 
82
            end
 
83
 
 
84
            % compare with UMFPACK
 
85
            try
 
86
                tic
 
87
                [L,U,p,q,R1] = lu (A, 'vector') ;
 
88
                t1 = max (1e-6, toc) ;
 
89
            catch
 
90
                % older version of MATLAB, which doesn't have 'vector' option
 
91
                tic
 
92
                [L,U,P,Q] = lu (A) ;
 
93
                t1 = max (1e-6, toc) ;
 
94
                [p ignore1 ignore2] = find (P') ;
 
95
                [q ignore1 ignore2] = find (Q) ;
 
96
                clear ignore1 ignore2 P Q
 
97
                R1 = speye (n) ;
 
98
            end
 
99
 
 
100
            if (Tlu (k) == -1)
 
101
                Tlu (k) = t1 ;
 
102
                LUnz (k) = nnz (L) + nnz (U) ;
 
103
            end
 
104
 
 
105
            % note that the scaling R1 and R2 are different with KLU and UMFPACK
 
106
            % UMFPACK:  L*U-P*(R1\A)*Q
 
107
            % KLU:      L*U-R2\(P*A*Q)
 
108
            %
 
109
            % R1 and R2 are related, via P, where R2 = P*R*P', or equivalently
 
110
            % R2 = R1 (p,p).
 
111
 
 
112
            rcond = min (abs (diag (U))) / max (abs (diag (U))) ;
 
113
            if (rcond < 1e-15)
 
114
                fprintf ('skip...\n') ;
 
115
                break ;
 
116
            end
 
117
 
 
118
            F.L = L ;
 
119
            F.U = U ;
 
120
            if (is64)
 
121
                F.p = int64(p) ;
 
122
                F.q = int64(q) ;
 
123
            else
 
124
                F.p = int32(p) ;
 
125
                F.q = int32(q) ;
 
126
            end
 
127
 
 
128
            F.R = R1(p,p) ;
 
129
            b = rand (n,1) ;
 
130
            x = klu (F, '\', b) ;
 
131
            y = klu (b', '/', F) ;
 
132
 
 
133
            fprintf ('solve with klu %g\n', ...
 
134
                norm (A*x-b,1)/norm(A,1)) ;
 
135
            fprintf ('solve with klu %g transpose\n', ...
 
136
                norm (y*A-b',1)/norm(A,1)) ;
 
137
 
 
138
 
 
139
            for nrhs = 1:10
 
140
                for do_b_imag = 0:1
 
141
                    b = rand (n, nrhs) ;
 
142
                    if (do_b_imag)
 
143
                        b = b + 1i * rand (n, nrhs) ;
 
144
                    end
 
145
 
 
146
                    % KLU backslash
 
147
                    tic ;
 
148
                    x = klu (A,'\',b) ;
 
149
                    t2 = max (1e-6, toc) ;
 
150
 
 
151
                    % KLU slash
 
152
                    xt = klu (b','/',A) ;
 
153
 
 
154
                    % KLU backslash with precomputed LU
 
155
                    tic
 
156
                    LU = klu (A) ;
 
157
                    z = klu (LU,'\',b) ;
 
158
                    t4 = max (1e-6, toc) ;
 
159
 
 
160
                    % KLU slash with precomputed LU
 
161
                    zt = klu (b','/',LU) ;
 
162
 
 
163
                    % UMFPACK
 
164
                    tic
 
165
                    rb = R1 \ b ;
 
166
                    y = U \ (L \ rb (p,:)) ;
 
167
                    y (q,:) = y ;
 
168
                    t3 = max (1e-6, toc) ;
 
169
 
 
170
                    yt = (L' \ (U' \ b (q,:))) ;
 
171
                    yt (p,:) = yt ;
 
172
                    yt = R1 \ yt ;
 
173
                    yt = yt' ;
 
174
 
 
175
                    if (Tklu (k) == -1)
 
176
                        Tlu (k) = Tlu (k) + t3 ;
 
177
                        Tklu (k) = t2 ;
 
178
                        Tklu2 (k) = t4 ;
 
179
                    end
 
180
 
 
181
                    err1 = max (err1, norm (A*x-b,1) / norm (A,1)) ;
 
182
                    err2 = max (err2, norm (A*y-b,1) / norm (A,1)) ;
 
183
                    err4 = max (err4, norm (A*z-b,1) / norm (A,1)) ;
 
184
 
 
185
                    terr1 = max (terr1, norm (xt*A-b',1) / norm (A,1)) ;
 
186
                    terr2 = max (terr2, norm (yt*A-b',1) / norm (A,1)) ;
 
187
                    terr4 = max (terr4, norm (zt*A-b',1) / norm (A,1)) ;
 
188
                end
 
189
            end
 
190
        end
 
191
 
 
192
        fprintf ('err %g %g %g\n', err1, err2, err4) ;
 
193
        if (err1 > 1e4*err2 | err4 > 1e4*err2)                              %#ok
 
194
            fprintf ('warning: KLU inaccurate!\n')
 
195
        end
 
196
 
 
197
        fprintf ('terr %g %g %g\n', terr1, terr2, terr4) ;
 
198
        if (terr1 > 1e4*terr2 | terr4 > 1e4*terr2)                          %#ok
 
199
            fprintf ('warning: KLU T inaccurate!\n')
 
200
        end
 
201
 
 
202
        lunzmax = max (LUnz (1:k)) ;
 
203
        loglog ( ...
 
204
            LUnz (1:k), Tklu (1:k) ./ Tlu (1:k), 'o', ...
 
205
            LUnz (1:k), Tklu2 (1:k) ./ Tlu (1:k), 'x', ...
 
206
            [10 lunzmax], [1 1], 'r-') ;
 
207
        drawnow
 
208
 
 
209
    end
 
210
 
 
211
catch
 
212
    % out-of-memory is OK, other errors are not
 
213
    disp (lasterr) ;
 
214
    if (isempty (strfind (lasterr, 'Out of memory')))
 
215
        error (lasterr) ;                                                   %#ok
 
216
    else
 
217
        fprintf ('test terminated early, but otherwise OK\n') ;
 
218
    end
192
219
end
 
220
 
 
221
close (h) ;