I use FUNCTION qromb ( from numerical recipes )in my program which
calls polint subroutine.
------------------------------------------------------
SUBROUTINE polint(xa,ya,x,y,dy)
USE nrtype; USE nrutil, ONLY : assert_eq,iminloc,nrerror
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya
REAL(SP), INTENT(IN) :: x
REAL(SP), INTENT(OUT) :: y,dy
!Given arrays xa and ya of length N, and given a value x, this
routine returns a value y,
!and an error estimate dy. If P(x) is the polynomial of degree N
- 1 such that P(xai) =
!yai, i = 1, . . . ,N, then the returned value y = P(x).
INTEGER(I4B) :: m,n,ns
REAL(SP), DIMENSION(size(xa)) :: c,d,den,ho
n=assert_eq(size(xa),size(ya),'polint')
c=ya !Initialize the tableau of c's and d's.
d=ya
ho=xa-x
ns=iminloc(abs(x-xa)) !Find index ns of closest table entry.
y=ya(ns) !This is the initial approximation to y.
ns=ns-1
do m=1,n-1 !For each column of the tableau,
den(1:n-m)=ho(1:n-m)-ho(1+m:n) !we loop over the current c's and
d's and up
if(any(den(1:n-m) == 0.0)) & !date them.
call nrerror('polint: calculation failure')
!This error can occur only if two input xa's are (to within
roundoff) identical.
den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
d(1:n-m)=ho(1+m:n)*den(1:n-m) !Here the c's and d's are updated.
c(1:n-m)=ho(1:n-m)*den(1:n-m)
if (2*ns < n-m) then !After each column in the tableau is
completed, we decide
!which correction, c or d, we want to add to our accumulating
!value of y, i.e., which path to take through
!the tableau--forking up or down. We do this in such a
!way as to take the most "straight line" route through the
!tableau to its apex, updating ns accordingly to keep track
!of where we are. This route keeps the partial approximations
!centered (insofar as possible) on the target x. The
!last dy added is thus the error indication.
dy=c(ns+1)
else
dy=d(ns)
ns=ns-1
end if
y=y+dy
end do
END SUBROUTINE polint
---------------------------------------------------------------
!FUNCTION qromb(func,a,b)*******
real(kind=8) FUNCTION qromb(a,b,yy)
USE nrtype; USE nrutil, ONLY :nrerror
!USE nr, ONLY :polint,trapzd
USE nr, ONLY :polint
IMPLICIT NONE
REAL(kind=8), INTENT(IN) :: a,b
REAL(kind=8) :: qromb
REAL(kind=8) ::func
REAL(kind=8), INTENT(IN) ::yy
!INTERFACE
!FUNCTION func(x)
!USE nrtype
!REAL(SP), DIMENSION(:), INTENT(IN) :: x
!REAL(SP), DIMENSION(size(x)) :: func
!END FUNCTION func
!END INTERFACE
INTEGER(I4B), PARAMETER :: JMAX=20,JMAXP=JMAX+1,K=5,KM=K-1
REAL(kind=8), PARAMETER :: EPS=1.0e-6_sp
!Returns the integral of the function func from a to b.
Integration is performed by Romberg's
!method of order 2K, where, e.g., K=2 is Simpson's rule.
!Parameters: EPS is the fractional accuracy desired, as
determined by the extrapolation error
!estimate; JMAX limits the total number of steps; K is the
number of points used in the
!extrapolation.
REAL(kind=8), DIMENSION(JMAXP) :: h,s !These store the
successive trapezoidal approximations
!and
REAL(SP) :: dqromb !their relative stepsizes.
INTEGER(I4B) :: j
h(1)=1.0
do j=1,JMAX
!call trapzd(func,a,b,s(j),j)
! call trapzd(a,b,s(j),j,yy)************************
call NewTrapz(a,b,s(j),j,yy)*********************
if (j >= K) then
call polint(h(j-KM:j),s(j-KM:j),0.0_sp,qromb,dqromb)
if (abs(dqromb) <= EPS*abs(qromb)) RETURN
end if
s(j+1)=s(j)
h(j+1)=0.25_sp*h(j) !This is a key step: The factor is 0.25 even
!though the stepsize is decreased by only
!0.5. This makes the extrapolation a polynomial
!in h2 as allowed by equation (4.2.1),
!not just a polynomial in h.
end do
call nrerror('qromb:too many steps')
END FUNCTION qromb
---------------------------------------------------------------
But I see some errors after compiling with gfortran (version 4.2.1
(SUSE Linux) ) such as the following :
new.f90:363.19:
call nrerror(\xE2\x80\x99qromb:too many steps\xE2\x80\x99)
1
Error: Syntax error in argument list at (1)
new.f90:353.18:
call polint(h(j-KM:j),s(j-KM:j),0.0_sp,qromb,dqromb)
1
Error: Type/rank mismatch in argument 'xa' at (1)
new.f90:377.36:
n=assert_eq(size(xa),size(ya),\xE2\x80\x99polint\xE2\x80\x99)
1
Error: Syntax error in argument list at (1)
new.f90:387.19:
call nrerror(\xE2\x80\x99polint: calculation failure
\xE2\x80\x99)
1
Error: Syntax error in argument list at (1)
----------------------------
Is there anyone can guide me what is the mistake which causes the
mentioned errors?
Thanks in advance & Happy the Starting of the Spring!
Fatemeh
A.
> Dear all ;
>
> I use FUNCTION qromb ( from numerical recipes )in my program
> which calls polint subroutine.
>
> ------------------------------------------------------
>
> SUBROUTINE polint(xa,ya,x,y,dy)
> USE nrtype; USE nrutil, ONLY : assert_eq,iminloc,nrerror
> IMPLICIT NONE
> REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya
<-- removed copyrighted code posted probably in violation of copyright -->
> REAL(kind=8), DIMENSION(JMAXP) :: h,s !These store the
> But I see some errors after compiling with gfortran (version
> 4.2.1 (SUSE Linux) ) such as the following :
>
> new.f90:363.19:
>
> call nrerror(\xE2\x80\x99qromb:too many steps\xE2\x80\x99)
> 1
> Error: Syntax error in argument list at (1)
The argument must be a proper character string.
> call polint(h(j-KM:j),s(j-KM:j),0.0_sp,qromb,dqromb)
> 1
> Error: Type/rank mismatch in argument 'xa' at (1)
The error is exactly what the error message tells you: The dummy argument is a real_sp array, the actual argument is a real_dp array. Learn the language and consult your compiler documentation.
> ----------------------------
> Is there anyone can guide me what is the mistake which causes
> the mentioned errors?
There are, but their existence does not relieve you from your obligation to try to understand the problem yourself; educate yourself before you ask for help here.
> Fatemeh
--
-- mecej4
\xhh looks like an "escape" mechanism to insert unusual characters
into text. But your argument to nrerror is NOT a character string. It
does not start with either of the accepted character string delimiters
- single quote / apostrophe OR double quote / quote.
In fact the sequence of bytes 0xE2,0x80,0x99 (sorry for the C
notation) looks like a Unicode symbol for a printer's single quote
mark.
call nrerror('negative argument to square root')
or
call nrerror("invalid argument to arc-sine here")
might be valid.
I suspect that the \xnn sequences are an artifact of OCR (optical
character recognition) - the source code might have been scanned in
from a printed copy.
> new.f90:353.18:
>
> call polint(h(j-KM:j),s(j-KM:j),0.0_sp,qromb,dqromb)
> 1
> Error: Type/rank mismatch in argument 'xa' at (1)
> new.f90:377.36:
>
> n=assert_eq(size(xa),size(ya),\xE2\x80\x99polint\xE2\x80\x99)
> 1
Sorry, I've seen this problematic code (qromb and polint) posted in
this newsgroup before by a different student. My own internal error
message is "NEC001 - not enough coffee to look at this now". I suspect
that the error message IS describing your actual problem. Either your
actual argument is the wrong size OR it is the wrong type (single vs
double precision).
[snip]
> ----------------------------
> Is there anyone can guide me what is the mistake which causes the
> mentioned errors?
>
> Thanks in advance & Happy the Starting of the Spring!
> Fatemeh
WARNING - I'm going to bop over to your last thread on this project
and make some comments on numerical matters. [smile]
-- e
Perhaps, but at least in the USA, one could argue that posting
portions of code for review and or criticism, and that DOES happen in
this newsgroup, constitutes "fair use". [I am not a lawyer, nor do I
play one on TV.]
>
> > REAL(kind=8), DIMENSION(JMAXP) :: h,s !These store the
> > But I see some errors after compiling with gfortran (version
> > 4.2.1 (SUSE Linux) ) such as the following :
>
> > new.f90:363.19:
>
> > call nrerror(\xE2\x80\x99qromb:too many steps\xE2\x80\x99)
> > 1
> > Error: Syntax error in argument list at (1)
>
> The argument must be a proper character string.
>
> > call polint(h(j-KM:j),s(j-KM:j),0.0_sp,qromb,dqromb)
> > 1
> > Error: Type/rank mismatch in argument 'xa' at (1)
[snip]
-- I do think _some_ levity is needed in this newsgroup. <smile>
--- e
-- e
Thanks to people who guided me.
I could solve my problem according to your guidelines.
> I suspect that the \xnn sequences are an artifact of OCR (optical
> character recognition) - the source code might have been scanned in
> from a printed copy.
Or cut&pasted from one of the many formats/viewers intended for 'real'
(English/etc.) text -- PDF, MSWord, RTF, HTML sometimes, etc.
I couldn't agree more. A bit of humour really lightens the load.:-)
--
Freedom - no pane, all gaiGN!
Code Art Now
http://codeartnow.com
Email: a...@codeartnow.com