0

I am trying to solve a linear equation in Fortran using the DGESV subroutine in LAPACK:

enter image description here

I have the following code:

program main
    implicit none (type, external)
    external         :: DGESV
    real(kind=8)     :: para_morse(3)     ! Vector b/x.
    real(kind=8)     :: pi, delta, Xtop, Ytop, Dtop
    real(kind=8)     :: Xhollow, Yhollow, Dhollow
    real(kind=8)     :: Xbridge, Ybridge, Dbridge
    real(kind=8)     :: a00,a01,a02,a10,a11,a12,a20,a21,a22
    INTEGER          :: N, NRHS, LDA, LDB, INFO
    INTEGER          :: IPIVOT(3)
    DOUBLE PRECISION :: A(3,3), B(3,1)

    pi    = 3.14159265359
    delta = 3.17*10**(-10)

    Xtop    = 0
    Ytop    = 0
    Dtop    = 4.9984788

    Xhollow = 0.5
    Yhollow = 0.5
    Dhollow = 7.60346595

    Xbridge = 0.5
    Ybridge = 0
    Dbridge = 5.56938086

    a00 = 1
    a10 = 1
    a20 = 1

    a01 = (cos((2*pi*Xbridge)/delta) + cos((2*pi*Ybridge)/delta))
    a11 = (cos((2*pi*Xtop)/delta)    + cos((2*pi*Ytop)/delta))
    a21 = (cos((2*pi*Xhollow)/delta) + cos((2*pi*Yhollow)/delta))

    a02 = (cos((2*pi*(Xbridge+Ybridge))/delta) + cos((2*pi*(Xbridge-Ybridge))/delta))
    a12 = (cos((2*pi*(Xtop+Ytop))/delta)       + cos((2*pi*(Xtop-Ytop))/delta))
    a22 = (cos((2*pi*(Xhollow+Yhollow))/delta) + cos((2*pi*(Xhollow-Yhollow))/delta))

    A   = reshape([ a00, a10, a20, a01, a11, a21, a02, a12, a22], [ 3, 3 ])
    B   = reshape([ Dbridge, Dtop, Dhollow], [ 3, 1 ])

    N       = 3
    NRHS    = 1
    IPIVOT  = 3
    LDA     = 3
    LDB     = 3


    CALL DGESV(N, NRHS, A, LDA, IPIVOT, B, LDB, INFO)

    if (info /= 0) then
        print '(a, i0)', 'Error: ', info
        stop
    end if

    print *, "Solution (x1, x2, x3): ", B
end program main

However, all my solutions are being displayed as NaN. I have solved the same equation using NumPy linalg.solve in python and I obtain the following result [ 7.6678536 -0.58501309 -31.07990139]. Can someone help me understand why I am getting NaN as the result?

francescalus
  • 30,576
  • 16
  • 61
  • 96
Jules
  • 35
  • 1
  • 7
  • 3
    What attempts have you made so far to debug your program, such as looking at intermediate calculations, using compiler flags – francescalus Feb 21 '21 at 13:11
  • 2
    It does not cause the main error, but be aware that `pi = 3.14159265359` is only a single-precision value. See https://stackoverflow.com/questions/6146005/is-there-a-better-double-precision-assignment-in-fortran-90 The same holds for other numerical values you use. – Vladimir F Героям слава Feb 21 '21 at 13:18
  • 1
    The error isn't exactly the same as in the linked question (`10**-10` is zero, rather than too large), but the misunderstanding is: `10**-10` is evaluated as an integer expression where you want it to be real (of kind 8 whatever that is): `10._8**-10`. – francescalus Feb 21 '21 at 13:22
  • @VladimirF Oh, thank you. I was not aware of this. – Jules Feb 21 '21 at 16:05
  • @francescalus yes! when I changed it to real, I am able to obtain some values. So, earlier I was getting a division by zero instead. Thank you so much for your help! – Jules Feb 21 '21 at 16:08
  • @VladimirF I am also able to obtain values closer to the required values by using your advice, Thank you. – Jules Feb 21 '21 at 16:10
  • 2
    Make sure you understand Valdimir's comment properly - almost all your maths in setting up the equations is in default, which is probably single, precision. Python by default uses (I believe) double precision. You will need to explicitly specify the kind you require for all those real constants if you want to get the same. – Ian Bush Feb 21 '21 at 16:12

0 Answers0