3
$\begingroup$

While trying to calculate the SmithDecomposition of an integer 293x329 matrix with less than 3000 non-zero elements in the range -4 to 3 several errors occur, none of them which seem to make any sense at all:

In[254]:= SmithDecomposition[tmat]

During evaluation of In[254]:= Part::partw: Part 296 of {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,<<245>>} does not exist.

During evaluation of In[254]:= Drop::drop: Cannot drop positions 1 through 296 in {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,<<245>>}.

During evaluation of In[254]:= Part::partw: Part 297 of {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,<<245>>} does not exist.

During evaluation of In[254]:= Part::partw: Part 298 of {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,<<245>>} does not exist.

During evaluation of In[254]:= General::stop: Further output of Part::partw will be suppressed during this calculation.

During evaluation of In[254]:= HermiteDecomposition::latm: Matrix contains an entry that is not rational.

During evaluation of In[254]:= Set::shape: Lists {System`SmithDecompositionDump`tmpu$84833,System`SmithDecompositionDump`hnf$84833} and HermiteDecomposition[{{1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,<<279>>},{0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,<<279>>},{0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,<<279>>},<<46>>,{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,<<279>>},<<2546>>}] are not the same shape.

During evaluation of In[254]:= HermiteDecomposition::latm: Matrix contains an entry that is not rational.

During evaluation of In[254]:= Set::shape: Lists {System`SmithDecompositionDump`tmpv$84833,System`SmithDecompositionDump`hnf$84833} and HermiteDecomposition[{{1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,<<2546>>},{0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,<<2546>>},{0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,<<2546>>},<<46>>,{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,<<2546>>},<<279>>}] are not the same shape.

During evaluation of In[254]:= HermiteDecomposition::latm: Matrix contains an entry that is not rational.

During evaluation of In[254]:= General::stop: Further output of HermiteDecomposition::latm will be suppressed during this calculation.

During evaluation of In[254]:= Set::shape: Lists {System`SmithDecompositionDump`tmpu$84833,System`SmithDecompositionDump`hnf$84833} and HermiteDecomposition[{{1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,<<2546>>},{0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,<<2546>>},{0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,<<2546>>},<<46>>,{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,<<2546>>},<<279>>}] are not the same shape.

During evaluation of In[254]:= General::stop: Further output of Set::shape will be suppressed during this calculation.

Out[254]= $Aborted

I checked that all elements of the matrix are integers and interestingly the HermiteDecomposition does work. Trying to calculate the SmithDecomposition of the Hermite upper diagonal form results in the same error though.

Does anybody know what is going on?

The matrix in compressed form (use Uncompress to get the actual matrix) is the following

