Program f : dgauss, gauss, gauss8, ordleg, ordleg8
program progf
implicit none
real roots(40)
real racp(20),pg(20),sia(20),rad(20),pgssin2(20),aerr,value
real sinm1(20),sinm2(20),sin2(20),maxerr,sum
integer i,nberr
real*8 value8,aerr8
real*8 racp8(20),pg8(20),sia8(20),rad8(20),pgssin28(20)
real*8 sinm18(20),sinm28(20),sin28(20),sum8
# first time with real*4
print *,'The roots are computed with dgauss and with gauss'
print *,'for comparison (real*4 are used)'
call dgauss(40,roots,0)
call gauss(20,racp,pg,sia,rad,pgssin2,sinm1,sinm2,sin2)
nberr=0
maxerr=0
aerr=0
do i=1,20
maxerr=max(maxerr,abs(1-roots(i)/racp(i)))
call ordleg(value,roots(i),40)
aerr=aerr+abs(value)
if(abs(1-roots(i)/racp(i)).gt.2e-6) then
nberr=nberr+1
endif
enddo
if (nberr.eq.0) then
print *,'The roots of the Legendre polynomial found by'
print *,' dgauss and by gauss are the same.'
else
print *,'There was',nberr,' difference(s) between the roots'
print *,' found by dgauss and by gauss.'
endif
print *,'The maximum relative error is:',maxerr
print *,'The sum of the absolute values of the'
print *,' Legendre polynomial at the roots is:',aerr
sum=0
do i=1,20
sum=sum+pg(i)
enddo
print *,'Sum of the Gauss weights is:',sum,'(abs. err.:',
* abs(sum-1),')'
print *,''
# this time with real*8
print *,'Now the roots will be computed by gauss8 using real*8'
call gauss8(20,racp8,pg8,sia8,rad8,pgssin28,sinm18,sinm28,sin28)
aerr8=0
sum8=0
do i=1,20
call ordleg8(value8,racp8(i),40)
aerr8=aerr8+abs(value8)
sum8=sum8+pg8(i)
enddo
print *,'Sum of the Gauss weights is:',sum8,'(abs. err.:',
* abs(sum8-1),')'
print *,'The sum of the absolute values of the'
print *,' Legendre polynomial at the roots is:',aerr8
stop
end