mycompsub_cmplex.f

parallel implementation of fft

c**********************************************************************c
c***************** form scalar product of 2 complex arrays ************c

        subroutine myscalconjin (c1,ps,nvecs)

        integer*4 nvecs,i
        real*8    c1(2*nvecs),ps(2)

        ps(1) = 0.0
        ps(2) = 0.0

c dir$l -novector

        do 100 i = 1, nvecs
           ps(1) = ps(1) + c1(2*i-1)*c1(2*i-1) + c1(2*i)*c1(2*i)
 100    continue

        return
        end

c**********************************************************************c
c***************** form scalar product of 2 complex arrays ************c

        subroutine mydiffrscal (c1,c2,tmpreal,nvecs)

        integer*4 nvecs,i
        real*8    c1(2*nvecs),c2(2*nvecs),tmpreal,tmpr,tmpc

        tmpreal = 0.0

c dir$l -novector

        do 100 i = 1, nvecs
	   tmpr    = c1(2*i-1)-c2(2*i-1)
	   tmpc    = c1(2*i)-c2(2*i)
           tmpreal = tmpreal + tmpr*tmpr + tmpc*tmpc 
 100    continue

        return
        end

c**********************************************************************c
c******* form the product of 2 complex arrays without conjugation *****c

        subroutine myscal (c1,c2,ps,nvecs)

        implicit  none
        integer*4 nvecs,i
        real*8    c1(2*nvecs),c2(2*nvecs),ps(2)

        ps(1) = 0.0
        ps(2) = 0.0

c dir$l -novector

        do 100 i = 1, nvecs
           ps(1) = ps(1) + c1(2*i-1)*c2(2*i-1) - c1(2*i)*c2(2*i)
           ps(2) = ps(2) + c1(2*i-1)*c2(2*i)   + c1(2*i)*c2(2*i-1)
 100    continue

        return
        end

c**********************************************************************c
c******************* normalize an array  ******************************c

        subroutine mynormal (c,nvecs,numall,fen_temp)

        implicit  none
        integer*4 nvecs,numall,i
        real*8    c(2*nvecs),p,fen_temp,x2,cmplx_scale

        call mypower (c,p,nvecs)

        p = p*fen_temp*fen_temp/float(numall)

        x2 = sqrt (p)
        cmplx_scale  = 1.0/(x2*x2)

        do 100 i = 1, nvecs
           c(2*i-1) = cmplx_scale*(c(2*i-1)*x2)
           c(2*i)   = cmplx_scale*(c(2*i)*x2)
 100    continue

        return
        end

c**********************************************************************c
c************** find the squared norm of an array *********************c
    
        subroutine mypower (c,pw,nvecs)

        implicit  none
C 7/05/01 by KSG the following include comm_info.h was added
        include 'comm_info.h'
        integer*4 nvecs,i
        real*8    c(2*nvecs),pw,pwtemp

        pw = 0.0

        do 100 i = 1, nvecs
           pw = pw + c(2*i-1)*c(2*i-1) + c(2*i)*c(2*i)
 100    continue

	call gdsum (pw,1,pwtemp)

        return
        end


c**********************************************************************c
c************** find the local squared norm of an array ***************c

        subroutine mypowerns (c,pw,nvecs)

        implicit  none
        integer*4 nvecs,i
C The following line was changed by KSG 8/29/01
C Because pwtemp is not used (it is no longer declared)
C        real*8    c(2*nvecs),pw,pwtemp
        real*8    c(2*nvecs),pw

        pw = 0.0

        do 100 i = 1, nvecs
           pw = pw + c(2*i-1)*c(2*i-1) + c(2*i)*c(2*i)
 100    continue

        return
        end

c**********************************************************************c
c************** find the squared norm of an array *********************c

        subroutine mypower2 (c1,pw1,c2,pw2,nvecs,adjust)

        implicit  none
C 7/05/01 by KSG the following include comm_info.h was added
        include "comm_info.h"
        integer*4 nvecs,i
        real*8    c1(2*nvecs),c2(2*nvecs),pw1,pw2
	real*8    tmp(2),work(2),adjust

	tmp(1) = 0.0
	tmp(2) = 0.0

        do 100 i = 1, nvecs
           tmp(1) = tmp(1) + c1(2*i-1)*c1(2*i-1) + c1(2*i)*c1(2*i)
 100    continue

        do 200 i = 1, nvecs
           tmp(2) = tmp(2) + c2(2*i-1)*c2(2*i-1) + c2(2*i)*c2(2*i)
 200    continue

        call gdsum (tmp,2,work)

	pw1 = adjust*tmp(1)
	pw2 = adjust*tmp(2)

        return
        end

