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(4, quiet);
0017
0018 ijs = [
0019 1 1 1205.63;
0020 4 1 -1205.63;
0021 25 1 17.3611;
0022 28 1 -17.3611;
0023 43 1 1;
0024 2 2 1024;
0025 8 2 -1024;
0026 26 2 16;
0027 32 2 -16;
0028 3 3 1164.84;
0029 6 3 -1164.84;
0030 27 3 17.0648;
0031 30 3 -17.0648;
0032 1 4 -1205.63;
0033 4 4 2201.34;
0034 5 4 -453.695;
0035 9 4 -542.009;
0036 13 4 0.776236;
0037 14 4 -0.451721;
0038 18 4 -0.324515;
0039 25 4 -17.3611;
0040 28 4 39.4759;
0041 29 4 -10.5107;
0042 33 4 -11.6041;
0043 37 4 -3.30738;
0044 38 4 1.94219;
0045 42 4 1.36519;
0046 4 5 -453.695;
0047 5 5 581.372;
0048 6 5 -127.677;
0049 13 5 -0.451721;
0050 14 5 0.568201;
0051 15 5 -0.11648;
0052 28 5 -10.5107;
0053 29 5 16.0989;
0054 30 5 -5.58824;
0055 37 5 1.94219;
0056 38 5 -3.2242;
0057 39 5 1.28201;
0058 3 6 -1164.84;
0059 5 6 -127.677;
0060 6 6 1676.74;
0061 7 6 -384.227;
0062 14 6 -0.11648;
0063 15 6 0.163059;
0064 16 6 -0.0465791;
0065 27 6 -17.0648;
0066 29 6 -5.58824;
0067 30 6 32.4374;
0068 31 6 -9.78427;
0069 38 6 1.28201;
0070 39 6 -2.4371;
0071 40 6 1.15509;
0072 6 7 -384.227;
0073 7 7 1141.16;
0074 8 7 -756.935;
0075 15 7 -0.0465791;
0076 16 7 0.371828;
0077 17 7 -0.325249;
0078 30 7 -9.78427;
0079 31 7 23.4822;
0080 32 7 -13.698;
0081 39 7 1.15509;
0082 40 7 -2.77221;
0083 41 7 1.61712;
0084 2 8 -1024;
0085 7 8 -756.935;
0086 8 8 1925.77;
0087 9 8 -144.836;
0088 16 8 -0.325249;
0089 17 8 0.844105;
0090 18 8 -0.518855;
0091 26 8 -16;
0092 31 8 -13.698;
0093 32 8 35.6731;
0094 33 8 -5.97513;
0095 40 8 1.61712;
0096 41 8 -2.80473;
0097 42 8 1.1876;
0098 4 9 -542.009;
0099 8 9 -144.836;
0100 9 9 686.845;
0101 13 9 -0.324515;
0102 17 9 -0.518855;
0103 18 9 0.84337;
0104 28 9 -11.6041;
0105 32 9 -5.97513;
0106 33 9 17.5792;
0107 37 9 1.36519;
0108 41 9 1.1876;
0109 42 9 -2.55279;
0110 10 10 1207.63;
0111 13 10 -1205.63;
0112 34 10 17.3611;
0113 37 10 -17.3611;
0114 11 11 1026;
0115 17 11 -1024;
0116 35 11 16;
0117 41 11 -16;
0118 12 12 1166.84;
0119 15 12 -1164.84;
0120 36 12 17.0648;
0121 39 12 -17.0648;
0122 4 13 0.776236;
0123 5 13 -0.451721;
0124 9 13 -0.324515;
0125 10 13 -1205.63;
0126 13 13 2190.83;
0127 14 13 -447.892;
0128 18 13 -535.137;
0129 28 13 3.30738;
0130 29 13 -1.94219;
0131 33 13 -1.36519;
0132 34 13 -17.3611;
0133 37 13 39.1419;
0134 38 13 -10.5107;
0135 42 13 -11.6041;
0136 4 14 -0.451721;
0137 5 14 0.568201;
0138 6 14 -0.11648;
0139 13 14 -447.892;
0140 14 14 573.221;
0141 15 14 -122.862;
0142 28 14 -1.94219;
0143 29 14 3.2242;
0144 30 14 -1.28201;
0145 37 14 -10.5107;
0146 38 14 15.5829;
0147 39 14 -5.58824;
0148 5 15 -0.11648;
0149 6 15 0.163059;
0150 7 15 -0.0465791;
0151 12 15 -1164.84;
0152 14 15 -122.862;
0153 15 15 1669.87;
0154 16 15 -379.651;
0155 29 15 -1.28201;
0156 30 15 2.4371;
0157 31 15 -1.15509;
0158 36 15 -17.0648;
0159 38 15 -5.58824;
0160 39 15 31.8704;
0161 40 15 -9.78427;
0162 6 16 -0.0465791;
0163 7 16 0.371828;
0164 8 16 -0.325249;
0165 15 16 -379.651;
0166 16 16 1131.92;
0167 17 16 -750.073;
0168 30 16 -1.15509;
0169 31 16 2.77221;
0170 32 16 -1.61712;
0171 39 16 -9.78427;
0172 40 16 23.1242;
0173 41 16 -13.698;
0174 7 17 -0.325249;
0175 8 17 0.844105;
0176 9 17 -0.518855;
0177 11 17 -1024;
0178 16 17 -750.073;
0179 17 17 1914.92;
0180 18 17 -138.499;
0181 31 17 -1.61712;
0182 32 17 2.80473;
0183 33 17 -1.1876;
0184 35 17 -16;
0185 40 17 -13.698;
0186 41 17 35.2181;
0187 42 17 -5.97513;
0188 4 18 -0.324515;
0189 8 18 -0.518855;
0190 9 18 0.84337;
0191 13 18 -535.137;
0192 17 18 -138.499;
0193 18 18 676.012;
0194 28 18 -1.36519;
0195 32 18 -1.1876;
0196 33 18 2.55279;
0197 37 18 -11.6041;
0198 41 18 -5.97513;
0199 42 18 17.0972;
0200 19 19 1.88667;
0201 25 19 -1;
0202 20 20 1.54931;
0203 26 20 -1;
0204 21 21 1.78346;
0205 27 21 -1;
0206 22 22 0.666667;
0207 34 22 -1;
0208 23 23 0.666667;
0209 35 23 -1;
0210 24 24 0.666667;
0211 36 24 -1;
0212 1 25 17.3611;
0213 4 25 -17.3611;
0214 19 25 -1;
0215 2 26 16;
0216 8 26 -16;
0217 20 26 -1;
0218 3 27 17.0648;
0219 6 27 -17.0648;
0220 21 27 -1;
0221 1 28 -17.3611;
0222 4 28 39.4759;
0223 5 28 -10.5107;
0224 9 28 -11.6041;
0225 13 28 3.30738;
0226 14 28 -1.94219;
0227 18 28 -1.36519;
0228 4 29 -10.5107;
0229 5 29 16.0989;
0230 6 29 -5.58824;
0231 13 29 -1.94219;
0232 14 29 3.2242;
0233 15 29 -1.28201;
0234 3 30 -17.0648;
0235 5 30 -5.58824;
0236 6 30 32.4374;
0237 7 30 -9.78427;
0238 14 30 -1.28201;
0239 15 30 2.4371;
0240 16 30 -1.15509;
0241 6 31 -9.78427;
0242 7 31 23.4822;
0243 8 31 -13.698;
0244 15 31 -1.15509;
0245 16 31 2.77221;
0246 17 31 -1.61712;
0247 2 32 -16;
0248 7 32 -13.698;
0249 8 32 35.6731;
0250 9 32 -5.97513;
0251 16 32 -1.61712;
0252 17 32 2.80473;
0253 18 32 -1.1876;
0254 4 33 -11.6041;
0255 8 33 -5.97513;
0256 9 33 17.5792;
0257 13 33 -1.36519;
0258 17 33 -1.1876;
0259 18 33 2.55279;
0260 10 34 17.3611;
0261 13 34 -17.3611;
0262 22 34 -1;
0263 11 35 16;
0264 17 35 -16;
0265 23 35 -1;
0266 12 36 17.0648;
0267 15 36 -17.0648;
0268 24 36 -1;
0269 4 37 -3.30738;
0270 5 37 1.94219;
0271 9 37 1.36519;
0272 10 37 -17.3611;
0273 13 37 39.1419;
0274 14 37 -10.5107;
0275 18 37 -11.6041;
0276 4 38 1.94219;
0277 5 38 -3.2242;
0278 6 38 1.28201;
0279 13 38 -10.5107;
0280 14 38 15.5829;
0281 15 38 -5.58824;
0282 5 39 1.28201;
0283 6 39 -2.4371;
0284 7 39 1.15509;
0285 12 39 -17.0648;
0286 14 39 -5.58824;
0287 15 39 31.8704;
0288 16 39 -9.78427;
0289 6 40 1.15509;
0290 7 40 -2.77221;
0291 8 40 1.61712;
0292 15 40 -9.78427;
0293 16 40 23.1242;
0294 17 40 -13.698;
0295 7 41 1.61712;
0296 8 41 -2.80473;
0297 9 41 1.1876;
0298 11 41 -16;
0299 16 41 -13.698;
0300 17 41 35.2181;
0301 18 41 -5.97513;
0302 4 42 1.36519;
0303 8 42 1.1876;
0304 9 42 -2.55279;
0305 13 42 -11.6041;
0306 17 42 -5.97513;
0307 18 42 17.0972;
0308 1 43 1;
0309 ];
0310 A = sparse(ijs(:, 1), ijs(:, 2), ijs(:, 3));
0311 b = [
0312 0;
0313 0;
0314 0;
0315 0;
0316 0;
0317 0;
0318 0;
0319 0;
0320 0;
0321 0;
0322 0;
0323 0;
0324 -0.00896054;
0325 -0.0617829;
0326 -0.0772931;
0327 -0.0230638;
0328 -0.0185934;
0329 -0.02;
0330 -0.336;
0331 -0.2755;
0332 -0.353;
0333 0;
0334 0;
0335 0;
0336 1.3;
0337 1.55;
0338 1.4;
0339 0;
0340 -0.9;
0341 0;
0342 -1;
0343 0;
0344 -1.25;
0345 0;
0346 0;
0347 0;
0348 0.167;
0349 -0.042;
0350 0.2835;
0351 -0.171;
0352 0.2275;
0353 -0.259;
0354 0;
0355 ];
0356 ex = [
0357 0;
0358 0.04612219791119214;
0359 0.0509334351882598;
0360 -0.05953571031927305;
0361 -0.09461814898578046;
0362 -0.008909854578010561;
0363 -0.05785829019394401;
0364 -0.02232729212460287;
0365 -0.1137760871247425;
0366 -0.03062777824802364;
0367 -0.01013282572376477;
0368 -0.005330939680091628;
0369 -0.02914165388753019;
0370 -0.03376073204420303;
0371 4.021341450281111e-05;
0372 -0.01727289094763518;
0373 -0.008382063634320435;
0374 -0.04854008812629265;
0375 -0.2663945795760685;
0376 -0.4548081594272799;
0377 -0.3787862287965494;
0378 -0.02580075363496279;
0379 -0.02801219343110931;
0380 -0.09165765332863518;
0381 -0.1665986614487812;
0382 -0.4291388294822789;
0383 -0.3225500876094941;
0384 3.967864472606993;
0385 5.577372159790927;
0386 3.762341481664236;
0387 5.858342308599034;
0388 3.951628532808602;
0389 6.471657726723339;
0390 -0.01720051102355974;
0391 -0.01867480495813735;
0392 -0.06110513277164123;
0393 -0.1239317474796463;
0394 -0.1848327901217308;
0395 -0.4283638085291988;
0396 0.07640167050820651;
0397 -0.1319901818980452;
0398 0.4366406661538687;
0399 0.0007894844305239289;
0400 ];
0401
0402 t = '''''';
0403 x = mplinsolve(A, b, '', []);
0404 t_is(x, ex, 12, t);
0405
0406 t = '\';
0407 x = mplinsolve(A, b, '\', []);
0408 t_is(x, ex, 12, t);
0409
0410 t = 'PARDISO';
0411 if exist('have_fcn', 'file') && have_fcn('pardiso') || have_pardiso()
0412 x = mplinsolve(A, b, 'PARDISO', struct('solver', 0));
0413 t_is(x, ex, 12, [t ' : direct']);
0414
0415 [x, info] = mplinsolve(A, b, 'PARDISO', struct('solver', 1));
0416 t_is(x, ex, 1, [t ' : iterative']);
0417 else
0418 t_skip(2, [t ' not available']);
0419 end
0420
0421 t_end;
0422
0423
0424 function TorF = have_pardiso()
0425 TorF = exist('pardisoinit', 'file') == 3 && ...
0426 exist('pardisoreorder', 'file') == 3 && ...
0427 exist('pardisofactor', 'file') == 3 && ...
0428 exist('pardisosolve', 'file') == 3 && ...
0429 exist('pardisofree', 'file') == 3;
0430 if TorF
0431 try
0432 A = sparse([1 2; 3 4]);
0433 b = [1;1];
0434 info = pardisoinit(11, 0);
0435 info = pardisoreorder(A, info, false);
0436 info = pardisofactor(A, info, false);
0437 [x, info] = pardisosolve(A, b, info, false);
0438 pardisofree(info);
0439 if any(x ~= [-1; 1])
0440 TorF = 0;
0441 end
0442 catch
0443 TorF = 0;
0444 end
0445 end