"1:eJztXcnO5LYRHi+XwE+RS24BWtRGHX00kFPyBE6QAeYQOJhxDnkXP2y6xarB/3d/\
bPFjiWxKzqXF1kLWzmJxqT/+/Ze/fvz+w4cPX364/vzt3z9//vLPHz9//vm/X/5w/f/\
jf3795V8///rpHx+/vb1ye+8vn778+tN314Jb+vXau+XT9fLh43dvX/n0ze3Wu68+/\
umbt3XcXrisv93663+6vd3168VdwkX+Deul78JlCpd5vQzhzWFcL6MLF/\
kX3pzCd1P4N4d/c3hzDrX48M/LvwCLX9bLMqx4dpdOrvp/DNfOyXUO1yvs4SrvOx+\
ufaBXN8jzQf6PUt8k9c1TuHr5L+07ac914blzoT3Xh/rd0MlVno9LuM7ynpf3lvC/\
7wIcfR++64fwfj8GfPpJrnNov/fC7yXUP1xCPYPAM0g9g+A1SD3DOMpV3p/lPR/+\
j0LP0YX3x17uy/fjFOofBY5RvxM4pktobxI4Jqln6kO70xjwnIS+\
k9B3WgIcs9B17kL9s3w/Cz6z8G0e5L7QaR4DPedZ7nt534f6Z61/Ce/7S3jPd/\
Jf8PTCLz92cpX/Aq8XvP0s3wkf/RLqW0QuF+HnIvUvbpCr/\
Bc5WYZRrvJf5GSZ5XvBY1nk/xX+\
71d5v3gpdHqnv0jhKntScFoYpDDqy1f5WQs3RQoFp3f0q27UO6O01c36std3rnz/\
PuiZFtxFC3qnl0bdKGC4WWp2yxgKfSc194pFP+\
ijUQuTFuZJC7MUvBSGi7Q1OC0o7oO2PszSxOClMF4EntENWpC2RiXmOAnKo9c7i7w8dfJo\
0kZvYh4Kk5Bu0kYnL8BPSrq5E+Bvwh0Ko1Q4T/pICX4T5LXgFWZ/EVB9J215ZYF3+\
nIv9Xglglca3iR3LSydPFoUi2XotaB3Jr0zS6OL10dLaN1duosWRin0emeYpTDpI6+\
PRBJcd5F3uk4LTgu9FoS8V+\
vrpCAS5TqvBeHOVR47LQiEzumdXguDFsavhUUKwgLnZi14LSxS6C8CT6+\
4904LvYDRD18LAkY/\
akHU6loQIvReXh6Ey1cABR4VbDf08vIw6J1J6hmUGsOsjxTmQdk0KllU+J0K/\
1VNtaAVjor7KNLiJtF3N3Xy+aRsmnq9IyrsplEfKTEnBUzVwU3Kr1lrvtn+UFC+zyo/\
s1Y4qyDNs36upJu1Qq8V+k4LvdTsBy2IkXFeFPbKZGnCK+m8snsRI+\
MWZcqiErUoERYVraX/\
ekdfVmFblD6LSt0iJvda0EfKgmX6ekcrVFCX2X387Ye3zt3NhXvv3HXgnru7FxyR+\
zfD3R7eHeDdEd11+O4E4MJtEXBF2ppBWxiDujTwAC7m+wV8j+h6j38r/L7Hn2/rvtY4v+\
9ptQ5zLkSl98CuFWAa7o/B2lhdtcWAQTJiXSTIeA/\
XehOCVYq28FWqMahkGF2zPibaubUxCME9G3kIMLpIPBzUskKMdHaNTGeku9fIEC1BNzEAh\
cDCmlPKUEArjiUB8zHdUGCW427ktV35ChjWsjJM7619WX9PsPVm3T4nXZR6u6+\
YbsH6us5ies/\
bW7zFEMUFr17HOfAuwvc6ljHTBjMCyzOsF1vXBzKEeHUhhAkGF0P4npkrwkinTo0wsrp7I\
FwKCdSnDRjcgyGBATuH6KEudEC95akRRgb9LN0H6i8H9OpZEEb95YCEzE3Y5YZGAHoYHuO\
LMUv3YQfYu3uslviu2a98GHCHu3AsEGEadLofxpU0caAnMEMV9lCiPaPuDHlxhKCQfx0Jv\
9Ri28oK6KNQrGAMER7E2cMyaWRY14BAwDDCFBJ2cLGFhEhAgz5jfxMbf2Y4bUZ4B5EmiAP\
1fcbWH4s0JBllc9KJzhl6zAq7JWIYlMaKZ11FejBphEQvxkzMtmYZVJcV0NDPGLB7VoRVe\
MkVcKbXbp8sliisLCQMcimLUQvhJ6YXe71WhCvzEgpkxHXHqNldd7tpsPAy7i+\
OUKSxjaUEHdtjexDCQpx1lXApp7emSE91HZqqkxQFTFlc/\
CeIsIcdm5mQbZABrYiYsOS0YPbSI2gT1GEioDSlRaTCXgPoH9R1bgsNyakYRLq7OUFfohj\
JsGLjOEjViEcp8hIxuLrkJexmzPWvayEJohPRtmaJzsk007HZQ5zprJiJOGDD5gVrhSVIS\
hMSxyJb6NrsQ+qqhMRxS2LAFRmjMkGtUuEr+\
0whcscjUoYH8GUUM0IyGACbcQTNzmEmCtEsL4noYBNGFgXAsIGrK5BY9LA3bRc9Jm5Tapm\
iWfSwL1yLaysI2DNsgT+\
FugSCP9hBiDgTdrehVBjDHqe1TuFg4jQx0mhhZocgJDPHXMrxsJtIAmHGn28hwmIfYZYJA\
WDUmlBBc4Sl0MAIh9fN5r+yAmFwmZkD+7qj9Jlcj2dyMbiFpjQocUJOlYceHAzQ+\
1Jr9Zl9QJaFs+vJROkxrbqL7jDFiq0yJHZeMkjgPv/\
1q1aLDbfMq0sxuFSYyrztj1oNhUwDpG5syQPBH9wY7uuYXdJ4mAElmohWVJ63RqsYFhzDw\
BO2EF88zm9hvtU+\
9LYPf8osG1xwZOP1MwG14x37r0hhWHG7212wg8CsG2emDfZ3TQMWdh8f0wy3hgdWxJCcUL\
fugntMRn5LBfjTesyABe4y627hx1RnsIB9MbRRe0i1fQP865eSdpf0cWOsc8LUKbV5vNb4\
LGBRZkF4uX4IYgGHjpUldQcsMGS4Bntr6RO3e7TW7vpMrFkYCzzmLzXkSIsrBchwHCzd1+\
jgMV948enRrBw+9aWDdCiHcU1bgvc+\
d1ADzo0x7N9KYVzgsJAVC1jv8bCAMBwOC9hvnkSH4KgD9pDnxjh9I8bx+kJkETvcQ+\
ITNLD/ALWFWnUbWVtkjcBxR2iU2qJOxMTg3Yj0tbA8fIdVzQyHfk8HDgR5MO/\
RiCwcrLtdwL5zo0z0luIFPEr03LrZgBbW1beH017Pr2+1Jk330kJi4UNlX4OhDoNx+\
iyMXXY6eKhvbGFJ3fW39j195nkReLJwG0sjzdTBcyj4KGRqoVG7K3DLzI13+\
DTlw0kJHBfig5MPJw+1FhwGmpl3zpi2XK0w4KOTm7DrDWgsPpO5Aa+C89fNS7TsHjQ+\
7PnUHvTr/RKOQ5gXDYxYqRX9mJINRJN2GLGmz0nEeEx48SfRQjuHGPrildHn0KFCy3+\
pbaR2DTAfctVEj16qb2EoSczuRE7XrDtSI9bO2A9/\
b3fUgU9TP8lMJhybwuPUO3w++Sno0EWOF2/BypU6BYmxiPbtZXZfLi1J4xP7Gzl1/\
Mw8tp9/1QLfWvDBS0Xda50OGShJ7LmnDsoptRmZwc3up7YrJczsQynZwd5rKdmp6pnjw+\
FP3TO0EN9pYO8kPvn+1JyvuzOvKjfhsf4dPic/\
st7oaEf2lOIQQ3Xi2KzIoXx2Hap5SmTlY1arHSYXuNlu7NVuEZlIQykrZ47IEikOOH2zzy\
gwo9D9T5UM1LHbHYIO5rOz2jhrtcxMW4ezPRyu5z1af0zQ99T9cQuerr0/\
LnXSTrv9cbM9L0wvgS1BZ8+\
g8P941FNemEcd7cajIlJy5qgEE49qYVaD4SZxODB1NjyzW6mBfUmleshIjokza8vRNABTv\
QUNOJqs23eIt+BVlFqbxFDSPEoqdio5gwXWllKpS1sYCUNvO3J+\
H34X860UzVqICmNtgadYRM5zhOdVQFnvmGQIlXteu1TbT1g0Wy4mc0ITq8+\
ZWJt9v3ShqBpO9XC4OIw92mxev74D3wqtdcf5MU6iWYQ9K6ZDDUT+\
D6dDdSPTdh0i4tVN9P7NjlCIxC9t5HUrtceG2KuBNbbuDvEdUsr87nb2R3LKtDD7W3N02+\
G0NIeTnapjf5xdp13ZiViNBs6HKdVvUpYAn/lSRqo7JnlQZW+l1FoqzAtzZBpnBDrcWp/\
X2zOG6i6S++fMI9Z2eXHmdeY14wQOZ2c6iVTXpWS7s+D2OWyGDma/utmxk8MZrU5xroTD+\
bNOceaAw1m1KI+/\
1JxlKflN73kdzg3WBHXq0qHuDH8DGONMZOfgJs5Pdmpu2qOh6XtsIpYLZ83B+dSKQZa+\
TszhHGm416u8Nv/1K/ZdJE/bsXp/fPKSwxnkDpcRCGKBc8WdBLdjZeByOKdbu/\
Di7G2Hkx0oJbgPOAluuG85mq2G2kLkdMPrJx3u3+\
pSJ5IX3Dy6xf0b9q7qxlao9WdQfnEfQHmOVjqU2/\
NqX8sKcWNy20VoVmpe2i4P6fsDYvTF1hOP/\
bE8lJISJq8ytr8QN5wFDMegK++jev2pRy6SPwv3LemrDFwk6xPjvTYQU8CZnHaYU2vhJK6\
6O9LK7PhzkexXpTj0e9vz6nD+rIj1PFjEEGfKKtX7R+\
wkXoVUyj4UOu2n2N5fyAuc0SrirTBzwnaf1j47yXg2xA7mYvsOGhj740xZkRloxkYRkXSc\
j6qynBE6hLM+nWLuK4axOYqCc/Tg2ZLK+\
WZa8CfN8ktFBKp6mfbcNA2suKu8m2aH3ZfmlZJMJpzKuwPPu4PZMTl6Gqa6OSpcaL+\
ew3mCWqDv0db/wvxDbfTSzUofzs5EnVxTZj+ki+RWMufyPdw+bGIcEMlVVJObLpJ3pwW+\
VR2TRbLm2OnQwr6v9Piki2TYqSqThXZfukh2GzuPmfjD62PQLpIJp2ovgleMRjK14N4/\
vZ8vdTKbg7lIOHjretCFdt7g7CuH83+b3XXIZF852r7xquetOCrzSSnbV8oaMfOxB9vbN+\
HzH6yZOUrNrjucoYSK8TewMw/\
n9jjcyk5IdZxzgfLwakpULMcq9uWIPuBgIzWODtgvKbUekditZJ7Ld5HcCK9fnxqTHUyHU\
usfzCsdIif6M6tvmh0BRk7TbyHKah+3mH3EBuaE28085SKn0+PVpXVXpRH+\
GT7pd4d8VHVP8sQ+TANeW6kdJ/h02XbziO2wrgL36GXW0LewysDBs/4iXjw8Ua+\
cn4qtBpR1aoyDNRavjyqFG9YszGN4F55TFrEE8BytHaQPnipVWbvheUwRDYBnFkV0CJ4BV\
HlMBs/q4eqFHILn5ERGgPDUmB6eocLuJPzth7d3b+904Pfb6++f7y8bL2fXFf1+q/\
kdH0QfU2BTmKA/m63Fa9irHor1GyTM5/\
kuYsDziaL1PpxPJFuOFh2onld82RoNMrRlJ7OzJYl7GmVe+\
N0mtVxarS7bUrlkeIy4si88v+Xe4u02YYtxq09j6vtm4mCktCxvDe/\
pvtUmJa6JHV2acrp3l7cv0NKSoUA51XCmJwe2XQn8/\
EIZ3UdVjnNvU082cdinyiRaP1GD5zrLSSrTjEPNxO1zDDRjY+\
4Zgx5gKknGLdEwSkBqA4kd03M671QxZbSp55nokVzYtKp7Qp7YMm96ckx/\
jsilmYNNCxGVTKCHqR3HVj2UiOxDh3yl5pDm5cUqm7tqHacPOdqT+ulmp1mMLmYUM+\
F7NsbY0K4NCm+\
PJAhVg2OYeJsGo7fbCxkmxvxltmDurI8GxS//qfV5hjU2vrtlnfKEqel6jvVla9Sj6slA/\
rnp2Ha5NpQr2ts+qxlA0EffS3Vbn8N276M+fIPhKEGQLUMFe+\
MqNHiGQo6gZxjuTVuepCEIm8Q6057ndAYUNvlGLx+\
0R22ws3DXl81fGvFIVaYcs5JoUHJUJkcNzJTYbNuKNodIYaFNnBeziGOOqjKSw9SSQcAgK\
VHKvR1o4uh0RCpgKBrI/3sI4jA/\
QMPaymilFAQkizI6sk3xNUO4P5wFJSsN9y27t7N8cYZmA45Hke3jFHysuAeP0EQ0ovNDa/\
HK4tDeTVHH6ZmG6KZmpdH7SdMZVjJuYzL0JyYW3OwrJTBpYhOT0dpiExWIfLEx6+\
cx7K3JtsWFPIN2qOkUoDcd8xe5Cs8b3cTwseczsHJbsPKFjRKzjfoZNHepNg3tRMYnMoFo\
hlsk8XT8lt5QHtUMhi4f+zhjKmHN2eACkkVgy6jak7FaTh/\
GG6Esp8v6SSKYVtSJmpl1Q6QAZYhhfseRpnvR+\
hPW2e1oSko0ajdcOUgZG3gmBDninabflOJtNsYZakYk3qHQ6yVqPVMrfVsd8EmxL7iX6TV\
QywyJQcqYt/HIETDkSSiAs4pU5/\
H6Cup88no8cyp4ajW3BXSnJ4kWNlU54qK9STzGYqeJT455eK7pyVNcacDyYrPbnqLnoO28\
jShGrrz1fmaxSm3OqDVZmGQBmm+\
xkpQmg173W3s2bEQ2TIw5ePs8ddlLmoRn7HmBtOM33YAPdw/\
uGJpjrHkKWVFdWXPdOVqyJVj5diutMep5Ri2J+mPFPLt+K/\
5pnxjf3VT7iA1M2beWaNlSlOap7uzTTqJWpVGFI8SGmU01ErwxZmQvn2Q7GYps1ieLf77d\
41Aw1J0hq49dFNc5gcUnibXJZXvVRbbs7SSanEnkmZS8ZTwKz86bzp+/u/fXLcIRrfV/\
RR/rPg=="