c**********************************************************************c
c************** find the squared norm of an array *********************c

        subroutine mypower3 (c1,pw1,c2,pw2,c3,pw3,nvecs,adjust)

        implicit  none
C 7/05/01 by KSG the following include comm_info.h was added
        include "comm_info.h"
        integer*4 nvecs,i
        real*8    c1(2*nvecs),c2(2*nvecs),c3(2*nvecs),pw1,pw2,pw3
        real*8    tmp(3),work(3),adjust

        tmp(1) = 0.0
        tmp(2) = 0.0
	tmp(3) = 0.0

        do 100 i = 1, nvecs
           tmp(1) = tmp(1) + c1(2*i-1)*c1(2*i-1) + c1(2*i)*c1(2*i)
 100    continue

        do 200 i = 1, nvecs
           tmp(2) = tmp(2) + c2(2*i-1)*c2(2*i-1) + c2(2*i)*c2(2*i)
 200    continue

        do 300 i = 1, nvecs
           tmp(3) = tmp(3) + c3(2*i-1)*c3(2*i-1) + c3(2*i)*c3(2*i)
 300    continue

        call gdsum (tmp,3,work)

        pw1 = adjust*tmp(1)
        pw2 = adjust*tmp(2)
        pw3 = adjust*tmp(3)

        return
        end

c**********************************************************************c
c**********************************************************************c
c*** compute rvec(i) = 1.0 - real(
c*** cvec1(i)*conj(cvec1(i))+cvec2(i)*conj(cvec2(i)))

	subroutine myreal (rvec,cvec1,cvec2,nvecs)

        implicit  none
	integer*4 nvecs,i
	real*8    rvec(nvecs),cvec1(2*nvecs),cvec2(2*nvecs)

	do 100 i = 1, nvecs
	   rvec(i) = 1.0-
     &        (cvec1(2*i-1)*cvec1(2*i-1)+cvec1(2*i)*cvec1(2*i))-
     &        (cvec2(2*i-1)*cvec2(2*i-1)+cvec2(2*i)*cvec2(2*i)) 
 100    continue

	return
	end 

