Write a Fortran program to implement the Jacobi algorithm.

[nitsmistry]

Hi,

I'm having a bit of a problem with a project i'm undertaking and i was wondering if anyone could help. The problem is:

Write a Fortran program to implement the Jacobi algorithm that includes the following steps.
(a) Find the off-diagonal element a[pq] of largest magnitude.
(b) Compute the Frobenius norm of the off-diagonal elements, off(A), and test whether convergence has been achieved, by setting the tolerance = 0.0001 and checking
(off(A)/frob(A)) < tolerance.
(c) Compute (theta), cos(theta) , sin(theta).
(d) Set A(n+1) = Rn(Transpose)ARn ; n = 1, 2, ....,N.
(e) When convergence is achieved, print out all the eigenvalues.
(f) Compute matrix P = n=1:N of Rn
(g) Print out all the eigenvectors corresponding to each eigenvalue.
Thanks for your time

 

[xwb]

What compiler/OS are you using
Which part are you having problems with?

 

[nitsmistry]

hi, im using plato 3 with the salford compiler, ive started writing the program but there are a few errors which occured. Also i think my code needs to be more efficient, my friend has helped me write some of it but i cant understand where a few bits have come from (for instance the block statement at the beginning)...
Heres what i have written so far neway..
(NB:my jacobi algorithm is supposed to compute the eigenvalues of a symmetric matrix)

*HPF$ DISTRIBUTE A (BLOCK, BLOCK)
*HPF$ ALIGN B(I,J) WITH A(I,J)
C arrays A and B with block distribution
PRINT *, '******** TEST_JACOBI_HPF ********'
C nest of two independent loops, iteration (i,j) will be executed
C on processor, which is owner of element A(i,j)
*HPF$ INDEPENDENT
DO 1 J = 1, K
*HPF$ INDEPENDENT
DO 1 I = 1, K
A(I,J) = 0.
IF(I.EQ.1 .OR. J.EQ.1 .OR. I.EQ.K .OR. J.EQ.K) THEN
B(I,J) = 0.
ELSE
B(I,J) = 1. + I + J
ENDIF
1 CONTINUE
DO 2 IT = 1, ITMAX
*HPF$ INDEPENDENT
DO 21 J = 2, K-1
*HPF$ INDEPENDENT
DO 21 I = 2, K-1
A(I,J) = B(I,J)
21 CONTINUE
*HPF$ INDEPENDENT
DO 22 J = 2, K-1
*HPF$ INDEPENDENT
DO 22 I = 2, K-1
B(I,J) = (A(I-1,J) + A(I,J-1) + A(I+1,J) + A(I,J+1)) / 4
22 CONTINUE
2 CONTINUE
3 OPEN (3, FILE='JACH.DAT', FORM='FORMATTED')
WRITE (3,*) B
CLOSE (3)
END

Once again, thank you for your time.

 

[xwb]

You have part of a program in F77 with lots of missing bits.

1) All arrays need to be declared before use
2) All variables need to be initialized before you refer to them.

What size does your array go to?

I don't know where the *HPF$ stuff comes from or what it means. In F77 any non numeric char in column 1 is a comment but this looks like a compiler directive of some sort.

 

[nitsmistry]

ok thank you, im goin to work on this program tonight and tomorrow. Hopefully il be a lot clearer by then. Thanks again!

 

[nitsmistry]

hello, i decided to start again...however im stuck on part (b) where i need to check for convergence...here is my code so far...

program JacobiAlgorithm
implicit none
real, allocatable,dimension ,:: A,D,invD,B
integer::n,p,q,r,s
real::tolerance,m
integer,dimension(2)::loc
print*,'What are the dimensions of your matrix'
read*,n
allocate(A(n,n))
allocate(D(n,n))
allocate(invD(n,n))
allocate(B(n,n))
print*,'The dimensions of this matrix is ',n, ' x' ,n,''
do p=1,n
do q=1,n
print*,'Input the entries of matrix A(',p,',',q,')'
read *,A(p,q)
end do
end do
do p=1,n
do q=1,n
if (A(p,q)/=A(q,p)) Then
print*,'This matrix is not symmetric'
stop
end if
end do
end do
print*,'Your matrix is'
do p=1,n
write(*,*)(A(p,q),q=1,n)
end do


do p=1,n
do q=1,n
if (p<q)then
B(p,q)=A(p,q)
else
B(p,q)=0
end if
end do
end do
m=maxval(abs(A),(p<q))
loc=maxloc(abs(B))
print*,'The maximum value where p<q is',m
print*,'The maximum value is in the location',loc
r=loc(1)
s=loc(2)
do r=1,n
do s=1,n

end program JacobiAlgorithm

 

[gummibaer]

Could you specify what you mean with convergence ?
Which entities of the algorithm should converge to which value ? Which result should be stable ? Where would this convergence depend on ?

Once we could define the problem in these terms, I think we would have made it.

Norbert