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