Solving a system of equations with numeric variables - Fortran

[leogt11]

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)

 

[SteamKing]

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).

 

[leogt11]

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:

 

[Mark44]

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:

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.

 

[SteamKing]

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.

 

[Mark44]

c is assigned a value of 10.25d0

I missed that line of the code...

 

[gsal]

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.