Linear solves with cached factorization sometimes gives different answers

[MarcBerliner]

When performing a linear solve, I occasionally get the incorrect answer only if I use a cached factorization. It is somewhat difficult to reliably reproduce this, but here's an example (the only minimal working example I could find is rather large...). I am using v0.2.3.

Thank you for the great work!

using KLU, LinearAlgebra, SparseArrays

A = sparse([1, 31, 61, 2, 32, 62, 3, 33, 63, 4, 34, 64, 5, 35, 65, 6, 36, 66, 7, 37, 67, 8, 38, 68, 9, 39, 69, 10, 40, 70, 11, 51, 71, 12, 52, 72, 13, 53, 73, 14, 54, 74, 15, 55, 75, 16, 56, 76, 17, 57, 77, 18, 58, 78, 19, 59, 79, 20, 80, 21, 51, 71, 22, 52, 72, 23, 53, 73, 24, 54, 74, 25, 55, 75, 26, 56, 76, 27, 57, 77, 28, 58, 78, 29, 59, 79, 30, 80, 1, 31, 32, 2, 31, 32, 33, 3, 32, 33, 34, 4, 33, 34, 35, 5, 34, 35, 36, 6, 35, 36, 37, 7, 36, 37, 38, 8, 37, 38, 39, 9, 38, 39, 40, 10, 39, 40, 41, 40, 41, 42, 41, 42, 43, 42, 43, 44, 43, 44, 45, 44, 45, 46, 45, 46, 47, 46, 47, 48, 47, 48, 49, 48, 49, 50, 49, 50, 51, 11, 21, 50, 51, 52, 12, 22, 51, 52, 53, 13, 23, 52, 53, 54, 14, 24, 53, 54, 55, 15, 25, 54, 55, 56, 16, 26, 55, 56, 57, 17, 27, 56, 57, 58, 18, 28, 57, 58, 59, 19, 29, 58, 59, 20, 30, 59, 60, 1, 61, 62, 81, 2, 61, 62, 63, 3, 62, 63, 64, 4, 63, 64, 65, 5, 64, 65, 66, 6, 65, 66, 67, 7, 66, 67, 68, 8, 67, 68, 69, 9, 68, 69, 70, 10, 69, 70, 11, 21, 71, 72, 12, 22, 71, 72, 73, 13, 23, 72, 73, 74, 14, 24, 73, 74, 75, 15, 25, 74, 75, 76, 16, 26, 75, 76, 77, 17, 27, 76, 77, 78, 18, 28, 77, 78, 79, 19, 29, 78, 79, 80, 20, 30, 79, 80, 81, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 61, 80], [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 27, 27, 27, 28, 28, 28, 29, 29, 29, 30, 30, 31, 31, 31, 32, 32, 32, 32, 33, 33, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 37, 38, 38, 38, 38, 39, 39, 39, 39, 40, 40, 40, 40, 41, 41, 41, 42, 42, 42, 43, 43, 43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 56, 56, 56, 56, 56, 57, 57, 57, 57, 57, 58, 58, 58, 58, 58, 59, 59, 59, 59, 60, 60, 60, 60, 61, 61, 61, 61, 62, 62, 62, 62, 63, 63, 63, 63, 64, 64, 64, 64, 65, 65, 65, 65, 66, 66, 66, 66, 67, 67, 67, 67, 68, 68, 68, 68, 69, 69, 69, 69, 70, 70, 70, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74, 74, 74, 74, 74, 75, 75, 75, 75, 75, 76, 76, 76, 76, 76, 77, 77, 77, 77, 77, 78, 78, 78, 78, 78, 79, 79, 79, 79, 79, 80, 80, 80, 80, 80, 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, 81, 81], [-1.0, -683116.1514329641, -0.09262591883836804, -1.0, -683116.1514329641, -0.09262591883836804, -1.0, -683116.1514329641, -0.09262591883836804, -1.0, -683116.1514329641, -0.09262591883836804, -1.0, -683116.1514329641, -0.09262591883836804, -1.0, -683116.1514329641, -0.09262591883836804, -1.0, -683116.1514329641, -0.09262591883836804, -1.0, -683116.1514329641, -0.09262591883836804, -1.0, -683116.1514329641, -0.09262591883836804, -1.0, -683116.1514329641, -0.09262591883836804, -2.000520644329273, -614387.719654895, -0.11207736179442528, -1.8318004565292179, -614387.719654895, -0.11207736179442528, -1.6943487912278106, -614387.719654895, -0.11207736179442528, -1.6138526171787806, -614387.719654895, -0.11207736179442528, -1.5621234147902863, -614387.719654895, -0.11207736179442528, -1.5266098779183053, -614387.719654895, -0.11207736179442528, -1.5019998539923227, -614387.719654895, -0.11207736179442528, -1.4857507347573913, -614387.719654895, -0.11207736179442528, -1.476764774286798, -614387.719654895, -0.11207736179442528, -1.4749552780455253, -0.11207736179442528, 1.000635908760994, -614387.719654895, -0.11207736179442528, 1.000341381245597, -614387.719654895, -0.11207736179442528, 1.0002177711736213, -614387.719654895, -0.11207736179442528, 1.0001537089822035, -614387.719654895, -0.11207736179442528, 1.000116496235308, -614387.719654895, -0.11207736179442528, 1.0000936546628076, -614387.719654895, -0.11207736179442528, 1.0000794799543697, -614387.719654895, -0.11207736179442528, 1.0000711443528876, -614387.719654895, -0.11207736179442528, 1.0000673012639298, -614387.719654895, -0.11207736179442528, 1.0000675900694103, -0.11207736179442528, -0.000536898323310069, 326.7883746460036, -326.7883746460036, -0.0005425929078755278, -326.7883746460036, 655.5469266773066, -328.7585520313031, -0.0005551714391717656, -328.7585520313031, 661.339844949825, -332.58129291852185, -0.00057752612470627, -332.58129291852185, 671.8404569229615, -339.2591640044397, -0.0006157970631492979, -339.2591640044397, 689.9827378435225, -350.72357383908286, -0.000682458509509591, -350.72357383908286, 721.4321060587135, -370.70853221963074, -0.0007934959998253688, -370.70853221963074, 776.2311579477577, -405.52262572812697, -0.0009092535963912138, -405.52262572812697, 865.2985073358996, -459.77588160777265, -0.0009187915012681606, -459.77588160777265, 983.4311390377767, -523.6552574300041, -0.0008320654200983155, -523.6552574300041, 1609.7717942986123, -1086.1165368686081, -1086.1165368686081, 24100.03231420114, -23013.91577733253, -23013.91577733253, 46045.17516614906, -23031.259388816532, -23031.259388816532, 46078.36731246605, -23047.10792364951, -23047.10792364951, 46108.562267269735, -23061.454343620226, -23061.454343620226, 46135.746133436354, -23074.29178981613, -23074.29178981613, 46159.90537132222, -23085.613581506088, -23085.613581506088, 46181.02679657345, -23095.413215067365, -23095.413215067365, 46199.097578011046, -23103.684362943684, -23103.684362943684, 46214.10523556625, -23110.42087262256, -23110.42087262256, 25610.206337430296, -2499.785464807739, -0.0010369665689314224, 6.590729834299752e-7, -2499.785464807739, 3814.9113911696522, -1315.125926361913, -0.000862100419729013, 3.538167270847678e-7, -1315.125926361913, 2615.785944757199, -1300.6600183952864, -0.0007196418087567145, 2.257039160431882e-7, -1300.6600183952864, 2584.5740703994124, -1283.914052004126, -0.0006362134040104998, 1.5930813367572092e-7, -1283.914052004126, 2551.281920012776, -1267.3678680086502, -0.0005825998641051258, 1.207398394611211e-7, -1267.3678680086502, 2519.4587907000673, -1252.0909226914173, -0.000545792677057987, 9.706621785144531e-8, -1252.0909226914173, 2490.863093168667, -1238.77217047725, -0.000520286184765981, 8.237516783233921e-8, -1238.77217047725, 2466.6882104223037, -1227.9160399450536, -0.0005034451595024182, 7.373592570084611e-8, -1227.9160399450536, 2447.8144639772117, -1219.8984240321581, -0.0004941318678034394, 6.975284468308397e-8, -1219.8984240321581, 2434.881524078571, -0.000492256457122771, 7.005217047272613e-8, -1214.9831000464135, 1.0, 0.000536898323310069, -1.0, 1.0, 1.0, 0.0005425929078755278, 1.0, -2.0, 1.0, 0.0005551714391717656, 1.0, -2.0, 1.0, 0.00057752612470627, 1.0, -2.0, 1.0, 0.0006157970631492979, 1.0, -2.0, 1.0, 0.000682458509509591, 1.0, -2.0, 1.0, 0.0007934959998253688, 1.0, -2.0, 1.0, 0.0009092535963912138, 1.0, -2.0, 1.0, 0.0009187915012681606, 1.0, -2.0, 1.0, 0.0008320654200983155, 1.0, -1.0, 0.0010369665689314224, -6.590729834299752e-7, -1.0, 1.0, 0.000862100419729013, -3.538167270847678e-7, 1.0, -2.0, 1.0, 0.0007196418087567145, -2.257039160431882e-7, 1.0, -2.0, 1.0, 0.0006362134040104998, -1.5930813367572092e-7, 1.0, -2.0, 1.0, 0.0005825998641051258, -1.207398394611211e-7, 1.0, -2.0, 1.0, 0.000545792677057987, -9.706621785144531e-8, 1.0, -2.0, 1.0, 0.000520286184765981, -8.237516783233921e-8, 1.0, -2.0, 1.0, 0.0005034451595024182, -7.373592570084611e-8, 1.0, -2.0, 1.0, 0.0004941318678034394, -6.975284468308397e-8, 1.0, -2.0, 1.0, 0.000492256457122771, -7.005217047272613e-8, 1.0, -1.0, -1.0, 3.5756529395451995e-8, 1.9195531556317828e-8, 1.2245056578758705e-8, 8.6429015291434e-9, 6.55046483853513e-9, 5.266105349743969e-9, 4.469075144039891e-9, 4.000372116416478e-9, 3.784278492271054e-9, 3.8005169499489525e-9, 3.963390812251965e-6, -5.3321737917734765e-6], 81, 81)
b = [1.9737985253393334e-17, 4.955672137029691e-17, 3.8330061984093687e-16, 4.097847566510897e-15, 5.091398240102173e-14, 7.201399273724842e-13, 1.0823055864343182e-11, 1.3538632618015052e-10, 1.1473095071239852e-9, 6.963765873676213e-9, -7.809576078499119e-8, -1.8724934918801236e-7, -5.443895888331384e-9, -2.382154674036239e-9, -1.5494741575172198e-9, -5.305814444704892e-11, -4.450552797594077e-11, -1.6371748976377388e-11, -4.0694323080834195e-12, -9.420450170643819e-9, -1.7063272063410017e-9, 3.293021719749384e-11, -5.74928909471849e-12, -4.348690790493546e-12, -1.7062597777500879e-12, 9.901863891959797e-13, 1.526891898033429e-12, 1.834890322431126e-12, 2.088784859107366e-12, 2.8255276799942745e-11, -5.0546372643012205e-14, -5.870304242705515e-15, -6.56627530126741e-14, 5.088984789125561e-14, -2.3869795029440866e-14, 7.593925488436071e-14, -3.441691376337985e-14, -3.6415315207705135e-14, -6.106226635438361e-14, 2.842170943040401e-14, 4.547473508864641e-13, -1.0378017889500768e-12, -1.0146605777805462e-12, 4.1424502716935763e-13, -8.197331702319843e-13, 1.355811296566145e-12, -3.6186331708876196e-13, 4.7004067305067565e-14, -9.705569681273118e-13, 4.547473508864641e-13, -2.842170943040401e-14, 3.375077994860476e-14, -1.454392162258955e-14, 7.993605777301127e-15, -1.765254609153999e-14, 4.052314039881821e-15, 1.0547118733938987e-15, -4.6074255521944e-15, -2.220446049250313e-16, 0.0, 0.0, 1.7763568394002505e-15, 8.881784197001252e-16, -8.881784197001252e-16, -8.881784197001252e-16, 0.0, -8.881784197001252e-16, 0.0, 0.0, 0.0, -2.7755575615628914e-17, 0.0, -5.551115123125783e-17, 0.0, 0.0, 2.7755575615628914e-17, 0.0, 2.7755575615628914e-17, -2.1006716490468147e-6, 2.9570568886860826e-6, 0.0];

