Fortran Program for Simulation of Coupled Harmonic Oscillator

program oscillator implicit none real(kind=8):: amass(2) !masses integer,parameter :: ndim=1 !no. of dimensions = 3 integer,parameter :: natoms = 2 !no. of atoms integer,parameter :: imax=400 integer :: tmax real(kind=8):: mxyz(ndim,natoms) !mean position coordinates real(kind=8):: rxyz(ndim,natoms) !instanteous position coordi real(kind=8):: qxyz(ndim,natoms) !coordi wrt to mean position real(kind=8):: vxyz(ndim,natoms) !velocity real(kind=8):: dt=0.1d0 !time interval real(kind=8):: k(2) !k1=m1*w^2,same for other real(kind=8):: kc !of coupling spring integer :: i,j,istep,e=0 real(kind=8):: fxyz(ndim,natoms) !past force real(kind=8):: ffxyz(ndim,natoms) !new force real(kind=8),allocatable :: v10(:,:),v20(:,:) ,v1(:,:),v2(:,:) !v0 is the velocity at some reference time usually t0 or t1 !v1,v2 to hold instantaneous velocity real(kind=8),allocatable :: v0xyz(:,:,:) real(kind=8),allocatable :: vtxyz(:,:,:) real(kind=8):: c !for VAF real(kind=8),allocatable :: vaf(:) ! vaf(imax) integer :: kk,n !kk,n used in line 216 in is index in fourier transform real(kind=8):: b !219 temp varaible to values real(kind=8), allocatable :: M1(:,:),M2(:,:) !matrix for real components of fourier transform real(kind=8), allocatable :: M(:,:) !matrix for imag compo of foureir transform complex(kind=8)::cpx !complex variable not used finally real(kind=8),allocatable :: g1(:) !fourier transform of VAF real(kind=8),allocatable :: g2(:) !fourier transform of VAF real(kind=8),allocatable :: g(:) !fourier transform of VAF real(kind=8)::ke,pe,te !energy symbols real(kind=8),allocatable :: vv(:,:) !pe matrix real(kind=8),allocatable :: vvv(:,:) real(kind=8),allocatable :: tt(:,:) !ke matrix real(kind=8),allocatable :: sqrtmassinv(:) real(kind=8),allocatable :: kkk(:) real(kind=8),allocatable :: lenn(:) real(kind=8)::x,y,z integer,allocatable :: point(:) !=========== !sagar real(kind=8):: v_sum!=====================!for finding the eigen values INTEGER LDA !leading dimension of A herea A= VV PARAMETER ( LDA = natoms*ndim ) INTEGER LWMAX PARAMETER ( LWMAX = 1000 ) INTEGER INFO, LWORK DOUBLE PRECISION A( LDA, natoms*ndim ), W( natoms*ndim ), WORK( LWMAX ) DOUBLE PRECISION ww(natoms*ndim) !eigen frequeniesinteger :: ITYPE = 1 !tell the type of problem to be solvedcharacter(len=1) :: JOBZ ='V' !tell that vectors will also be calculatedcharacter(len=1) :: UPLO = 'U' !upper part of A and B are storedinteger,parameter :: LDB=natoms*ndim!leading dim of B!double precision BB( LDB, natoms ) !INOUT!===================!for solving linear equations !real(kind=8),allocatable :: BBB(:,:) !required by subroutine to store values while making result real(kind=8),allocatable :: TTT(:,:) ! temperory arry to store A INTEGER LB !real(kind=8),allocatable :: QQQ(:,:) ! Normal Coordinates, these will not be req. now real(kind=8),allocatable :: XXX(:,:) ! arr. to hold qxyz !AQ=q,AX=B, our aim is finding Q , q=XXX, Q=QQQ integer :: NRHS = 1 integer :: IPIV(natoms*ndim) !allocate (BBB(natoms,natoms)) allocate (TTT(natoms*ndim,natoms*ndim)) !BBB=0.d0 !initialization TTT=0.d0!-----------------

