About function acos in IMSL FNL

[Zhiyong_b_]

hello, i was using intel fortran with IMSL FNL and faced some problems.

I used the special function by using 'ACOS_INT'. when calculating cos(x) for x<1 and x>1 separately it works fun, but it did not work when I try do x from a1(<1) to a2(>1). Here follows my fortran code:

program main
    include'link_fnl_shared.h'
      USE ACOS_INT
      IMPLICIT   NONE

      REAL*8, PARAMETER:: T = 1
      REAL*8, PARAMETER:: MUS = -1.95
      REAL*8, PARAMETER:: MUN = -1.95
      REAL*8, PARAMETER:: W   = 2.0
      REAL*8, PARAMETER:: DELTAL = 0.02
      REAL*8, PARAMETER:: DELTAR = 0.02
      REAL*8, PARAMETER:: NY  = 30
      REAL*8, PARAMETER:: NX  = 60
      REAL*8, PARAMETER:: NUMR = 5
      REAL*8, PARAMETER:: KBT = 0.001 
      REAL*8, PARAMETER:: PI = 4.d0*datan(1.d0)
      INTEGER    N1
      COMPLEX*8::    VALUE1, Z
      real*8:: omega=0.5d0, omegal, omegar
      complex*8, dimension(NY+1):: PL1, PL2, PR1, PR2
      real*8, dimension(NY+1):: Q, L1, L2, R1, R2
      integer:: i

      OMEGAR = SQRT(OMEGA**2 + DELTAR**2)
      OMEGAL = SQRT(OMEGA**2 + DELTAL**2)

      DO N1 = 1, NY+1
      Q(N1) = N1*PI/(2+NY)
      L1(N1) = -(MUN + 2*T*COS(Q(N1)+OMEGAL))/(2*T)
      L2(N1) = -(MUN + 2*T*COS(Q(N1)-OMEGAL))/(2*T)
      R1(N1) = -(MUN + 2*T*COS(Q(N1)+OMEGAR))/(2*T)
      R2(N1) = -(MUN + 2*T*COS(Q(N1)-OMEGAR))/(2*T)
      END DO

      PL1 = ACOS(L1)
      PL2 = ACOS(L2)
      PR1 = ACOS(R1)
      PR2 = ACOS(R2)

      open(11, file='data1.txt')
      do i =1, ny
      
      WRITE (11,*) L1(i),PR2(i)
99999 FORMAT (' ACOS((', F6.3, ',', F6.3, ')) = (', &
         F6.3, ',', F6.3, ')')
      end do
END

If anyone don't think I described it clear, I will be more clear.

Thanks a lot.

[jimdempseyatthecove]

My guess is that some element(s) of L1, L2, R1, R2 are (slightly) outside the range of -1.0 to 1.0. This can happed due to round off errors in the expression that produced the values. Try

L1(N1) = MIN(MAX(-(MUN + 2*T*COS(Q(N1)+OMEGAL))/(2*T), -1.0), 1.0)
...

Jim Dempsey
 

[mecej4 · Accepted solution]

The IMSL routine ACOS is defined to take a scalar argument. When you attempt to pass an array, the type matching rules probably cause the compiler to generate a call to the intrinsic function ACOS, which takes a real argument. This intrinsic function will not work for arguments whose absolute value is > 1.0.

To test the conjecture, I replaced the USE statement by "USE ACOS_INT, ACOSIMSL=>ACOS", changed "ACOS" to "ACOSIMSL" in your assignment statements, and changed the type of the arrays L1, L2, R1, R2 to COMPLEX*8. The compiler now says:

xacos.f90(35): error #6362: The data types of the argument(s) are invalid.   [ACOSIMSL]
      PL1 = ACOSIMSL(L1)
---------------------^

Therefore, one solution is to compute PL1, etc., inside the DO loop, passing one array element at a time to the IMSL function ACOS.

[Zhiyong_b_]

jimdempseyatthecove wrote:

My guess is that some element(s) of L1, L2, R1, R2 are (slightly) outside the range of -1.0 to 1.0. This can happed due to round off errors in the expression that produced the values. Try

