Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solving a system of equations with numeric variables - Fortran

  1. Nov 4, 2013 #1
    Hello,

    I've been trying to solve a system of equations but I'm getting a lot of troubles when I tried to insert inside a matrix a numeric variable. This is my code. I've tried both schemes, i.e., (1) introducing all elements of the matrix by hand (real numbers) and (2) introducing numeric variables inside the matrix that are the result of other calculations. I've called this numeric variable c. Under the first scheme (1) the code is working but when I replace this by numeric variable it doesn't compile anymore. I would really appreciate if somebody could help me with this. Thank you

    program DD


    parameter (n=6, ia=6)
    integer :: i, j
    !integer :: n, ai
    dimension a(ia,n),b(n),r(n),l(n),s(n),x(n)
    dimension aa(ia,n),e(n)
    real :: c


    !Matrix that contains a numeric varibable
    c=10.25d0
    a(1,:) = (/c, 0d0, 0d0, 120d0, 49.4d0, 37d0/)
    a(2,:) = (/0.429, 0.273, 0.0, 0.0, 0.775, 0.885/)
    a(3,:) = (/0.0, 0.0, 0.111, 1.0, 0.225, 0.075/)
    a(4,:) = (/0.571, 0.727, 0.889, 0.0, 0.0, 0.04/)
    a(5,:) = (/1.0, -1.9, 0.0, 0.0, 0.0, 0.0/)
    a(6,:)=(/0.0, -0.61, 0.0, 0.0, 1.0, 0.0/)
    b =(/14.75, 0.375, 0.075, 0.55, 0.0, 0.0/)

    !Matrix constructed by introducing all values by hand (Code working under this scheme)
    !a(1,:) = (/10.25d0, 0d0, 0d0, 120d0, 49.4d0, 37d0/)
    !a(2,:) = (/0.429, 0.273, 0.0, 0.0, 0.775, 0.885/)
    !a(3,:) = (/0.0, 0.0, 0.111, 1.0, 0.225, 0.075/)
    !a(4,:) = (/0.571, 0.727, 0.889, 0.0, 0.0, 0.04/)
    !a(5,:) = (/1.0, -1.9, 0.0, 0.0, 0.0, 0.0/)
    !a(6,:)=(/0.0, -0.61, 0.0, 0.0, 1.0, 0.0/)
    !b =(/14.75, 0.375, 0.075, 0.55, 0.0, 0.0/)


    !save a

    do 2 i=1,n
    do 2 j=1,n
    aa(i,j) = a(i,j)
    2 continue

    call gauss(a,ia,n,l,s)
    call solve(a,ia,n,l,b,x)
    print *,' Initial solution is:'
    print 5,x
    print *

    do 4 k=1,5
    call residual(aa,ia,n,b,x,r)
    print *,' At iteration',k
    print *,' Residual vector is:'
    print 5,r
    call solve(a,ia,n,l,r,e)
    print *,' Error vector is:'
    print 5,e

    do 3 i=1,n
    x(i) = x(i) + e(i)

    3 continue

    print *,' Solution is:'
    print 5,x
    write(*,'(/a,3f12.2)')' x1= ',(x(1))
    write(*,'(/a,3f12.2)')' x2 = ',(x(2))
    write(*,'(/a,3f12.2)')' x3 = ',(x(3))
    write(*,'(/a,3f12.2)')' x4 = ',(x(4))
    write(*,'(/a,3f12.2)')' x5 = ',(x(5))
    write(*,'(/a,3f12.2)')' x6 = ',(x(6))

    print *
    4 continue
    5 format (3x,4(e13.6,2x))
    stop

    end

    subroutine residual(a,ia,n,b,x,r)

    !Residual vector r = b-Ax

    dimension a(ia,n),b(n),r(n),x(n)
    double precision sum
    do 3 i=1,n
    sum = 0.0
    do 2 j=1,n
    sum = sum + a(i,j)*x(j)
    2 continue
    r(i) = b(i) - sum

    3 continue

    return

    end

    subroutine gauss(a,ia,n,l,s)


    dimension a(ia,n),l(n),s(n)
    real max
    do 2 i=1,n
    l(i) = i
    s(i) = 0.0
    do 2 j=1,n
    s(i) = max(s(i),abs(a(i,j)))
    2 continue
    do 4 k=1,n-1
    rmax = 0.0
    do 3 i=k,n
    r = abs(a(l(i),k))/s(l(i))
    if (r .gt. rmax) then
    j = i
    rmax = r
    endif
    3 continue

    lk = l(j)
    l(j) = l(k)
    l(k) = lk
    do 4 i = k+1,n
    xmult = a(l(i),k)/a(lk,k)
    a(l(i),k) = xmult
    do 4 j = k+1,n
    a(l(i),j) = a(l(i),j) - xmult*a(lk,j)
    4 continue

    return
    end

    subroutine solve(a,ia,n,l,b,x)

    !LUx = b

    dimension a(ia,n),l(n),b(n),x(n),z(4)

    do 3 i=1,n
    sum = b(l(i))
    do 2 j=1,i-1
    sum = sum - a(l(i),j)*z(j)
    2 continue
    z(i) = sum
    3 continue
    do 5 i=n,1,-1
    sum = z(i)
    do 4 j=i+1,n
    sum = sum - a(l(i),j)*x(j)
    4 continue
    x(i) = sum/a(l(i),i)
    5 continue

    return

    end



    This are the errors that I'm getting

    Compilation error time: 0 memory: 4032 signal:0

    prog.f95:35.19:

    do 2 j=1,n
    1
    Warning: Obsolescent feature: Shared DO termination label 2 at (1)
    prog.f95:103.19:

    do 2 j=1,n
    1
    Warning: Obsolescent feature: Shared DO termination label 2 at (1)
    prog.f95:119.23:

    do 4 i = k+1,n
    1
    Warning: Obsolescent feature: Shared DO termination label 4 at (1)
    prog.f95:122.26:

    do 4 j = k+1,n
    1
    Warning: Obsolescent feature: Shared DO termination label 4 at (1)
    prog.f95:14.20:

    a(1,:) = (/c, 0d0, 0d0, 120d0, 49.4d0, 37d0/)
    1
    Error: Element in REAL(4) array constructor at (1) is REAL(8)
    prog.f95:86.16:

    r(i) = b(i) - sum
    1
    Warning: Possible change of value in conversion from REAL(8) to REAL(4) at (1)
     
  2. jcsd
  3. Nov 4, 2013 #2

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    You have REAL(8) constants in the first row of matrix a and REAL(4) constants elsewhere.

    Try making the first row of matrix a all of type REAL (4). Or, if you are worried about the precision of the result, make everything REAL (8).
     
  4. Nov 4, 2013 #3
    How do I make the first row of matrix a all of type REAL (4). You mean like this real*4 :: c

    prog.f95:9.12:

    real*4 :: c
    1
    Error: GNU Extension: Nonstandard type declaration REAL*4 at (1)
    prog.f95:43.19:
     
  5. Nov 4, 2013 #4

    Mark44

    Staff: Mentor

    You're missing the parentheses.

    real(4) :: c

    Also, and this is very important, you cannot do what you're trying to do. You have a variable named c in your array. You didn't say, but I would guess that you are trying to write a program that solves for c. Fortran doesn't work this way.

    Since you haven't given c a value by either initialization or assignment, c will have whatever value is lying around in the memory that is set aside for c. Any calculation you do will be with this "garbage" value.
     
  6. Nov 4, 2013 #5

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    c is assigned a value of 10.25d0

    IDK why the constants in the first row of the a matrix are assigned real(8) values, but these should probably be expressed as real(4) constants.
     
  7. Nov 4, 2013 #6

    Mark44

    Staff: Mentor

    I missed that line of the code...
     
  8. Nov 4, 2013 #7
    Please be more diligent; you are mixing various programming practices and not even good ones.

    You are mixing F77 and F90 practices.
    Please do not use " dimension " alone to declare variables without using a type...these variables end up being the size of the default REAL..whatever that is...could be 4 or 8 bytes depending on compilation flags. The lowercase "L" ends up as integer, of course.
    The use of " d " at the end of a number is for double precision, but you did not define any "double precision " variables, all you defined were REAL
    real*4 is correct to define a scalar of 4 bytes long; real(4) is an array of reals of the default size.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Solving a system of equations with numeric variables - Fortran
Loading...