!----------------- allocate(M1(imax,imax)) allocate(M2(imax,imax)) allocate(M(imax,imax)) allocate(vaf(imax)) allocate(g1(imax)) allocate(g2(imax)) allocate(g(imax)) allocate(vv(natoms*ndim,natoms*ndim)) allocate(vvv(natoms*ndim,natoms*ndim)) allocate(tt(natoms,natoms)) allocate (XXX(natoms*ndim,imax)) !temp arr. to hold generl. coordi. ! allocate (QQQ(natoms,imax)) ! normal coordi. allocate (sqrtmassinv(natoms*ndim)) allocate (kkk(natoms*ndim)) allocate (lenn(natoms)) allocate ( v10(ndim,imax)) allocate ( v20(ndim,imax)) allocate ( v1(ndim,imax)) allocate ( v2(ndim,imax)) tmax = imax/10 != 40 allocate (v0xyz(ndim,tmax,natoms)) allocate (vtxyz(ndim,imax,natoms)) allocate (point(imax)) point=0 do i =1,tmax!=40 point(i*10)= 10*i end dog1=0.d0!initializationg2=0.d0g=0.d0 read (*,*) amass(1:natoms)! amass(1:2)= (/1.d0,1.d0/) do i=1,natoms do j=1,ndim sqrtmassinv(ndim*(i-1)+j)=1.d0/dsqrt(amass(i))! massinv(ndim*(i-1)+2)=1.d0/dsqrt(amass(i))! massinv(ndim*(i-1)+3)=1.d0/dsqrt(amass(i)) end do end do print*,"sqrtmassinv",sqrtmassinv(1:natoms*ndim) print*,"mass1 and mass2", amass(1:2)! natoms = 2 !dt = 0.1d0! kc = 2.0d0! k(1:2) = (/7.d0,1.d0/) read (*,*) kc,k(1:natoms) print*,"kc,k1 and k1",kc,k(1:2) do i=1,natoms do j=1, ndim kkk(ndim*(i-1)+j)=k(i) !kkk(ndim*(i-1)+2)=k(i) !kkk(ndim*(i-1)+3)=k(i) end do end do print*,"kkk:",KKK(1:ndim*natoms)rxyz=0.d0vxyz=0.d0mxyz=0.d0fxyz=0.d0ffxyz=0.d0qxyz=0.d0 read (*,*) rxyz(1:ndim,1)! rxyz(1:ndim,1)=(/2.2d0/)!,0.d0/)!,0.d0/) print*,"pos1",rxyz(1:ndim,1) read (*,*) rxyz(1:ndim,2)

! rxyz(1:ndim,1)=(/2.2d0/)!,0.d0/)!,0.d0/) print*,"pos1",rxyz(1:ndim,1) read (*,*) rxyz(1:ndim,2) ! rxyz(1:ndim,2)=(/5.2d0/)!,0.d0/)!,0.d0/) print*,"pos2",rxyz(1:ndim,2) read (*,*) vxyz(1:ndim,1) ! vxyz(1:ndim,1)=(/0.d0/)!,0.d0/)!,0.d0/) print*,"vel1",vxyz(1:ndim,1) read (*,*) vxyz(1:ndim,2) ! vxyz(1:ndim,2)=(/0.d0/)!,0.d0/)!,0.d0/) print*,"vel2",vxyz(1:ndim,2) read (*,*) mxyz(1:ndim,1) !mxyz(1:ndim,1)=(/2.d0/)!,0.d0/)!,0.d0/) print*,"mean position:",mxyz(1:ndim,1) read (*,*) mxyz(1:ndim,2) !mxyz(1:ndim,2)=(/5.d0/)!,0.d0/)!,0.d0/) print*,"mean position:",mxyz(1:ndim,2) do i=1,natoms do j=1,ndim qxyz(j,i) = rxyz(j,i) - mxyz(j,i) !qxyz(2,i) = rxyz(2,i) - mxyz(2,i) !qxyz(3,i) = rxyz(3,i) - mxyz(3,i) end do end do print*,"qxyz1",qxyz(1:ndim,1) print*,"qxyz2",qxyz(1:ndim,2) lenn=0.d0 do i=1,natoms do j=1,ndim lenn(i) = lenn(i) + qxyz(j,i)**2 !print*,"len",i,dsqrt(qxyz(1,i)**2+qxyz(2,i)**2+qxyz(2,i)**2) !print*,"len2",dsqrt(qxyz(1,2)**2+qxyz(2,2)**2+qxyz(2,2)**2) end do end do do i=1,natoms print*,"len",i,dsqrt(lenn(i)) end do open (11,file="trajectory.xyz",status='unknown') open (12,file="pos",status='unknown') open (13,file="VVAF",status='unknown') open (14,file="ene",status='unknown')do istep =1,imax !fortran 90 print*,"--------------" do i=1, 2 do j=1,ndim fxyz(j,i) = -1*k(i)*qxyz(j,i) + ((-1)**i)* kc*( qxyz(j,1)-qxyz(j,2) ) !fxyz(1,i) = -1*k(i)*qxyz(1,i) + ((-1)**i)* kc*( qxyz(1,1)-qxyz(1,2) ) !fxyz(2,i) = -1*k(i)*qxyz(2,i) + ((-1)**i)* kc*( qxyz(2,1)-qxyz(2,2) ) !fxyz(3,i) = -1*k(i)*qxyz(3,i) + ((-1)**i)* kc*( qxyz(3,1)-qxyz(3,2) ) end do end do print*, print*,"fxyz" print*,fxyz(1:ndim,1) print*,fxyz(1:ndim,2)