Update

According to the customer service the bug should be resolved in the current release of version 12.0 (and was already resolved for 12.1):

Hello Gert Vercleyen,

In September 2019 you reported an issue with Mathematica wherein SmithDecomposition problem with rank deficient case. We believe that the issue has been resolved in the current 12.0 release of Mathematica.

Thank you for your report and we look forward to a continued, productive relationship with you.

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ Hm. With the matrix that you posted in compressed form, I cannot reproduce the problem (running version 12). $\endgroup$
    – Henrik Schumacher
    Commented Jun 5, 2019 at 21:06
  • 1
    $\begingroup$ Works here too... maybe try restarting kernel to rule out lingering definitions. $\endgroup$
    – bill s
    Commented Jun 5, 2019 at 23:21

1 Answer 1

Reset to default
6
$\begingroup$

You have run into a bug in SmithDecomposition, but the example you've posted doesn't trigger it.

First, the bug. It only occurs in very rare corner cases, so much so that if you sit there and feed random matrices to SmithDecomposition you'll never see it fail. This presumably explains how it got through testing. The effect of the bug is that SmithDecomposition prints the error messages you described, then hangs indefinitely in a busy loop.

Your example matrix however works fine for me, so I'm guessing you pasted the wrong one or it got corrupted somehow. The simplest matrix I know of which triggers the bug is the following, which is a 12x12 integer matrix:

