Commit 9b9473ce by Daniel Wortmann

### Merge branch 'aixcel' of https://iffgit.fz-juelich.de/fleur/fleur into aixcel

```Conflicts:
eigen/hsmt_sph.F90```
parents bb1a6d24 3d41caec
 ... ... @@ -70,15 +70,15 @@ CONTAINS !\$OMP MASTER IF ((atoms%invsat(na).EQ.0) .OR. (atoms%invsat(na).EQ.1)) THEN IF ((atoms%invsat(na) == 0) .OR. (atoms%invsat(na) == 1)) THEN !---> if this atom is the first of two atoms related by inversion, !---> the contributions to the overlap matrix of both atoms are added !---> at once. where it is made use of the fact, that the sum of !---> these contributions is twice the real part of the contribution !---> of each atom. note, that in this case there are twice as many !---> (2*(2*l+1)) k-vectors (compare abccoflo and comments there). IF (atoms%invsat(na).EQ.0) invsfct = 1 IF (atoms%invsat(na).EQ.1) invsfct = 2 IF (atoms%invsat(na) == 0) invsfct = 1 IF (atoms%invsat(na) == 1) invsfct = 2 ! DO lo = 1,atoms%nlo(ntyp) ... ... @@ -127,7 +127,7 @@ CONTAINS !+t3e DO nkvec = 1,invsfct* (2*l+1) locol= lapw%nv(jintsp)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix IF (MOD(locol-1,mpi%n_size).EQ.mpi%n_rank) THEN IF (MOD(locol-1,mpi%n_size) == mpi%n_rank) THEN locol=(locol-1)/mpi%n_size+1 !this is the column in local storage !-t3e IF (hmat%l_real) THEN ... ... @@ -174,7 +174,7 @@ CONTAINS !---> local orbitals at the same atom and with itself DO nkvec = 1,invsfct* (2*l+1) locol = lapw%nv(jintsp)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix IF (MOD(locol-1,mpi%n_size).EQ.mpi%n_rank) THEN IF (MOD(locol-1,mpi%n_size) == mpi%n_rank) THEN locol=(locol-1)/mpi%n_size+1 !this is the column in local storage !-t3e !---> calculate the hamiltonian matrix elements with other ... ...
 ... ... @@ -313,41 +313,41 @@ SUBROUTINE hsmt_sph_gpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_s END SUBROUTINE hsmt_sph_gpu #endif SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_shift,usdus,fj,gj,smat,hmat) USE m_constants, ONLY : fpi_const,tpi_const USE m_types IMPLICIT NONE TYPE(t_input),INTENT(IN) :: input TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_noco),INTENT(IN) :: noco TYPE(t_atoms),INTENT(IN) :: atoms TYPE(t_lapw),INTENT(IN) :: lapw TYPE(t_usdus),INTENT(IN) :: usdus CLASS(t_mat),INTENT(INOUT) :: smat,hmat ! .. ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: n,isp,iintsp,jintsp COMPLEX, INTENT(IN) :: chi ! .. ! .. Array Arguments .. REAL, INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,input%jspins) REAL, INTENT (IN) :: e_shift!(atoms%ntype,input%jspins) REAL, INTENT (IN) :: fj(:,0:,:),gj(:,0:,:) ! .. ! .. Local Scalars .. REAL tnn(3), elall,fct,fct2,fjkiln,gjkiln,ddnln,ski(3) REAL apw_lo1,apw_lo2,apw1,w1 SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_shift,usdus,fj,gj,smat,hmat) USE m_constants, ONLY : fpi_const,tpi_const USE m_types IMPLICIT NONE TYPE(t_input),INTENT(IN) :: input TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_noco),INTENT(IN) :: noco TYPE(t_atoms),INTENT(IN) :: atoms TYPE(t_lapw),INTENT(IN) :: lapw TYPE(t_usdus),INTENT(IN) :: usdus CLASS(t_mat),INTENT(INOUT) :: smat,hmat ! .. ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: n,isp,iintsp,jintsp COMPLEX, INTENT(IN) :: chi ! .. ! .. Array Arguments .. REAL, INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,input%jspins) REAL, INTENT (IN) :: e_shift!(atoms%ntype,input%jspins) REAL, INTENT (IN) :: fj(:,0:,:),gj(:,0:,:) ! .. ! .. Local Scalars .. REAL tnn(3), elall,fct,fct2,fjkiln,gjkiln,ddnln,ski(3) REAL apw_lo1,apw_lo2,apw1,w1 COMPLEX capw1 INTEGER kii,ki,kj,l,nn INTEGER kj_BlockSize, kj_BlockNum, kjb, kj_start, kj_end ! .. ! .. Local Arrays .. REAL fleg1(0:atoms%lmaxd),fleg2(0:atoms%lmaxd),fl2p1(0:atoms%lmaxd) REAL fl2p1bt(0:atoms%lmaxd) REAL qssbti(3),qssbtj(3) REAL, ALLOCATABLE :: gdot(:) REAL, ALLOCATABLE :: VecHelpS(:),VecHelpH(:) COMPLEX capw1 INTEGER kii,ki,kj,l,nn,kj_end ! .. ! .. Local Arrays .. REAL fleg1(0:atoms%lmaxd),fleg2(0:atoms%lmaxd),fl2p1(0:atoms%lmaxd) REAL fl2p1bt(0:atoms%lmaxd) REAL qssbti(3),qssbtj(3) REAL, ALLOCATABLE :: plegend(:,:) REAL, ALLOCATABLE :: VecHelpS(:),VecHelpH(:) COMPLEX, ALLOCATABLE :: cph(:),cp_h(:,:) LOGICAL apw(0:atoms%lmaxd) ... ... @@ -364,112 +364,118 @@ SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_s CALL lapw%phase_factors(1,atoms%taual(:,nn),noco%qss,cp_h(:,nn)) ENDDO !\$OMP PARALLEL DEFAULT(NONE)& !\$OMP SHARED(lapw,atoms,noco,mpi,input,usdus,smat,hmat,cp_h)& !\$OMP SHARED(jintsp,iintsp,n,fleg1,fleg2,fj,gj,isp,fl2p1,el,e_shift,fl2p1bt,chi)& !\$OMP PRIVATE(kii,ki,ski,kj,plegend,l,kj_end,qssbti,qssbtj,fct2)& !\$OMP PRIVATE(cph,nn,tnn,fjkiln,gjkiln)& !\$OMP PRIVATE(w1,apw_lo1,apw_lo2,ddnln,elall,fct,apw1)& !\$OMP PRIVATE(capw1,VecHelpS,VecHelpH) ALLOCATE(cph(MAXVAL(lapw%nv))) ALLOCATE(plegend(MAXVAL(lapw%nv),0:atoms%lmaxd)) ALLOCATE(VecHelpS(MAXVAL(lapw%nv)),VecHelpH(MAXVAL(lapw%nv))) plegend=0.0 plegend(:,0)=1.0 qssbti=MERGE(- noco%qss/2,+ noco%qss/2,jintsp.EQ.1) qssbtj=MERGE(- noco%qss/2,+ noco%qss/2,iintsp.EQ.1) !\$OMP DO SCHEDULE(DYNAMIC,1) DO ki = mpi%n_rank+1, lapw%nv(jintsp), mpi%n_size kii=(ki-1)/mpi%n_size+1 ski = lapw%gvec(:,ki,jintsp) + qssbti !---> legendre polynomials DO kj = 1,ki plegend(kj,1) = DOT_PRODUCT(lapw%gk(:,kj,iintsp),lapw%gk(:,ki,jintsp)) END DO DO l = 1,atoms%lmax(n) - 1 plegend(:ki,l+1) = fleg1(l)*plegend(:ki,1)*plegend(:ki,l) - fleg2(l)*plegend(:ki,l-1) END DO !---> set up phase factors cph = 0.0 DO nn = SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n)) tnn = tpi_const*atoms%taual(:,nn) DO kj = 1,ki cph(kj) = cph(kj) +& cp_h(ki,nn)/cp_h(kj,nn) ! CMPLX(COS(DOT_PRODUCT(ski-lapw%gvec(:,kj,iintsp)-qssbtj,tnn)),& ! SIN(DOT_PRODUCT(lapw%gvec(:,kj,iintsp)+qssbtj-ski,tnn))) ! IF (iintsp.NE.jintsp) cph(kj)=CONJG(cph(kj)) END DO END DO !---> update overlap and l-diagonal hamiltonian matrix kj_end = MIN(ki,lapw%nv(iintsp)) VecHelpS = 0.d0 VecHelpH = 0.d0 DO l = 0,atoms%lmax(n) fjkiln = fj(ki,l,jintsp) gjkiln = gj(ki,l,jintsp) ! IF (input%l_useapw) THEN w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + & usdus%us(l,n,isp)*usdus%duds(l,n,isp) ) apw_lo1 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( gjkiln * w1 +& fjkiln * usdus%us(l,n,isp) * usdus%dus(l,n,isp) ) apw_lo2 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( fjkiln * w1 +& gjkiln * usdus%uds(l,n,isp) * usdus%duds(l,n,isp) ) ! ENDIF ddnln = usdus%ddn(l,n,isp) elall = el(l,n,isp) IF (l<=atoms%lnonsph(n)) elall=elall-e_shift!(isp) IF (smat%l_real) THEN DO kj = 1,ki fct = plegend(kj,l)*fl2p1(l)*& ( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln ) fct2 = plegend(kj,l)*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) ) smat%data_r(kj,kii)=smat%data_r(kj,kii)+REAL(cph(kj))*fct hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + REAL(cph(kj)) * ( fct * elall + fct2) !+APW IF (input%l_useapw) THEN apw1 = REAL(cph(kj)) * plegend(kj,l) * & kj_BlockSize = 1000 !\$OMP PARALLEL DEFAULT(NONE)& !\$OMP SHARED(cp_h,kj_BlockSize,lapw,atoms,noco,mpi,input,usdus,smat,hmat)& !\$OMP SHARED(jintsp,iintsp,n,fleg1,fleg2,fj,gj,isp,fl2p1,el,e_shift,fl2p1bt,chi)& !\$OMP PRIVATE(kii,ki,ski,kj,gdot,l,qssbti,qssbtj,fct2,plegend)& !\$OMP PRIVATE(cph,nn,tnn,fjkiln,gjkiln)& !\$OMP PRIVATE(w1,apw_lo1,apw_lo2,ddnln,elall,fct,apw1)& !\$OMP PRIVATE(capw1,VecHelpS,VecHelpH)& !\$OMP PRIVATE(kj_start,kj_end,kj_BlockNum, kjb) ALLOCATE(cph(MAXVAL(lapw%nv))) ALLOCATE(gdot(0:atoms%lmaxd)) ALLOCATE(VecHelpS(MAXVAL(lapw%nv)),VecHelpH(MAXVAL(lapw%nv))) gdot=0.0 qssbti=MERGE(- noco%qss/2,+ noco%qss/2,jintsp.EQ.1) qssbtj=MERGE(- noco%qss/2,+ noco%qss/2,iintsp.EQ.1) !\$OMP DO SCHEDULE(DYNAMIC,1) DO ki = mpi%n_rank+1, lapw%nv(jintsp), mpi%n_size kii=(ki-1)/mpi%n_size+1 ski = lapw%gvec(:,ki,jintsp) + qssbti kj_BlockNum = ki/kj_BlockSize + 1 DO kjb = 1, kj_BlockNum kj_start = (kjb-1) * kj_BlockSize + 1 kj_end = MIN(ki,kjb * kj_BlockSize,lapw%nv(iintsp)) !---> legendre polynomials DO kj = kj_start,kj_end gdot(kj,1) = DOT_PRODUCT(lapw%gk(:,kj,iintsp),lapw%gk(:,ki,jintsp)) END DO !---> set up phase factors cph(kj_start:kj_end) = 0.0 DO nn = SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n)) tnn = tpi_const*atoms%taual(:,nn) DO kj = kj_start,kj_end cph(kj) = cph(kj) +cp_h(ki,nn)/cp_h(kj,nn) ! CMPLX(COS(DOT_PRODUCT(ski-lapw%gvec(:,kj,iintsp)-qssbtj,tnn)),& ! SIN(DOT_PRODUCT(lapw%gvec(:,kj,iintsp)+qssbtj-ski,tnn))) ! IF (iintsp.NE.jintsp) cph(kj)=CONJG(cph(kj)) END DO END DO !---> update overlap and l-diagonal hamiltonian matrix VecHelpS = 0.d0 VecHelpH = 0.d0 DO l = 0,atoms%lmax(n) fjkiln = fj(ki,l,jintsp) gjkiln = gj(ki,l,jintsp) IF (input%l_useapw) THEN w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + & usdus%us(l,n,isp)*usdus%duds(l,n,isp) ) apw_lo1 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( gjkiln * w1 +& fjkiln * usdus%us(l,n,isp) * usdus%dus(l,n,isp) ) apw_lo2 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( fjkiln * w1 +& gjkiln * usdus%uds(l,n,isp) * usdus%duds(l,n,isp) ) ENDIF ddnln = usdus%ddn(l,n,isp) elall = el(l,n,isp) IF (l<=atoms%lnonsph(n)) elall=elall-e_shift!(isp) IF (smat%l_real) THEN DO kj = kj_start,kj_end plegend = LegendrePoly(l,gdot(kj)) fct = plegend*fl2p1(l)*& ( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln ) fct2 = plegend*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) ) smat%data_r(kj,kii)=smat%data_r(kj,kii)+REAL(cph(kj))*fct hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + REAL(cph(kj)) * ( fct * elall + fct2) !+APW IF (input%l_useapw) THEN apw1 = REAL(cph(kj)) * plegend * & ( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) ) hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + apw1 ENDIF !-APW ENDDO ELSE DO kj = 1,kj_end fct = plegend(kj,l)*fl2p1(l)*& hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + apw1 ENDIF !-APW ENDDO ELSE plegend = LegendrePoly(l,gdot(kj)) DO kj = kj_start,kj_end fct = plegend*fl2p1(l)*& ( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln ) fct2 = plegend(kj,l)*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) ) VecHelpS(kj) = VecHelpS(kj) + fct VecHelpH(kj) = VecHelpH(kj) + fct*elall + fct2 IF (input%l_useapw) THEN capw1 = cph(kj)*plegend(kj,l)& * ( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) ) hmat%data_c(kj,kii)=hmat%data_c(kj,kii) + capw1 ENDIF END DO ENDIF !---> end loop over l ENDDO IF (.not.smat%l_real) THEN smat%data_c(:kj_end,kii)=smat%data_c(:kj_end,kii) + chi*cph(:kj_end) * VecHelpS(:kj_end) hmat%data_c(:kj_end,kii)=hmat%data_c(:kj_end,kii) + chi*cph(:kj_end) * VecHelpH(:kj_end) ENDIF !---> end loop over ki ENDDO !\$OMP END DO DEALLOCATE(plegend) DEALLOCATE(cph) DEALLOCATE(VecHelpS,VecHelpH) !\$OMP END PARALLEL CALL timestop("spherical setup") RETURN END SUBROUTINE hsmt_sph_cpu fct2 = plegend*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) ) VecHelpS(kj) = VecHelpS(kj) + fct VecHelpH(kj) = VecHelpH(kj) + fct*elall + fct2 IF (input%l_useapw) THEN capw1 = cph(kj)*plegend *& ( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) ) hmat%data_c(kj,kii)=hmat%data_c(kj,kii) + capw1 ENDIF END DO ENDIF !---> end loop over l ENDDO IF (.not.smat%l_real) THEN smat%data_c(kj_start:kj_end,kii)=smat%data_c(kj_start:kj_end,kii) + chi*cph(kj_start:kj_end) * VecHelpS(kj_start:kj_end) hmat%data_c(kj_start:kj_end,kii)=hmat%data_c(kj_start:kj_end,kii) + chi*cph(kj_start:kj_end) * VecHelpH(kj_start:kj_end) ENDIF !---> end loop over kj-blocks ENDDO !---> end loop over ki ENDDO !\$OMP END DO DEALLOCATE(gdot) DEALLOCATE(cph) DEALLOCATE(VecHelpS,VecHelpH) !\$OMP END PARALLEL CALL timestop("spherical setup") RETURN END SUBROUTINE hsmt_sph_cpu END MODULE m_hsmt_sph
 ... ... @@ -58,15 +58,15 @@ CONTAINS CALL lapw%phase_factors(i,atoms%taual(:,na),noco%qss,cph(:,i)) ENDDO IF ((atoms%invsat(na).EQ.0) .OR. (atoms%invsat(na).EQ.1)) THEN IF ((atoms%invsat(na) == 0) .OR. (atoms%invsat(na) == 1)) THEN !---> if this atom is the first of two atoms related by inversion, !---> the contributions to the overlap matrix of both atoms are added !---> at once. where it is made use of the fact, that the sum of !---> these contributions is twice the real part of the contribution !---> of each atom. note, that in this case there are twice as many !---> (2*(2*l+1)) k-vectors (compare abccoflo and comments there). IF (atoms%invsat(na).EQ.0) invsfct = 1 IF (atoms%invsat(na).EQ.1) invsfct = 2 IF (atoms%invsat(na) == 0) invsfct = 1 IF (atoms%invsat(na) == 1) invsfct = 2 con = fpi_const/SQRT(cell%omtil)* ((atoms%rmt(ntyp))**2)/2.0 ... ... @@ -82,7 +82,7 @@ CONTAINS DO nkvec = 1,invsfct* (2*l+1) !Each LO can have several functions !+t3e locol = lapw%nv(jintsp)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix IF (MOD(locol-1,mpi%n_size).EQ.mpi%n_rank) THEN IF (MOD(locol-1,mpi%n_size) == mpi%n_rank) THEN locol=(locol-1)/mpi%n_size+1 !this is the column in local storage !-t3e k = lapw%kvec(nkvec,lo,na) ... ... @@ -107,7 +107,7 @@ CONTAINS !---> orbitals at the same atom, if they have the same l DO lop = 1, (lo-1) lp = atoms%llo(lop,ntyp) IF (l.EQ.lp) THEN IF (l == lp) THEN fact3 = con**2 * fl2p1 * (& alo1(lop)*(alo1(lo) + & clo1(lo)*ud%uulon(lo,ntyp,isp))+& ... ...
 ... ... @@ -13,7 +13,7 @@ MODULE m_fleur_arguments CHARACTER(len=200) :: values END TYPE t_fleur_param INTEGER,PARAMETER:: no_params=6 INTEGER,PARAMETER:: no_params=7 TYPE(t_fleur_param) :: fleur_param(no_params)=(/& t_fleur_param(0,"-old","Generate input file for old fleur versions",""),& ... ... @@ -21,6 +21,7 @@ MODULE m_fleur_arguments t_fleur_param(0,"-genEnpara","generate an 'enpara' file",""),& t_fleur_param(0,"-electronConfig","explicitely write the electron configuration into inp.xml",""),& t_fleur_param(0,"-fast_defaults","generate more aggressive (and less stable) input parameters for faster calculations",""),& t_fleur_param(0,"-kpts_gw","add alternative k point set for GW",""),& t_fleur_param(0,"-h","print this help message","")& /) ... ...
 ... ... @@ -32,6 +32,7 @@ CONTAINS CALL print_argument("-explicit") CALL print_argument("-electronConfig") CALL print_argument("-fast_defaults") CALL print_argument("-kpts_gw") CALL print_argument("-h") WRITE(*,'(a)')"" WRITE(*,'(a)')"Please check the documentation on www.flapw.de for more details" ... ...
 ... ... @@ -279,6 +279,13 @@ SUBROUTINE w_inpXML(& 212 FORMAT(' ') WRITE (fileNum,212) kpts%kPointDensity(1),kpts%kPointDensity(2),kpts%kPointDensity(3),kptGamma END IF IF(juDFT_was_argument("-kpts_gw")) THEN WRITE(fileNum,'(a)') ' ' WRITE(fileNum,'(a)') ' ' WRITE(fileNum,'(a)') ' ' END IF WRITE (fileNum,'(a)') ' ' ! ... ...
 ... ... @@ -13,7 +13,7 @@ MODULE m_fleur_arguments CHARACTER(len=200) :: values END TYPE t_fleur_param INTEGER,PARAMETER:: no_params=22 INTEGER,PARAMETER:: no_params=23 TYPE(t_fleur_param) :: fleur_param(no_params)=(/& !Input options t_fleur_param(0,"-toXML","Convert an old 'inp' file into the new XML format",""),& ... ... @@ -60,6 +60,7 @@ MODULE m_fleur_arguments !Output t_fleur_param(0,"-no_out","Do not open the 'out' file but write to stdout",""),& t_fleur_param(0,"-genEnpara","Generate an 'enpara' file for the energy parameters",""),& t_fleur_param(0,"-kpts_gw","add alternative k point set for GW in all outputs for the XML input file",""),& t_fleur_param(0,"-h","Print this message",""),& t_fleur_param(0,"-no_send","Do not send usage data","")& !HDF density ... ...
 ... ... @@ -28,6 +28,7 @@ USE m_setupMPI USE m_cdn_io USE m_fleur_info USE m_broyd_io USE m_checks USE m_prpqfftmap USE m_writeOutHeader ... ... @@ -523,6 +524,7 @@ CALL results%init(dimension,input,atoms,kpts,noco) IF (mpi%irank.EQ.0) THEN IF(input%gw.NE.0) CALL resetBroydenHistory() CALL setStartingDensity(noco%l_noco) END IF ... ...
 ... ... @@ -32,6 +32,7 @@ math/DoubleFactorial.f90 math/ExpSave.f90 math/intgr.F90 math/ylm4.F90 math/LegendrePoly.f90 ) if (FLEUR_USE_FFTMKL) set(fleur_F90 \${fleur_F90} math/mkl_dfti.f90) ... ...
 module m_LegendrePoly implicit none contains pure function LegendrePoly(l,x) result(p) implicit none integer, intent(in) :: l real, intent(in) :: x(:) real :: p(size(x)) select case(l) case(0) p = 1.0 case(1) p = x case(2) p = x**2 - 1./3. case(3) p = x**3 - 3./5. * x case(4) p = x**4 - 6./7. * x**2 + 3./35. case(5) p = x**5 - 10./9. * x**3 + 5./21. * x case(6) p = x**6 - 15./11. * x**4 + 5./11. * x**2 - 5./231 case(7) p = x**7 - 21./13. * x**5 + 105./143. * x**3 - 35./429. * x case(8) p = x**8 - 28./15. * x**6 + 14./13. * x**4 - 28./143. * x**2 + 7./1287. case(9) p = x**9 - 36./17. * x**7 + 126./85. * x**5 - 84./221. * x**3 + 17./656. case(10) p = x**10 - 45./19. * x**8 + 630./323. * x**6 - 210./323. * x**4 + 106./1413. * x**2 - 1./733. end select end function LegendrePoly pure function LegendrePoly_scalar(l,x) result(p) implicit none integer, intent(in) :: l real, intent(in) :: x real :: p select case(l) case(0) p = 1.0 case(1) p = x case(2) p = x**2 - 1./3. case(3) p = x**3 - 3./5. * x case(4) p = x**4 - 6./7. * x**2 + 3./35. case(5) p = x**5 - 10./9. * x**3 + 5./21. * x case(6) p = x**6 - 15./11. * x**4 + 5./11. * x**2 - 5./231 case(7) p = x**7 - 21./13. * x**5 + 105./143. * x**3 - 35./429. * x case(8) p = x**8 - 28./15. * x**6 + 14./13. * x**4 - 28./143. * x**2 + 7./1287. case(9) p = x**9 - 36./17. * x**7 + 126./85. * x**5 - 84./221. * x**3 + 17./656. case(10) p = x**10 - 45./19. * x**8 + 630./323. * x**6 - 210./323. * x**4 + 106./1413. * x**2 - 1./733. end select end function LegendrePoly_scalar end module m_LegendrePoly
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!