Home     |     .Net Programming    |     cSharp Home    |     Sql Server Home    |     Javascript / Client Side Development     |     Ajax Programming

Ruby on Rails Development     |     Perl Programming     |     C Programming Language     |     C++ Programming     |     IT Jobs

Python Programming Language     |     Laptop Suggestions?    |     TCL Scripting     |     Fortran Programming     |     Scheme Programming Language


 
 
Cervo Technologies
The Right Source to Outsource

MS Dynamics CRM 3.0

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.

if it was so simple :)
I have run serial version on same matrices and get real number.

On May 18, 4:42 pm, Carramba <u@example.net> wrote:

Hi,
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