# Factorizing the sparsity pattern of A
A_sparsity = deepcopy(A)
A_sparsity.nzval .= ones(length(A_sparsity.nzval))

factor = klu(A_sparsity)

# updating the factorization with the real values of A
klu!(factor,A)

# I get incorrect results only with the cached factorization
[factor\b klu(A)\b lu(A)\b]
81×3 Matrix{Float64}:
  6.1124e-12   -8.33057e-12  -8.33057e-12
  9.45051e-12  -1.2095e-11   -1.2095e-11
  2.38783e-11  -3.00322e-11  -3.00322e-11
  6.77289e-11  -8.49574e-11  -8.49574e-11
  2.00907e-10  -2.51996e-10  -2.51996e-10
  6.2692e-10   -7.87668e-10  -7.87668e-10
  2.08081e-9   -2.63565e-9   -2.63565e-9
  6.8747e-9    -8.97056e-9   -8.97056e-9
  1.94325e-8   -2.75142e-8   -2.75142e-8
  4.34986e-8   -7.45746e-8   -7.45746e-8
  ⋮                          
 -0.00298534   -0.00762745   -0.00762745
 -0.0029854    -0.00762748   -0.00762748
 -0.00298545   -0.00762754   -0.00762754
 -0.00298551   -0.00762764   -0.00762764
 -0.00298559   -0.0076278    -0.0076278
 -0.0029857    -0.00762804   -0.00762804
 -0.00298586   -0.00762841   -0.00762841
 -0.00324806   -0.00763106   -0.00763106
  0.00170187   -0.00268456   -0.00268456

[rayegun]

I'll take a look at this, thanks!

[rayegun]

@DrTimothyAldenDavis do you have any thoughts here? Is this just a failed refactor? I'm not sure how low rgrowth would be in a failed refactor, but this refactor gives an rgrowth of 6.21e-13 which seems really low. Should I be checking this for the user?

I'm not immediately seeing any errors in the wrapper, but I'll keep looking.

[DrTimothyAldenDavis]

That is a low number. You're likely seeing roundoff errors amplified into garbage because the matrix is too ill conditioned. Check the condition # of the matrix.

[rayegun]

Using the newly wrapped KLU.condest we do indeed see a very large condition number for that matrix. So I'll close this and recommend checking KLU.condest or KLU.rgrowth or KLU.rcond in the future.

[oscardssmith]

I think this should likely be re-opened. why should klu_l_refactor struggle when klu_l_factorworks fine? The condest is around 1e11 for both the out of place and the in place version, so it doesn't make sense (to me at least) that the condition of the matrix is responsible.

[DrTimothyAldenDavis]

klu_factor does numerical partial pivoting. klu_refactor does not. If the change in values requires a revision to the numerical pivoting, then klu_refactor should not be used.