Need help pls.

azy71

10-01-2008, 10:15 AM

Dear Dave,

Thanks for your comment.

No is not a home work, I need to write a code that will solve a system of equations whose matrix of coefficeints is banded.

I have checked Numerical Recipes Section 2.4 you referred me to, the section discussed matrix inversion, which I am afraid I can not make out how to employ that for the purpose I envisaged.

Thanks.

davekw7x

10-01-2008, 12:57 PM

No is not a home work I hope that you didn't think I was accusing you of something; I just wanted to know what kind of help you would like (and whether there were any particular things you were supposed to use or not use).

Here's an example program, that I compiled with GNU gfortran.

!

! Illustration of use of bandec and bandbks to

! solve a system of equations whose matrix of

! coefficients is a banded matrix

!

! This particular example is for a "pentadiagonal"

! matrix. That is: m1 = 2 and m2 = 2

!

! davekw7x

! September, 2008

!

PROGRAM xbandec

USE nrtype

USE nr

IMPLICIT NONE

INTEGER(I4B) :: m1 = 2, m2 = 2

REAL(SP) :: d

REAL(SP), DIMENSION(8) :: x

Real(SP), DIMENSION(8) :: b

REAL(SP), DIMENSION(8,2) :: al

REAL(SP), DIMENSION(8,5) :: a ! m1+m2+1 = 5

INTEGER(I4B), DIMENSION(8) :: ix

INTEGER(I4B) :: i, j, idx

! The full-sized matrix

REAL(SP), DIMENSION(8,8) :: biggie

biggie(1,:)=(/ 1.0, 2.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)

biggie(2,:)=(/ 4.0, 5.0, 6.0, 7.0, 0.0, 0.0, 0.0, 0.0 /)

biggie(3,:)=(/ 8.0, 9.0, 1.0, 2.0, 3.0, 0.0, 0.0, 0.0 /)

biggie(4,:)=(/ 0.0, 4.0, 5.0, 6.0, 7.0, 8.0, 0.0, 0.0 /)

biggie(5,:)=(/ 0.0, 0.0, 9.0, 1.0, 2.0, 3.0, 4.0, 0.0 /)

biggie(6,:)=(/ 0.0, 0.0, 0.0, 5.0, 6.0, 7.0, 8.0, 9.0 /)

biggie(7,:)=(/ 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 4.0 /)

biggie(8,:)=(/ 0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 6.0, 7.0 /)

x =(/ 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 /)

write(*,'(" Here is the original complete matrix:")')

do i=1,8

write(*,'(10I4)') (int(biggie(i,j)),j=1,8)

end do

write(*,*)

! Form the banded matrix from the original

do i = 1,8

do j=1,m1+m2+1

idx = i+j-3

if ((idx .LT. 1) .OR. (idx .GT. 8)) then

a(i,j) = 0

else

a(i,j)=biggie(i,idx)

end if

end do

end do

write(*,'(" This is the banded matrix [a]:")')

do i=1,8

write(*, '(10I4)') (int(a(i,j)),j=1,m1+m2+1)

end do

write(*,*)

write(*,'(" Vector [x]:")')

write(*,'(F9.1)') (x(i), i=1,8)

write(*,*)

call banmul(a(1:8,1:m1+m2+1),2,2,x,b)

write(*,'(" After multiplying [a] times [x] this is [b]:")')

write(*,'(F9.2)') (b(i), i=1,8)

write(*,*)

call bandec(a(1:8,1:m1+m2+1),m1,m2,al(1:8,1:2),ix(1:8), d)

call banbks(a(1:8,1:m1+m2+1),m1,m2,al(1:8,1:2),ix(1:8), b)

write(*,'(" After solving [a][x] = [b] for [x]:")')

write(*,'(" Original Solved")')

write(*,'(" x x")')

write(*,'(" ------- -------")')

write(*,'(F12.5,F12.5)') (x(i), b(i),i=1,8)

!

! After bandec:

! The determinant is "d" times the product of the first column of [a]

!

do i=1,8

d = d * a(i,1)

end do

write (*, '(/" Determinant of [a] = ", 1PE12.5/)') d

!

! Just for comparison, use ludcmp on the original

! matrix and calculate its determinant

!

!call ludcmp(biggie(1:8,1:8),ix(1:8),d)

!do i=1,8

!d = d*biggie(i,i)

!end do

!write(*,'(/" From ludcmp, determinant of biggie = ",1PE12.5)') d

END PROGRAM xbandec

I defined a system of 8 equations and 8 unknowns. The matrix of coefficients is a banded matrix with non-zero elements on the main diagonal and the two subdiagonals on either side of the main diagonal.

I showed how to get a banded matrix, a from the original matrix so that the Numerical Recipes functions could be applied.

I defined an 8-element vector, x and showed how to use banmul to multiply the banded matrix times the vector. The result is in vector b.

Then I showed how to use bandec and banbks to solve the sytstem a x = b for x

Finally I showed the solution from the system alongside the original vector.

Output:

Here is the original complete matrix:

1 2 3 0 0 0 0 0

4 5 6 7 0 0 0 0

8 9 1 2 3 0 0 0

0 4 5 6 7 8 0 0

0 0 9 1 2 3 4 0

0 0 0 5 6 7 8 9

0 0 0 0 1 2 3 4

0 0 0 0 0 5 6 7

This is the banded matrix [a]:

0 0 1 2 3

0 4 5 6 7

8 9 1 2 3

4 5 6 7 8

9 1 2 3 4

5 6 7 8 9

1 2 3 4 0

5 6 7 0 0

Vector [x]:

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

After multiplying [a] times [x] this is [b]:

1.40

6.00

5.20

13.00

8.70

22.00

7.00

12.80

After solving [a][x] = [b] for [x]:

Original Solved

x x

------- -------

0.10000 0.10000

0.20000 0.20000

0.30000 0.30000

0.40000 0.40000

0.50000 0.50000

0.60000 0.60000

0.70000 0.70000

0.80000 0.80000

Determinant of [a] = 6.18900E+04

I threw in the calculations for the determinant just for kicks.

Regards,

Dave

davekw7x

10-02-2008, 09:39 AM

...need to use the resources from Numerical Recipes...

In addition to the main program file, which I called xband.f90, you need the following:

nrtype.f90

nr.f90

nrutil.f90

bandec.f90

banbks.f90

banmul.f90

Source files nr.f90, nrtype.f90 and nrutil.f90 are for modules; the others are the source code for the subprograms used by the test program.

For my compiler (GNU gfortran), I compile the modules separately.

Assuming all are in the same directory and you are using GNU gfortran, you can build the application by:

gfortran -c nrtype.f90

gfortran -c nr.f90

gfortran -c nrutil.f90

gfortran -c bandec.f90

gfortran -c banmul.f90

gfortran -c banbks.f90

gfortran -c xband.f90

gfortran xband.o bandec.o banmul.o banbks.o nrutil.o -oxband

In fact I create a Makefile that uses object and module outputs in another directory where I have previously compiled them into a library, but for a one-off application like this, the above sequence works on my Linux system.

...would you advice me typing them out? or is there a faster way ...

The Fortran versions are no longer supported by the authors but the latest published release, Numerical Recipes Source Code CD-ROM 3rd Edition, has legacy code for the Fortran functions and modules. The code is copyrighted.

Regards,

Dave