DO i=1,natoms do j=1, ndim rxyz(j,i) =rxyz(j,i)+vxyz(j,i)*dt +fxyz(j,i)*dt*dt*0.5d0/amass(i) !rxyz(1,i) =rxyz(1,i)+vxyz(1,i)*dt +fxyz(1,i)*dt*dt*0.5d0/amass(i) !rxyz(2,i) =rxyz(2,i)+vxyz(2,i)*dt +fxyz(2,i)*dt*dt*0.5d0/amass(i) !rxyz(3,i) =rxyz(3,i)+vxyz(3,i)*dt +fxyz(3,i)*dt*dt*0.5d0/amass(i) end do END DO print*,"pos1",rxyz(1:ndim,1) print*,"pos2",rxyz(1:ndim,2) do i=1,natoms do j=1,ndim qxyz(j,i) = rxyz(j,i) - mxyz(j,i) !qxyz(2,i) = rxyz(2,i) - mxyz(2,i) !qxyz(3,i) = rxyz(3,i) - mxyz(3,i) end do end do print*,"qxyz1",qxyz(1:ndim,1) print*,"qxyz2",qxyz(1:ndim,2) print*,"qxyz2",qxyz(1:ndim,2) lenn=0.d0 do i=1,natoms do j=1,ndim lenn(i) = lenn(i) + qxyz(j,i)**2 !print*,"len",i,dsqrt(qxyz(1,i)**2+qxyz(2,i)**2+qxyz(2,i)**2) !print*,"len2",dsqrt(qxyz(1,2)**2+qxyz(2,2)**2+qxyz(2,2)**2) end do end do do i=1,natoms print*,"len",i,dsqrt(lenn(i)) end do !print*,"len1",dsqrt(qxyz(1,1)**2+qxyz(2,1)**2+qxyz(2,1)**2) !print*,"len2",dsqrt(qxyz(1,2)**2+qxyz(2,2)**2+qxyz(2,2)**2)

do i=1, 2 do j=1,ndim ffxyz(j,i) = -1*k(i) *qxyz(j,i)+((-1)**i) *kc*( qxyz(j,1)-qxyz(j,2) ) !ffxyz(1,i) = -1*k(i) *qxyz(1,i)+((-1)**i) *kc*( qxyz(1,1)-qxyz(1,2) ) !ffxyz(2,i) = -1*k(i) *qxyz(2,i)+((-1)**i) *kc*( qxyz(2,1)-qxyz(2,2) ) !ffxyz(3,i) = -1*k(i) *qxyz(3,i)+((-1)**i) *kc*( qxyz(3,1)-qxyz(3,2) ) end do end do print*, print*,"ffxyz" print*,ffxyz(1:ndim,1) print*,ffxyz(1:ndim,2) DO i=1,natoms do j=1,ndim vxyz(j,i)=vxyz(j,i)+0.5D0*dt*( fxyz(j,i) + ffxyz(j,i) )/amass(i) !vxyz(1,i)=vxyz(1,i)+0.5D0*dt*( fxyz(1,i) + ffxyz(1,i) )/amass(i) !vxyz(2,i)=vxyz(2,i)+0.5D0*dt*( fxyz(2,i) + ffxyz(2,i) )/amass(i) !vxyz(3,i)=vxyz(3,i)+0.5D0*dt*( fxyz(3,i) + ffxyz(3,i) )/amass(i) end do END DO print*,"vel1",vxyz(1:ndim,1) print*,"vel2",vxyz(1:ndim,2) if (istep==point(istep) .or. istep==1) then e = e + 1! vxyz(ndim,tmax,natoms) finally e is tmax=40 print*,"e",e v0xyz(1:ndim,e,1)=vxyz(1:ndim,1) v0xyz(1:ndim,e,2)=vxyz(1:ndim,2) !v10(1:ndim,e) = vxyz(1:ndim,1) !Ref vel !v20(1:ndim,e) = vxyz(1:ndim,2) !print*,"v10",v10(1:ndim) !print*,"v10",v20(1:ndim) end if !=========writing the trajectories=======

