Anonymous user
Talk:Multi-dimensional array: Difference between revisions
no edit summary
No edit summary |
|||
Line 65:
:I agree that I should not have mumbled that mathematics uses "row major" since in the Platonic space the issue of storage doesn't arise, but I was trying to explain the difference between the expectations engendered by mathematical notation and layout, and the "formula translated" with an eye to the many who have been disconcerted by the results of WRITE(...) A after a READ(...) A followed by some calculations. And indeed there are computations that may be expressed row-wise or column-wise and there may be mathematical reasons to prefer one over the other or to regard both as equivalent. But when staring at the realisation of one in Fortran, where when the mathematics calls for (i,j) and (j,i) one may find (j,i) and (i,j) perhaps with DO J...; DO I... when DO I...; DO J... were expected, but that after some inspection seems to be consistent and indeed correctly implementing the algorithm, so, why? Aha, someone has discovered "thrashing" and been driven to taking steps... [[User:Dinosaur|Dinosaur]] ([[User talk:Dinosaur|talk]]) 01:19, 26 October 2016 (UTC)
::I was probably too unkind, sorry for this. Patience is not my primary quality. Anyway, let me show you why all this discussion is really pointless, on an example: gaussian pivoting, to solve the equation AX=B, where A is (n,n), X is (n,p) and B is (n,p). A first attempt could be this (I removed row swapping to shorten the code, as it won't change the idea):
::<lang fortran>do i = 1, n - 1
do j = i + 1, n
s = a(j, i) / a(i, i)
do k = i + 1, n
a(j, k) = a(j, k) - s * a(i, k)
end do
do k = 1, p
b(j, k) = b(j, k) - s * b(i, k)
end do
end do
end do
do i = n, 1, -1
do j = i + 1, n
s = a(i, j)
do k = 1, p
b(i, k) = b(i, k) - s * b(j, k)
end do
end do
do k = 1, p
b(i, k) = b(i, k) / a(i, i)
end do
end do</lang>
::It works, but as you can see, all inner loops run on rows: it would be good in C or another "row-order language". It's very easy to write the same thing with column-order inner loops:
::<lang fortran>do i = 1, n - 1
do k = i + 1, n
s = a(i, k) / a(i, i)
a(i, k) = s
do j = i + 1, n
a(j, k) = a(j, k) - s * a(j, i)
end do
end do
end do
do k = 1, p
do i = 1, n
s = b(i, k) / a(i, i); b(i, k) = s
do j = i + 1, n
b(j, k) = b(j, k) - s * a(j, i)
end do
end do
do i = n, 1, -1
s = b(i, k)
do j = i - 1, 1, -1
b(j, k) = b(j, k) - s * a(j, i)
end do
end do
end do</lang>
::There is absolutely nothing awkward in this alternate code. The number of loops is identical. However, with n=p=1000, to invert the matrix A (with B=eye(1000)), it's 10 times faster. Not 10 percent, but 10 times. So, basically, you need to do this correctly. The problem would be identical in C: one of these algorithms should be 10 times faster, but it will be the first. The problem does not lie in itroducing risk or complication in the program: they are equally complicated. There is no confusion either. You don't have to actually know how the matrix is stored to understand the algorithm. You have to know this, however, to understand why one is much slower than the other.<br/>
::I would like to add: when writing programs, I almost never enter a matrix with READ/WRITE, but they are built from data that is read from a file (or passed by a DLL call from a higher level language such as Excel/VBA or R).
::There is no source of confusion: in A(i,j), i is the row index as one would expect, and j the column index. You will be confused only if you think the storage dictates the interpretation. But it does not. You don't have to worry about th storage except in those situations:
::*Mixing languages that do not use the same storage order (it happens if I call Fortran code from SAS/IML for instance).
::*Optimization: cache misses cost '''much'''. Always manage to get your inner loops to work on contiguous data, whenever it's possible.
::*Disk storage, although usually it's not a problem: either you store quick and dirty memory dumps, and you don't care how it's stored, or you store data in a widespread format for interchange, and either you would use a library, or you will have to take care of many details beyond storage order, anyway.<br/>
::[[User:Arbautjc|Arbautjc]] ([[User talk:Arbautjc|talk]]) 16:44, 26 October 2016 (UTC)
| |||