"1:eJxTTMoPSuNhYGAoZgESPpnFJai8zP9AkAnkMGAl8MqimQQWFwYRIPHMDyCtAiAWB1y\
X\\
NIj4AjcTrO4HkIvNJBuQshIQMRXEZQGxguFaz/\\
wDspJArPsgMSEQkY3DTUogghdE3AYpYwOxPsFN0gIRT0HcjyAC7DpcbtoH1JZ5E+\\
izzFygRGbObyDL5wfMpG9AAzLbgNoyk74BWdY/\\
gcRsZuxuEgIRCiBiI8jCQhArBx5OfiDiCUjiJYiYCOIewXQTlePuMEjDLhABthocMG5wXa\
\\
ogwgJEXIOr24UjnE4AnZ358DmQMAaFzgs+\\
oNhFCZhJFsC4yIwEhkBmPig8U0DqeJ9hN0lwMTdDZubhlf8zOw6H/8+\\
MjxFhyPz4UpkBSXbtpk6GzAXLu4AxI332f6armB5W38U6rfif+\\
eBSGNA4q7MMmVLpCf8ze7V+/IPI7k32/p/5oTzlf2Ylrz5j5q11X/9nxi6xwWrSk/LpQN+\
\\
ZRTJkzmicxpC588nv/5nP5B2gblJaFcuQ+Vpk8v/MKzzVQIOt3v/PbPk19z8AJZFAWg=="

