Strange results using dgesv (lapack) via fortran 90

[cjm2176]

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

 

[Mark44]

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.

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]

You are solving a linear matrix equation Ax=b. Have you verified that A and b have the same values in both cases?

 

[cjm2176]

Yes they are the same, would it help if I posted these arrays?

 

[hotvette]

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.

 

[cjm2176]

the arrays are attached

 

[hotvette]

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]

DGESV worked perfectly for me. There must be some error in your program.

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