Ordine per cicli nidificati e trasposizione matrici

Credevate che fosse finita la telenovela del Fortran, invece no!
Di ritorno dalle ferie si riprende con due puntate: * 2 puntate 2 *, ma brevissime. Questa è la prima.

La questione è un classico ma siccome salta fuori continuamente lo ripeto ancora una volta. L’ultima! (A meno che…).

Supponiamo di avere una matrice bidimensionale, per esempio punti X, Y del piano cartesiano, come conviene gestirli?
Il Fortran, a differenza degli altri linguaggi (C/C++, Basic, Python, tutti) ordina i dati per colonna. Fintanto che i dati sono pochi la differenza tra la versione corretta e quella usuale è oggi, con le macchine attuali, trascurabile. Diventa però sensibile e poi non è bello.

Ecco un esempio, minimo, didascalico. Talmente didascalico che lo riciclerò tutte le volte che la questione salta fuori. (Capita sovente dopo le ferie).

Le due versioni seguenti producono lo stesso risultato: rempiono la parte sopra la diagonale principale (compresa) con la somma degli indici di riga e colonna. Con una differenza: la versione per colonne è 3.2 volte più veloce!

program righe
implicit none

integer, parameter :: nmax = 10000
integer :: r, c
integer :: mat(nmax, nmax)

do r = 1, nmax
    do c = 1, nmax
        if (r >= c) then
            mat(r, c) = r * c
        else
            mat(r, c) = 0
        end if
    end do
end do

! per controllo
print *, mat(nmax, 1), mat(nmax, nmax), mat(nmax, 1)

end

program colonne
implicit none

integer, parameter :: nmax = 10000
integer :: r, c
integer :: mat(nmax, nmax)

do c = 1, nmax
    do r = 1, nmax
        if (r >= c) then
            mat(r, c) = r * c
        else
            mat(r, c) = 0
        end if
    end do
end do

! per controllo
print *, mat(nmax, 1), mat(nmax, nmax), mat(1, nmax)

end

OK, si tratta di secondi ma voglio solo vedere quella giusta. Avvisati nèh!

Però metti che io la mia matrice la voglio usare in un programma scritto, per esempio, in C++; devo leggere i dati in modo innaturale?
No! il Fortran ha pensato anche a te e a questo tuo caso: esiste la funzione transpose(). Visualizzarne il funzionamento con una matrice così grossa sarebbe problematico, diciamo che la riduciamo a 10 * 10, ecco:

program inverti
implicit none

integer, parameter :: nmax = 10
integer :: r, c
integer :: mat(nmax, nmax), inv(nmax, nmax)

do c = 1, nmax
    do r = 1, nmax
        if (r >= c) then
            mat(r, c) = r * c
        else
            mat(r, c) = 0
        end if
    end do
end do

inv = transpose(mat)

print '(10i4)' , mat
print *
print '(10i4)' , inv

end