0001 function t_mplinsolve(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
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
0026
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
0580 if have_feature('pardiso')
0581 if have_feature('pardiso_object')
0582 tols = [6 5 12 12 6 5];
0583 else
0584 tols = [13 13 12 12 1 2];
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;