split a program in two subroutine

roddur · ‎11-14-2009

Hello friends,
I have this code in hand and my job is to split this in two subroutine, from atompp. (Sorry for putting the long code, but hope this may be helpful)

[cpp]c#define LM
      program lm
!    subroutine lm()
      implicit none
! Common blocks:
! --- dynamical memory allocation:
      integer wksize
      parameter(wksize=8000000)
      integer w(wksize)
      common/wkarr/w
! Controls for IO
      integer ishow,ncnt1,ncnt2,mcount
      parameter (ncnt1=1,ncnt2=2,mcount=20)
      character*4 cquit,err,output
      character*8 prgnam
! Parameters for options, structure, basis
      integer nbas,nclass,nl,nlrho,nsp,mmixat
      integer nit,nitat,ndel(3),ltmax(3),ifmt3d
      double precision bound(3,6),origin(3),plat(9,6),alat(2)
      logical switch(31)
!
! for iterations
      integer niter,mmixpq
      double precision beta,rmsdel,rms2,toleto,wc
! Class parameters
      integer mxclas
      parameter (mxclas=200)
      integer lmx(mxclas),lmxrho(mxclas),nrclas(mxclas),
     .        konfig(0:3,mxclas)
      double precision dq(mxclas),qc(mxclas),qtot(2*mxclas),
     .               ves(2*mxclas),vrmax(2*mxclas),wsr(mxclas),z(mxclas)
      character*4 clabl(mxclas)
! Symmetry group
      integer nsymop,ngen
      double precision symopm(9,48),symopv(3,48)
      character*72 symops(48),gens(48)
! BZ - integration / empty-sphere finder
      integer nkxyz(3),norder,nopts,npts
      double precision emin,eminc,emax,emaxc,range,tolef,width
! Empty-sphere finder
      double precision rmaxes,rmines
      integer nrxyz(3)
! Scaling wsr
      double precision asavol,facvol,gamma,ommax1(3),ommax2(3)
! Structure constant parameters
      integer ialpha,itrans,lmaxw,ndimin,jbasdn,ldn,mdn,nsmesh
      double precision dele,deltr,kap2,kfit,rmaxs,rmaxs2,eps
! Parameters relating to Ewald sums
      integer nkdmx
      double precision as,awald,tolews,vol
! for rhofit
      integer rfout
! for SYML
      integer mxline
      parameter (mxline=99)
      integer nq(mxline),nline
      double precision q1(3,mxline),q2(3,mxline)
      character*1 lab1(mxline),lab2(mxline)
! Local variables, functions and offsets to dynamic arrays
      integer fopn,i,ibas,ic,icalc,iprint,isym(10),j,lunit,len2,mxkp,
     .               nkd,nkg,nkp,npp,nrxc,obeta,oadot0,oalph0,obas,obk,
     .        odlat,oenu,oevrl,oglat,oibk,oiclas,oidmod,oidxdn,onrcls,
     .        opnu,opp,oqnu,osigma,osigm2,otrtab,owk,owtkp,og
      parameter(npp=6)
      double precision avw,avwsr,dovl1(3),dovl2(3),efermi,
     .       qb(9),vmtz
      logical begmom,convrg,lsmall,t,f,lpot,lsymst,lsdm1,nowrit
      character*30 csym(7)
      character*72 name,tit
      parameter(t=.true.,f=.false.)
