Fortran 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

Attachments

• 3 KB Views: 239
Related Programming and Computer Science News on Phys.org

Mark44

Mentor
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
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

Homework Helper
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:

cjm2176

the arrays are attached

Attachments

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

hotvette

Homework Helper
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
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')
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

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