|
|
 |
 |
 |
 |
Fortran Programming Language
|
 |
 |
 |
 |
|
|
|
|
Matrix multiplication
Hi!
Im really newbie in Fortran, so I have started easy and written some
simple program for serial matrix multiplication, and some for parallel
matrix vector multiplication. So I moved on to matrix*matrix
multiplication, but then printing it out I get only zeros as answer..
but serial computation shows otherwise. Maybe you can point me out were
in code I have lost my soul? Thank you in advance!
PROGRAM parallel_mv
! global arrays
real, dimension(16,16) :: a
real, dimension(16,16) :: b,y
! variables for BLACS initialization
!and processor grid creation
integer iam,nprocs,ictxt,nprow,npcol,myrow,mycol
!variables needed for distributing global
!arrays across the proc grid
integer desca(9), descb(9),descy(9),m,n,mb,nb,rsrc,csrc
integer llda,lldb,info
! local arrays
real, dimension(8,8) :: la
! local array for matrix b
real, dimension(8,8) :: lb
! local array for matrix y
real, dimension(8,8) ::ly
!Initializing the BLACS library (STEP 2)
call blacs_pinfo (iam,nprocs)
call blacs_get(-1,0,ictxt)
!Creating and using the processor grid (STEP 3)
nprow=2; npcol=2
call blacs_gridinit(ictxt,'r',nprow,npcol)
call blacs_gridinfo(ictxt,nprow,npcol,myrow,mycol)
! Making the array descriptor vectors (STEP 4)
m=16; n=16
mb=8; nb=8
rsrc=0; csrc=0
llda=8
call descinit(desca,m,n,mb,nb,rsrc,csrc,ictxt,llda,info)
m=16; n=16
mb=8; nb=8
rsrc=0; csrc=0
lldb=8
call descinit(descb,m,n,mb,nb,rsrc,csrc,ictxt,lldb,info)
call descinit(descy,m,n,mb,nb,rsrc,csrc,ictxt,lldb,info)
! Filling the global arrays A,b
open(unit=12,file="a.data")
read(12,*) a
open(unit=13,file="b.data")
read(13,*) b
! Each processor fills in its
!local arrays with correct elements
! from the global arrays (STEP 5)
print *,"filling local matrix"
if(myrow.eq.0.and.mycol.eq.0) then
do i_loc=1,8
do j_loc=1,8
!fill a
la(i_loc,j_loc)=a(i_loc,j_loc)
!fill b
la(i_loc,j_loc)=b(i_loc,j_loc)
end do
end do
end if
if(myrow.eq.1.and.mycol.eq.0) then
do i_loc=1,8
do j_loc=1,8
!fill a
la(i_loc,j_loc)=a(i_loc,j_loc)
!fill b
la(i_loc,j_loc)=b(i_loc,j_loc)
end do
end do
end if
if(myrow.eq.0.and.mycol.eq.1) then
do i_loc=1,8
do j_loc=1,8
!fill a
la(i_loc,j_loc)=a(i_loc,j_loc)
!fill b
la(i_loc,j_loc)=b(i_loc,j_loc)
end do
end do
end if
if(myrow.eq.1.and.mycol.eq.1) then
do i_loc=1,8
do j_loc=1,8
!fill a
la(i_loc,j_loc)=a(i_loc,j_loc)
!fill b
la(i_loc,j_loc)=b(i_loc,j_loc)
end do
end do
end if
! Call the ScaLAPACK routine (STEP 6)
n=4
print *,"caling psgemm"
call psgemm('n','n',n,n,n,1.,la,1,1,desca(1),lb,1,1,descb(1),
& 0.,c,1,1,descy(1))
! Each processor prints out its part of the product
!vector y (STEP 7)
if(myrow.eq.0.and.mycol.eq.0) then
do i=1,8
do j=1,8
print *,ly(i), ly(j)
end do
end do
end if
if(myrow.eq.1.and.mycol.eq.0) then
do i=1,8
do j=1,8
print *,ly(i), ly(j)
end do
end do
end if
! Release the proc grid and BLACS library (STEP 8)
call blacs_gridexit(ictxt)
call blacs_exit(0)
STOP
END
In a previous article, Carramba <u@example.net> wrote:
>Hi!
>Im really newbie in Fortran, so I have started easy and written some
>simple program for serial matrix multiplication, and some for parallel
>matrix vector multiplication. So I moved on to matrix*matrix
>multiplication, but then printing it out I get only zeros as answer..
>but serial computation shows otherwise. Maybe you can point me out were
>in code I have lost my soul? Thank you in advance!
maybe all your numbers are all less than 1.0.
Chris
m@skyway.usask.ca skrev:
> In a previous article, Carramba <u @example.net> wrote: >> Hi! >> Im really newbie in Fortran, so I have started easy and written some >> simple program for serial matrix multiplication, and some for parallel >> matrix vector multiplication. So I moved on to matrix*matrix >> multiplication, but then printing it out I get only zeros as answer.. >> but serial computation shows otherwise. Maybe you can point me out were >> in code I have lost my soul? Thank you in advance! > maybe all your numbers are all less than 1.0.
I have run serial version on same matrices and get real number.
> Chris
On May 18, 4:42 pm, Carramba <u@example.net> wrote:
> Hi! > Im really newbie in Fortran, so I have started easy and written some > simple program for serial matrix multiplication, and some for parallel > matrix vector multiplication. So I moved on to matrix*matrix > multiplication, but then printing it out I get only zeros as answer.. > but serial computation shows otherwise. Maybe you can point me out were > in code I have lost my soul? Thank you in advance! > PROGRAM parallel_mv
> ! global arrays
> real, dimension(16,16) :: a
> real, dimension(16,16) :: b,y
> ! variables for BLACS initialization
> !and processor grid creation
> integer iam,nprocs,ictxt,nprow,npcol,myrow,mycol
> !variables needed for distributing global
> !arrays across the proc grid
> integer desca(9), descb(9),descy(9),m,n,mb,nb,rsrc,csrc
> integer llda,lldb,info
> ! local arrays
> real, dimension(8,8) :: la
> ! local array for matrix b
> real, dimension(8,8) :: lb
> ! local array for matrix y
> real, dimension(8,8) ::ly
> !Initializing the BLACS library (STEP 2)
> call blacs_pinfo (iam,nprocs)
> call blacs_get(-1,0,ictxt)
> !Creating and using the processor grid (STEP 3)
> nprow=2; npcol=2
> call blacs_gridinit(ictxt,'r',nprow,npcol)
> call blacs_gridinfo(ictxt,nprow,npcol,myrow,mycol)
> ! Making the array descriptor vectors (STEP 4)
> m=16; n=16
> mb=8; nb=8
> rsrc=0; csrc=0
> llda=8
> call descinit(desca,m,n,mb,nb,rsrc,csrc,ictxt,llda,info)
> m=16; n=16
> mb=8; nb=8
> rsrc=0; csrc=0
> lldb=8
> call descinit(descb,m,n,mb,nb,rsrc,csrc,ictxt,lldb,info)
> call descinit(descy,m,n,mb,nb,rsrc,csrc,ictxt,lldb,info)
> ! Filling the global arrays A,b
> open(unit=12,file="a.data")
> read(12,*) a
> open(unit=13,file="b.data")
> read(13,*) b
> ! Each processor fills in its
> !local arrays with correct elements
> ! from the global arrays (STEP 5)
> print *,"filling local matrix"
> if(myrow.eq.0.and.mycol.eq.0) then
> do i_loc=1,8
> do j_loc=1,8
> !fill a
> la(i_loc,j_loc)=a(i_loc,j_loc)
> !fill b
> la(i_loc,j_loc)=b(i_loc,j_loc)
> end do
> end do
> end if
> if(myrow.eq.1.and.mycol.eq.0) then
> do i_loc=1,8
> do j_loc=1,8
> !fill a
> la(i_loc,j_loc)=a(i_loc,j_loc)
> !fill b
> la(i_loc,j_loc)=b(i_loc,j_loc)
> end do
> end do
> end if
> if(myrow.eq.0.and.mycol.eq.1) then
> do i_loc=1,8
> do j_loc=1,8
> !fill a
> la(i_loc,j_loc)=a(i_loc,j_loc)
> !fill b
> la(i_loc,j_loc)=b(i_loc,j_loc)
> end do
> end do
> end if
> if(myrow.eq.1.and.mycol.eq.1) then
> do i_loc=1,8
> do j_loc=1,8
> !fill a
> la(i_loc,j_loc)=a(i_loc,j_loc)
> !fill b
> la(i_loc,j_loc)=b(i_loc,j_loc)
> end do
> end do
> end if
> ! Call the ScaLAPACK routine (STEP 6)
> n=4
> print *,"caling psgemm"
> call psgemm('n','n',n,n,n,1.,la,1,1,desca(1),lb,1,1,descb(1),
> & 0.,c,1,1,descy(1))
> ! Each processor prints out its part of the product
> !vector y (STEP 7)
> if(myrow.eq.0.and.mycol.eq.0) then
> do i=1,8
> do j=1,8
> print *,ly(i), ly(j)
> end do
> end do
> end if
> if(myrow.eq.1.and.mycol.eq.0) then
> do i=1,8
> do j=1,8
> print *,ly(i), ly(j)
> end do
> end do
> end if
> ! Release the proc grid and BLACS library (STEP 8)
> call blacs_gridexit(ictxt)
> call blacs_exit(0)
> STOP
> END
I think that the problem is in the way you read the matrix from file.
I don't think you can open and read files like that in MPI-1. Since
I'm not 100% sure (and I don't have time to check it now) I suggest to
print the values of your matrix before calling the product. I suggest
you to use the PDLAREAD subroutine that reads the matrix on the root
processor and distributes it. PDLAREAD is not included in ScaLAPACK
but you can find it here:
http://netlib.org/scalapack/examples/pdlaread.f
Another option could be to use MPI-2 that has parallel-IO features.
Hope it helps
Alfredo Buttari
----------
http://alfredobuttari.wordpress.com
|
|
|
|
|
|
Add to del.icio.us |
Digg this |
Stumble it |
Powered by Megasolutions Inc