LU factorization

[Tamiris_Crepalde]

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

[gib]

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.