Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Problem with fortran and lapack

  1. Nov 28, 2016 #1
    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:

    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: Nov 29, 2016
  2. jcsd
  3. Nov 28, 2016 #2

    Mark44

    Staff: Mentor

  4. Nov 28, 2016 #3
    Hello Mark. Thank you very much for your feedback. I'm not really sure about that, because its the first time I'm attempting to use lapack. I have set those parameters while following this example: http://physics.oregonstate.edu/~landaur/nacphy/lapack/codes/linear-f.html , which works well, but uses sgesv (I have assumed that for dgesv should be the same, and I had took a look at the lapack documentation, but as I said, I'm not an expert).

    I have tried as you said to put that parameter equal to zero, but I get this error when I run the program:

    Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

    Backtrace for this error:
    #0 0x7F9354A59E08
    #1 0x7F9354A58F90
    #2 0x7F93546AA4AF
    #3 0x7F9354F11F18
    #4 0x400F2D in MAIN__ at matrix.f:?
    Segmentation fault (core dumped)

    That parameter is associated with an output from lapack.
     
  5. Nov 28, 2016 #4
    Here are my sugggestions:

    1) The Lapack routine dgesv solves the linear system AX=B, where A and B are matrices. So, if you want to solve for the case where the right-hand side is a vector you should declare B as a two-dimensional array of size (n,1) where n is the number of rows.
    2) Before you call dgesv you need to compute the LU factorization of A by using the routine DGETRF. This routine produces an array which contains the pivot indices. In your case it is the variable called pivot. But, I suppose that for the moment this array only contains zeros because you don't call DGETRF.
    3) As Mark44 said please print out and post the value of the variable "ok".
     
  6. Nov 29, 2016 #5

    Mark44

    Staff: Mentor

    No, the 6th parameter is an output variable. You don't set it -- the function sets it, with 0 indicating a successful completion. If the value is negative or positive, the value gives useful information, as explained in the documentation page I linked to.

    From the page in the link:
    SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
    Of these six parameters, A and B are input/output, which means they must be set when the subroutine is called, but can be changed by the subroutine.
    IPIV and INFO are output parameters, which means they can be anything when the subroutine is called, but get set by the subroutine.
     
  7. Nov 29, 2016 #6

    DrClaude

    User Avatar

    Staff: Mentor

    This is not necessary. Since Fortran is "pass by reference," the address of the first element is passed to the routine, and the shape of b doesn't matter.

    That's incorrect. dgsev calls dgetrf itself on line 167.
     
  8. Nov 29, 2016 #7

    DrClaude

    User Avatar

    Staff: Mentor

    There is a first problem here: pivot has to be the same size as the actual matrix A, so you should declare it as size mxpts.

    The dimensions here are all wrong.
    SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )

    The first 3 parameters are fine, but LDA = 3 is incorrect. LDA is the distance in memory between the the first element in the column and the first element in the column next to it. In other words, it is the first dimension of A, so mxpts.

    As I said above, IPIV is of size N, which is not how you declared pivot.

    As for LDA, LDB should be the first dimension of B, so also mxpts.

    This is all explained in the comments at the top of the LAPACK routine. You can also check this information on netlib.
     
  9. Nov 29, 2016 #8
    I've found the problem, the thing seems to be that I have specified the matrix to be from a size which doesn't correspond to the size I was trying to solve. I have changed the value of the parameter mxpts, and then used mxpts as the input to call dgesv.

    Thank you all.
     
  10. Nov 29, 2016 #9
    You are totally right. Thank you very much.
     
  11. Nov 30, 2016 #10
    I don't know if I should start another topic for this. The thing is that fortran after running dgesv to solve ##A\vec{x}=b## gives me this message:


    user@user-desktop:~/Desktop/radtransf$ ./matrix.x
    rmax= 10.000 npoints= 5000 Dr= 2.00000E-03
    Rows,hrows= 10000 5000
    Info= 0
    Note: The following floating-point exceptions are signalling: IEEE_DENORMAL

    What could be the problem?
     
  12. Nov 30, 2016 #11

    DrClaude

    User Avatar

    Staff: Mentor

    That seems to be some kind of underflow. What compiler are you using?

    Since you get back info = 0, I wouldn't worry too much. You can always check the result: Does Ax give you back b?
     
    Last edited: Nov 30, 2016
  13. Nov 30, 2016 #12
    I'm using gfortran. I'll try to do Ax to get b. Now one of the ##x_i## is wrong, and I know that because it is set to zero when constructing the matrix, its a boundary condition, but when lapack does its work it seems to change its value to some small but non zero number.
     
  14. Nov 30, 2016 #13

    DrClaude

    User Avatar

    Staff: Mentor

    From what I could find on the web, there was a problem with gfortran throwing floating point exceptions. I really wouldn't worry too much.


    You have to have some tolerance for numerical errors. What kind of numbers are you getting instead of 0?
     
  15. Nov 30, 2016 #14
    4.6338717745555646E-015

    With typical values beeing like: E-04 to E-06 the algorithm doesn't seems to be working. I'm trying to solve an integrodifferential equation, and I get nusty things.

    BTW: Have you ever worked with the ##S_N## approximation to the Boltzmann equation?

    After checking, its clearly not working.

    This is the code, i'll attach it with the input if you want to check.

    Code (Fortran):

            program radtransf                            !Nov 29 2016
     
     
    !...... This program calculates the solution for the radiative transfer
    !...... equation in the Sn approximation, with n=2. The program starts
    !...... by building a Matrix for the discretized system of coupled
    !...... differential equations, obtained by the Diamond difference scheme.


    !...... Version 0.0.1

            implicit none
     
            real*8 Dr,rmax
            integer npoints
            common/bkinput/ Dr,rmax,npoints


    !....... Input file
            open (unit=10,file='radtrans.inp',status='old')


    !....... call input
            call readinput

    !....... Matrix generation and resolution
            call matrix


    !...... Close files

            close (unit=10)
            close (unit=17)
            close (unit=18)
            close (unit=19)
            close (unit=20)
     
            stop
            end
     
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            subroutine readinput
            implicit none
     
    !...... Input parameters

     
            real*8 Dr,rmax
            integer npoints
            common/bkinput/ Dr,rmax,npoints
     
            real*8 rzero, one, two, five
            data rzero,one,two,five/0.0d0,1.0d0,2.0d0,5.0d0/
     
     
            namelist /griddata/Dr,npoints,rmax
         
            Dr=rzero
            npoints=0
            rmax=rzero
     
    !...... read griddata

            read(10,griddata)
            if (rmax.eq.rzero) rmax=Dr*npoints
            if (Dr.eq.0) Dr=rmax/npoints
            write(*,111) rmax,npoints,Dr
    111     format(/5x,'rmax=',f9.3,5x,'npoints=',i8,5x,'Dr=',1pg12.5)


            return
            end

    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


          subroutine matrix

    !.... This program constructs a matrix A, a vector b, and then solves
    !.... the system Ax=b using the lapack package.

     
          implicit none


          real*8 Dr,rmax
          integer npoints
          common/bkinput/ Dr,rmax,npoints
     
          integer mxpts
          parameter (mxpts=10000)
          integer i,j,rows,hrows,cols,k,aux,pivot(mxpts),info

     
          real*8 A(mxpts,mxpts),b(mxpts),v(mxpts),c,x,mu1,mu2


          real*8 one, two, four, hundred
          data one, two, four, hundred/ 1.0d0, 2.0d0, 4.0d0, 100.0d0/


    !..... Output files

    !      open (unit=16,file='rhs.dat',status='unknown')
          open (unit=17,file='matrix.dat',status='unknown')
          open (unit=18,file='Psi1.dat',status='unknown')
          open (unit=19,file='Psi2.dat',status='unknown')
          open (unit=20,file='check.dat',status='unknown')

    !..... Inputs
     
          rows=2*npoints
          hrows=npoints
          cols=rows !square matrix

          print*,'Rows,hrows=',rows,hrows

    !..... Coefficients
    !..... Initialization

            mu1=dsqrt(1.0d0/3.0d0)
            mu2=-1.0d0*dsqrt(1.0d0/3.0d0)

          do i=1,rows
            b=0.0d0
            do j=1,cols
              A(i,j)=0.0d0
           enddo
          enddo

    !..... Build matrix A

          do i=1,hrows
            do j=1,cols
            if(i.eq.1 .and. i.ne.j) then
              A(i,j)=0.0d0
            else if(i.eq.1 .and. i.eq.j) then
              A(i,j)=1.0d0
            else if(i.ne.1 .and. i.eq.j) then
              A(i,j)=two*mu1
            else if(i.ne. 1 .and. i-1 .eq. j) then
              A(i,j)=-two*mu1
            else if(i.ne.1 .and. i+j.eq.rows+1) then
              A(i,j)=-hundred*Dr
            else if(i.ne.1) then
              A(i,j)=0.0d0
            endif
           enddo
          enddo

          do i=hrows+1,rows
            do j=1,cols
            if(i.eq.hrows+1 .and. i.ne.j) then
              A(i,j)=0.0d0
            else if(i.eq.hrows+1 .and. i.eq.j) then
              A(i,j)=1.0d0
            else if(i.ne.hrows+1 .and. i.eq.j) then
              A(i,j)=-two*mu2
            else if(i.ne.hrows+1 .and. i-1.eq.j) then
              A(i,j)=two*mu2
            else if(i.ne.hrows+1 .and. i+j.eq.rows+1) then
              A(i,j)=-hundred*Dr
            else if(i.ne.hrows+1) then
              A(i,j)=0.0d0
            endif
           enddo
          enddo

    !        do i=1,rows
    !           do j=1,cols
    !                print 200,i,j,A(i,j)
    !           enddo
    !       enddo

    !..... Vector preparation for the right hand side of Ax=b.

            do i=1,rows
                     b(i)=0.01d0*Dr
            enddo

           b(1)=0.0d0
           b(hrows+1)=0.0d0


    !    do i=1,rows
    !      print*,b(i)
    !    end do


    !200   format(4x,"A(",i5,",",i5,")=",1pg15.8)  !,f10.8)


    !..... Print Matrix

            do i=1,rows
                    do j=1,cols
                        write (17,200)i,j,A(i,j)
                    enddo
            enddo

    200   format(4x,"A(",i5,",",i5,")=",1pg15.8)  !,f10.8)

    !      do i=1,rows
    !        print*,( A(i,j), j=1,cols)
    !      enddo

    !..... Find the solution x=A-1*b

            call dgesv (rows, 1, A, mxpts, pivot, b, mxpts, info)

            print*,"Info=",info


        do j=1,hrows+1
               x=(j-1)*Dr
               write(18,*)x,b(j)
        end do

     
        do i=hrows,rows
               x=(i-hrows)*Dr
               write(19,*)x,b(hrows+1+rows-i)
        end do

    !..... Verification for Ax=b

          do i=1,rows
            do j=1,cols
              v(i)=A(i,j)*b(j)
            write(20,*)v(i)
           enddo
          enddo

    !      do i=1,rows
    !        b(i)=0.01d0*Dr
    !        write(20,*)v(i),b(i)
    !      enddo

           return
           end

    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
    Ok, now I see that I was thinking too much symbolically.

    I think now the matrix product is ok:

    Code (Fortran):

    !..... Verification for Ax=b
    !..... Initialization

          do i=1,rows
            v(i)=0.0d0
          enddo

          rsum=0.0d0
          sumaux=0.0d0

    !..... Verification for Ax=b

          do i=1,rows
            do j=1,cols
              sumaux=A(i,j)*b(j)
              rsum=rsum+sumaux
              v(i)=rsum
           enddo
           write(20,*)v(i)
          enddo

     
    Anyway, only the first element of the array gives the correct value.
     

    Attached Files:

    Last edited: Nov 30, 2016
  16. Nov 30, 2016 #15
    Please, notice that the matrix A is changed by dgesv. So, to do the test you need to make a copy of A before calling that subroutine.
     
  17. Dec 1, 2016 #16

    Mark44

    Staff: Mentor

    @Telemachus, here's what I said in post #5:
    When you're using the lapack routines (or, in fact, the routines in any library), it's very important to read and understand the documenation for the routines.
     
  18. Dec 1, 2016 #17
    Thanks. And sorry for not giving the proper attention to those important recommendations, I am too anxious trying to make this work, so sorry for that, and thank you all for your patience, I'll try to not commit those mistakes again, and to be more lucid in the future.

    I have rebuilt the matrix after calling dgesv. Still not giving the correct right hand side after matrix multiplication. Could this be a consequence of the matrix A being ill conditioned? I haven't proven it is actually ill conditioned, but I think it could be the case.
     
    Last edited: Dec 1, 2016
  19. Dec 6, 2016 #18
    Hi again. I'm having additional problems with this code. When I try to build a bigger matrix, so I can get a finer grid, I have an upper bound. When I try to set the parameter mxpts>20000, I get this error when I try to compile it:


    :~/Desktop/radtransf$ gfortran radtrans.f -llapack -o radtrans.x
    /tmp/ccNJ2Cfs.o: In function `matrix_':
    radtrans.f:(.text+0x198): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
    radtrans.f:(.text+0x1a3): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
    radtrans.f:(.text+0x494): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
    radtrans.f:(.text+0x6d0): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
    radtrans.f:(.text+0x77f): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
    radtrans.f:(.text+0x922): relocation truncated to fit: R_X86_64_32 against `.bss'
    radtrans.f:(.text+0x9eb): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
    radtrans.f:(.text+0xacc): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
    radtrans.f:(.text+0xda8): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
    radtrans.f:(.text+0xfe4): relocation truncated to fit: R_X86_64_PC32 against symbol `bkinput_' defined in COMMON section in /tmp/ccNJ2Cfs.o
    radtrans.f:(.text+0x1098): additional relocation overflows omitted from the output

    If the parameter (which is used to store the arrays for matrix and vectors) is smaller than 20000, then the program compiles without any problem.
     
  20. Dec 6, 2016 #19
    Maybe, you don't have enough RAM memory. How much memory do you have on your computer?
    With the dimension 20 000, you will have 4*10^8 matrix elements, where each of them is a double-precision number. And with several arrays/matrices it gets even worse.
     
  21. Dec 6, 2016 #20
    The compilller knows that? I'm short of ram, 3.8gb.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Problem with fortran and lapack
  1. Fortran problem (Replies: 1)

  2. Fortran LAPACK DYSEV (Replies: 1)

Loading...