- #1

Telemachus

- 835

- 30

Hi there. I'm trying to solve a linear system which I have constructed with the aim of learning how to do this in fortran. The idea is to solve:

##A\vec{x}=\vec{b}##

The thing is that when I call the subroutine degsv from lapack to solve the linear system, and then write the solution, it only prints the vector b, and I can't find what I am doing wrong. It should solve the system by calling dgesv, and then print x in b.

This is the code:

<<Moderator's note: please use code tags around your code>>

##A\vec{x}=\vec{b}##

The thing is that when I call the subroutine degsv from lapack to solve the linear system, and then write the solution, it only prints the vector b, and I can't find what I am doing wrong. It should solve the system by calling dgesv, and then print x in b.

This is the code:

Fortran:

```
program matrix
!... This program constructs a matrix A, a vector b, and then solves
!... the system Ax=b using the lapack package.
implicit none integer mxpts
parameter (mxpts=5000)
integer i,j,rows,cols,k,aux,pivot(3),ok
real*8 A(mxpts,mxpts),b(mxpts),c,x(mxpts)
open (unit=10,file='matrix.dat',status='unknown')
open (unit=20,file='rhs.dat',status='unknown')
!... Inputs
rows=3
cols=rows !square matrix
!... Coefficients
!... Initialization do i=1,rows
b=0.0d0
do j=1,cols
A(i,j)=0.0d0
enddo
enddo
!... Matrix preparation
do i=1,rows
do j=1,cols
k=i+j
aux=mod(k,3)
c=1.0d0
if(aux.eq.1) c=-1.0d0
! print*,k,aux,c
A(i,j)=(i+j)*j*c
print 200,i,j,A(i,j)
enddo
enddo
!... Vector preparation for the right hand side of A*x=b.
do i=1,rows
b(i)=i
print*,i,b(i)
enddo
200 format(4x,"A(",i1,",",i1,")=",1pg15.8) !,f10.8)
!... 4x: 4 espacios en blanco, i1 imprime entero y le asigna 1 espacio
!... 1pg10.5 imprime real, le asigna 10 espacios donde 5 son para
!... decimales.
!... Print Matrix
do i=1,rows
do j=1,cols
write (10,*)A(i,j)
enddo
enddo
! do i=1,rows
! print*,( A(i,j), j=1,cols)
! enddo
!... Find the solution x=A-1*b
call dgesv (rows, 1, A, 3, pivot, b, 3, ok)
c
c parameters in the order as they appear in the function call
c order of matrix A, number of right hand sides (b), matrix A,
c leading dimension of A, array that records pivoting,
c result vector b on entry, x on exit, leading dimension of b
c return value
c
c print the vector x
do i=1,rows
print*,b(i)
write (20,*) b(i)
end do
stop
end
```

<<Moderator's note: please use code tags around your code>>

Last edited by a moderator: