Fortran 90 - Prodotto tra matrice e vettore

G.D.5
Buona sera a tutti. Sto lavorando in Fortran 90 (per l'università) e sto cercando di scrivere un programma che utilizzi una subroutine per il calcolo del prodotto tra una matrice a m righe e n colonne per un vettore a n righe.
Il problema è che il compilatore non mi da errori ma poi quando mando in esecuzione il programma mi fornisce dei risulatati che non mi tornano.
Se invece scrivo il programma per intero senza subroutine, il risultato restituito è corretto.

Il programma senza subroutine è questo:

program matrice_vettore
real::a(100,100),b(100),c(100)
integer::m,n,i,j
print*,"Inserire il numero delle righe e delle colonne della matrice"
read*,m,n
print*,"Inserire gli elementi della matrice"
do i=1,m
   do j=1,n
      read*,a(i,j)
   enddo
enddo
print*,"Inserire gli elementi del vettore"
do i=1,n
   read*,b(i)
enddo
do i=1,m
   c(i)=0
   do j=1,n
      c(i)=c(i)+a(i,j)*b(j)
   enddo
enddo
print*,"Il vettore risultato è"
do i=1,m
print*,c(i)
enddo
end


e funziona.

Con la subroutine invece mi viene così:

program matrice_vettore
real::a(100,100),b(100),c(100)
integer::m,n,i,j
print*,"Inserire il numero delle righe e delle colonne della matrice"
read*,m,n
print*,"Inserire gli elementi della matrice"
do i=1,m
   do j=1,n
      read*,a(i,j)
   enddo
enddo
print*,"Inserire gli elementi del vettore"
do i=1,n
   read*,b(i)
enddo
call mat_vet(m,n,a,b,c)
print*,"Il vettore risultato è"
do i=1,m
print*,c(i)
enddo
end


subroutine mat_vet(m,n,a,b,c)
real::a(m,n),b(n),c(m)
integer::m,n,i,j
do i=1,m
   c(i)=0
   do j=1,n
      c(i)=c(i)+a(i,j)*b(j)
   enddo
enddo
end


e non funzione: se per esempio uso la matrice a 2 righe e 3 colonne con 1,2,3 come prima riga e 4,5,6 come seconda riga e poi l vettore a 3 righe con 1,2,3 in prima, seconda e terza riga rispettivamente, anziché darmi 14,32 come risultato mi da 1,4.

Qualcuno saprebbe dirmi dove sbaglio nella subroutine?

Risposte
GIBI1
Prova a modificare la subroutine:

subroutine mat_vet(m,n,a,b,c)
real::a(1,1),b(1),c(1)
integer::m,n,i,j
do i=1,m
c(i)=0
do j=1,n
c(i)=c(i)+a(i,j)*b(j)
enddo
enddo
end

Umby2
a me sembra tutto corretto....

L'unica cosa un po' strana è la dimensione della matrice e dei vettori:

nel programma "main" li hai definiti a 100 (in modo fisso)
nel "sub" sono definiti di grandezza m e n

G.D.5
Sì, ma la subroutine deve essere realizzata in maniera tale da potere essere adeoperata per qualunque programma principale che preveda l'utilizzo del prodotto matrice-vettore. Se fisso le dimensioni della matrice e/o del vettore, limito l'utilizzo della subroutine.

Umby2
E' evidente che devi dimensionarla in modo opportuno (semmai per un valore che tu ritieni possa essere il massimo). Potresti semmai anche controllare questi valori (m, n) non superano il tuo "limite".

In generale, tra il main ed il sub le variabili di "link" dovrebbero essere SEMPRE dimensionate allo stesso modo. (Almeno cosi' è in Cobol, ma penso che il concetto possa estendersi anche ad altri linguaggi).

Umby2
la matrice b dovrebbe essere (1:n)

penso, un semplice errore di trascrizione

G.D.5
Il mio problema è che questa soluzione non è didatticamente accettabile. L'unico strumento che posso usare per manipolare le matrici è il comando ld.
Però non capisco se in questo caso possa essere di qualche utilità.

GIBI1
Scusa Sergio, hai provato a sostituire nella subroutine

real::a(m,n),b(n),c(m)

con

real::a(1,1),b(1),c(1)


lasciando tutto il resto invariato.

G.D.5
Innanzitutto vi ringrazio per l'enorme aiuto che mi state dando.
Non vi nascondo che tutto quello che mi state dicendo non lo sapevo, quindi mi state insegnando molto.

Permane però un problema: stamane ho consultato telefonicamente un collega e mi ha confermato che la soluzione "didattica" al mio pasticcio sta nell'uso del comando lda (leading dimension matrice a), il problema è che questo mio collega ha spudoratamente scopiazzato il codice (non so da chi e da dove) e io non ho voglia di fare un copia-incolla da un file eventualmente inviatomi.
Ho provato a dichiarare lda=m nel main program dopo l'istruzione read*,m,n, quindi ho aggiunto lda ai parametri di passaggio tra main program e subroutine, ma non gira.

program matrice_vettore
real::a(100,100),b(100),c(100)
integer::m,n,i,j
print*,"Inserire il numero delle righe e delle colonne della matrice"
read*,m,n
lda=m
print*,"Inserire gli elementi della matrice"
do i=1,m
   do j=1,n
      read*,a(i,j)
   enddo
enddo
print*,"Inserire gli elementi del vettore"
do i=1,n
   read*,b(i)
enddo
call mat_vet(m,n,a,b,c,lda)
print*,"Il vettore risultato è"
do i=1,m
print*,c(i)
enddo
end


subroutine mat_vet(m,n,a,b,c,lda)
real::a(lda,n),b(n),c(m)
integer::m,n,i,j
do i=1,m
   c(i)=0
   do j=1,n
      c(i)=c(i)+a(i,j)*b(j)
   enddo
enddo
end


Qualcuno ha idee?

GIBI1
Per quanto ne so, ed è molto poco:

con i vettori non c'è alcun problema, qualunque cosa fai va sempre bene: b(n),b(1), b(100) ecc.

con le matrici bidimensionali accorre la prima dimensione corretta: a(100,n), a(k,1); con k=100.

Magari ci sono altri artifici (ad es. io utilizzo il “common”), che solo un informatico ti può dare.

lorven
L'uso di array ad allocazione dinamica risolve il problema.
program matrice_vettore
real, allocatable :: a(:,:), b(:), c(:)
integer :: m,n,i,j
print*,"inserire il numero delle righe e delle colonne della matrice"
read*,m,n
allocate (a(1:m,1:n),b(1:n),c(1:m))
print*,"inserire gli elementi della matrice"
do i=1,m
   do j=1,n
      read*,a(i,j)
   enddo
enddo

print*,"inserire gli elementi del vettore"
do i=1,n
   read*,b(i)
enddo
call mat_vet(m,n,a,b,c)
print*,"il vettore risultato è"
do i=1,m
   print*,c(i)
enddo
deallocate(a,b,c)
read*
end

subroutine mat_vet(m,n,a,b,c)
   integer :: m, n, i, j
   real :: a(1:m,1:n), b(1:n), c(1:m)
   
do i=1,m
   c(i)=0
   do j=1,n
      c(i)=c(i)+a(i,j)*b(j)
   enddo
enddo
 
end

:D

G.D.5
Eh, ho capito che l'allocazione dinamica risolve il problema, purtroppo non è argomento di corso, quindi sono felice di averla imparata così quando lavoro per i fatti miei uso quella, ma all'uni non posso.

G.D.5
"Sergio":

a) il tuo prof ti "costringe" ad usare tecniche di programmazione ormai superate col Fortran 90, tipiche del Fortran 77;


Purtroppo sì.
Sergio sei un grande: sei riuscito a spiegare in quattro righe quello che io riesco a speigare solo in venti. E' esattamente questo il problema.
Lunedì prendo anche il codice di cui mi parlava stamane il mio collega, verifico che funzioni, lo posto sul forum e, se ti e/o vi va, lo commentiamo.
Frattanto, buon fine settimana.

G.D.5
"Sergio":

E.... la birra?


Intendi l'emotico? Non l'ho messa perché non bevo: non reggo alcun tipo di bevanda alcolica :partyman:
Però ti faccio compagnia con la Coca-Cola :drinkers:

G.D.5
Intendi il cioccolatino al liquore riepieno di grappa?

G.D.5
"Sergio":

No no: bicchierino di grappa. Aiuta ;-)


