- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page