1function test3 (nmat) 2%TEST3 test for BTF 3% Requires UFget 4% Example: 5% test3 6% See also btf, maxtrans, strongcomp, dmperm, UFget, 7% test1, test2, test3, test4, test5. 8 9% Copyright 2007, Timothy A. Davis, http://www.suitesparse.com 10 11doplot = 1 ; 12dopause = 0 ; 13dostrong = 1 ; 14 15index = UFget ; 16f = find (index.nrows == index.ncols) ; 17[ignore i] = sort (index.nnz (f)) ; 18f = f (i) ; 19clear i 20 21% short test set: seg faults, lots of blocks, lots of work, and so on: 22nasty = [ 23 % --- various test matrices (no seg fault, quick run time) 24 -(1:8)' % generated matrices 25 904 % vanHeukelum/cage3 (5-by-5) 26 819 % Simon/raefsky6 (permuted triangular matrix) 27 % 28 % --- older seg faults: 29 264 % HB/west0156, causes older strongcomp_recursive to fail 30 824 % TOKAMAK/utm300 (300-by-300), causes older code to fail 31 868 % Pothen/bodyy4 32 % 33 % --- seg faults in old MATLAB dmperm 34 290 % Averous/epb3 35 983 % Sanghavi/ecl32 36 885 % Pothen/tandem_dual 37 879 % Pothen/onera_dual 38 955 % Schenk_IBMSDS/2D_54019_highK 39 957 % Schenk_IBMSDS/3D_51448_3D 40 958 % Schenk_IBMSDS/ibm_matrix_2 41 912 % vanHeukelum/cage11 42 924 % Andrews/Andrews 43 960 % Schenk_IBMSDS/matrix-new_3 44 862 % Kim/kim1 45 544 % Hamm/scircuit 46 897 % Norris/torso2 47 801 % Ronis/xenon1 48 53 % HB/bcsstk31 49 958 % Schenk_IBMSDS/matrix_9 50 844 % Cunningham/qa8fk 51 845 % Cunningham/qa8fk 52 821 % Simon/venkat25 53 822 % Simon/venkat50 54 820 % Simon/venkat01 55 812 % Simon/bbmat 56 804 % Rothberg/cfd1 57 54 % HB/bcsstk32 58 913 % vanHeukelum/cage12 59 846 % Boeing/bcsstk39 60 972 % Schenk_IBMSDS/para-10 61 974 % Schenk_IBMSDS/para-5 62 975 % Schenk_IBMSDS/para-6 63 976 % Schenk_IBMSDS/para-7 64 977 % Schenk_IBMSDS/para-8 65 978 % Schenk_IBMSDS/para-9 66 961 % Schenk_ISEI/barrier2-10 67 962 % Schenk_ISEI/barrier2-11 68 963 % Schenk_ISEI/barrier2-12 69 964 % Schenk_ISEI/barrier2-1 70 965 % Schenk_ISEI/barrier2-2 71 966 % Schenk_ISEI/barrier2-3 72 967 % Schenk_ISEI/barrier2-4 73 968 % Schenk_ISEI/barrier2-9 74 851 % Chen/pkustk05 75 979 % Kamvar/Stanford 76 374 % Bova/rma10 77 % 78 % --- lots of time: 79 395 % DRIVCAV/cavity16 80 396 % DRIVCAV/cavity17 81 397 % DRIVCAV/cavity18 82 398 % DRIVCAV/cavity19 83 399 % DRIVCAV/cavity20 84 400 % DRIVCAV/cavity21 85 401 % DRIVCAV/cavity22 86 402 % DRIVCAV/cavity23 87 403 % DRIVCAV/cavity24 88 404 % DRIVCAV/cavity25 89 405 % DRIVCAV/cavity26 90 1109 % Sandia/mult_dcop_01 91 1110 % Sandia/mult_dcop_02 92 1111 % Sandia/mult_dcop_03 93 376 % Brethour/coater2 94 284 % ATandT/onetone2 95 588 % Hollinger/mark3jac100 96 589 % Hollinger/mark3jac100sc 97 452 % Grund/bayer01 98 920 % Hohn/sinc12 99 590 % Hollinger/mark3jac120 100 591 % Hollinger/mark3jac120sc 101 809 % Shyy/shyy161 102 448 % Graham/graham1 103 283 % ATandT/onetone1 104 445 % Garon/garon2 105 541 % Hamm/bcircuit 106 592 % Hollinger/mark3jac140 107 593 % Hollinger/mark3jac140sc 108 435 % FIDAP/ex40 109 912 % Hohn/sinc15 110 894 % Norris/lung2 111 542 % Hamm/hcircuit 112 752 % Mulvey/finan512 113 753 % Mulvey/pfinan512 114 564 % Hollinger/g7jac180 115 565 % Hollinger/g7jac180sc 116 566 % Hollinger/g7jac200 117 567 % Hollinger/g7jac200sc 118 748 % Mallya/lhr34 119 749 % Mallya/lhr34c 120 922 % Hohn/sinc18 121 447 % Goodwin/rim 122 807 % Rothberg/struct3 123 286 % ATandT/twotone 124 982 % Tromble/language 125 953 % Schenk_IBMNA/c-73 126 890 % Norris/heart1 127 750 % Mallya/lhr71 128 751 % Mallya/lhr71c 129 925 % FEMLAB/ns3Da 130 827 % Vavasis/av41092 131 931 % FEMLAB/sme3Db 132 1297 % GHS_index/boyd2 133 1301 % GHS_indef/cont-300 134 % 135 % --- lots of time, and seg faults: 136 285 % ATandT/pre2 137 % --- huge matrix, turn off plotting 138 940 % Shenk/af_shell1, memory leak in plot, after call to btf, once. 139 % ---- 140]' ; 141 142% maxtrans_recursive causes a seg fault on these matrices, because of 143% stack overflow (this is expected) 144skip_list_maxtrans_recursive = 285 ; 145 146% p = dmperm (A) in MATLAB 7.4 causes a seg fault on these matrices: 147skip_list_dmperm = [285 1301 1231 1251 1232 1241] ; 148 149% [p,q,r] = dmperm (A) in MATLAB 7.4 causes a seg fault on these matrices: 150skip_list_dmperm_btf = ... 151[ 285 879 885 290 955 957 958 924 960 897 959 844 845 ... 152 821 822 820 804 913 846 972 974:978 961:968 979 940 ... 153 1422 1513 1412 1510 1301 1231 1251 1434 1213 1232 1241 1357 1579 1431 1281] ; 154% length(skip_list_dmperm_btf) 155 156% time intensive 157skip_costly = [1514 1297 1876 1301] ; 158 159% strongcomp (recursive) causes a seg fault on these matrices because of 160% stack overflow (this is expected). 161skip_list_strongcomp_recursive = ... 162[983 285 879 885 290 955 957 958 912 924 960 862 544 897 801 53 959 844 845 ... 163 821 822 820 812 804 54 913 846 972 974:978 961:968 851 374 940] ; 164skip_list_strongcomp_recursive = ... 165[ skip_list_strongcomp_recursive 592 593 752 753 807 286 982 855 566 567 ] ; 166 167% matrices with the largest # of nonzeros in the set (untested) 168toobig = [ 169928 853 852 356 761 368 973 895 805 849 932 ... 170803 854 936 802 850 537 856 898 857 859 971 937 ... 171914 858 980 896 806 538 863 369 938 860 941 942 ... 172943 944 945 946 947 948 915 939 916 ] ; 173 174f = [ -(1:8) f ] ; 175% f = nasty ; 176 177h = waitbar (0, 'BTF test 3 of 6') ; 178 179if (nargin < 1) 180 nmat = 1000 ; 181end 182nmat = min (nmat, length (f)) ; 183f = f (1:nmat) ; 184 185try 186 187 for matnum = 1:nmat 188 189 waitbar (matnum/nmat, h) ; 190 191 j = f (matnum) ; 192 193 if (any (j == toobig) | any (j == skip_costly)) %#ok 194 fprintf ('\n%4d: %3d %s/%s too big\n', ... 195 matnum, j, index.Group{j}, index.Name{j}) ; 196 continue ; 197 end 198 199 rand ('state', 0) ; 200 201 % clear all unused variables. 202 % nothing here is left that is proportional to the matrix size 203 clear A p1 p2 p3 q3 r3 match1 match2 match4 pa ra sa qa B C pb rb pc rc 204 clear jumble B11 B12 B13 B21 B22 B23 B31 B32 B33 pjumble qjumble ans 205 clear c kbad kgood 206 % whos 207 % pause 208 209 if (j > 0) 210 Problem = UFget (j, index) ; 211 name = Problem.name ; 212 A = Problem.A ; 213 clear Problem 214 else 215 % construct the jth test matrix 216 j = -j ; 217 if (j == 1 | j == 2) %#ok 218 B11 = UFget ('Grund/b1_ss') ; % 7-by-7 diagonal block 219 B11 = B11.A ; 220 B12 = sparse (zeros (7,2)) ; 221 B12 (3,2) = 1 ; 222 B13 = sparse (ones (7,5)) ; 223 B21 = sparse (zeros (2,7)) ; 224 B22 = sparse (ones (2,2)) ; % 2-by-2 diagonal block 225 B23 = sparse (ones (2,5)) ; 226 B31 = sparse (zeros (5,7)) ; 227 B32 = sparse (zeros (5,2)) ; 228 B33 = UFget ('vanHeukelum/cage3') ; % 5-by-5 diagonal block 229 B33 = B33.A ; 230 A = [ B11 B12 B13 ; B21 B22 B23 ; B31 B32 B33 ] ; 231 name = '(j=1 test matrix)' ; 232 end 233 if (j == 2) 234 pjumble = [ 10 7 11 1 13 12 8 2 5 14 9 6 4 3 ] ; 235 qjumble = [ 3 14 2 11 1 8 5 7 10 12 4 13 9 6 ] ; 236 A = A (pjumble, qjumble) ; 237 name = '(j=2 test matrix)' ; 238 elseif (j == 3) 239 A = sparse (1) ; 240 elseif (j == 4) 241 A = sparse (0) ; 242 elseif (j == 5) 243 A = sparse (ones (2)) ; 244 elseif (j == 6) 245 A = sparse (2,2) ; 246 elseif (j == 7) 247 A = speye (2) ; 248 elseif (j == 8) 249 A = sparse (2,2) ; 250 A (2,1) = 1 ; 251 end 252 if (j > 2) 253 full (A) 254 end 255 end 256 257 [m n] = size (A) ; 258 if (m ~= n) 259 continue ; 260 end 261 fprintf ('\n%4d: ', matnum) ; 262 fprintf (' =========================== Matrix: %3d %s\n', j, name) ; 263 fprintf ('n: %d nz: %d\n', n, nnz (A)) ; 264 265 if (nnz (A) > 6e6) 266 doplot = 0 ; 267 end 268 269 %----------------------------------------------------------------------- 270 % now try maxtrans 271 tic 272 match1 = maxtrans (A) ; 273 t = toc ; 274 s1 = sum (match1 > 0) ; 275 fprintf ('n-sprank: %d\n', n-s1) ; 276 fprintf ('maxtrans: %8.2f seconds\n', t) ; 277 singular = s1 < n ; 278 279 if (doplot) 280 clf 281 subplot (2,4,1) 282 spy (A) 283 title (name) ; 284 end 285 286 p1 = match1 ; 287 if (any (p1 <= 0)) 288 % complete the permutation 289 badrow = find (p1 <= 0) ; 290 291 badcol = ones (1,n) ; 292 badcol (p1 (p1 > 0)) = 0 ; 293 badcol = find (badcol) ; 294 295 p1 (badrow) = badcol ; 296 297 % construct the older form of match1 298 match1 (badrow) = -p1 (badrow) ; 299 end 300 if (any (sort (p1) ~= 1:n)) 301 error ('!!') ; 302 end 303 304 B = A (:,p1) ; 305 306 if (doplot) 307 subplot (2,4,2) 308 hold off 309 spy (B) 310 hold on 311 badcol = find (match1 < 0) ; 312 Junk = sparse (badcol, badcol, ones (length (badcol), 1), n, n) ; 313 % if (~isempty (A)) 314 % spy (Junk, 'ro') ; 315 % end 316 title ('maxtrans') ; 317 end 318 319 d = nnz (diag (B)) ; 320 if (d ~= s1) 321 error ('bad sprank') ; 322 end 323 clear B 324 325 %----------------------------------------------------------------------- 326 % try p = dmperm(A) 327 skip_dmperm = any (j == skip_list_dmperm) ; 328 329 if (~skip_dmperm) 330 tic 331 match4 = dmperm (A) ; 332 t = toc ; 333 fprintf ('p=dmperm(A): %8.2f seconds\n', t) ; 334 s4 = sum (match4 > 0) ; 335 singular4 = (s4 < n) ; 336 337 if (doplot) 338 if (~singular4) 339 subplot (2,4,3) 340 spy (A (match4,:)) 341 title ('dmperm') ; 342 end 343 end 344 if (singular ~= singular4) 345 error ('s4?') ; 346 end 347 if (s1 ~= s4) 348 error ('bad sprank') ; 349 end 350 else 351 fprintf ('p=dmperm(A): skip\n') ; 352 end 353 354 %----------------------------------------------------------------------- 355 nblocks = -1 ; 356 skip_dmperm_btf = any (j == skip_list_dmperm_btf) ; 357 if (~skip_dmperm_btf) 358 % get btf form 359 tic 360 [pa,qa,ra,sa] = dmperm (A) ; 361 t = toc ; 362 fprintf ('[p,q,r,s]=dmperm(A): %8.2f seconds\n', t) ; 363 nblocks = length (ra) - 1 ; 364 fprintf ('nblocks: %d\n', nblocks) ; 365 if (~singular4) 366 checkbtf (A, pa, qa, ra) ; 367 if (doplot) 368 subplot (2,4,4) 369 drawbtf (A, pa, qa, ra) 370 title ('dmperm blocks') 371 end 372 end 373 else 374 fprintf ('[p,q,r,s]=dmperm(A): skip\n') ; 375 end 376 377 jumble = randperm (n) ; 378 379 %----------------------------------------------------------------------- 380 % try strongcomp, non-recursive version 381 382 %------------------------------------------------------------------- 383 % try strongcomp on original matrix 384 B = A (:,p1) ; 385 tic ; 386 [pb,rb] = strongcomp (B) ; 387 t = toc ; 388 fprintf ('strongcomp %8.2f seconds\n', t) ; 389 if (~singular & ~skip_dmperm_btf & (length (rb) ~= nblocks+1)) %#ok 390 error ('BTF:invalid (rb)') ; 391 end 392 checkbtf (B, pb, pb, rb) ; 393 if (doplot) 394 subplot (2,4,5) 395 drawbtf (B, pb, pb, rb) ; 396 title ('strongcomp') ; 397 end 398 399 %------------------------------------------------------------------- 400 % try btf on original matrix 401 tic ; 402 [pw,qw,rw] = btf (A) ; 403 t = toc ; 404 fprintf ('btf %8.2f seconds nblocks %d\n', ... 405 t, length (rw)-1) ; 406 407 if (any (pw ~= pb)) 408 error ('pw') ; 409 end 410 if (any (rw ~= rb)) 411 error ('rw') ; 412 end 413 if (any (abs (qw) ~= p1 (pw))) 414 error ('qw') ; 415 end 416 c = diag (A (pw,abs (qw))) ; 417 if (~singular & ~skip_dmperm_btf & (length (rw) ~= nblocks+1)) %#ok 418 error ('BTF:invalid (rw)') ; 419 end 420 checkbtf (A, pw, abs (qw), rw) ; 421 422 kbad = find (qw < 0) ; 423 kgood = find (qw > 0) ; 424 if (any (c (kbad) ~= 0)) 425 error ('kbad') ; 426 end 427 if (any (c (kgood) == 0)) %#ok 428 error ('kgood') ; 429 end 430 431 if (doplot) 432 subplot (2,4,6) 433 drawbtf (A, pw, abs (qw), rw) ; 434 if (n < 500) 435 for k = kbad 436 plot ([k (k+1) (k+1) k k]-.5, ... 437 [k k (k+1) (k+1) k]-.5, 'r') ; 438 end 439 end 440 title ('btf') ; 441 end 442 443 %------------------------------------------------------------------- 444 % try [p,q,r] = strongcomp (A, qin) form 445 tic 446 [pz,qz,rz] = strongcomp (A, match1) ; 447 t = toc ; 448 fprintf ('[p,q,r]=strongcomp(A,qin)%8.2f seconds\n', t) ; 449 if (any (pz ~= pb)) 450 error ('pz') ; 451 end 452 if (any (rz ~= rb)) 453 error ('rz') ; 454 end 455 if (any (abs (qz) ~= p1 (pz))) 456 error ('qz') ; 457 end 458 c = diag (A (pz,abs (qz))) ; 459 if (~singular & ~skip_dmperm_btf & (length (rz) ~= nblocks+1)) %#ok 460 error ('BTF:invalid (rz)') ; 461 end 462 checkbtf (A, pz, abs (qz), rz) ; 463 464 kbad = find (qz < 0) ; 465 kgood = find (qz > 0) ; 466 if (any (c (kbad) ~= 0)) 467 error ('kbad') ; 468 end 469 if (any (c (kgood) == 0)) %#ok 470 error ('kgood') ; 471 end 472 473 if (doplot) 474 subplot (2,4,7) 475 drawbtf (A, pz, abs (qz), rz) ; 476 if (n < 500) 477 for k = kbad 478 plot ([k (k+1) (k+1) k k]-.5, ... 479 [k k (k+1) (k+1) k]-.5, 'r') ; 480 end 481 end 482 title ('strongcomp(A,qin)') ; 483 end 484 485 %------------------------------------------------------------------- 486 % try strongcomp again, on a randomly jumbled matrix 487 C = sparse (B (jumble, jumble)) ; 488 tic ; 489 [pc,rc] = strongcomp (C) ; 490 t = toc ; 491 fprintf ('strongcomp (rand) %8.2f seconds\n', t) ; 492 if (~singular & ~skip_dmperm_btf & (length (rc) ~= nblocks+1)) %#ok 493 error ('BTF:invalid (rc)') ; 494 end 495 checkbtf (C, pc, pc, rc) ; 496 if (doplot) 497 subplot (2,4,8) 498 drawbtf (C, pc, pc, rc) ; 499 title ('strongcomp(rand)') ; 500 end 501 502 if (length (rc) ~= length (rb)) 503 error ('strongcomp random mismatch') ; 504 end 505 506 %----------------------------------------------------------------------- 507 if (doplot) 508 drawnow 509 end 510 511 if (matnum ~= nmat & dopause) %#ok 512 input ('Hit enter: ') ; 513 end 514 515 end 516 517catch 518 % out-of-memory is OK, other errors are not 519 disp (lasterr) ; 520 if (isempty (strfind (lasterr, 'Out of memory'))) 521 error (lasterr) ; %#ok 522 else 523 fprintf ('test terminated early, but otherwise OK\n') ; 524 end 525end 526 527close (h) ; 528