Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!

LU factorization

Tamiris_Crepalde
Beginner
105 Views

Hello,

I don't know if that kind of subject is allowed here, if not allowed can you indicate me some forum that sover these subject?

I'm doing a code for LU factorization but my code has some problems with large matrix. When I use for small matrices it works fine but when I try use for a matrix 118x118 the answer is wrong. Someone can help me?

 

      program LU
      

      implicit none

      character*256 :: filename
      character*256 :: imsg
      real          :: start,finish,t
      integer       :: ul,erro,dim,i,j,k,n,s
      real*8,dimension(:,:),allocatable :: a,l,u
      real*8,dimension(:),allocatable   :: x,y,b

      filename    = ''
      ul          = 3
      s           = 10
      erro        = 0
      dim         = 0

      !matriz A
      write(*,'(A)') 'Informe o nome do arquivo que contem a matriz A:'
      read(*,'(A)') filename

      write(*,'(A)') 'Informe a dimensao da matriz:'
      read(*,'(I)') dim

      allocate(a(dim,dim))
      allocate(l(dim,dim))
      allocate(u(dim,dim))
      a = 0.0
      l = 0.0
      u = 0.0


      open(ul,file=filename,status='old',action='read',iostat=erro,
     +form='formatted',iomsg=imsg)
      if(erro.gt.0) write(*,10)erro
      if(erro.ne.'') write(*,'(A)')imsg
      read(ul,*) a


      !matriz b
      write(*,'(A)') 'Informe o nome do arquivo que contem a matriz b:'
      read(*,'(A)') filename

      allocate(x(dim))
      allocate(y(dim))
      allocate(b(dim))
      x = 0.0
      y = 0.0
      b = 0.0


      open(ul,file=filename,status='old',action='read',iostat=erro,
     +form='formatted')
      if(erro.gt.0) write(*,10)erro
      if(erro.ne.'') write(*,'(A)')imsg
      read(ul,*) b
      
      u=a
            
      do i=1,dim
          l(i,i)=1
      enddo
             
      do k=1,dim
          do i=k+1,dim
              do j=k,dim
                  u(j,i)=u(j,i)-l(k,i)*u(j,k)
                  l(k,i) = u(k,i)/u(k,k)
              enddo
          enddo
      enddo

      do k=1,dim
          y(k)=b(k)
          do i=1,(dim-1)
              if(i.ne.k) y(k)=y(k)-l(i,k)*y(i)
          enddo
      enddo
      do k=dim,1,-1
          x(k)=y(k)
          do i=k+1,dim
              x(k)=x(k)-u(i,k)*x(i)
          enddo
          x(k)=x(k)/u(k,k)
      enddo
      
      
      write(*,20) x

      
      pause
      
      deallocate(A)
      deallocate(l)
      deallocate(u)
      deallocate(x)
      deallocate(y)
      deallocate(b)

   10 format(1X,'Ocorreu um erro ao abrir o arquivo. IOSTAT=', I6/)
   20 format(20000(F18.14,2X))
   40 format(F18.14 /)
   50 format(/ 3(F8.6,2X))    

      end program

 

thanks

0 Kudos
1 Reply
gib
New Contributor II
105 Views

I made a small change to your program, generating a(:,:) and b(:) as uniformly distributed random numbers (generated with random_number(R)), instead of reading them from files.

If I check the resulting solution x by computing y = ax like this:

y = matmul(a,x)

and comparing y with b, they are completely different.

 But if I do this:
    do i = 1,dim
        y(i) = 0
        do j = 1,dim
            y(i) = y(i) + a(j,i)*x(j)   ! not a(i,j)*x(j)
        enddo
        write(*,'(i4,2f12.8)') i,b(i),y(i)
    enddo

I get good agreement between b and y with dim ranging up to 1000.  It seems that you have the row and column indices of 'a' swapped.  In Fortran A(i,j) is the entry in row i, column j of the matrix A. 

 

Reply