Fortran Strange results using dgesv (lapack) via fortran 90

  • Thread starter cjm2176
  • Start date
8
0
Hello all,

I am trying to learn fortran 90 by rewriting some simple matlab codes I have in fortran.

I tried to rewrite a linear, 1D finite element code for an elliptic equation and my fortran and matlab codes both end up assembling the same system matrixes (K and f), but the solution to the linear system I am getting by calling dgsev in fortran makes no sense.

heres what I get using dgsev
5.84805285111992524E-295
-5.75492768302598656E+018
-2.50280951573243828E-214
2.50555212078780657E-292
-0.0000000000000000
-0.0000000000000000
-9.67185004345343066E+025
-9.67140655691703340E+025
-0.0000000000000000

and from matlab
5.0000000000000000
4.3750000000000000
3.7499999999999991
3.1249999999999991
2.5000000000000000
1.8750000000000000
1.2500000000000000
0.62500000000000000
0.0000000000000000

I am lost on this one, any ideas?

Thanks!
cjm2176
 

Attachments

32,766
4,474
Your code is very difficult to read, at least using Notepad in Windows 7. I have pasted it into [ code] and [ /code] tags (without the leading spaces) below.
Code:
program FEMelliptic
implicit none

!declarations
real(kind = kind(0.0d0)), parameter :: T = 1.0d0, force = 0.0d0, L = 1.0d0, uL = 5.0d0, uR = 0.0d0
integer, parameter :: nnp = 9
integer, parameter :: nel = nnp - 1
integer :: node, I, el, j, info

real(kind = kind(0.0d0)), dimension(nnp,1) :: f, d, n_bc, P, node_x, e_bc, bc
real(kind = kind(0.0d0)), dimension(2,1) :: fe
real(kind = kind(0.0d0)), dimension(2,2) :: ke
real, dimension(nnp-1,2) :: connec  
real, dimension(nnp,nnp) :: K
real :: ans
real(kind = kind(0.0d0)), external :: ddot, dscal
integer, dimension(2,1) :: connec_el
integer, dimension(nnp) :: ipiv

!set bcs
bc(1,1) = 2.0d0
e_bc(1,1) = uL

bc(nnp,1) = 2.0d0
e_bc(nnp,1) = uR

!initialize system matrices
K(:,:) = 0.0d0
f(:,1) = 0.0d0
d(:,1) = 0.0d0

!Mesh information
!node coords
do node = 1, nnp
	node_x(node,1) = L*real(node-1)/(nnp-1)
end do
!connectivity
do I = 1, nel
	connec(I,1) = I
	connec(I,2) = I+1
end do

!assemble FE matrixes
do el = 1, nel
	connec_el(1,1) = connec(el,1)
	connec_el(2,1) = connec(el,2)
	call element_matrix(el, ke, fe)

	do j = 1, 2
		do i = 1, 2
			K(connec_el(i,1),connec_el(j,1)) = K(connec_el(i,1),connec_el(j,1)) + ke(i,j)
		end do
	end do
	
	do j = 1, 2
		f(connec_el(j,1),1) = f(connec_el(j,1),1) + fe(j,1)
	end do
end do

!apply bcs
do i = 1, nnp
	if (bc(I,1) == 2.0d0) then
		f(:,1) = f(:,1) - K(:,i)*e_bc(i,1)
		K(i,:) = 0.0d0
		K(:,i) = 0.0d0
		K(i,i) = 1.0d0
		f(i,1) = e_bc(i,1)
	end if
end do


!solve linear system
call dgesv(nnp, 1, K, nnp, ipiv, f, nnp, info)

write(*,*) f
 


contains
!compute element matrixes
subroutine element_matrix(el, ke, fe)
implicit none

integer, intent(in) :: el
real(kind = kind(0.0d0)), dimension(2,2), intent(out) :: ke
real(kind = kind(0.0d0)), dimension(2,1), intent(out) :: fe
real(kind = kind(0.0d0)), dimension(2,1) :: w, gp, xe
real(kind = kind(0.0d0)), dimension(1,2) :: B, N, ff
real(kind = kind(0.0d0)), dimension(2,2) :: c
real(kind = kind(0.0d0)) :: he, Jac, alpha
integer :: i, j
ke(1:2,1:2) = 0.0d0
fe(1:2,1) = 0.0d0

xe(1,1) = node_x(connec_el(1,1),1)
xe(2,1) = node_x(connec_el(2,1),1)
he = abs(xe(1,1)-xe(2,1))
Jac = he/2

w(1,1) = 1.0d0
w(2,1) = 1.0d0

gp(1,1) = -0.577350269189626d0
gp(2,1) = 0.577350269189626d0

