Demo entry 6643638

Doolittle

   

Submitted by anonymous on Sep 28, 2017 at 11:53
Language: Fortran. Code size: 2.5 kB.

program doolittle

implicit none
integer::ndim
real*8 ,dimension(:,:),allocatable :: amat,lmat,umat
real*8 ,dimension(:),allocatable :: bvec,xvec,yvec
integer::irow,m,n
ndim=9
allocate(amat(ndim,ndim),lmat(ndim,ndim),umat(ndim,ndim),bvec(ndim),xvec(ndim),yvec(ndim))
call load()
call calclu()
call calcy()
call calcx()

bvec=MATMUL(amat,xvec)
print"(9f8.3)",(bvec(irow),irow=1,9)
print *,"---------------recalculate vecter b done------------------------------"

contains
    subroutine load()
        open(file="./amat.txt",unit=10)
        read(10,*)amat
        close(10)
        open(file="./bvec.txt",unit=10)
        read(10,*)bvec
        close(10)
        amat=TRANSPOSE(amat)
        print"(9f8.3)",(amat(irow,:),irow=1,9)
        print *,"-----------------load matrix A done----------------------------"
        print"(9f8.3)",(bvec(irow),irow=1,9)
        print *,"-----------------load vecter b done----------------------------"
    end subroutine
    
    subroutine calclu()
        umat(1,:)=amat(1,:)
        lmat(:,1)=amat(:,1)/umat(1,1)
        do m=2,ndim
            do n=m,ndim
                umat(m,n)=amat(m,n)-DOT_PRODUCT(lmat(m,1:m-1),umat(1:m-1,n))
                lmat(n,m)=amat(n,m)-DOT_PRODUCT(lmat(n,1:m-1),umat(1:m-1,m))
                lmat(n,m)=lmat(n,m)/umat(m,m)
            end do
        end do
        print"(9f8.3)",(lmat(irow,:),irow=1,9)
        print *,"--------------calculate lower matrix done-------------------------------"
        print"(9f8.3)",(umat(irow,:),irow=1,9)
        print *,"--------------calculate upper matrix done-------------------------------"
    end subroutine
    
    subroutine calcy()
        yvec(1)=bvec(1)/lmat(1,1)
        do irow=2,ndim
            yvec(irow)=bvec(irow)-DOT_PRODUCT(lmat(irow,1:irow-1),yvec(1:irow-1))
            yvec(irow)=yvec(irow)/lmat(irow,irow)
        end do
        print"(9f8.3)",(yvec(irow),irow=1,9)
        print *,"---------------calculate vecter y done------------------------------"
    end subroutine
    
    subroutine calcx()
        xvec(ndim)=yvec(ndim)/umat(ndim,ndim)
        do m=1,ndim-1
            irow=ndim-m
            xvec(irow)=yvec(irow)-DOT_PRODUCT(umat(irow,irow+1:ndim),xvec(irow+1:ndim))
            xvec(irow)=xvec(irow)/umat(irow,irow)
        end do
        print"(9f8.3)",(xvec(irow),irow=1,9)
        print *,"---------------calculate vecter x done------------------------------"
    end subroutine

end program

This snippet took 0.00 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).