program gauss_main2
integer n, np, m, mp, i
parameter (np=10, mp=20)
real a(np,np), b(np,mp)
character*40 filename1, filename2write (*,*)'Please tell me the dimension n of the matrix A(n,n):'
read (*,*) n
if (n.gt.np) stop 'n can not be larger than np'
write (*,*) 'Please give the size m of the augmented b(n,m):'
read (*,*) m
if (m.gt.mp) stop 'm can not be larger than mp'write (*,*) 'Please give the file name that contains A :'
read (*,*) filename1
write (*,*) 'Please give the file name that contains b :'
read (*,*) filename2open (unit=20, file=filename1)
do i=1,n
read (20,*) ( a(i,j), j=1,n )
enddo
close(20)open (unit=20, file=filename2)
do i=1,n
read (20,*) ( b(i,j), j=1,m )
enddo
close(20)
call gaussj(a,n,np,b,m,mp)
write(*,*) 'The inverse matrix is :'
do i=1,n
write (*,*) ( a(i,j), j=1,n )
enddo
write (*,*)
write (*,*) 'and the solution matrix is :'
do i=1,n
write (*,*) ( b(i,j), j=1,m )
enddoend
program lu_main
integer n, np, m, i
parameter (np=10)
real a(np,np), b(np), d, det
integer indx(np)
character*40 filename1, filename2write (*,*)'Please tell me the dimension n of the matrix A(n,n):'
read (*,*) nwrite (*,*) 'Please give the file name that contains A :'
read (*,*) filename1
open (unit=20, file=filename1)
do i=1,n
read (20,*) ( a(i,j), j=1,n )
enddo
close(20)write (*,*) 'Please give the file name that contains b :'
read (*,*) filename2
open (unit=20, file=filename2)
do i=1,n
read (20,*) b(i)
enddo
close(20)
call ludcmp(a,n,np,indx,d)call lubksb(a,n,np,indx,b)
write (*,*) 'The solution matrix is :'
do i=1,n
write (*,*) b(i)
enddodet = 1.0
do i=1,n
det = det * a(i,i)
end do
det = d * det
write(*,*) 'The determinant is', detend
program polint_main
c In subroutine polint, array declaration is used for xa(n) and ya(n),
c for this we need to reclare larger or equal size xa(n) and ya(n) in
c this main program. As for NMAX and c(NAMX) and d(NAMX) in that routine,
c they are taken care locally in that routine alreay with specific size
c (parameter) given when declaring array, therefore we don't need to deal
c with that.
implicit none
integer max_tabu, ngrids_max, n_tabu, m, i, j, k, choice
integer k1, k2, more_points, ngrids
parameter (max_tabu = 100, ngrids_max = 2000)
real xa(max_tabu), ya(max_tabu), x, y, dy, delta_x
real x_plot(ngrids_max), y_plot(ngrids_max),dy_plot(ngrids_max)
real x_min, x_max
character*40 filename1
c Ask filename and the number of tabulated points in the file.
write(*,*) 'Type in the file name of tabulated values (x,y)'
read (*,*) filename1
write(*,*) 'How many points are there in the tabulated set ?'
read (*,*) n_tabu
if (n_tabu.gt.max_tabu) stop 'Too many tabulated values,
& please recompile the main program with a larger max_tabu.'
c Open the file and read in the tabulated values, and write to screen
c for double check.
open (unit = 20, file = filename1)
do i = 1, n_tabu
read (20,*) xa(i), ya(i)
write(*,*) xa(i), ya(i)
end doclose (20)
c Ask the order of polynomial to be usedfor interpolation, not higher
c than 10.
1000 continue
write(*,*) 'How many points do you want to form polynomial ?'
write(*,*) 'For eample, 2 points mean linear interpolation.'
read (*,*) mif (m.gt.10) then
write(*,*) 'Waring : Neumarial Recipe default max is 10.'
write(*,*) 'Please decress your value or set a bigger NMAX'
write(*,*) 'in the subroutine polint.'
write(*,*) 'Your choice (0 : exit) or (1: input again).'
read (*,*) choice
if (choice.eq.0) stop 'Terminaed by user.'
if (choice.eq.1) goto 1000
goto 1000
endif
c One can choose to type in x one by one or to generte a continuous set
c for graph plotting.
write(*,*) 'What do you want to do ?'
write(*,*)'0:input x, find y ; 1:generate a set of points (x,y)'
read (*,*) choice
c This part get x one by one, using "locate" to find the right interval.
if (choice.eq.0) then
1002 write(*,*) 'Type in your x:'
read(*,*) x
call locate (xa, n_tabu, x, j)
k = min( max(j-(m-1)/2,1) , n_tabu+1-m )
call polint(xa(k), ya(k), m, x, y, dy)
write(*,*) 'Interpolated value is ', y, ', error is ', dy
write(*,*) 'One more point (yes=1;no=0) ?'
read (*,*) more_points
if (more_points.eq.1) goto 1002
endif
c This part generate a set of x, using "hunt" to find the right interval.
c The resulting (x,y) is written out into a file named by user.
if (choice.eq.1) then
write(*,*) 'Type your range x_initial and x-final:'
read (*,*) x_min, x_max
write(*,*) 'How many points do you want to plot ?'
read (*,*) ngrids
if(ngrids.le.1) stop 'Need at least two grids.'delta_x = ( x_max - x_min ) / (ngrids - 1)
j = 1
do i=1, ngrids
x = x_min + (i-1)*delta_x
call hunt(xa, n_tabu, x, j)
k = min( max(j-(m-1)/2,1) , n_tabu+1-m )
x_plot(i) = x
write(*,*) 'x is ',x_plot(i),' j is ',j,' k is',k
call polint(xa(k),ya(k),m,x_plot(i),y_plot(i),dy_plot(i))
end dowrite(*,*) 'Plrease give a name to your output file:'
read (*,*) filename1
open (unit=30, file=filename1)
do i=1,ngrids
write(30,*) x_plot(i), y_plot(i)
end do
close(30)write(*,*) 'Info: Error estimate will be written into file
&error.dat'
open (unit=40, file='error.dat')
do i=1,ngrids
write(40,*) x_plot(i), dy_plot(i)
end do
close(40)end if
end
program spline_main
implicit none
integer max_tabu, ngrids_max, n_tabu, i, choice
integer more_points, ngrids
parameter (max_tabu = 100, ngrids_max = 2000)
real xa(max_tabu), ya(max_tabu), y2a(max_tabu), x, y, delta_x
real x_plot(ngrids_max), y_plot(ngrids_max), yp1, ypn
real x_min, x_max
character*40 filename1
c Ask filename and the number of tabulated points in the file.
write(*,*) 'Type in the file name of tabulated values (x,y)'
read (*,*) filename1
write(*,*) 'How many points are there in the tabulated set ?'
read (*,*) n_tabu
if (n_tabu.gt.max_tabu) stop 'Too many tabulated values,
& please recompile the main program with a larger max_tabu.'
c Open the file and read in the tabulated values, and write to screen
c for double check. open (unit = 20, file = filename1)
do i = 1, n_tabu
read (20,*) xa(i), ya(i)
write(*,*) xa(i), ya(i)
end doclose (20)
c Ask slopes at two ends. (If > 10E30 then natural spline is set.)
write(*,*) 'Please type in slopes at bothend (> 10E30 -> y2a=0.'
read (*,*) yp1, ypn
c Call subroutine spline to solve for y2a array.
call spline(xa, ya, n_tabu, yp1, ypn, y2a)
c One can choose to type in x one by one or to generte a continuous set
c for graph plotting.
write(*,*) 'What do you want to do ?'
write(*,*)'0:input x, find y ; 1:generate a set of points (x,y)'
read (*,*) choice
c This part get x one by one.
if (choice.eq.0) then
1002 write(*,*) 'Type in your x:'
read(*,*) x
call splint(xa, ya, y2a, n_tabu, x, y)
write(*,*) 'Interpolated value is ', y
write(*,*) 'One more point (yes=1;no=0) ?'
read (*,*) more_points
if (more_points.eq.1) goto 1002
endif
c This part generate a set of x.
c The resulting (x,y) is written out into a file named by user.
if (choice.eq.1) then
write(*,*) 'Type your range x_initial and x-final:'
read (*,*) x_min, x_max
write(*,*) 'How many points do you want to plot ?'
read (*,*) ngrids
if (ngrids.le.1) stop 'Need at least two grids.'delta_x = ( x_max - x_min ) / (ngrids - 1)
do i=1, ngrids
x_plot(i) = x_min + (i-1)*delta_x
call splint(xa,ya,y2a,n_tabu,x_plot(i),y_plot(i))
end dowrite(*,*) 'Plrease give a name to your output file:'
read (*,*) filename1
open (unit=30, file=filename1)
do i=1,ngrids
write(30,*) x_plot(i), y_plot(i)
end do
close(30)end if
end
program polint_pgplot
integer ntab_max, nplot_max
parameter (ntab_max = 20, nplot_max = 200)
real x_tab(ntab_max), y_tab(ntab_max), x_min, x_max, delta_x
real x_plot(nplot_max), y_plot(nplot_max), dy_plot(nplot_max)
integer ntab, i, nplot, m, j, pgopen
character*40 filename1write(*,*)'Please tell me the number of tabulated values:'
read(*,*) ntabwrite(*,*) 'Which filename contains the tabulated value ?'
read(*,*) filename1open(unit=10, file=filename1)
do i=1,ntab
read(10,*) x_tab(i), y_tab(i)
end do
close(10)write(*,*) 'Type in min and max of x:'
read(*,*) x_min, x_maxwrite(*,*) 'How many points to draw ?'
read(*,*) nplotwrite(*,*) 'How many point to form polynomial ?'
read(*,*) mdelta_x = (x_max - x_min) / (nplot - 1)
do i=1, nplot
x_plot(i) = x_min + (i-1) * delta_x
call hunt(x_tab,ntab,x_plot(i),j)
k = min( max(j-(m-1)/2,1) , ntab+1-m )
call polint(x_tab(k), y_tab(k), m, x_plot(i), y_plot(i),
& dy_plot(i))
end do
do i=1,nplot,10
write(*,*) 'x_plot and y_plot,', x_plot(i), y_plot(i)
end doif ( pgopen('?') .le. 0 ) stop
call pgenv(2.0, 6.0, 0.0, 40.0, 0, 1)
call pgline(nplot, x_plot, y_plot)
call pgpt(ntab, x_tab, y_tab, 9)
call pgclosend
program lu_inverse_main
implicit none
integer n, np, m, i, j, k
parameter (np=10)
real a(np,np), b(np), d, det
real a_inv(np,np), a_ori(np,np), aainv(np,np)
integer indx(np)
character*40 filename1, filename2write (*,*)'Please tell me the dimension n of the matrix A(n,n):'
read (*,*) nwrite (*,*) 'Please give the file name that contains A :'
read (*,*) filename1open (unit=20, file=filename1)
do i=1,n
read (20,*) ( a(i,j), j=1,n )
enddo
close(20)write (*,*) 'Please give the file name that contains b :'
read (*,*) filename2open (unit=20, file=filename2)
do i=1,n
read (20,*) b(i)
enddo
close(20)
call ludcmp(a,n,np,indx,d)call lubksb(a,n,np,indx,b)
write (*,*) 'The solution matrix is :'
do i=1,n
write (*,*) b(i)
enddodet = 1.0
do i=1,n
det = det * a(i,i)
end do
det = d*det
write(*,*) 'The determinant is', det
c Calculate inverse of the matrix
do j=1,n
do i=1,n
b(i) = 0.0
end do
b(j) = 1.0
call lubksb(a,n,np,indx,b)
do i=1,n
a_inv(i,j) = b(i)
end do
end do
write(*,*)
write(*,*) 'The inverse is'
do i=1,n
write(*,*) ( a_inv(i,j), j=1,n )
end do
c Double check whether A*A_inv=I.
open (unit=20, file=filename1)
do i=1,n
read (20,*) ( a_ori(i,j), j=1,n )
enddo
close(20)
do j=1,n
do i=1,n
aainv(i,j) = 0.0
do k=1,n
aainv(i,j) = aainv(i,j) + a_ori(i,k)*a_inv(k,j)
end do
end do
end do
write(*,*)
write(*,*) 'Check A*A_inv result ...'
do i=1,n
write(*,*) ( aainv(i,j), j=1,n )
end do
end
program spline_pgplot_main
implicit none
integer max_tabu, ngrids_max, n_tabu, i, choice
integer more_points, ngrids, pgopen
parameter (max_tabu = 100, ngrids_max = 2000)
real xa(max_tabu), ya(max_tabu), y2a(max_tabu), x, y, delta_x
real x_plot(ngrids_max), y_plot(ngrids_max), yp1, ypn
real x_min, x_max
character*40 filename1
c Ask filename and the number of tabulated points in the file.
write(*,*) 'Type in the file name of tabulated values (x,y)'
read (*,*) filename1
write(*,*) 'How many points are there in the tabulated set ?'
read (*,*) n_tabu
if (n_tabu.gt.max_tabu) stop 'Too many tabulated values,
& please recompile the main program with a larger max_tabu.'
c Open the file and read in the tabulated values, and write to screen
c for double check.
open (unit = 20, file = filename1)
do i = 1, n_tabu
read (20,*) xa(i), ya(i)
write(*,*) xa(i), ya(i)
end doclose (20)
c Ask slopes at two ends. (If > 10E30 then natural spline is set.)
write(*,*) 'Please type in slopes at bothend (> 10E30 -> y2a=0.'
read (*,*) yp1, ypn
c Call subroutine spline to solve for y2a array.
call spline(xa, ya, n_tabu, yp1, ypn, y2a)
c One can choose to type in x one by one or to generte a continuous set
c for graph plotting.
write(*,*) 'What do you want to do ?'
write(*,*)'0:input x, find y ; 1:generate a set of points (x,y)'
read (*,*) choice
c This part get x one by one.
if (choice.eq.0) then
1002 write(*,*) 'Type in your x:'
read(*,*) x
call splint(xa, ya, y2a, n_tabu, x, y)
write(*,*) 'Interpolated value is ', y
write(*,*) 'One more point (yes=1;no=0) ?'
read (*,*) more_points
if (more_points.eq.1) goto 1002
endif
c This part generate a set of x.
c The resulting (x,y) is written out into a file named by user.
if (choice.eq.1) then
write(*,*) 'Type your range x_initial and x-final:'
read (*,*) x_min, x_max
write(*,*) 'How many points do you want to plot ?'
read (*,*) ngrids
if(ngrids.le.1) stop 'Need at least two grids.'delta_x = ( x_max - x_min ) / (ngrids - 1)
do i=1, ngrids
x_plot(i) = x_min + (i-1)*delta_x
call splint(xa,ya,y2a,n_tabu,x_plot(i),y_plot(i))
end doif ( pgopen('?') .le. 0 ) stop
call pgpap(8.0,0.75)
call pgenv(0.0, 4.0, 0.0, 20.0, 0, 1)
call pgline(ngrids, x_plot, y_plot)
call pgpt(n_tabu, xa, ya, 9)
call pgclos
c write(*,*) 'Plrease give a name to your output file:'
c read (*,*) filename1
c open (unit=30, file=filename1)
c do i=1,ngrids
c write(30,*) x_plot(i), y_plot(i)
c end do
c close(30)
end if
end