do i = 1, 2
	call Nmatrix_lin (gp(i,1), N)
	call Bmatrix_lin (he, B)
	
	alpha = w(i,1)*T*Jac
	c = w(i,1)*T*Jac*matmul(transpose(B),B)
	ke = ke + c
	
	alpha = w(i,1)*force*Jac
	do j = 1, 2
		ff(1,j) = N(1,j)*alpha
	end do
	
	fe = fe + transpose(ff)
end do

end subroutine element_matrix
!shape functions
subroutine Nmatrix_lin (psi, N)
implicit none
real(kind = kind(0.0d0)), intent(in) :: psi
real(kind = kind(0.0d0)), dimension(1,2), intent(out) :: N

N(1,1) = .5d0*(1.0d0-psi)
N(1,2) = .5d0*(1.0d0+psi)

end subroutine Nmatrix_lin
!derivatives of shape functions
subroutine Bmatrix_lin (he, B)
implicit none
real(kind = kind(0.0d0)), intent(in) :: he
real(kind = kind(0.0d0)), dimension(1,2), intent(out) :: B

B(1,1) = -1.0d0/he
B(1,2) = 1.0d0/he

end subroutine Bmatrix_lin

end program FEMelliptic
 

hotvette

Homework Helper
982
3
You are solving a linear matrix equation Ax=b. Have you verified that A and b have the same values in both cases?
 
8
0
Yes they are the same, would it help if I posted these arrays?
 

hotvette

Homework Helper
982
3
Yes they are the same, would it help if I posted these arrays?
Sure, I have LAPACK routines and could check it. Tab or space delimited text file would be great.
 
Last edited:
8
0
the arrays are attached
 

Attachments

  • 1 KB Views: 239
  • 125 bytes Views: 224

hotvette

Homework Helper
982
3
Looks pretty straightforward (K is tri-diagonal). Just for kicks I used Excel matrix functions and got the expected answer. I'll try it tonight w/ LAPACK.
 

hotvette

Homework Helper
982
3
DGESV worked perfectly for me. There must be some error in your program.

Code:
PROGRAM PF
IMPLICIT NONE
INTEGER          NIN, NOUT, NMAX, LDA, I, IFAIL, INFO, J, N
PARAMETER        (NIN=5,NOUT=6, NMAX=10)
PARAMETER        (LDA=NMAX)
INTEGER          IPIV(NMAX)        
REAL (KIND=8)    A(LDA,NMAX), B(NMAX)
OPEN (NIN,file='PF.TXT',STATUS='OLD')
OPEN (NOUT,file='OUT.TXT',STATUS='OLD')
READ (NIN,*) N
READ (NIN,*) ((A(I,J),J=1,N),I=1,N)
READ (NIN,*) (B(I),I=1,N)
WRITE (NOUT,*) 'MATRIX A'
DO I=1,N
WRITE (NOUT,*) (A(I,J),J=1,N)
END DO
WRITE (NOUT,*) 'VECTOR B'
WRITE (NOUT,*) (B(I),I=1,N)
CALL DGESV(N,1,A,LDA,IPIV,B,N,INFO)
WRITE (NOUT,*) 'RESULT'
WRITE (NOUT,*) (B(I),I=1,N)
CLOSE (UNIT=NOUT)
CLOSE (UNIT=NIN)
!N =  9
!MATRIX A
!1. 0. 0. 0. 0. 0. 0. 0. 0.
!0. 16. -8. 0. 0. 0. 0. 0. 0.
!0. -8. 16. -8. 0. 0. 0. 0. 0.
!0. 0. -8. 16. -8. 0. 0. 0. 0.
!0. 0. 0. -8. 16. -8. 0. 0. 0.
!0. 0. 0. 0. -8. 16. -8. 0. 0.
!0. 0. 0. 0. 0. -8. 16. -8. 0.
!0. 0. 0. 0. 0. 0. -8. 16. 0.
!0. 0. 0. 0. 0. 0. 0. 0. 1.
!VECTOR B
!5. 40. 0. 0. 0. 0. 0. 0. 0.
!RESULT
!5. 4.375 3.7499999999999996 3.124999999999999 2.5 1.875 1.25 0.625 0.
END
 

Related Threads for: Strange results using dgesv (lapack) via fortran 90

Replies
22
Views
2K
  • Posted
Replies
1
Views
6K
Replies
3
Views
1K
Replies
5
Views
5K
Replies
4
Views
4K
  • Posted
Replies
1
Views
2K
  • Posted
Replies
7
Views
4K
  • Posted
Replies
1
Views
1K

Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving
Top