y=0.d0;z=0.d0 write(11,*) 2 write(11,*) write(11,*) "o", rxyz(1:ndim,1),y,z!0.d0,0.d0 write(11,*) "o", rxyz(1:ndim,2),y,z!0.d0,0.d0 write(12,*) istep, rxyz(1:ndim,1),rxyz(1:ndim,2) !----------------------------------------- v1(1:ndim,istep)=vxyz(1:ndim,1) !print*,"v1",v1(1:ndim) v2(1:ndim,istep)=vxyz(1:ndim,2) !print*,"v2",v2(1:ndim) !vtxyz(ndim,imax,natoms) vtxyz(1:ndim,istep,1) = vxyz(1:ndim,1) vtxyz(1:ndim,istep,2) = vxyz(1:ndim,2) !velocity auto corelation function !http://www.ccp5.ac.uk/DL_POLY/Democritus/Theory/vaf.html ! vaf(istep) = dot_product(v1,v10)/dot_product(v10,v10)/natoms/natoms + dot_product(v2,v20)/dot_product(v20,v20)/natoms/natoms !VAF !vaf(istep) = dot_product(vxyz(:,1),v10(:,1))/natoms + dot_product(vxyz(:,2),v20(:,1))/natoms !VAF !vaf(istep) = dot_product(v1(:,istep),v10(:,1))/natoms + dot_product(v2(:,istep),v20(:,1))/natoms !VAF ! vaf(istep) = dot_product(v1(:,istep),v0xyz(:,1,1))/natoms + dot_product(v2(:,istep),v0xyz(:,1,2))/natoms !VAF !v0xyz(ndim,tmax,matoms) !vtxyz(ndim,imax,natoms) ! do i=1, natoms ! vaf(istep) = vaf(istep)+ dot_product(vxyz(:,i),v0xyz(:,1,i))/natoms ! + dot_product(vxyz(:,2),v0xyz(:,1,2))/natoms !VAF ! end do

!========= writing the energies ========== ke= 0.5d0*amass(1)*vxyz(1,1)**2 + 0.5d0*amass(2)*vxyz(1,2)**2 pe= 0.5d0*k(1)*qxyz(1,1)**2 + 0.5d0*k(2)*qxyz(1,2)**2 + 0.5d0*kc*(qxyz(1,1)-qxyz(1,2))**2 te= ke+pe write(14,*) istep*dt ,ke,pe,te !---------------------------------------- do i= 1, natoms do j =1,ndim XXX(ndim*(i-1)+j,istep) = qxyz(j,i) !AX=B; AQ=q; A=A,X=Q,B=q=XXX ! XXX(i,istep) = dsqrt (qxyz(1,i)**2+qxyz(2,i)**2+qxyz(3,i)**2) end do enddo print*,"xxx(1),xxx(2)",XXX(1,istep),XXX(2,istep) end do vaf(:)=0.d0 do istep=1,imax do i=1, natoms do j=1,tmax vaf(istep) = vaf(istep)+ dot_product(vtxyz(:,istep,i),v0xyz(:,j,i))/natoms/tmax/dot_product(v0xyz(:,j,i),v0xyz(:,j,i)) ! + dot_product(vxyz(:,2),v0xyz(:,1,2))/natoms !VAF !vaf(istep) = vaf(istep)+ vtxyz(1,istep,i)*v0xyz(1,j,i) /dot_product(v0xyz(:,j,i),v0xyz(:,j,i)) end do end do end do !v_sum=(vtxyz(1,1,1)*vtxyz(1,1,1))+(vtxyz(1,1,2)*vtxyz(1,1,2)) !============================ !writing the VAF in the file do i=1, imax write(13,*) i,vaf(i)!/v_sum end do !----------------------------

