t_mplinsolve

T_MPLINSOLVE  Tests of MIPS/MATPOWER linear solvers.

0001 function t_mplinsolve(quiet)
0002 %T_MPLINSOLVE  Tests of MIPS/MATPOWER linear solvers.
0003 
0004 %   MIPS
0005 %   Copyright (c) 2015-2018, Power Systems Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   This file is part of MIPS.
0009 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0010 %   See https://github.com/MATPOWER/mips for more info.
0011 
0012 if nargin < 1
0013     quiet = 0;
0014 end
0015 
0016 t_begin(66, quiet);
0017 
0018 isoctave = exist('OCTAVE_VERSION', 'builtin') == 5;
0019 if isoctave
0020     lu_warning_id = 'Octave:lu:sparse_input';
0021     s = warning('query', lu_warning_id);
0022     warning('off', lu_warning_id);
0023 end
0024 
0025 %% non-sparse A matrix crashes some MATLAB lu() calls on certain
0026 %% combinations of MATLAB and macOS versions
0027 mlv = have_feature('matlab', 'vnum');
0028 skipcrash = ~isempty(mlv) && mlv < 8.003 && mlv > 7.013 && ...
0029     strcmp(computer, 'MACI64');
0030 
0031 ijs = [
0032     1 1 1205.63;
0033     4 1 -1205.63;
0034     25 1 17.3611;
0035     28 1 -17.3611;
0036     43 1 1;
0037     2 2 1024;
0038     8 2 -1024;
0039     26 2 16;
0040     32 2 -16;
0041     3 3 1164.84;
0042     6 3 -1164.84;
0043     27 3 17.0648;
0044     30 3 -17.0648;
0045     1 4 -1205.63;
0046     4 4 2201.34;
0047     5 4 -453.695;
0048     9 4 -542.009;
0049     13 4 0.776236;
0050     14 4 -0.451721;
0051     18 4 -0.324515;
0052     25 4 -17.3611;
0053     28 4 39.4759;
0054     29 4 -10.5107;
0055     33 4 -11.6041;
0056     37 4 -3.30738;
0057     38 4 1.94219;
0058     42 4 1.36519;
0059     4 5 -453.695;
0060     5 5 581.372;
0061     6 5 -127.677;
0062     13 5 -0.451721;
0063     14 5 0.568201;
0064     15 5 -0.11648;
0065     28 5 -10.5107;
0066     29 5 16.0989;
0067     30 5 -5.58824;
0068     37 5 1.94219;
0069     38 5 -3.2242;
0070     39 5 1.28201;
0071     3 6 -1164.84;
0072     5 6 -127.677;
0073     6 6 1676.74;
0074     7 6 -384.227;
0075     14 6 -0.11648;
0076     15 6 0.163059;
0077     16 6 -0.0465791;
0078     27 6 -17.0648;
0079     29 6 -5.58824;
0080     30 6 32.4374;
0081     31 6 -9.78427;
0082     38 6 1.28201;
0083     39 6 -2.4371;
0084     40 6 1.15509;
0085     6 7 -384.227;
0086     7 7 1141.16;
0087     8 7 -756.935;
0088     15 7 -0.0465791;
0089     16 7 0.371828;
0090     17 7 -0.325249;
0091     30 7 -9.78427;
0092     31 7 23.4822;
0093     32 7 -13.698;
0094     39 7 1.15509;
0095     40 7 -2.77221;
0096     41 7 1.61712;
0097     2 8 -1024;
0098     7 8 -756.935;
0099     8 8 1925.77;
0100     9 8 -144.836;
0101     16 8 -0.325249;
0102     17 8 0.844105;
0103     18 8 -0.518855;
0104     26 8 -16;
0105     31 8 -13.698;
0106     32 8 35.6731;
0107     33 8 -5.97513;
0108     40 8 1.61712;
0109     41 8 -2.80473;
0110     42 8 1.1876;
0111     4 9 -542.009;
0112     8 9 -144.836;
0113     9 9 686.845;
0114     13 9 -0.324515;
0115     17 9 -0.518855;
0116     18 9 0.84337;
0117     28 9 -11.6041;
0118     32 9 -5.97513;
0119     33 9 17.5792;
0120     37 9 1.36519;
0121     41 9 1.1876;
0122     42 9 -2.55279;
0123     10 10 1207.63;
0124     13 10 -1205.63;
0125     34 10 17.3611;
0126     37 10 -17.3611;
0127     11 11 1026;
0128     17 11 -1024;
0129     35 11 16;
0130     41 11 -16;
0131     12 12 1166.84;
0132     15 12 -1164.84;
0133     36 12 17.0648;
0134     39 12 -17.0648;
0135     4 13 0.776236;
0136     5 13 -0.451721;
0137     9 13 -0.324515;
0138     10 13 -1205.63;
0139     13 13 2190.83;
0140     14 13 -447.892;
0141     18 13 -535.137;
0142     28 13 3.30738;
0143     29 13 -1.94219;
0144     33 13 -1.36519;
0145     34 13 -17.3611;
0146     37 13 39.1419;
0147     38 13 -10.5107;
0148     42 13 -11.6041;
0149     4 14 -0.451721;
0150     5 14 0.568201;
0151     6 14 -0.11648;
0152     13 14 -447.892;
0153     14 14 573.221;
0154     15 14 -122.862;
0155     28 14 -1.94219;
0156     29 14 3.2242;
0157     30 14 -1.28201;
0158     37 14 -10.5107;
0159     38 14 15.5829;
0160     39 14 -5.58824;
0161     5 15 -0.11648;
0162     6 15 0.163059;
0163     7 15 -0.0465791;
0164     12 15 -1164.84;
0165     14 15 -122.862;
0166     15 15 1669.87;
0167     16 15 -379.651;
0168     29 15 -1.28201;
0169     30 15 2.4371;
0170     31 15 -1.15509;
0171     36 15 -17.0648;
0172     38 15 -5.58824;
0173     39 15 31.8704;
0174     40 15 -9.78427;
0175     6 16 -0.0465791;
0176     7 16 0.371828;
0177     8 16 -0.325249;
0178     15 16 -379.651;
0179     16 16 1131.92;
0180     17 16 -750.073;
0181     30 16 -1.15509;
0182     31 16 2.77221;
0183     32 16 -1.61712;
0184     39 16 -9.78427;
0185     40 16 23.1242;
0186     41 16 -13.698;
0187     7 17 -0.325249;
0188     8 17 0.844105;
0189     9 17 -0.518855;
0190     11 17 -1024;
0191     16 17 -750.073;
0192     17 17 1914.92;
0193     18 17 -138.499;
0194     31 17 -1.61712;
0195     32 17 2.80473;
0196     33 17 -1.1876;
0197     35 17 -16;
0198     40 17 -13.698;
0199     41 17 35.2181;
0200     42 17 -5.97513;
0201     4 18 -0.324515;
0202     8 18 -0.518855;
0203     9 18 0.84337;
0204     13 18 -535.137;
0205     17 18 -138.499;
0206     18 18 676.012;
0207     28 18 -1.36519;
0208     32 18 -1.1876;
0209     33 18 2.55279;
0210     37 18 -11.6041;
0211     41 18 -5.97513;
0212     42 18 17.0972;
0213     19 19 1.88667;
0214     25 19 -1;
0215     20 20 1.54931;
0216     26 20 -1;
0217     21 21 1.78346;
0218     27 21 -1;
0219     22 22 0.666667;
0220     34 22 -1;
0221     23 23 0.666667;
0222     35 23 -1;
0223     24 24 0.666667;
0224     36 24 -1;
0225     1 25 17.3611;
0226     4 25 -17.3611;
0227     19 25 -1;
0228     2 26 16;
0229     8 26 -16;
0230     20 26 -1;
0231     3 27 17.0648;
0232     6 27 -17.0648;
0233     21 27 -1;
0234     1 28 -17.3611;
0235     4 28 39.4759;
0236     5 28 -10.5107;
0237     9 28 -11.6041;
0238     13 28 3.30738;
0239     14 28 -1.94219;
0240     18 28 -1.36519;
0241     4 29 -10.5107;
0242     5 29 16.0989;
0243     6 29 -5.58824;
0244     13 29 -1.94219;
0245     14 29 3.2242;
0246     15 29 -1.28201;
0247     3 30 -17.0648;
0248     5 30 -5.58824;
0249     6 30 32.4374;
0250     7 30 -9.78427;
0251     14 30 -1.28201;
0252     15 30 2.4371;
0253     16 30 -1.15509;
0254     6 31 -9.78427;
0255     7 31 23.4822;
0256     8 31 -13.698;
0257     15 31 -1.15509;
0258     16 31 2.77221;
0259     17 31 -1.61712;
0260     2 32 -16;
0261     7 32 -13.698;
0262     8 32 35.6731;
0263     9 32 -5.97513;
0264     16 32 -1.61712;
0265     17 32 2.80473;
0266     18 32 -1.1876;
0267     4 33 -11.6041;
0268     8 33 -5.97513;
0269     9 33 17.5792;
0270     13 33 -1.36519;
0271     17 33 -1.1876;
0272     18 33 2.55279;
0273     10 34 17.3611;
0274     13 34 -17.3611;
0275     22 34 -1;
0276     11 35 16;
0277     17 35 -16;
0278     23 35 -1;
0279     12 36 17.0648;
0280     15 36 -17.0648;
0281     24 36 -1;
0282     4 37 -3.30738;
0283     5 37 1.94219;
0284     9 37 1.36519;
0285     10 37 -17.3611;
0286     13 37 39.1419;
0287     14 37 -10.5107;
0288     18 37 -11.6041;
0289     4 38 1.94219;
0290     5 38 -3.2242;
0291     6 38 1.28201;
0292     13 38 -10.5107;
0293     14 38 15.5829;
0294     15 38 -5.58824;
0295     5 39 1.28201;
0296     6 39 -2.4371;
0297     7 39 1.15509;
0298     12 39 -17.0648;
0299     14 39 -5.58824;
0300     15 39 31.8704;
0301     16 39 -9.78427;
0302     6 40 1.15509;
0303     7 40 -2.77221;
0304     8 40 1.61712;
0305     15 40 -9.78427;
0306     16 40 23.1242;
0307     17 40 -13.698;
0308     7 41 1.61712;
0309     8 41 -2.80473;
0310     9 41 1.1876;
0311     11 41 -16;
0312     16 41 -13.698;
0313     17 41 35.2181;
0314     18 41 -5.97513;
0315     4 42 1.36519;
0316     8 42 1.1876;
0317     9 42 -2.55279;
0318     13 42 -11.6041;
0319     17 42 -5.97513;
0320     18 42 17.0972;
0321     1 43 1;
0322 ];
0323 A = sparse(ijs(:, 1), ijs(:, 2), ijs(:, 3));
0324 b = [
0325     0;
0326     0;
0327     0;
0328     0;
0329     0;
0330     0;
0331     0;
0332     0;
0333     0;
0334     0;
0335     0;
0336     0;
0337     -0.00896054;
0338     -0.0617829;
0339     -0.0772931;
0340     -0.0230638;
0341     -0.0185934;
0342     -0.02;
0343     -0.336;
0344     -0.2755;
0345     -0.353;
0346     0;
0347     0;
0348     0;
0349     1.3;
0350     1.55;
0351     1.4;
0352     0;
0353     -0.9;
0354     0;
0355     -1;
0356     0;
0357     -1.25;
0358     0;
0359     0;
0360     0;
0361     0.167;
0362     -0.042;
0363     0.2835;
0364     -0.171;
0365     0.2275;
0366     -0.259;
0367     0;
0368 ];
0369 ex = [
0370 0;
0371 0.04612219791119214;
0372 0.0509334351882598;
0373 -0.05953571031927305;
0374 -0.09461814898578046;
0375 -0.008909854578010561;
0376 -0.05785829019394401;
0377 -0.02232729212460287;
0378 -0.1137760871247425;
0379 -0.03062777824802364;
0380 -0.01013282572376477;
0381 -0.005330939680091628;
0382 -0.02914165388753019;
0383 -0.03376073204420303;
0384 4.021341450281111e-05;
0385 -0.01727289094763518;
0386 -0.008382063634320435;
0387 -0.04854008812629265;
0388 -0.2663945795760685;
0389 -0.4548081594272799;
0390 -0.3787862287965494;
0391 -0.02580075363496279;
0392 -0.02801219343110931;
0393 -0.09165765332863518;
0394 -0.1665986614487812;
0395 -0.4291388294822789;
0396 -0.3225500876094941;
0397 3.967864472606993;
0398 5.577372159790927;
0399 3.762341481664236;
0400 5.858342308599034;
0401 3.951628532808602;
0402 6.471657726723339;
0403 -0.01720051102355974;
0404 -0.01867480495813735;
0405 -0.06110513277164123;
0406 -0.1239317474796463;
0407 -0.1848327901217308;
0408 -0.4283638085291988;
0409 0.07640167050820651;
0410 -0.1319901818980452;
0411 0.4366406661538687;
0412 0.0007894844305239289;
0413 ];
0414 
0415 t = ''''' : ';
0416 x = mplinsolve(A, b, '');
0417 t_is(x, ex, 12, [t 'x']);
0418 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0419 x = mplinsolve(full(A), b, '');
0420 t_is(x, ex, 12, [t 'x (full A)']);
0421 t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0422 
0423 t = '\ : ';
0424 x = mplinsolve(A, b, '\');
0425 t_is(x, ex, 12, [t 'x']);
0426 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0427 x = mplinsolve(full(A), b, '\');
0428 t_is(x, ex, 12, [t 'x (full A)']);
0429 t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0430 
0431 t = 'LU : ';
0432 x = mplinsolve(A, b, 'LU');
0433 t_is(x, ex, 12, [t 'x']);
0434 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0435 if skipcrash
0436     t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
0437 else
0438     x = mplinsolve(full(A), b, 'LU');
0439     t_is(x, ex, 12, [t 'x (full A)']);
0440     t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0441 end
0442 
0443 t = 'LU3 : ';
0444 x = mplinsolve(A, b, 'LU3');
0445 t_is(x, ex, 12, [t 'x']);
0446 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0447 if skipcrash
0448     t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
0449 else
0450     x = mplinsolve(full(A), b, 'LU3');
0451     t_is(x, ex, 12, [t 'x (full A)']);
0452     t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0453 end
0454 t = 'LU, nout = 3, vec = 1, thresh = 1 : ';
0455 opt = struct('nout', 3, 'vec', 1, 'thresh', 1);
0456 x = mplinsolve(A, b, 'LU', opt);
0457 t_is(x, ex, 12, [t 'x']);
0458 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0459 if skipcrash
0460     t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
0461 else
0462     x = mplinsolve(full(A), b, 'LU', opt);
0463     t_is(x, ex, 12, [t 'x (full A)']);
0464     t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0465 end
0466 
0467 t = 'LU3a : ';
0468 x = mplinsolve(A, b, 'LU3a');
0469 t_is(x, ex, 12, [t 'x']);
0470 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0471 if skipcrash
0472     t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
0473 else
0474     x = mplinsolve(full(A), b, 'LU3a');
0475     t_is(x, ex, 12, [t 'x (full A)']);
0476     t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0477 end
0478 t = 'LU, nout = 3, vec = 1 : ';
0479 opt = struct('nout', 3, 'vec', 1);
0480 x = mplinsolve(A, b, 'LU', opt);
0481 t_is(x, ex, 12, [t 'x']);
0482 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0483 if skipcrash
0484     t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
0485 else
0486     x = mplinsolve(full(A), b, 'LU', opt);
0487     t_is(x, ex, 12, [t 'x (full A)']);
0488     t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0489 end
0490 
0491 t = 'LU4 : ';
0492 x = mplinsolve(A, b, 'LU4');
0493 t_is(x, ex, 12, [t 'x']);
0494 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0495 t = 'LU, nout = 4, vec = 1 : ';
0496 opt = struct('nout', 4, 'vec', 1);
0497 x = mplinsolve(A, b, 'LU', opt);
0498 t_is(x, ex, 12, [t 'x']);
0499 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0500 
0501 t = 'LU5 : ';
0502 x = mplinsolve(A, b, 'LU5');
0503 t_is(x, ex, 12, [t 'x']);
0504 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0505 t = 'LU, nout = 5, vec = 1 : ';
0506 opt = struct('nout', 5, 'vec', 1);
0507 x = mplinsolve(A, b, 'LU', opt);
0508 t_is(x, ex, 12, [t 'x']);
0509 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0510 
0511 t = 'LU3m : ';
0512 x = mplinsolve(A, b, 'LU3m');
0513 t_is(x, ex, 12, [t 'x']);
0514 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0515 if skipcrash
0516     t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
0517 else
0518     x = mplinsolve(full(A), b, 'LU3m');
0519     t_is(x, ex, 12, [t 'x (full A)']);
0520     t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0521 end
0522 t = 'LU, nout = 3, vec = 0, thresh = 1 : ';
0523 opt = struct('nout', 3, 'vec', 0, 'thresh', 1);
0524 x = mplinsolve(A, b, 'LU', opt);
0525 t_is(x, ex, 12, [t 'x']);
0526 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0527 if skipcrash
0528     t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
0529 else
0530     x = mplinsolve(full(A), b, 'LU', opt);
0531     t_is(x, ex, 12, [t 'x (full A)']);
0532     t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0533 end
0534 
0535 t = 'LU3am : ';
0536 x = mplinsolve(A, b, 'LU3am');
0537 t_is(x, ex, 12, [t 'x']);
0538 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0539 if skipcrash
0540     t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
0541 else
0542     x = mplinsolve(full(A), b, 'LU3am');
0543     t_is(x, ex, 12, [t 'x (full A)']);
0544     t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0545 end
0546 t = 'LU, nout = 3, vec = 0 : ';
0547 opt = struct('nout', 3, 'vec', 0);
0548 x = mplinsolve(A, b, 'LU', opt);
0549 t_is(x, ex, 12, [t 'x']);
0550 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0551 if skipcrash
0552     t_skip(2, [t 'potential MATLAB crash with non-sparse A']);
0553 else
0554     x = mplinsolve(full(A), b, 'LU', opt);
0555     t_is(x, ex, 12, [t 'x (full A)']);
0556     t_is(norm(b - A*x), 0, 12, [t '||b - A*x|| (full A)']);
0557 end
0558 
0559 t = 'LU4m : ';
0560 x = mplinsolve(A, b, 'LU4m');
0561 t_is(x, ex, 12, [t 'x']);
0562 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0563 t = 'LU, nout = 4, vec = 0 : ';
0564 opt = struct('nout', 4, 'vec', 0);
0565 x = mplinsolve(A, b, 'LU', opt);
0566 t_is(x, ex, 12, [t 'x']);
0567 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0568 
0569 t = 'LU5m : ';
0570 x = mplinsolve(A, b, 'LU5m');
0571 t_is(x, ex, 12, [t 'x']);
0572 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0573 t = 'LU, nout = 5, vec = 0 : ';
0574 opt = struct('nout', 5, 'vec', 0);
0575 x = mplinsolve(A, b, 'LU', opt);
0576 t_is(x, ex, 12, [t 'x']);
0577 t_is(norm(b - A*x), 0, 12, [t '||b - A*x||']);
0578 
0579 %% PARDISO
0580 if have_feature('pardiso')
0581     if have_feature('pardiso_object')
0582         tols = [6 5 12 12 6 5];     %% tolerances for PARDISO v6
0583     else
0584         tols = [13 13 12 12 1 2];   %% tolerances for PARDISO v5
0585     end
0586     vb = false;
0587 
0588     t = 'PARDISO (direct) : ';
0589     opt = struct('pardiso', struct('solver', 0, 'verbose', vb));
0590     x = mplinsolve(A, b, 'PARDISO', opt);
0591     t_is(x, ex, tols(1), [t 'x']);
0592     t_is(norm(b - A*x), 0, tols(2), [t '||b - A*x||']);
0593 
0594     t = 'PARDISO (direct, symmetric indefinite) : ';
0595     opt = struct('pardiso', struct('solver', 0, 'mtype', -2, 'verbose', vb));
0596     x = mplinsolve(A, b, 'PARDISO', opt);
0597     t_is(x, ex, tols(3), [t 'x']);
0598     t_is(norm(b - A*x), 0, tols(4), [t '||b - A*x||']);
0599 
0600     t = 'PARDISO (iterative) : ';
0601     opt = struct('pardiso', struct('solver', 1, 'verbose', vb));
0602     [x, info] = mplinsolve(A, b, 'PARDISO', opt);
0603     t_is(x, ex, tols(5), [t 'x']);
0604     t_is(norm(b - A*x), 0, tols(6), [t '||b - A*x||']);
0605 else
0606     t_skip(6, [t ' not available']);
0607 end
0608 
0609 if isoctave
0610     warning(s.state, lu_warning_id);
0611 end
0612 
0613 t_end;