Ah, ecco :-D

G.D.5
Rieccomi qui. Ieri ho avuto il codice dal mio collega, ma purtroppo manco funzionava. Ho preso però spunto da quel codice ho apportato alcune modifiche e adesso funziona. Onestamente non ho capito nemmeno io come faccia a funzionare.

Quello che segue è il codice del main program e quello della subroutine.

program prodotto_matrice_vettore
real::a(100,100),b(100),c(100)
integer::m,n,i,j
print*,"Programma per il calcolo del prodotto di una matrice per un vettore"
print*,"Inserire le dimensioni della matrice (righe,colonne)"
read*,m,n
print*,"Inserire gli elementi della matrice"
do i=1,m
   do j=1,n
      read*,a(i,j)
   enddo
enddo
print*,"Inserire gli elementi del vettore di lunghezza pari al numero di colonne della matrice"
do i=1,n
   read*,b(i)
enddo
call prod_mat_vet(m,n,a,b,c,100)
print*,"Il vettore risultato è"
do i=1,m
   print*,c(i)
enddo
end


subroutine prod_mat_vet(m,n,a,b,c,lda)
real::a(lda,n),b(n),c(m)
integer::m,n,i,j
do i=1,m
c(i)=0
do j=1,n
c(i)=c(i)+a(i,j)*b(j)
enddo
enddo
return
end


