Es posible declarar un arreglo local en un subprograma en Fortran, pero es muy poco usado. Por lo general, todos los arreglos son declarados (y dimensionados) en el programa principal y entonces son pasados a los subprogramas conforme se van necesitando.
y := alpha*x + y
donde alpha es un escalar y x e y son vectores. Se muestra una subrutina
simple para esto:
subroutine saxpy (n, alpha, x, y)
integer n
real alpha, x(*), y(*)
c
c Saxpy: Calcula y := alpha*x + y,
c donde x e y son vectores de longitud n (al menos).
c
c Variables Locales
integer i
c
do 10 i = 1, n
y(i) = alpha*x(i) + y(i)
10 continue
c
return
end
La única característica nueva es el uso del asterisco en las declaraciones de x(*) e y(*). Con la notación anterior se indica que x e y son arreglos de longitud arbitraria. La ventaja es que se puede usar la misma subrutina para vectores de cualquier longitud. Se debe recordar que como Fortran esta basado en llamadas por referencia, no se requiere espacio adicional, ya que la subrutina trabaja directamente en el arreglo de elementos de la rutina o programa que la llamo. Es la responsabilidad del programador asegurarse que los vectores x e y han sido realmente declarados para tener longitud n o más en algún lugar del programa. Un error común en Fortran 77 sucede cuando se intenta acceder arreglos del elemento fuera de los límites.
Se pudieron también haber declarado los arreglos como:
real x(n), y(n)
Muchos programadores prefieren usar la notación asterisco para recalcar que la
"longitud verdadera del arreglo" es desconocida.
Algunos programas viejos de Fortran 77 podrían declarar arreglos de longitud variable de la siguiente forma:
real x(1), y(1)
La sintaxis anterior es válida, aunque los arreglos sean mayores que uno, pero
es un estilo pobre de programación y no se recomienda hacerlo.
subroutine matvec (m, n, A, lda, x, y)
integer m, n, lda
real x(*), y(*), A(lda,*)
c
c Calcular y = A*x, donde A es m por n y guardado en un arreglo
c con dimensión principal lda.
c
c Variables locales
integer i, j
c Inicializar y
do 10 i = 1, m
y(i) = 0.0
10 continue
c Producto Matriz-vector por saxpy en las columnas de A.
c Observar que la longitud de cada columna de A es m, y no n
do 20 j = 1, n
call saxpy (m, x(j), A(1,j), y)
20 continue
return
end
Hay varias cosas importantes de comentar. Primero, se debe observar que a
pesar de que se pretendió escribir el código tan general como fuera posible
para permitir dimensiones arbitrarias de m y n, se necesita
todavía especificar la dimensión principal de la matriz A. La declaración de
longitud variable (*) puede ser solamente usado para la última dimensión de un arreglo. La razón de lo anterior es la forma como Fortran 77 guarda un arreglo multidimensional.
Cuando se calcula y=A*x, se necesita acceder las columnas de A. La j-ésima columna de A es A(1:m,j). Sin embargo, en Fortran 77 no se puede usar una sintaxis de subarreglo (pero se puede hacer en Fortran 90). En vez de eso se da un apuntador al primer elemento en la columna, el cual es A(1,j) (no es realmente un apuntador, pero puede ser útil pensarlo como si fuera). Se sabe que las siguientes localidades de memoria contendrán los siguientes elementos del arreglo en esta columna. La subrutina saxpy tratará a A(1,j) como el primer elemento de un vector, y no sabe nada de que este vector es una columna de una matriz.
Finalmente, se debe observar que se ha tomado por convención que las matrices tienen m renglones y n columnas. El índice i es usado como un índice de renglones (de 1 a m), mientras el índice j es usado como un índice de columnas (de 1 a n). Muchos programas de Fortran que manejan álgebra lineal usan esta notación, lo cual facilita mucho la lectura del código.
Se muestra un ejemplo muy simple. Otra operación básica con vectores es el escalamiento, por ejemplo, multiplicar cada elemento de un vector por la misma constante. La subrutina es la siguiente:
subroutine escala(n, alpha, x)
integer n
real alpha, x(*)
c
c Variables Locales
integer i
do 10 i = 1, n
x(i) = alpha * x(i)
10 continue
return
end
Supongamos que ahora se tiene una matriz de m por n que se
quiere escalar. En vez de escribir una nueva subrutina para lo anterior, se
puede hacer tratando la matriz como un vector y usar la subrutina escala.
Una versión sencilla se da a continuación:
integer m, n
parameter (m=10, n=20)
real alpha, A(m,n)
c Algunas sentencias definen A ...
c Ahora se quiere escalar A
call escala(m*n, alpha, A)
Observar que este ejemplo trabaja porque se asume que la dimensión declarada
de x es igual a la dimensión actual de la matriz guardada en A. Esto no sucede
siempre. Por lo general la dimensión principal lda es diferente de la
dimensión actual de m, y se debe tener mucho cuidado para hacerlo
correctamente. Se muestra un subroutina más robusta para escalar una matriz
que usa la subrutina escala
subroutine mescala(m, n, alpha, A, lda)
integer m, n, lda
real alpha, A(lda,*)
c
c Variables Locales
integer j
do 10 j = 1, n
call escala(m, alpha, A(1,j) )
10 continue
return
end
Esta versión trabaja aún cuando m no sea igual a lda ya que
se escala una columna cada vez y solamente se procesan los m primeros
elementos de cada columna (dejando los otros sin modificar).
integer nmax
parameter (nmax=40)
real A(nmax, nmax)
Declarar apropiadamente los vectores x e y e inicializarlos.
m=10, n=20, A(i,j) = i+j-2 para
1<=i<=m y 1<=j<=n, x(j) = 1 para 1<=j<=n.
Calcular y = A*x llamando la subrutina matvec dada previamente.
Mostrar el resultado de y.