Issue with Out of Bound array for RDF code
by Anupama Sharma from LinuxQuestions.org on (#5HJ4M)
Hello,
I am fairly new to fortran coding and tried to write a RDF (radial distribution function) code for the analysis where during compilation I am getting out of bound array error for line 168 in my code.
Can anyone please check my code and help me to resolve the issue. I will be greatfull for the kindness.
Code:!RADIAL DISTRIBUTION FUNCTION CODE (WORK FOR BOTH NVT & NPT)
program gofr
!===============================================================================
! Gromacs trajectroy reading utility available on github
!===============================================================================
use gmxfort_trajectory
use gmxfort_utils
implicit none
type(Trajectory) :: trj
!=================================================================================
! This section is for the trajectory file input variables
!=================================================================================
integer :: ninput,ntot,inp,nsites,nframes ! input for gofr
integer,parameter :: maxinp=150 ! maximum number of trajectory files as an input
character(len=100) :: infile(maxinp)
integer :: ibegin(maxinp),iend(maxinp),nskip,nread,nused,intv
integer :: i1,j1,u1,n1
integer :: i,j,k,ia,im,ja,jm
integer :: total_atom,natms,nstep
integer :: nmol_cta,nmol_br,nmol_water,natom_cta,natom_br,natom_water
integer :: atom_cbw,ifile
integer :: total_frame
integer :: iatom_cta,iatom_br,iatom_water,imat,jmat
integer :: imol_cta,imol_br,imol_water
integer :: no_atname1
integer :: ibin,nprtot,npr(1000),nbins
integer :: iflg,keytrj,imcon,nmol_diff,new_file
real :: myatom(3)
real(kind=8),allocatable :: chge(:),weight(:)!,yyy(:),zzz(:)
real(kind=8),allocatable :: vxx(:),vyy(:),vzz(:),fxx(:),fyy(:),fzz(:)
real(kind=8):: xx1(1),yy1(1),zz1(1),xxx(1),yyy(1),zzz(1)
real(kind=8),allocatable :: rw(:,:,:),rcta(:,:,:),rbr(:,:,:)
real(kind=8) :: boxl(3,3),dr(3),rmsx,rmsy,rmsz
real(kind=8) :: timstp
real(kind=8),parameter :: au=0.5291770
!================================================================
real(kind=8),parameter :: drg=0.05 ! spacing in g(r)
!================================================================
real(kind=8) :: cell(3,3),rcell(9),det,dt0
!================================================================
real(kind=8) :: rmax,rdr,rsq,r,rr,rrs,cn,g,gavg,fact,cutoff_vol
real(kind=8) :: pi
real(kind=8),allocatable :: rs(:),gs(:)
real :: start,finish
character(len=10) :: nchar
character(len=100) :: outfile1,outfile2,a
character(len=8),allocatable :: atnm_rw(:,:),atnm_rcta(:,:),atnm_rbr(:,:),x(:,:)
character(len=8) :: atname1,atname2
character(len=8),allocatable :: atname(:)
character(len=80) :: title
! write (*,*) "type the trjectory file name"
call cpu_time(start)
!=======================================================================
open(unit=2,file='br_br_gofr.inp',action='read')
read(2,*)outfile1
read(2,*)outfile2
open(unit=3,form='formatted',file=outfile1,&
action='write')
open(unit=4,form='formatted',file=outfile2,&
action='write')
!=======================================================================
!=======================================================================
! Begining of the code for gofr calculation
!=======================================================================
write(*,*)'Enter the no of trajectory file you want to analyze'
read(2,*)ninput ! (maxinp=150)
write(*,*)'Along with this also enter the number of frames you',&
' want to work with.'
read(2,*)nframes ! number of time-steps used for binning of atoms
write(*,*)'Enter the number of different molecules types'
read(2,*)nmol_diff ! (1)
write(*,*)'Enter the no of molecules of all types'
read(2,*)nmol_water ! (884)
write(*,*)'Enter the no of atom in each of the diff molecules'
read(2,*) natom_water ! (3)
if(ninput>maxinp)then
write(*,*)'Too many inputs'
stop
endif
write(*,*)'Pair correlation atom pair X-Y: X'
read(2,*)atname1
write(*,*)'Pair correlation atom pair X-Y: Y'
read(2,*)atname2
write(*,*)'Enter the number of X atom'
read(2,*)no_atname1
!=======================================================================
atom_cbw=0
pi=datan(1.0d0)*4
intv=10
write(*,*)'pi====>',pi
total_atom=nmol_water*natom_water
natms=total_atom
!=======================================================================
ntot=0 ! total number of inputs in trajectory
!=======================================================================
write(*,*)'enter the skip value'
read(2,*)nskip
read(2,*)rmax ! half the box length
!write(4,*)'nskip,rmax=======>',nskip,rmax
if(nskip/=0)nskip=1
!ntot=ntot/nskip
!total_frame=ntot
write(*,*)'Finished reading inputs'
!=======================================================================
!GOFR CALCULATION
!=======================================================================
rdr=1.0d0/drg
nbins=nint(rmax*rdr)+1 !number of bins for calculation has been evaluated
if(real(nbins)*drg>rmax)nbins=nbins-1
do j=1,nbins
npr(j)=0
enddo
nprtot=0
!write(4,*)'rmax =',rmax,'nbins =',nbins
!=======================================================================
! from this line of code below trajectory file name was taken directly
! from the file
!=======================================================================
fread: do inp=1,ninput,1
write(*,*) 'Enter the name of trajectory file',inp
read (2,*) infile(inp)
read (2,*) ibegin(inp),iend(inp)
!=======================================================================
!Number of atoms was going to be read and passed to a variable name nsites
!=======================================================================
call trj%open(infile(inp))
nsites = trj%natoms()
!-----------------------------------------------------------------------
allocate(rw(nmol_water,natom_water,3))!,rcta(nmol_cta,natom_cta,3),rbr(nmol_br,natom_br,3))
allocate(atnm_rw(nmol_water,natom_water))!,atnm_rcta(nmol_cta,natom_cta),atnm_rbr(nmol_br,natom_br))
! allocate(atname(nsites))
!allocate(xx1(nsites),yy1(nsites),zz1(nsites))
!allocate(xxx(nsites),yyy(nsites),zzz(nsites))
!-----------------------------------------------------------------------
n1 = trj%read_next(nframes) !trj%read_next lib allow us to read frame wise data
!=======================================================================
!write(4,*)'number of sites = ',nsites
!write(4,*)'Total atom=',total_atom
! allocate(xxx(nsites),yyy(nsites),zzz(nsites))
! allocate(xx1(nsites),yy1(nsites),zz1(nsites))
!=======================================================================
!here the code will loop over the provided number of frames you want to
!work with as an input trajectory
!=======================================================================
nframe: do i1=1, n1 ! for nframes frames as mentioned
! write (3,*)"frame=========",i1
! write(*,*) i1,"frame============"
cell= trj%box(i1)
call invert (cell,rcell,det)
site: do j1=1,nsites
myatom = trj%x(i1,j1)
xxx = trj%frameArray(i1)%xyz(1,j1)
yyy = trj%frameArray(i1)%xyz(2,j1)
zzz = trj%frameArray(i1)%xyz(3,j1)
write(3,*) j1,"atom-----------------------------"
write(3,*) xxx , yyy, zzz
!write(*,*)'h1',j1,xxx(j1)
xx1=xxx(j1)*rcell(1)+yyy(j1)*rcell(4)+zzz(j1)*rcell(7) !
yy1=xxx(j1)*rcell(2)+yyy(j1)*rcell(5)+zzz(j1)*rcell(8) !
zz1=xxx(j1)*rcell(3)+yyy(j1)*rcell(6)+zzz(j1)*rcell(9)
enddo site
!write(3,*) j1,a
! write(4,*) j1,"atom "," frame",i1
!
!write(4,*) j1," atom ",xx1,yy1,zz1
!100 format (i5,2x,) !
wt: do imol_water=1,nmol_water
do iatom_water=1,natom_water
atom_cbw=atom_cbw+1
rw(imol_water,iatom_water,1)=xx1(atom_cbw)
rw(imol_water,iatom_water,2)=yy1(atom_cbw)
rw(imol_water,iatom_water,3)=zz1(atom_cbw)
atnm_rw(imol_water,iatom_water)=atname(atom_cbw)
!write(*,*) rw
enddo
enddo wt
do imat=1,3
do jmat=1,3
boxl(imat,jmat)=cell(imat,jmat)
enddo
enddo
br1: do im=1,nmol_water-1
do ia=1,natom_water
if(atnm_rw(im,ia)==atname1)then !pair function atm1 (cta)
br2: do jm=(im+1),nmol_water
do ja=1,natom_water
if(atnm_rw(jm,ja)==atname2)then !pair function atm2 (cta)
dr(:)=rw(im,ia,:)-rw(jm,ja,:)
dr(:)=dr(:)-anint(dr(:)) !Box length L=1
!CALCULATE CORRECTED REAL SPACE COORDINATES
rmsx=dr(1)*cell(1,1)+dr(2)*cell(2,1)+dr(3)*cell(3,1)
rmsy=dr(1)*cell(1,2)+dr(2)*cell(2,2)+dr(3)*cell(3,2)
rmsz=dr(1)*cell(1,3)+dr(2)*cell(2,3)+dr(3)*cell(3,3)
dr(1)=rmsx;dr(2)=rmsy;dr(3)=rmsz
rsq=dot_product(dr,dr)
r=(sqrt(rsq))
ibin=int(r*rdr)+1
if(ibin<=nbins)then
npr(ibin)=npr(ibin)+1
nprtot=nprtot+1
endif
endif
enddo
enddo br2
endif
enddo
enddo br1
enddo nframe
!
call trj%close()
enddo fread
100 continue
!-----------------------------------------------------------------------
cutoff_vol=((4.0d0/3.0d0)*pi*(rmax**3))
!-----------------------------------------------------------------------
cn=0.0
fact=(cutoff_vol)/(4.0d0*pi*(drg**3)*real(nprtot))
allocate(rs(nbins),gs(nbins))
do i=1,nbins
r=(real(i)-0.5d0)*drg
rr=((real(i)+0.5d0)**2+(1.d0/12.d0))
g=fact*(real(npr(i))/(rr))
!--------------------------------------------------------
cn=cn+real(npr(i))/real(nused*(no_atname1-1))
!--------------------------------------------------------
rs(i)=(real(i)-0.5d0)*drg
rrs=((real(i)+0.5d0)**2+(1.d0/12.d0))
gs(i)=fact*(real(npr(i))/(rrs))
write(3,2)r,g,cn
enddo
!=========================================================
2 format(2x,f10.5,f10.5,f10.5)
3 format(2x,f10.5,f10.5)
!-----------------------------------------
if(a=='smooth')then
write(4,3)rs(1),gs(1)
write(4,3)rs(2),gs(2)
do i=3,nbins-2
gavg=0.0d0
do j=i-2,i+2
gavg=gavg+gs(j)
enddo
gavg=(gavg/5.0)
write(4,3)rs(i),gavg
enddo
write(4,3)rs(nbins-1),gs(nbins-1)
write(4,3)rs(nbins),gs(nbins)
endif !For 5 point smoothing
!-----------------------------------------------------------------------
! deallocate(xxx,yyy,zzz,xx1,yy1,zz1)
deallocate(rw,atnm_rw,atname)
call cpu_time(finish)
print'("Time= ",f6.3,"seconds")',finish-start,finish,start
end program gofr
!---------------------------------------------------------------
subroutine invert(box,b,d) !subroutine sname (dummy arugments)
!---------------------------------------------------------------
! dl_poly subroutine to invert a 3*3 matrix using cofactors
! author - w. smith april 1992
!---------------------------------------------------------------
real(kind=8):: box(3,3) ! data passed to subroutine that cann't be changed
real(kind=8)::b(9),d
real(kind=8):: r ! r is a temporary variable, not accessible to main program
write(4,*) "in subroutine"
do i=1,3
write(4,*) (box(i,j),j=1,3)
enddo
! calculate adjoint matrix cofactors with simultaneous doing transpose of mtx
b(1)=box(2,2)*box(3,3)-(box(2,3)*box(3,2))
b(2)=box(3,2)*box(1,3)-(box(1,2)*box(3,3))
b(3)=box(1,2)*box(2,3)-(box(2,2)*box(1,3))
b(4)=box(3,1)*box(2,3)-(box(2,1)*box(3,3))
b(5)=box(1,1)*box(3,3)-(box(3,1)*box(1,3))
b(6)=box(2,1)*box(1,3)-(box(1,1)*box(2,3))
b(7)=box(2,1)*box(3,2)-(box(2,2)*box(1,3))
b(8)=box(1,2)*box(3,1)-(box(1,1)*box(3,2))
b(9)=box(1,1)*box(2,2)-(box(2,1)*box(1,2))
! calculate determinant
d=box(1,1)*b(1)+box(1,2)*b(4)+box(1,3)*b(7)
r=0.d0
if(abs(d)>0.d0)r=1.d0/d
! complete inverse matrix with all elements of matrix calculation
b(1)=r*b(1)
b(2)=r*b(2)
b(3)=r*b(3)
b(4)=r*b(4)
b(5)=r*b(5)
b(6)=r*b(6)
b(7)=r*b(7)
b(8)=r*b(8)
b(9)=r*b(9)
do i=1,9
write(4,*)i," element of inverse matrix ",b(i)
enddo
write(4,*)"determinant",d
return ! to return to the main program
end subroutine invert ! subroutine program was closed
!-g -fcheck=all -Wall -fbacktrace command to check error
!----------------------------------------------------------------------
I am fairly new to fortran coding and tried to write a RDF (radial distribution function) code for the analysis where during compilation I am getting out of bound array error for line 168 in my code.
Can anyone please check my code and help me to resolve the issue. I will be greatfull for the kindness.
Code:!RADIAL DISTRIBUTION FUNCTION CODE (WORK FOR BOTH NVT & NPT)
program gofr
!===============================================================================
! Gromacs trajectroy reading utility available on github
!===============================================================================
use gmxfort_trajectory
use gmxfort_utils
implicit none
type(Trajectory) :: trj
!=================================================================================
! This section is for the trajectory file input variables
!=================================================================================
integer :: ninput,ntot,inp,nsites,nframes ! input for gofr
integer,parameter :: maxinp=150 ! maximum number of trajectory files as an input
character(len=100) :: infile(maxinp)
integer :: ibegin(maxinp),iend(maxinp),nskip,nread,nused,intv
integer :: i1,j1,u1,n1
integer :: i,j,k,ia,im,ja,jm
integer :: total_atom,natms,nstep
integer :: nmol_cta,nmol_br,nmol_water,natom_cta,natom_br,natom_water
integer :: atom_cbw,ifile
integer :: total_frame
integer :: iatom_cta,iatom_br,iatom_water,imat,jmat
integer :: imol_cta,imol_br,imol_water
integer :: no_atname1
integer :: ibin,nprtot,npr(1000),nbins
integer :: iflg,keytrj,imcon,nmol_diff,new_file
real :: myatom(3)
real(kind=8),allocatable :: chge(:),weight(:)!,yyy(:),zzz(:)
real(kind=8),allocatable :: vxx(:),vyy(:),vzz(:),fxx(:),fyy(:),fzz(:)
real(kind=8):: xx1(1),yy1(1),zz1(1),xxx(1),yyy(1),zzz(1)
real(kind=8),allocatable :: rw(:,:,:),rcta(:,:,:),rbr(:,:,:)
real(kind=8) :: boxl(3,3),dr(3),rmsx,rmsy,rmsz
real(kind=8) :: timstp
real(kind=8),parameter :: au=0.5291770
!================================================================
real(kind=8),parameter :: drg=0.05 ! spacing in g(r)
!================================================================
real(kind=8) :: cell(3,3),rcell(9),det,dt0
!================================================================
real(kind=8) :: rmax,rdr,rsq,r,rr,rrs,cn,g,gavg,fact,cutoff_vol
real(kind=8) :: pi
real(kind=8),allocatable :: rs(:),gs(:)
real :: start,finish
character(len=10) :: nchar
character(len=100) :: outfile1,outfile2,a
character(len=8),allocatable :: atnm_rw(:,:),atnm_rcta(:,:),atnm_rbr(:,:),x(:,:)
character(len=8) :: atname1,atname2
character(len=8),allocatable :: atname(:)
character(len=80) :: title
! write (*,*) "type the trjectory file name"
call cpu_time(start)
!=======================================================================
open(unit=2,file='br_br_gofr.inp',action='read')
read(2,*)outfile1
read(2,*)outfile2
open(unit=3,form='formatted',file=outfile1,&
action='write')
open(unit=4,form='formatted',file=outfile2,&
action='write')
!=======================================================================
!=======================================================================
! Begining of the code for gofr calculation
!=======================================================================
write(*,*)'Enter the no of trajectory file you want to analyze'
read(2,*)ninput ! (maxinp=150)
write(*,*)'Along with this also enter the number of frames you',&
' want to work with.'
read(2,*)nframes ! number of time-steps used for binning of atoms
write(*,*)'Enter the number of different molecules types'
read(2,*)nmol_diff ! (1)
write(*,*)'Enter the no of molecules of all types'
read(2,*)nmol_water ! (884)
write(*,*)'Enter the no of atom in each of the diff molecules'
read(2,*) natom_water ! (3)
if(ninput>maxinp)then
write(*,*)'Too many inputs'
stop
endif
write(*,*)'Pair correlation atom pair X-Y: X'
read(2,*)atname1
write(*,*)'Pair correlation atom pair X-Y: Y'
read(2,*)atname2
write(*,*)'Enter the number of X atom'
read(2,*)no_atname1
!=======================================================================
atom_cbw=0
pi=datan(1.0d0)*4
intv=10
write(*,*)'pi====>',pi
total_atom=nmol_water*natom_water
natms=total_atom
!=======================================================================
ntot=0 ! total number of inputs in trajectory
!=======================================================================
write(*,*)'enter the skip value'
read(2,*)nskip
read(2,*)rmax ! half the box length
!write(4,*)'nskip,rmax=======>',nskip,rmax
if(nskip/=0)nskip=1
!ntot=ntot/nskip
!total_frame=ntot
write(*,*)'Finished reading inputs'
!=======================================================================
!GOFR CALCULATION
!=======================================================================
rdr=1.0d0/drg
nbins=nint(rmax*rdr)+1 !number of bins for calculation has been evaluated
if(real(nbins)*drg>rmax)nbins=nbins-1
do j=1,nbins
npr(j)=0
enddo
nprtot=0
!write(4,*)'rmax =',rmax,'nbins =',nbins
!=======================================================================
! from this line of code below trajectory file name was taken directly
! from the file
!=======================================================================
fread: do inp=1,ninput,1
write(*,*) 'Enter the name of trajectory file',inp
read (2,*) infile(inp)
read (2,*) ibegin(inp),iend(inp)
!=======================================================================
!Number of atoms was going to be read and passed to a variable name nsites
!=======================================================================
call trj%open(infile(inp))
nsites = trj%natoms()
!-----------------------------------------------------------------------
allocate(rw(nmol_water,natom_water,3))!,rcta(nmol_cta,natom_cta,3),rbr(nmol_br,natom_br,3))
allocate(atnm_rw(nmol_water,natom_water))!,atnm_rcta(nmol_cta,natom_cta),atnm_rbr(nmol_br,natom_br))
! allocate(atname(nsites))
!allocate(xx1(nsites),yy1(nsites),zz1(nsites))
!allocate(xxx(nsites),yyy(nsites),zzz(nsites))
!-----------------------------------------------------------------------
n1 = trj%read_next(nframes) !trj%read_next lib allow us to read frame wise data
!=======================================================================
!write(4,*)'number of sites = ',nsites
!write(4,*)'Total atom=',total_atom
! allocate(xxx(nsites),yyy(nsites),zzz(nsites))
! allocate(xx1(nsites),yy1(nsites),zz1(nsites))
!=======================================================================
!here the code will loop over the provided number of frames you want to
!work with as an input trajectory
!=======================================================================
nframe: do i1=1, n1 ! for nframes frames as mentioned
! write (3,*)"frame=========",i1
! write(*,*) i1,"frame============"
cell= trj%box(i1)
call invert (cell,rcell,det)
site: do j1=1,nsites
myatom = trj%x(i1,j1)
xxx = trj%frameArray(i1)%xyz(1,j1)
yyy = trj%frameArray(i1)%xyz(2,j1)
zzz = trj%frameArray(i1)%xyz(3,j1)
write(3,*) j1,"atom-----------------------------"
write(3,*) xxx , yyy, zzz
!write(*,*)'h1',j1,xxx(j1)
xx1=xxx(j1)*rcell(1)+yyy(j1)*rcell(4)+zzz(j1)*rcell(7) !
yy1=xxx(j1)*rcell(2)+yyy(j1)*rcell(5)+zzz(j1)*rcell(8) !
zz1=xxx(j1)*rcell(3)+yyy(j1)*rcell(6)+zzz(j1)*rcell(9)
enddo site
!write(3,*) j1,a
! write(4,*) j1,"atom "," frame",i1
!
!write(4,*) j1," atom ",xx1,yy1,zz1
!100 format (i5,2x,) !
wt: do imol_water=1,nmol_water
do iatom_water=1,natom_water
atom_cbw=atom_cbw+1
rw(imol_water,iatom_water,1)=xx1(atom_cbw)
rw(imol_water,iatom_water,2)=yy1(atom_cbw)
rw(imol_water,iatom_water,3)=zz1(atom_cbw)
atnm_rw(imol_water,iatom_water)=atname(atom_cbw)
!write(*,*) rw
enddo
enddo wt
do imat=1,3
do jmat=1,3
boxl(imat,jmat)=cell(imat,jmat)
enddo
enddo
br1: do im=1,nmol_water-1
do ia=1,natom_water
if(atnm_rw(im,ia)==atname1)then !pair function atm1 (cta)
br2: do jm=(im+1),nmol_water
do ja=1,natom_water
if(atnm_rw(jm,ja)==atname2)then !pair function atm2 (cta)
dr(:)=rw(im,ia,:)-rw(jm,ja,:)
dr(:)=dr(:)-anint(dr(:)) !Box length L=1
!CALCULATE CORRECTED REAL SPACE COORDINATES
rmsx=dr(1)*cell(1,1)+dr(2)*cell(2,1)+dr(3)*cell(3,1)
rmsy=dr(1)*cell(1,2)+dr(2)*cell(2,2)+dr(3)*cell(3,2)
rmsz=dr(1)*cell(1,3)+dr(2)*cell(2,3)+dr(3)*cell(3,3)
dr(1)=rmsx;dr(2)=rmsy;dr(3)=rmsz
rsq=dot_product(dr,dr)
r=(sqrt(rsq))
ibin=int(r*rdr)+1
if(ibin<=nbins)then
npr(ibin)=npr(ibin)+1
nprtot=nprtot+1
endif
endif
enddo
enddo br2
endif
enddo
enddo br1
enddo nframe
!
call trj%close()
enddo fread
100 continue
!-----------------------------------------------------------------------
cutoff_vol=((4.0d0/3.0d0)*pi*(rmax**3))
!-----------------------------------------------------------------------
cn=0.0
fact=(cutoff_vol)/(4.0d0*pi*(drg**3)*real(nprtot))
allocate(rs(nbins),gs(nbins))
do i=1,nbins
r=(real(i)-0.5d0)*drg
rr=((real(i)+0.5d0)**2+(1.d0/12.d0))
g=fact*(real(npr(i))/(rr))
!--------------------------------------------------------
cn=cn+real(npr(i))/real(nused*(no_atname1-1))
!--------------------------------------------------------
rs(i)=(real(i)-0.5d0)*drg
rrs=((real(i)+0.5d0)**2+(1.d0/12.d0))
gs(i)=fact*(real(npr(i))/(rrs))
write(3,2)r,g,cn
enddo
!=========================================================
2 format(2x,f10.5,f10.5,f10.5)
3 format(2x,f10.5,f10.5)
!-----------------------------------------
if(a=='smooth')then
write(4,3)rs(1),gs(1)
write(4,3)rs(2),gs(2)
do i=3,nbins-2
gavg=0.0d0
do j=i-2,i+2
gavg=gavg+gs(j)
enddo
gavg=(gavg/5.0)
write(4,3)rs(i),gavg
enddo
write(4,3)rs(nbins-1),gs(nbins-1)
write(4,3)rs(nbins),gs(nbins)
endif !For 5 point smoothing
!-----------------------------------------------------------------------
! deallocate(xxx,yyy,zzz,xx1,yy1,zz1)
deallocate(rw,atnm_rw,atname)
call cpu_time(finish)
print'("Time= ",f6.3,"seconds")',finish-start,finish,start
end program gofr
!---------------------------------------------------------------
subroutine invert(box,b,d) !subroutine sname (dummy arugments)
!---------------------------------------------------------------
! dl_poly subroutine to invert a 3*3 matrix using cofactors
! author - w. smith april 1992
!---------------------------------------------------------------
real(kind=8):: box(3,3) ! data passed to subroutine that cann't be changed
real(kind=8)::b(9),d
real(kind=8):: r ! r is a temporary variable, not accessible to main program
write(4,*) "in subroutine"
do i=1,3
write(4,*) (box(i,j),j=1,3)
enddo
! calculate adjoint matrix cofactors with simultaneous doing transpose of mtx
b(1)=box(2,2)*box(3,3)-(box(2,3)*box(3,2))
b(2)=box(3,2)*box(1,3)-(box(1,2)*box(3,3))
b(3)=box(1,2)*box(2,3)-(box(2,2)*box(1,3))
b(4)=box(3,1)*box(2,3)-(box(2,1)*box(3,3))
b(5)=box(1,1)*box(3,3)-(box(3,1)*box(1,3))
b(6)=box(2,1)*box(1,3)-(box(1,1)*box(2,3))
b(7)=box(2,1)*box(3,2)-(box(2,2)*box(1,3))
b(8)=box(1,2)*box(3,1)-(box(1,1)*box(3,2))
b(9)=box(1,1)*box(2,2)-(box(2,1)*box(1,2))
! calculate determinant
d=box(1,1)*b(1)+box(1,2)*b(4)+box(1,3)*b(7)
r=0.d0
if(abs(d)>0.d0)r=1.d0/d
! complete inverse matrix with all elements of matrix calculation
b(1)=r*b(1)
b(2)=r*b(2)
b(3)=r*b(3)
b(4)=r*b(4)
b(5)=r*b(5)
b(6)=r*b(6)
b(7)=r*b(7)
b(8)=r*b(8)
b(9)=r*b(9)
do i=1,9
write(4,*)i," element of inverse matrix ",b(i)
enddo
write(4,*)"determinant",d
return ! to return to the main program
end subroutine invert ! subroutine program was closed
!-g -fcheck=all -Wall -fbacktrace command to check error
!----------------------------------------------------------------------