Community
cancel
Showing results for
Did you mean:
Highlighted
Beginner
23 Views

## LU factorization

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:'

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

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

+form='formatted',iomsg=imsg)
if(erro.gt.0) write(*,10)erro
if(erro.ne.'') write(*,'(A)')imsg

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

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

+form='formatted')
if(erro.gt.0) write(*,10)erro
if(erro.ne.'') write(*,'(A)')imsg

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

Tags (1)
Highlighted
New Contributor II
23 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.