Strange results using dgesv (lapack) via fortran 90

Click For Summary

Discussion Overview

The discussion revolves around discrepancies observed when solving a linear system using the LAPACK routine dgesv in Fortran 90, compared to results obtained from MATLAB. The context involves rewriting a finite element code for an elliptic equation, focusing on the assembly of system matrices and the subsequent solution of the linear system.

Discussion Character

  • Technical explanation
  • Debate/contested
  • Mathematical reasoning

Main Points Raised

  • One participant reports unexpected results from the dgesv routine, providing both the Fortran and MATLAB output for comparison.
  • Another participant suggests verifying that the matrices A and b are identical in both implementations.
  • A later reply indicates that the system matrix K is tri-diagonal and mentions successful results using Excel matrix functions.
  • One participant claims that dgesv worked perfectly for them, implying there may be an error in the original program.
  • Another participant offers to check the arrays if they are provided in a suitable format.
  • One participant shares a simplified program that successfully uses dgesv, including sample input and output data for clarity.

Areas of Agreement / Disagreement

Participants express differing experiences with the dgesv routine, with some achieving expected results while others encounter significant discrepancies. The discussion remains unresolved as participants explore potential errors in the original code.

Contextual Notes

There are indications of potential issues related to the readability of the code, the format of the input data, and the handling of boundary conditions, but these aspects have not been conclusively addressed.

cjm2176
Messages
8
Reaction score
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

Technology news on Phys.org
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
 
You are solving a linear matrix equation Ax=b. Have you verified that A and b have the same values in both cases?
 
Yes they are the same, would it help if I posted these arrays?
 
cjm2176 said:
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:
the arrays are attached
 

Attachments

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

Similar threads

Fortran Problem with fortran and lapack
  • · Replies 22 ·
Replies
22
Views
5K
Fortran Creating a distribution with specific mean and variance in FORTRAN 90
  • · Replies 3 ·
Replies
3
Views
4K
MATLAB Fortran 90 vs Matlab: Comparing Program Speed
  • · Replies 5 ·
Replies
5
Views
11K
Undergrad FORTRAN 90 - DVERK? - initial condition solver
  • · Replies 3 ·
Replies
3
Views
6K
MATLAB Using Matlab functions in Fortran -
  • · Replies 1 ·
Replies
1
Views
4K