! External calls
      external avwsr,chkcls,chkdis,chksym,chkwsr,defdr,defi,dscal,fexit,
     .         finits,fopn,icopy,iprint,lunit,ioctrl,allocm,wkinit
      integer isum,ldim,nderiv,neighm,nmesh,ntet,
     .        oadot,oadoto,oalpha,oalpho,ocg,ocy,odmad,odmat,
     .        oeband,oemom,oiax,oidn,oidtet,oindcg,oirr,oixirr,
     .        ojcg,onpr,opold,oqold,orwats,os,osdot,osdm1,osdm1s,
     .        otrad,otral
      double precision detot,ebmax,ebmin,etot,etoto,sumdq,zval
      character*4 cnit
      data prgnam/'LM'/
      data ialpha/0/,ifmt3d/1/,itrans/0/,jbasdn/1/,ldn/0/,lmaxw/8/,
     .     ltmax/2,2,2/, mdn/0/, mmixat/2/,
     .     mmixpq/1/, nbas/1/, nclass/1/, ndimin/350/,ngen/0/,nit/1/,
     .     nitat/30/,niter/30/,nkdmx/250/,nkxyz/4,4,4/,nl/3/,
     .     nopts /8001/, norder/0/,npts/1001/,nrxc/1/,nrxyz/0,0,0/,
     .     nsp/1/
      data as/.2d1/, beta/.5d0/,dele/.05d0/,deltr/.1d0/, efermi/-.25d0/,
     .     kap2/.0d0/,kfit/0.d0/ ,emax/.2d1/ ,emaxc/.2d1/ ,emin/-.2d1/,
     .     eminc/-.2d1/,eps/.0d0/,facvol/.1d1/,gamma/.0d0/,
     .     ommax1/.16d0,.18d0,.20d0/,ommax2/.40d0,.45d0,.50d0/,
     .     range/.5d1/  ,rmaxes/.40d1/,
     .     rmaxs/.32d1/,rmaxs2/.35d1/,rmines/.9d0/,rms2/99.9999999d0/,
     .     rmsdel/.1d-4/,tolef/.1d-5/,toleto/.1d-4/,tolews/.1d-5/,
     .     vmtz/-.75d0/,wc/-.1d1/,width/.5d-2/
      data clabl/mxclas*'E'/, cquit/'----'/, csym/7*'undef'/,
     .     err/'ERR '/, gens/'E',47*' '/, output/'LM  '/
      data convrg/f/, switch/t,t,f,t,f,f,f,f,f,t,  t,t,t,f,t,f,f,t,f,f,
     .                       f,f,f,f,t,t,f,f,f,t,  f/

!------------------------------------------
    write(*,*) "LM called"
!------------------------------------------

! --- time measure
      call cpusrt('TOTL',t)
      call cpusrt('PRWK',t)
! --- Set up dynamic memory allocation and other initialization ---
      call wkinit(wksize)
     
! --- Set up the file logical units ---
      call finits
! --- read CTRL file
      do ishow = ncnt1,ncnt2
        do i = 0, mcount
          call ioctrl(i,prgnam,w(oadot0),alat,w(oalph0),as,avw,
     .                w(obas),beta,bound,clabl,cquit,csym,dele,deltr,
     .                efermi,emax,emaxc,emin,eminc,eps,err,facvol,
     .                gamma,gens,ialpha,icalc,w(oiclas),w(oidmod),
     .                w(oidxdn),ifmt3d,ishow,isym,itrans,jbasdn,kap2,
     .                kfit,konfig,lab1,lab2,ldn,lmaxw,lmx,lmxrho,
     .                lsmall,ltmax,mcount,mdn,mmixat,mmixpq,nbas,
     .                nclass,ndel,ndimin,ngen,nitat,niter,nkdmx,nkp,
     .                nkxyz,nl,nline,nlrho,norder,nopts,npts,nq,nrclas,
     .                nrxc,nrxyz,nsp,nsymop,ommax1,ommax2,origin,output,
     .             plat,w(opnu),w(opp),q1,q2,w(oqnu),range,rfout,rmaxes,
     .                rmaxs,rmaxs2,rmines,rmsdel,w(osigma),w(osigm2),
     .                switch,tit,tolef,toleto,tolews,ves,vmtz,vrmax,wc,
     .                width,wsr,z)
          if (i.eq.6.and.ishow.eq.1) then
! --------- Allocate memory for pp's, moments, basis ....
            call allocm(mxclas,nbas,nclass,nl,nlrho,npp,nsp,oadot0,
     .                  oalph0,obas,oenu,oevrl,oiclas,oidmod,oidxdn,
     .                  opnu,opp,oqnu,osigma,osigm2)
          elseif (i.eq.8.and.ishow.eq.1) then
! --------- Check the symmetry and eventually complete the basis
            call chksym(alat,clabl,csym,gens,isym,switch(30),switch(10),
     .                lsmall,nbas,nclass,ngen,nrclas,nsymop,obas,oiclas,
     .                otrtab,plat,prgnam,symopm,symops,symopv,switch(3))
