Problemas del Curso Física 4

Instructor: Darío Mitnik

Introducción Básica de Fortran - Parte 2

Subrutinas

Hemos visto como utilizar funciones en un programa Fortran. Ahora veremos como utilizar subrutinas. Una subrutina es algo asi como un programa dentro de un programa. Separar un programa en subrutinas sería análogo a separar un libro en capítulos. Si un programa está bien escrito, debería contar con unas pocas líneas, en las cuales se llama a ejecutar subrutinas (las cuales a su vez, llamarían a ejecutar otras subrutinas...). Veamos un ejemplo de un programa no muy bien escrito. El programa se llama sumatorynosub.f y es así: 

	program sumatory

c....... This program calculates the sumatory S=Sum(1/i)

c	**** NO use of subroutines ****

c....... variables declaration
	implicit real*4(a-h,o-z)

	parameter (mxpoints=10000)         !! npoints: points in x

	dimension sumasc(0:mxpoints+1),sumdesc(0:mxpoints+1)

	print*,'Sumatory of 1/i  from 1 to Npoints'
50	print*,' '
	print*,' give the number of points Npoints: (le.0 to exit)'
	read*,npoints
c......  test for quit
	if (npoints.le.0) go to 500
c......  test for dimensions
	if (npoints.gt.mxpoints) then
		print*,' the maximum number of points should be '//
     +                 ' smaller than ',mxpoints
     		go to 50
	endif

c....... initialization
	do 100 i=0,npoints+1
		sumasc(i)  = 0.0
		sumdesc(i) = 0.0
100	continue

c.......  sumatory in ascendent order
	  do 200 i=1,npoints
	       sumasc(i) = sumasc(i-1) + 1.0/i
200       continue

c.......  sumatory in descendent order
	  do 300 i=npoints,1,-1
	       sumdesc(i) = sumdesc(i+1) + 1.0/i
300       continue

c.......  print results
 	  open(unit=10,file='sumatoryns.dat',status='unknown')
	  do 400 i=1,npoints
	    write(10,425) i,sumasc(i),sumdesc(npoints-i+1)
400	  continue
	  close(unit=10)
425	  format(5x,i6,2(5x,f15.6))	
	  print*,'Sum(1/i) for ',npoints,':   ascendent=:',
     +            sumasc(npoints),'  descendent=:',sumdesc(1)	  

	  go to 50
	  
500	stop
	end

Ejercicio:

Qué hace este programa? Los resultados "ascendentes" y "descendentes" son diferentes? Deberían ser iguales? Cómo se puede mejorar los resultados? (Probar que mejoran!)



El programa sumatory podría estar mejor escrito si usamos subrutinas. En este caso, podriamos separar algunas etapas del mismo, como hicimos en sumatorysub.f. Allí separamos los siguientes pasos: 

	program sumatory

c....... This program calculates the sumatory S=Sum(1/i)

c	**** WITH subroutines ****

c....... variables declaration
	implicit real*4(a-h,o-z)

	parameter (mxpoints=10000)         !! npoints: points in x
	common/blocksum/sumasc(0:mxpoints+1),sumdesc(0:mxpoints+1),
     +                  npoints

        ......
	......

c.......  initialization
	call initialization

c.......  sumatory in ascendent order
	call sumascendent

c.......  sumatory in descendent order
	call sumdescendent

c.......  print results
	call printresults

c.......  print totals
          print*,'Sum(1/i) for ',npoints,':   ascendent=:',
     +            totasc,'  descendent=:',totdesc


	  go to 50
	  
500	stop
	end
c
c---------------------------------------------------------------------
c
	subroutine initialization
	implicit real*4(a-h,o-z)

c....... initialization of variables

	parameter (mxpoints=10000)         
	common/blocksum/sumasc(0:mxpoints+1),sumdesc(0:mxpoints+1),
     +                  npoints

	do 100 i=0,npoints+1
		sumasc(i)  = 0.0
		sumdesc(i) = 0.0
100	continue

	return
	end
	
c
c---------------------------------------------------------------------
c	
	subroutine sumascendent
	implicit real*4(a-h,o-z)

c.......  sumatory in ascendent order

	parameter (mxpoints=10000)         
	common/blocksum/sumasc(0:mxpoints+1),sumdesc(0:mxpoints+1),
     +                  npoints
     
c.......  sumatory in ascendent order
	do 200 i=1,npoints
	       sumasc(i) = sumasc(i-1) + 1.0/i
200     continue

	return
	end
c
c---------------------------------------------------------------------
c
	subroutine sumdescendent
	implicit real*4(a-h,o-z)

c.......  sumatory in descendent order

	parameter (mxpoints=10000)         
	common/blocksum/sumasc(0:mxpoints+1),sumdesc(0:mxpoints+1),
     +                  npoints

	do 300 i=npoints,1,-1
	       sumdesc(i) = sumdesc(i+1) + 1.0/i
300     continue

	return
	end
c
c---------------------------------------------------------------------
c
	subroutine printresults
	implicit real*4(a-h,o-z)

c.......  print results on 'sumatorys.dat'

	parameter (mxpoints=10000)         
	common/blocksum/sumasc(0:mxpoints+1),sumdesc(0:mxpoints+1),
     +                  npoints

	
	open(unit=10,file='sumatorys.dat',status='unknown')
	  
	do 400 i=1,npoints
	    write(10,425) i,sumasc(i),sumdesc(npoints-i+1)
400	continue
	close(unit=10)
425	format(5x,i6,2(5x,f15.6))	
c.......  print totals
          totasc = sumasc(npoints)
          totdesc = sumdesc(1)

	
	return
	end
c
c---------------------------------------------------------------------
c


Importante!

Las variables en una subrutina son locales. Esto implica que su valor sólo es conocido por la subrutina y no por el resto del programa. Si se desea pasar una variable de una subrutina a otra (o al programa principal), se la debe poner en common blocks, como hicimos por ejemplo escribiendo: 


        common/blocksum/sumasc(0:mxpoints+1),sumdesc(0:mxpoints+1),
     +                  npoints

En este caso vemos como los arrays sumasc, sumdesc, y el entero npoints son llevados de un lugar del programa al otro.

Ejercicio:

Revisar los resultados del nuevo programa!! (Esto es obvio, cuando se escribe un nuevo programa se debe chequear que haga lo que le pedimos). Graficar los resultados del archivo de salida y compararlos con los del programa anterior. Si hay algún problema, deberíamos ser capaces de solucionarlo...


Ejercicio:

Una versión aún más mejorada del programa, utilizaría un solo loop, en lugar de 4. Qué modificación sugiere para que el programa corra mas rápido? Modifique y compruebe. (Solución: sumatory.f).

Ejercicio:

Agregar al programa una subrutina que reste 1/i del total hallado. Si el programa está bien escrito, el resultado final tiene que ser 0.

Ejercicio:

El programa que utilizamos, también sirve como aproximación para el cálculo de la integral de la función f(x)=1/x. Escribir un programa que genere un archivo X-Y de f(x). Utilizando el programa de gráfica xmgrace, graficar la función y su integral. Graficar el resultado analítico y comparar. Graficar el resultado obtenido por nuestro programa y comparar. Explicar el orígen de las diferencias.