Dear all,
I am relatively new with Fortran, but I am making good progress.
I am stuck with a problem on how to solve a pentadiagonal matrix. Could someone assist me on how to go about this please?

Dear all,
...how to solve a pentadiagonal matrix... What do you mean "solve a ... matrix?"

Exactly what kind of help would you like?

Are you asking how to solve a system of equations whose matrix of coefficients is a banded matrix?

If so, then look at Numerical Recipes Section 2.4: You call bandec and banbks with m1 = 2 and m2 = 2.

If there is something that you don't understand (how to set up the problem or whatever) then ask specific questions.

Is this a homework problem?

What resources are you allowed to use?

Is it that the problem that you are "stuck with" is to find the solution to a particular set of equations?

Or is it that you are supposed to write a program that does the deed for a general problem of this type? (Using Numerical Recipes Fortran functions?)

Or is it that you just want an example using Numerical Recipes Fortran functions?

Or what?


Regards,

Dave

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.

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

Dear Dave,

Well I have not feel in anyway that you are accusing me of "something". In fact I fully comprehend why you posed those questions and where you are aiming at.

Thanks for your reply and the sample code.
But advice me please if I need to use the resources from Numerical Recipes.

For instance if I need to compile your sample code so as to see how it works, how do I have access to the subroutines you called? would you advice me typing them out? or is there a faster way to obtain them electronically?

Best regards.

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