! --------- Average Wigner Seitz radius ---
            avw = avwsr(plat,alat,vol,nbas)
            if (.not.switch(3).and..not.lsmall)
     .       call chkcls(w(obas),clabl,w(oiclas),w(oidmod),w(oidxdn),
     .                   w(otrtab),konfig,lmx,nbas,nclass,nl,nrclas,nsp,
     .                   nsymop,w(opnu),w(opp),prgnam,w(oqnu),wsr,z)
          elseif (i.eq.9.and.ishow.eq.1) then
           call sclwsr(alat,asavol,w(obas),clabl,dovl1,dovl2,facvol,
     .                 gamma,w(oiclas),switch(25),nbas,nclass,nrclas,
     .                 ommax1,ommax2,plat,wsr,z)
           call chkwsr(alat,clabl,nclass,nrclas,plat,wsr,z)
           if(dabs(asavol-facvol).gt.1.d-5)
     .        call errmsg
     .     (' Not space filling: increase number of ES or OMMAX $',4)
          endif
        enddo
        if (ishow.eq.0) then
          call fexit(0,' HELP=T encountered')
        elseif (ishow.eq.1) then
          call chkdis(alat,w(obas),clabl,w(oiclas),nbas,nclass,plat,z)
          call gtname(len2,nclass,nrclas,z,name)
          tit = name(1:len2)//', '//csym(2)
         if (switch(17).or.switch(8).or.prgnam.eq.'LMHART') then
! -------- Get irreducible mesh points
           nmesh=(ndel(1)+1)*(ndel(2)+1)*(ndel(3)+1)
           call defi(oirr  ,  nmesh)
           call defi(oixirr,4*nmesh)
           call makirr(bound,w(oirr),w(oixirr),ndel,nsymop,origin,
     .                 plat,symopm,symopv)
           call rlse(oixirr)
         endif        
! --- Set up BZ mesh ---
      mxkp = nkxyz(1)*nkxyz(2)*nkxyz(3)
      call defi (oibk ,mxkp)
      call defdr(owtkp,mxkp)
      call defdr(obk  ,mxkp*3)
      call bzmesh(w(obk),w(oibk),nkxyz(1),nkxyz(2),nkxyz(3),nkp,nsymop,
     .            plat,qb,symopm,w(owtkp))
      call rlse(obk)
      call defdr(obk,3*nkp)
      if (switch(11).and.(switch(12).or.prgnam.eq.'LMDOS')) then
! ----- find inequivalent tetrahedra ---
        call defi(oidtet,mxkp*30)
        call tetirr(w(oibk),w(oidtet),nkxyz(1),nkxyz(2),nkxyz(3),
     .              ntet,qb)
        call rlse(oidtet)
        call defi(oidtet,ntet*5)
      endif
       call cpustp('PRWK',f,f)
! --- Calculate localized structure constants ---
      call defdr(orwats,-nclass)
      nderiv=1
      nsmesh=1
      lpot=.false.
      lsymst=.true.
      lsdm1 =.false.
      nowrit=.false.
      call strrs(w(oadot0),alat,w(oalph0),avw,w(obas),clabl,dele,deltr,
     .         switch(16),ialpha,icalc,w(oiclas),itrans,kap2,lmaxw,lmx,
     .        lpot,lsdm1,lsymst,nbas,nclass,nderiv,ndimin,neighm,nl,
     .         nowrit,nrclas,nsmesh,oadot,oadoto,oalpha,oalpho,ocg,
     .         ocy,oiax,oindcg,ojcg,onpr,os,osdot,osdm1,osdm1s,
     .         otrad,otral,plat,rmaxs,w(orwats),w(osigma),wsr)
! --- Generate Madelung matrix ---
      call defdr(odmad,nbas*nbas)
      call defdr(odlat,3*nkdmx)
      call defdr(oglat,3*nkdmx)
      call lattc(alat,as,awald,w(odlat),w(oglat),nkd,nkdmx,nkg,
     .           plat,tolews)
      call madmat(alat,awald,w(obas),w(odlat),w(odmad),w(oglat),nbas,
     .            nkd,nkg,plat)
      call rlse(odlat)
      if (cquit.eq.'MAD ') call fexit(0,' LM: Q=MAD encountered')
      call cpusrt('ITET',t)
      do nit=0,niter
        begmom = switch(4).or.(nit.gt.0)
        write(cnit,'(a1,i3)')'I',nit
        call cpusrt(cnit,t)