c**********************************************************************c
c**********************************************************************c
c*** compute rvec(i) = real(cvec1(i)*conj(cvec1(i))

        subroutine dreal_conj (rvec,cvec1,nvecs)

        implicit  none
        integer*4 nvecs,i
        real*8    rvec(nvecs),cvec1(2*nvecs)

        do 100 i = 1, nvecs
           rvec(i) = cvec1(2*i-1)*cvec1(2*i-1)+cvec1(2*i)*cvec1(2*i) 
 100    continue

        return
        end

c**********************************************************************c
c**********************************************************************c
c*** function mycdexp computes and returns a complex number mycdexp with
c*** mycdexp = dexp(zz) = dexp (x,y) = exp(x)*(dcos(y),dsin(y))

	function mycdexp (x,y)

	implicit   none
	complex*16 mycdexp
	real*8     exptm,x,y

	exptm   = dexp(x)
	mycdexp = dcmplx(exptm*dcos(y),exptm*dsin(y))

	return
	end

c**********************************************************************c
c**********************************************************************c
c*** function mycdexp_cmplx computes and returns a complex number mycdexp
c*** where mycdexp = cdexp(zz) = cdexp (0.0,imag(zz))
c*** mycdexp = (dcos(imag(zz),dsin(imag(zz))

        function mycdexp_cmplx (y)

	implicit   none
        complex*16 mycdexp_cmplx
        real*8     y 

        mycdexp_cmplx = dcmplx (dcos(y),dsin(y))

        return
        end

c**********************************************************************c
c**********************************************************************c
c*** function mycdexp_real computes and returns a complex number mycdexp 
c*** where mycdexp = cdexp(zz) = cdexp (real(zz),0.0)
c*** mycdexp = (dexp(real(zz),0.0)

        function mycdexp_real (x)
	
	implicit   none
        complex*16 mycdexp_real
	real*8     x

        mycdexp_real = dcmplx (dexp(x),0.0)

        return
        end

c**********************************************************************c
c**********************************************************************c
c*** compute cvec(i) = exp(z(i)) where z(i)=(0.0,cmplx_imag(i))

	subroutine cdexp_cmplx (cvec,cmplx_imag,nvecs)

	implicit  none
	integer*4 nvecs,i
	real*8    cvec(2*nvecs),cmplx_imag(nvecs),exptmp

        exptmp = 1.0

	do 100 i = 1, nvecs
           cvec(2*i-1) = dcos(cmplx_imag(i))
	   cvec(2*i)   = dsin(cmplx_imag(i))
 100    continue

        return
        end

c**********************************************************************c
c**********************************************************************c
c*** compute cvec(i) = exp(z(i)) where z(i)=(cmplx_real(i),0.0)

        subroutine cdexp_real (cvec,cmplx_real,nvecs)

        implicit  none
        integer*4 nvecs,i
        real*8    cvec(2*nvecs),cmplx_real(nvecs)

        do 100 i = 1, nvecs
           cvec(2*i-1) = dexp (cmplx_real(i)) 
           cvec(2*i)   = 0.0
 100    continue

        return
        end

c**********************************************************************c
c**********************************************************************c
c*** compute cvec(i) = exp(z(i)) 

        subroutine cdexp_all (cvec,zz,nvecs)

        implicit  none
        integer*4 nvecs,i
        real*8    cvec(2*nvecs),zz(2*nvecs),exptmp

        do 100 i = 1, nvecs
	   exptmp      = dexp(zz(2*i-1)) 
           cvec(2*i-1) = exptmp*dcos(zz(2*i))
           cvec(2*i)   = exptmp*dsin(zz(2*i))
 100    continue

        return
        end

c**********************************************************************c
c**********************************************************************c
c*** compute mydreal = zz*conj(zz) 

        function mydreal (zz)

        implicit none
        real*8   mydreal,zz(2)

	mydreal = zz(1)*zz(1)+zz(2)*zz(2)

        return
        end

c**********************************************************************c
c**********************************************************************c
c*** compute mydreal = zz*conj(zz)

        function mydiv (zz1,zz2)

        implicit   none
        complex*16 mydiv,zz1,zz2
	real*8     cmterm,x1,y1,x2,y2

	x1 = dreal(zz1)
	y1 = dimag(zz1)

        x2 = dreal(zz2)
        y2 = dimag(zz2) 

	cmterm = 1.0/(x2*x2+y2*y2)
	mydiv = dcmplx (cmterm*(x1*x2+y1*y2),cmterm*(x2*y1-x1*y2))

	return
	end

c**********************************************************************c
c**********************************************************************c

C******************* find the s.d. of two arrays **********************

        subroutine myecart3 (c1,c2,sigma1,c3,c4,sigma2,c5,c6,sigma3)

        include 'timeheader.h'
C 7/05/01 by KSG the following include comm_info.h was added
        include "comm_info.h"
	complex*16 c1(n,nz),c2(n,nz),c3(n,nz),c4(n,nz)
	complex*16 c5(n,nz),c6(n,nz)

	real*8     fen,adjust
	real*8     sigma1,sigma2,sigma3,tmpr(3),work(3)

	common/norm/fen

        call mydiffrscal (c1,c2,tmpr(1),nmnz)
        call mydiffrscal (c3,c4,tmpr(2),nmnz)
        call mydiffrscal (c5,c6,tmpr(3),nmnz)

        call gdsum (tmpr,3,work)

	adjust = (FEN*FEN)/DFLOAT(nsq)

	sigma1 = dsqrt (adjust*tmpr(1))
        sigma2 = dsqrt (adjust*tmpr(2))
        sigma3 = dsqrt (adjust*tmpr(3))

	return
	end

c**********************************************************************c
c**********************************************************************c

	subroutine adjustall (ctemp,ntemp)

	include 'timeheader.h'
C 7/05/01 by KSG the following include comm_info.h was added
        include "comm_info.h"
        integer*4 i,ntemp
	real*8    ctemp(2*ntemp),work(24)
        real*8    fen,adjust

        common/norm/fen

	adjust = (FEN*FEN)/DFLOAT(nsq)

	call gdsum (ctemp,2*ntemp,work)

	do 100 i = 1, 2*ntemp
	   ctemp(i) = adjust*ctemp(i)
 100    continue

	return
	end

c**********************************************************************c
c**********************************************************************c

        subroutine myscal2 (c1,c2,ps1,c3,c4,ps2,nvecs,adjust)
C 7/05/01 by KSG the following include comm_info.h was added
        include "comm_info.h"
        integer*4  nvecs
	real*8     adjust
	complex*16 ps1,ps2,work(2),tmp(2)
	complex*16 c1(nvecs),c2(nvecs),c3(nvecs),c4(nvecs) 
	
        call myscalconj (c1,c2,tmp(1),nvecs)
	call myscalconj (c3,c4,tmp(2),nvecs)

        call gdsum (tmp,4,work)

        ps1 = adjust*tmp(1)
	ps2 = adjust*tmp(2)

	return
	end

c**********************************************************************c
c**********************************************************************c