I reported this bug almost a year ago, and the developers responded fairly quickly. However, it doesn't look like it will actually be fixed until next release. If you really need this function to work properly, then this answer lists some possible replacements. Alternatively, here's a hacky workaround:

SmithDecomposition[{{0}}];
Unprotect[SmithDecomposition];
Clear[SmithDecomposition];
SmithDecomposition[mat_] := 
 Module[{uu, dd, vv, diags, col = 0, tmpu, tmpv, ord},
  {uu, dd, vv} = System`SmithDecompositionDump`diagonalize[mat];
  diags = Select[Diagonal[dd], # != 0 &];
  While[col + 1 < Length[diags],
   col++;
   If[System`SmithDecompositionDump`divides[diags[[col]], 
     GCD[Sequence @@ Drop[diags, col]]], Continue[]];
   vv = Transpose[vv];
   Do[
    dd[[j, col]] = diags[[j]];
    vv[[col]] += vv[[j]]
    , {j, col + 1, Length[diags]}
    ];
   vv = Transpose[vv];
   {tmpu, dd, tmpv} = System`SmithDecompositionDump`diagonalize[dd];
   uu = tmpu.uu;
   vv = vv.tmpv;
   diags = Select[Diagonal[dd], # != 0 &];
   ];
  ord = OrderingBy[Diagonal[dd], # /. (0 -> Infinity) &];
  If[! OrderedQ[ord], {uu, dd, vv} = {uu[[ord]], DiagonalMatrix[Diagonal[dd][[ord]]], vv[[;; , ord]]}];
  {uu, dd, vv}
  ]
Protect[SmithDecomposition];

This is the internal code, extracted with this method, and empirically and badly tweaked to make it work.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I have reported the bug as well and there's a big chance it will be fixed by next release. The matrix I posted triggered it on my computer at that time but meanwhile I have found many other examples that trigger it. It seems that most of the cases for which I needed it (e.g. finding cocycles of groups) consisted matrices of which the decomposition is very hard compared to the general case... $\endgroup$
    – Gert
    Commented Feb 22, 2020 at 15:29
  • $\begingroup$ I'm not sure exactly when it was fixed, but it's working now in 12.1. $\endgroup$
    – PBS
    Commented Dec 31, 2020 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.