La cosa strana è che passando $m$ in luogo di $100$ alla subroutine dal main program (mi riferisco alla stringa "subroutine prod_mat_vet(m,n,a,b,c,100)" il programma torna a non funzionare, mentre per quanto Sergio aveva argomentato nell'ultimo post sulla questione, avrebbe dovuto fungere.
Cosa mi sto perdendo?
Boh... Vado a meditare!

G.D.5
Ah, ecco. Sergio sei un grande.

svuotapista
Scusate, ho un problema molto simile che mi sta coinvolgendo in questi giorni.
Devo eseguire un algoritmo che implementi il prodotto tra matrice A e vettore u generico.
La matrice A è simmetrica ed è memorizzata (in spagnolo è 'almacenada' - sono in interscambio) come skyline (en perfil). Cioè in Fortran è memorizzata come vettore w che riporta tutti gli elementi partendo dalla prima colonna in poi esclusi gli elementi nulli prima di un elemento non nullo. Mi spiego meglio, data ad esempio:
1 0 4 0 0
x 2 0 0 1
x x 3 1 0
x x x 4 1
x x x x 5
la si memorizza con un vettore w = [1,2,4,0,3,1,4,1,0,1,5] e con un altro vettore di "punteros" (vettore posizione) l, dove l è la posizione dell'elemento che sta sulla diagonale principale.
In questo caso ho: l=[1,2,5,7,11].
E' 5 giorni che ci sbatto la testa sopra ma francamente mi mancano due spunti, ovvero come richiamare il generico elemento della matrice corrispondente a quello che devo moltiplicare e come moltiplicarlo.
In pratica tutto.
-.-
Scusate, grazie per chi mi volesse aiutare.

apatriarca
Ciao svuotapista, benvenuto/a nel forum.

Non conosco bene fortran (ne so qualcosa ma non ricordo molto bene la sintassi), per cui ti risponderò in termini matematici. Gli elementi non nulli, sopra la diagonale, della colonna i-esima sono contenuti negli indici da l(i-1) a l(i) (dove suppongo che l(0) = 0 in modo da far funzionare questa formula). Supponiamo quindi di voler ottenere il valore all'indice (i,j) della matrice. Per prima cosa ci assicuriamo che sia \( i \le j \). Se non fosse così è necessario/sufficiente scambiare i due indici. A questo punto calcoliamo \( d = l(i) - (j - i). \) Se \( d > l(i) \) allora w(d) contiene il valore cercato, altrimenti questo è nullo. Ti ricordo che, siccome la matrice è simmetrica, puoi anche pensare di scambiare il ruolo di righe e colonne se fosse più comodo pensare in termini di righe.

Rispondi
Per rispondere a questa discussione devi prima effettuare il login.