! ----  This is the beginning of zeroth iteration ---
! ----- Calculate the charges and ves for coming iteration. Necessary
! ----- that this precede ATOMSC since it requires electrostatic shifts
        call getq(clabl,dq,konfig,lmx,nclass,nl,nrclas,nsp,w(opnu),qc,
     .            w(oqnu),qtot,sumdq,z,zval)
        call madpot(clabl,w(odmad),w(oiclas),nbas,nclass,qtot,ves,wsr)
        call atomsc(begmom,clabl,switch(14),detot,dq,w(oenu),
     .              etot,w(oevrl),kap2,konfig,switch(17).or.switch(8),
     .              switch(29),switch(21),switch(6),switch(18),lmx,
     .              switch(28),switch(1),mmixat,nclass,nitat,nl,
     .              switch(9),nrclas,nrxc,nsp,w(opnu),qc,w(oqnu),qtot,
     .              ves,vmtz,vrmax,wsr,z)
        if (begmom) call enutpp(w(oenu),w(oidmod),nclass,nl,npp,nsp,
     .                          w(opp))
        call atompp(avw,clabl,kap2,switch(6),lmx,switch(1),nclass,
     .              nl,nsp,w(opp),w(osigma),ves,vmtz,wsr,z)
        if (nit.ne.0) detot = etot - etoto
        etoto = etot
        if (begmom.and.iprint().ge.0)
     .    write(lunit(1),301)nit,niter,sumdq,etot,rms2,detot
        if (cquit.eq.'ATOM') call fexit(0,' LM: Q=ATOM encountered')
        if (switch(6)) call fexit(0,' LM: FREE=T encountered')
! ----- Check for convergence ---
        convrg=(rms2.lt.rmsdel.and.dabs(detot).lt.toleto)
        if (.not.convrg.and.nit.ne.niter) then
! ------- Make copy of old p's and q's and downfolding indices ---
          call defdr(opold,nl*nsp*nclass)
          call defdr(oqold,nl*nsp*nclass*3)
          call defi (oidn ,nl*nclass)
          call icopy(nl*nclass,w(oidxdn),1,w(oidn),1)
          call dcopy(nl*nsp*nclass  ,w(opnu),1,w(opold),1)
          call dcopy(nl*nsp*nclass*3,w(oqnu),1,w(oqold),1)
          if(switch(8)) call defdr(oemom,-3*nl**4*nsp*nbas)
*===============================================*
c    call writepp(w(opp),nl,nsp,nclass)
*===============================================*
! ------- Band calculation ---
          call bndasa(alat,avw,w(obas),w(obk),clabl,dele,
     .                ebmax,ebmin,efermi,emaxc,eminc,w(oemom),
     .                w(oevrl),w(oiclas),w(oidtet),w(oidn),kap2,
     .                ldim,lmx,nbas,nclass,neighm,nkp,nkxyz,nl,
     .                norder,npts,nrclas,nsp,nsymop,ntet,oadot,oalpha,
     .                odmat,oeband,oiax,onpr,os,osdot,plat,w(opp),
     .                w(oqnu),range,w(osigma),switch,symopm,symopv,
     .                tolef,vmtz,width,wsr,w(owtkp),zval)
          if (switch(23)) then
            call bzlabs(lab1,lab2,nline,q1,q2)
            call wrlmfs(alat,w(obk),w(oeband),efermi,w(oibk),isym(1),
     .                  ldim,nkxyz(1),nkxyz(2),nkxyz(3),nkp,nsp,plat)
          endif
          if (switch(17))
     .      call dens(alat,avw,bound,clabl,ifmt3d,w(oirr),kap2,
     .                switch(29),switch(21),lmx,switch(24),nbas,nclass,
     .                ndel,neighm,nl,nsp,obas,odmat,oiax,oiclas,oidn,
     .                onpr,opp,origin,os,plat,wsr,z)
          if (switch(8))
     .      call rhofit(alat,avw,w(obas),bound,clabl,deltr,w(oemom),
     .                  w(oiclas),w(oidn),ifmt3d,w(oirr),kfit,lmaxw,
     .                  lmx,lmxrho,nbas,nclass,ndel,nl,nlrho,nrclas,
     .                  nsp,origin,plat,rfout,rmaxs2,w(orwats),
     .                  w(osigm2),wsr,z,zval)
! ------- Transform to the orthogonal representation
          call ortrep(clabl,w(oidn),nclass,nl,nsp,w(opp),w(oqnu))
!-------- Shift enu to center of gravity and calculate new moments
!-------- and backup old qnu's and pnu's
          call enutcg(avw,clabl,w(oevrl),w(oidmod),kap2,lmx,
     .                nclass,nl,nsp,w(opnu),w(opp),w(oqnu),wsr)
! ------- Mixing moments and pnu's
          call mixpq(beta,switch(13),clabl,lmx,mmixpq,nclass,nl,nsp,
     .               w(opnu),w(opold),w(oqnu),w(oqold),rms2,wc)
          call rlse(opold)
          if (cquit.eq.'/home/rudra/Recursion/ASR/BAND_A') call fexit(0,' LM: Q=BAND encountered')
        endif
        call cpustp(cnit,f,f)
        if (convrg.or.switch(17).or.switch(8)) goto 1
      enddo