!===================================== !finding fourier tranf for real compo. do kk=1,imax !k for freq do n=1,imax !n for no of VAF M1(kk,n)=0.d0 b=0.d0 !b=-2.d0*3.141d0*((-imax/2+kk)-1)*((-imax/2+n)-1)/imax b=-2.d0*3.141d0*((-imax/2+kk)-1)*(n-1)/imax M1(kk,n)=cos(b) end do end do call dgemm('n','n' ,imax,1,imax,& 1.d0 ,M1,imax ,VAF,imax,& 1.d0 ,g1,imax) !g1 is the fourier tranform for real part !-------------------------------------- !===================================== !finding fourier tranf for imaginary part do kk=1,imax !k for freq do n=1,imax !n for no of VAF M2(kk,n)=0.d0 b=0.d0 !b=-2.d0*3.141d0*((-imax/2+kk)-1)*((-imax/2+n)-1)/imax b=-2.d0*3.141d0*((-imax/2+kk)-1)*(n-1)/imax M2(kk,n)=sin(b) end do end do call dgemm('n','n' ,imax,1,imax,& 1.d0 ,M2,imax ,VAF,imax,& 1.d0 ,g2,imax) !g2 is the fourier transform for imag part !-------------------------------------

!================================ ! in this part i have taken magnitude of complex number ! this is just for watching what is the outcome do i=1 , imax do j=1, imax M(i,j)=dsqrt( M1(i,j)**2 + M2(i,j)**2 ) end do end do call dgemm('n','n' ,imax,1,imax,& 1.d0 ,M,imax ,VAF,imax,& 1.d0 ,g,imax) !g is the F.T. !-------------------------------- !================================= !writing the FT in the files open(111,file="g1",status="unknown") do kk=0,imax-1 !imax = imax write(111,*) ((-imax/2+kk))*2.0d0*3.141d0/imax/dt,g1(kk+1)**2 end do print*,"radian freq interval :", 2.0d0*3.141d0/imax/dt open(112,file="g2",status="unknown") do kk=0,imax-1 write(112,*) ((-imax/2+kk))*2.0d0*3.141d0/imax/dt,g2(kk+1)**2 end do open(113,file="g",status="unknown") do kk=0,imax-1 write(113,*) ((-imax/2+kk)-1)*kk*2.0d0*3.141d0/imax/dt,g(kk+1)**2 end do !=========== hessian matrix =============

!=========== hessian matrix ============= !do i=1,natoms*ndim do i=1,natoms*ndim do j=1, natoms*ndim if(i==j) then vvv(i,j)= (kc + kkk(i)) * sqrtmassinv(i) * sqrtmassinv(j) else vvv(i,j)= (-1)*kc * sqrtmassinv(i) * sqrtmassinv(j) end if end do end do print*, "Mass Weighted Hessian matrix" do i=1,natoms*ndim print*,vvv(i,1:natoms*ndim) end do !vvv=0.d0 !...................................... !============ Mass Weighted Hessian ====================! do i=1,natoms! do j=1, natoms! if(i==j) then ! vv(i,j)= (kc + k(i)) / (dsqrt(amass(i)) * dsqrt(amass(j)))! else! vv(i,j)= (-1)*kc / (dsqrt(amass(i)) * dsqrt(amass(j)))! end if ! end do! end do! print*, "pe matrix"! do i=1,natoms*ndim! print*,vv(i,1:natoms*ndim)! end do !======= mass matrix matrix========! do i=1,natoms! do j=1, natoms! if(i==j) then! tt(i,j)= amass(i)! else! tt(i,j)= 0.d0! end if! end do! end do!! print*, "ke matrix"! do i=1,natoms!! print*,tt(i,1:natoms)! end do

