1function [C,R,E,B,X, err] = spqr_gpu3 (ordering, A)
2% info = spqr_gpu2 (ordering,A)
3%   ordering: 1 colamd
4%   ordering: 2 metis
5
6% write the matrix to a file
7mwrite ('A.mtx', A) ;
8
9if (exist ('R.mtx','file'))
10    delete ('R.mtx') ;
11end
12if (exist ('C.mtx','file'))
13    delete ('C.mtx') ;
14end
15if (exist ('E.txt','file'))
16    delete ('E.txt') ;
17end
18
19setenv('LD_LIBRARY_PATH', '/usr/local/cuda/lib64:/usr/lib:/lib/x86_64-linux-gnu:/usr/lib/x86_64-linux-gnu:/lib64')
20if (ordering == 1)
21    system ('tcsh demo_colamd3.sh') ;
22else
23    system ('tcsh demo_metis3.sh') ;
24end
25
26atanorm = norm (A'*A,1) ;
27
28[m n] = size (A) ;
29R = mread ('R.mtx') ;
30C = mread ('C.mtx') ;
31B = ones (m,1) ;
32E = load ('E.txt') ;
33S = A (:,E) ;
34err (1) = norm (R'*R - S'*S, 1) / atanorm  ;
35X = R\C ;
36X (E) = X ;
37
38[C2, R2] = qr (S, B, 0) ;
39whos
40err (2) = norm (R2'*R2 - S'*S, 1) / atanorm ;
41X2 = R2\C2 ;
42X2 (E) = X2 ;
43err (3) = norm (X-X2) / norm (X2) ;
44
45x = A\B ;
46err (4) = norm (A'*(A*x-B))  / atanorm ;
47err (5) = norm (A'*(A*X-B))  / atanorm ;
48err (6) = norm (A'*(A*X2-B)) / atanorm ;
49
50
51% normalize the diagonal of R
52e = spdiags (sign (full (diag (R))), 0, n, n) ;
53R = e*R ;
54err (7) = norm (R'*R - S'*S, 1) / atanorm ;
55
56e = spdiags (sign (full (diag (R2))), 0, n, n) ;
57R2 = e*R2 ;
58err (8) = norm (R2'*R2 - S'*S, 1) / atanorm ;
59
60
61% err
62% log10(err)
63
64fprintf ('1: GPU    norm(R''R-S''S)               %12.4e\n', err (1)) ;
65fprintf ('2: matlab norm(R''R-S''S)               %12.4e\n', err (2)) ;
66fprintf ('3: GPU vs matlab norm (X-X2)/norm(X)  %12.4e\n', err (3)) ;
67fprintf ('4: x=A\\b norm(A''(Ax-b))/norm(A''A)     %12.4e\n', err (4)) ;
68fprintf ('5: GPU    norm(A''(Ax-b))/norm(A''A)    %12.4e\n', err (5)) ;
69fprintf ('6: matlab norm(A''(Ax-b))/norm(A''A)    %12.4e\n', err (6)) ;
70fprintf ('7: GPU    norm(R''R-S''S) nonneg        %12.4e\n', err (7)) ;
71fprintf ('8: matlab norm(R''R-S''S) nonneg        %12.4e\n', err (8)) ;
72fprintf ('\n-------------------------------------------------------\n') ;
73
74if (err (2) > 1e-10 || err (5) > max (1e-12, 1e3 * err (4)))
75    warning ('!')
76end
77