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