L1(N1) = MIN(MAX(-(MUN + 2*T*COS(Q(N1)+OMEGAL))/(2*T), -1.0), 1.0)

...

Jim Dempsey

 

yes, you are right. Some of L1, L2, R1, R2's elements are out of the range of [-1, 1], my problem here, for example: if L1 is ranging from -1 to 1, the ACOS function works fine, and when all L1 elements are bigger than 1, it works fine too. But unfortunately, in here like L1,  its range is from 0 to 1.9 and it did not give a result for L1(i) > 1.   That's the reason I post here. And of course, for x > 1, ACOS(x) is a complex number which I am pretty sure.

[Zhiyong_b_]

mecej4 wrote:

The IMSL routine ACOS is defined to take a scalar argument. When you attempt to pass an array, the type matching rules probably cause the compiler to generate a call to the intrinsic function ACOS, which takes a real argument. This intrinsic function will not work for arguments whose absolute value is > 1.0.

To test the conjecture, I replaced the USE statement by "USE ACOS_INT, ACOSIMSL=>ACOS", and changed "ACOS" to "ACOSIMSL" in your assignment statements, and changed the type of the arrays L1, L2, R1, R2 to COMPLEX*8. The compiler now says:

xacos.f90(35): error #6362: The data types of the argument(s) are invalid.   [ACOSIMSL]
      PL1 = ACOSIMSL(L1)
---------------------^

Therefore, one solution is to compute PL1, etc., inside the DO loop, passing one array element at a time to the IMSL function ACOS.

In fact, I changed the 27th line of my code into:

DO N1 = 1, NY+1
   TQ     =  N1*PI / (2+NY)
   TL1   = -( MUN + 2*T*COS( TQ + OMEGAL ) ) / (2*T)
   TL2   = -( MUN + 2*T*COS( TQ -  OMEGAL ) ) / (2*T)
   TR1  = -( MUN + 2*T*COS( TQ + OMEGAR ) ) / (2*T)
   TR2  = -( MUN + 2*T*COS( TQ -  OMEGAR ) ) / (2*T)
   Q(N1)  =  TQ
   L1(N1) =  TL1
   L2(N1) =  TL2
   R1(N1) = TR1
   R2(N1) = TR2
END DO

and it works fine now

[jimdempseyatthecove]

And is your former coding of:

      PL1 = ACOS(L1)
      PL2 = ACOS(L2)
      PR1 = ACOS(R1)
      PR2 = ACOS(R2)

the same?

If so, the use of local temps inside did not "fix" the problem, it made the symptoms go away. mecej4 is suggesting

DO N1 = 1, NY+1
  Q(N1)   = N1*PI/(2+NY)
  L1(N1)  = -(MUN + 2*T*COS(Q(N1)+OMEGAL))/(2*T)
  L2(N1)  = -(MUN + 2*T*COS(Q(N1)-OMEGAL))/(2*T)
  R1(N1)  = -(MUN + 2*T*COS(Q(N1)+OMEGAR))/(2*T)
  R2(N1)  = -(MUN + 2*T*COS(Q(N1)-OMEGAR))/(2*T)
  PL1(N1) = ACOS(L1(N1))
  PL2(N1) = ACOS(L2(N1))
  PR1(N1) = ACOS(R1(N1))
  PR2(N1) = ACOS(R2(N1))
END DO

His point being, in your former use PL1 = ACOS(L1) you were performing an array statement... that was calling the standard Fortran library function ACOS that requires argument to be in -1.000000000000000000.... to +1.0000000000000000000....

The difference with the above loop is PL1(N1) = ACOS(L1(N1)) is now the ACOS is being called with a scalar (single real*8) and because you are also using IMSL ACOS_INT, the use of scalar with ACOS calls the ACOS_INT function which performs a fixup on the input arg.

Check to see if the above loop vectorizes. If it does not then insert the MIN(MAX fixup as this may have a better chance of vectorization.

Jim Dempsey