Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations Mike Lewis on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Fortran Matrix Help

Status
Not open for further replies.

AdYrL

Technical User
Oct 18, 2008
5
SG
Hi,this is a program i have for solving a 3x3 matrix. The thing is, that it won't solve. I suspect that there's somthing wrong with my definition, it'll appear as 0 errors but like 7 warnings. Can anybody help me with this?

program testing

double precision A(3,3),B(3),Z(1000)
real N,lud
! Initialising the Matrix A

A(1,1)=1
A(1,2)=2
A(1,3)=3

A(2,1)=4
A(2,2)=2
A(2,3)=2

A(3,1)=1
A(3,2)=3
A(3,3)=5

!Initialising the matrix B
B(1)=2
B(2)=3
B(3)=2

call SOLVER(A,B,N,lud,Z)

print*,"The values are:",Z(1),Z(2),Z(3)

end program testing

!*************************************************
subroutine SOLVER(A,B,N,lud,Z)

real lda,N,ipvt(1000),info,lud,IDAMAX,j,k,kp1,l,nm1,kb

double precision A(1000,1000),B(1000),Z(1000),t,AMD(1000,1000)

common /ludcmp/ipvt,AMD

nm1=N-1

do 5 i=1,N
Z(i)=B(i)
5 continue

if (lud.eq.0) goto 99

do 6 i=1,N
do 6 j=1,N
AMD(i,j)=A(i,j)
6 continue

!This part for decomposing A is taken from DGEFA.

info=0

if (nm1.lt.1) go to 70
do 60 k=1,nm1
kp1=k+1
l=IDAMAX(N-k+1,AMD(k,k),1)+k-1
ipvt(k)=l
if (AMD(l,k).eq.0.0d0) goto 40
if (l.eq.k) goto 10
t=AMD(l,k)
AMD(l,k)=AMD(k,k)
AMD(k,k)=t
10 continue
t=-1.0d0/AMD(k,k)
call DSCAL(N-k,t,AMD(k+1,k),1)
do 30 j=kp1,N
t=AMD(l,j)
if (l.eq.k) go to 20
AMD(l,j)=AMD(k,j)
AMD(k,j)=t
20 continue
call DAXPY(N-k,t,AMD(k+1,k),1,AMD(k+1,j),1)
30 continue
go to 50
40 continue
info=k
50 continue
60 continue
70 continue
ipvt(N)=N
if (AMD(N,N).eq.0.0d0) info=N
if (info.ne.0) pause 'Division by zero in SOLVER!'

99 continue

!This part for finding solution of AX=B is taken from DGESL.

if (nm1.lt.1) goto 130
do 120 k=1,nm1
l=ipvt(k)
t=Z(l)
if (l.eq.k) goto 110
Z(l)=Z(k)
Z(k)=t

....
....
...
....

Thanks!
 
Before you use any subroutine you should know what it should do.
Here you should solve the system of 3 equations
A*Z = B.
If A is 3x3 matrix, then Z should have 3 elements too. So deklare first
Code:
double precision A(3,3),B(3),Z(3)

Then calling the SOLVER subroutine has this form
Code:
call SOLVER(A,B,N,lud,Z)

So, before the above call you have to set these arguments too:
N - I guess it's number of equations, i.e. 3 in your case
lud - I guess it's flag for LU-decomposition of the matrix. If you don't want to do LU-decomposition, then lud=0 (see your source)

Try this
Code:
...
N=3
lud=0
call SOLVER(A,B,N,lud,Z)
...
and your program will probably work
 
Hi,thanks for your help but i've tried out what you have suggested, but the same problem has still occured. It still showed 0 errors and 5 warnings when compiling. Adding the n=3 and lud =0 before the calling of my solver program resulted in a windows error box. Do u know what's the problem still? Thanks!
 
If you want to help, you should post the 5 warnings exactly.
 
The warning message actually says:

Warning: In the call to IDAMAX, actual argument #1 does not match the type and kind of the corresponding dummy argument.
l=IDAMAX(N-k+1,AMD(k,k),1)+k-1

Warning: In the call to DSCAL, actual argument #1 does not match the type and kind of the corresponding dummy argument.
call DSCAL(N-k,t,AMD(k+1,k),1)

Warning: In the call to DAXPY, actual argument #1 does not match the type and kind of the corresponding dummy argument.
call DAXPY(N-k,t,AMD(k+1,k),1,AMD(k+1,j),1)

somthing to do with N-k

Thanks
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top