1     call cpustp('ITET',f,f)
! --- End of the self-consistency loop

        if(convrg) call atompp1(avw,clabl,kap2,switch(6),lmx,
     .              switch(1),nclass,
     .              nl,nsp,w(opp),w(osigma),ves,vmtz,wsr,z)

      if (convrg) write(lunit(1),302) rms2,dabs(detot)
      emax=(idnint(ebmax*10)+1)*0.1d0
      emin=(idnint(ebmin*10)-1)*0.1d0
 301  format(/72('#')//'ITER',i3,' OUT OF',i3,':   MAG MOM=',f12.8,
     .             3x,'ETOT=',f17.8,/22x,'RMS DQ=',f12.8,
     .            2x,'DETOT=',f17.8//72('#')/)
 302  format(/' Jolly good show!  You converged to rms DQ=',f10.8,
     .       /' Jolly good show!  You converged to  DETOT=',f10.8)
        endif
      enddo
      call timres(clabl,nbas,nclass,nit)
!    write(*,*) "IMHERE"
      call fexit(0,prgnam)
      End
[/cpp]

What I am tryong to do is make another variable name for each of the variable present here in the first half of the code, add them as a argument etc. like

[cpp]	d_ngen=ngen
	d_nit=nit
	d_nitat=nitat
	d_niter=niter
    d_ialpha=ialpha
[/cpp]

and so on and put these as argument on the subroutine!!
what I am confused is is there any other, better way of transfering this value, that are got updated at the end of the 1st subroutine to the second one as well as to the calling program?

jimdempseyatthecove · ‎11-14-2009

Roddur

1) You may be able to declare those arguments to the subroutine as pass by value instead of pass by reference.

subroutine foo(ngen, nit, nitat, niter, ...)
integer, value :: ngen, nit, nitat, niter, ...

2) Place all those arguments into a user defined type in your main subroutine. Then declare a second instance of that type and perform the copy of that type and then pass the reference to the subroutine

type yourContext
sequence
integer, value :: ngen, nit, nitat, niter, ...
end type yourContext

type(yourContext) :: mine,nextLevel
! initialize mine for use in current subroutine
mine%ngen = ...
mine%nit = ...
...
!use mine in current subroutine
...
nextLevel = mine ! make copyfor next level
call yourNextLevel(nextLevel)

-----------------------

subroutine yourNextLevel(mine)
! ^^ use minehere
! declare the referenct to incomming context
type(yourContext) :: mine
! declare local copy for next level
type(yourContext) :: theirs

...
!use mine in current subroutine
...
nextLevel = mine ! make copyfor next level
call yourNextLevel(nextLevel) ! could be recursive or other subroutine

Jim Dempsey

jimdempseyatthecove · ‎11-14-2009

You may need to declare an interface too

Jim

roddur · ‎11-14-2009

Dear Jim,
As your first option, I cannot do that because i have to pass many arrays as well!
I mean, how can I pass the arrays by argument?

TimP · ‎11-14-2009

Why not a module for the shared data?

roddur · ‎11-15-2009

Quoting - tim18

Why not a module for the shared data?

Ummm....I am not sure, but what in my mind for the module is that, if I use the module for declaring variables, it will reinitialize the variable every time. But my problem is I have to carry forth the value of the variables generated by the first part.

Is there any way to check via some compiler option which values are modified?
If there is only one code, I can do that using idb and print the value. But how here?

jimdempseyatthecove · ‎11-15-2009

Roddur,

Your original post indicated that you needed to make copies of the variables (as opposed to passing references to the variables).

Do you require copies of the arrays, or can the called subroutine simply use a reference to the array?

If reference for the array but copies for the scalars, then in

type yourType
sequence
integer : nat
...
real, pointer :: Array(:)
...
end type yourType

real :: someArray(5000)

type(yourType) :: mine, nextLevel
...
! do once
mine%Array => someArray

...
! your code can modifie either
someArray(I) = xxx
! or
mine%Array(I) = xxx
! update your level's copy of local data
mine%nat = n

! make copy of your level's data for next level
! note, this makes a copy of the pointer to Array and not the data of Array
nextLevel = mine
call yourSub(nextLevel)

-------------------------
subroutine yourSub(mine)
type(yourType) :: mine
! and if this Level calls self or next level
type(yourType) :: nextLevel

Jim Dempsey