!====== finding normal frequencies and e. vectors ========= A=vvv ! BB=tt LWORK = -1 !call dsygv( ITYPE,'V','U',natoms, A,LDA, BB,LDB,W,WORK,LWORK,INFO) ! CALL DSYEV( 'V', 'U', N1, A, LDA, W, WORK, LWORK, INFO ) CALL DSYEV( 'V', 'U', natoms*ndim, A, LDA, W, WORK, LWORK, INFO ) LWORK = MIN( LWMAX, INT( WORK( 1 ) ) ) ! Solve eigenproblem. ! !call dsygv( ITYPE,'V','U',natoms, A,LDA, BB,LDB,W,WORK,LWORK,INFO) ! CALL DSYEV( 'V', 'U', N1, A, LDA, W, WORK, LWORK, INFO ) CALL DSYEV( 'V', 'U', natoms*ndim, A, LDA, W, WORK, LWORK, INFO ) ! ! Check for convergence. IF( INFO.GT.0 ) THEN WRITE(*,*)'The algorithm failed to compute eigenvalues.' STOP END IF !Print*," eigenvalues" CALL PRINT_MATRIX( 'Eigenvalues', 1, natoms, W, 1 ) print*,"normal frequencies are:" do i=1,natoms*ndim ww(i)=dsqrt(abs( w(i))) ! ww(2)=dsqrt(abs( w(2))) ! ww(1)=sqrt(abs( w(1))/amass(1) ) ! ww(2)=sqrt(abs( w(2))/amass(2) ) end do print*,ww(1:natoms*ndim) print*,"---------------" ! Print eigenvectors. A= -A CALL PRINT_MATRIX( 'Eigenvectors (stored columnwise)', natoms*ndim, natoms*ndim, A, & LDA ) !============================================ !computing the normal coodinates !AQ=q, XXX= q, TTT = A, solving AX=B do i=1, imax TTT = A call DGESV(natoms*ndim,NRHS,TTT,LDA,IPIV,XXX(:,i),LDB, INFO) !on exit XXX will hold Q, thus a new array XXX of normal coord. is prepared end do IF( INFO.GT.0 ) THEN WRITE(*,*)'The diagonal element of the triangular factor of A,' WRITE(*,*)'U(',INFO,',',INFO,') is zero, so that' WRITE(*,*)'A is singular; the solution could not be computed.' STOP END IF !============================================ ! finding normal coordinates! TTT = A !call dcopy(natoms*natoms,A,1,TTT,1) ! T = AAA ! !---find inverse of A matrix---! IF(natoms.EQ.1) THEN! TTT(1,1)=1.D0/TTT(1,1)! ELSE! LB=(natoms-1)*natoms !Compute an LU factorization! CALL DGETRF(natoms,natoms,TTT,natoms,BBB(1,1),INFO)! IF(INFO.EQ.0) THEN! !Compute the inverse of a matrix using the LU factorization! CALL DGETRI(natoms,TTT,natoms,BBB(1,1),BBB(1,2),LB,INFO)! ENDIF! ENDIF

! print *, "--------inverse of A matrix------------" ! do i=1, natoms ! print*,"TTT(",i,",1:",natoms,")",TTT(i,1:natoms) ! end do !-----solving linear equations-------------- ! do i=1,imax ! QQQ(:,i) = matmul(TTT,XXX(:,i)) ! end do ! call dgemm('n','n' ,natoms,1,natoms,& ! 1.d0 ,TTT,natoms ,XXX,natoms,& ! 1.d0 ,QQQ,natoms) !=================printting norm cor. in file==================== !if both the particle vibrate in same phase (mode1) than Q2=zero, but Q=not zero but oscilatory !if both the particle vibrate in opp. phase (mode2) the Q1=0 but Q2=not zero !at any instant the displacement of pendulum is obtained by superposition of harmonic osci.of w1 and w2 !where w1 and w2 are normal frequencies !=====writing to file ===================== open (15,file="norml_cor",status='unknown') do i=1, imax !imax = no.of time samples=no. of md steps write(15,*) i, XXX(1:ndim*natoms,i) !will write all the nor. coordi. for istep one by one end do !========================================== print*,"kkk:",KKK(1:ndim*natoms) print*,"sqrtmassinv",sqrtmassinv(1:natoms*ndim) do i=1,natoms do j=1,ndim print*,"sqrtmassinv" ,ndim*(i-1)+j,":",sqrtmassinv(ndim*(i-1)+j)

! massinv(ndim*(i-1)+2)=1.d0/dsqrt(amass(i)) ! massinv(ndim*(i-1)+3)=1.d0/dsqrt(amass(i)) end do end do do i =1,imax!=400 ! print *,"point",i, point(i) end do end program !========================================== ! Auxiliary routine: printing a matrix. SUBROUTINE PRINT_MATRIX( DESC, M, N, A, LDA ) CHARACTER*(*) DESC INTEGER M, N, LDA DOUBLE PRECISION A( LDA, * ) INTEGER I, J WRITE(*,*) WRITE(*,*) DESC DO I = 1, M WRITE(*,9998) ( A( I, J ), J = 1, N ) END DO 9998 FORMAT( 11(:,1X,F6.2) ) RETURN END !============================================ !! For any query call me 9009035593

!! Email satyends@iitk.ac.in