Re: LAPACK library usage in fortran 90
- From: kiwanuka <robert.kiwanuka@xxxxxxxxx>
- Date: Sat, 24 Oct 2009 07:37:44 -0700 (PDT)
On Oct 24, 1:45 pm, jwm <jwmwal...@xxxxxxxxx> wrote:
On Oct 24, 3:49 am, kiwanuka <robert.kiwan...@xxxxxxxxx> wrote:Two reasons:
Hi all,
Has anyone managed to successfully use LAPACK libraries? I have
compiles them as a normal user for use with computations such as the
singular value decomposition (http://www.netlib.org/lapack/lug/
node53.html) but for some unclear reason I keep getting zeroes! Here
is a more complete account of what I've done:
1) Successfully compiled lapack-3.1.1.tgz fromhttp://www.netlib.org/lapack/index.html
2) Set up an example matrix whose SVD I found using python (see below)
3) Attempted to use the SVD subroutines, description of use in
lapack-3.1.1/html/dgesdd.f.html and lapack-3.1.1/html/dgesvd.f.html,
to find the SVD for the same matrix as in (2). I'm using fortran 90,
although I don't mind getting it working in F77.
4) I couldn't compile with the normal link options e.g.
$gortran -o exe file.f90 -L/path/to/libraries/ -llappack -lblas
I also tried making soft links to /usr/local/bin for the libraries so
that it was one of the paths I could use, but it didn't help. So I
tried first to copy the libraries to my working directory, which
didn't work! and second to provide full library paths without link
options e.g
$gortran -o exe file.f90 ../lapack-3.1.1/lapack_LINUX.a ../
lapack-3.1.1/blas_LINUX.a
These paths lead to the where I compiled the files from. This managed
to compile without errors.
The trouble is that the Fortran routines return zeroes for all the
matrices U,S, and VT which shouldn't be the case. I'm not sure if the
problem is with my provision of variables to or the libraries
themselves.
I've attached an example a file compiled with the command:
$gfortran -o out svd.f90 /home/robert/libraries/lapack-3.1.1/
lapack_LINUX.a /home/robert/libraries/lapack-3.1.1/blas_LINUX.a
For python comparison:
$python svd_test.py
+++++++++++++++++++
Files used: I've copied the contents below as I haven't yet worked out
how to attach them for this group.
+++++++++++++++++++
svd.f90
++++++++
program prog
INCLUDE 'ABA_PARAM.INC'
DIMENSION xmat(9,9), XDirtemp(9,9)
dimension :: U(9,9), S(9,9), VT(9,9)
dimension :: IWORK(72)
include "mydata.inc"
print *,'xmat'
write(*,730)((xmat(i,j),j=1,ubound(xmat,2)),i=1,ubound(xmat,1))
XDirtemp = xmat
M = ubound(xmat,1); N = ubound(xmat,2)
LDA = max(1,M); LDU = M; LDVT = N
LWORK = -1 !3*9**2 + 360
call DGESVD( 'A','A', 9, 9, XDirtemp, 9, S, U, 9, VT, 9, WORK, -1,
IWORK, INFO )
print *,'INFO',INFO,'M,N',M,N
print *,'U',ubound(U,1),ubound(U,2)
write(*,730) ((U(i,j),j=1,ubound(U,2)),i=1,ubound(U,1))
print *,'S',ubound(xmat,1),ubound(xmat,2)
write(*,100) ((S(i,j),j=1,ubound(S,2)),i=1,ubound(S,1))
print *,'VT',ubound(VT,1),ubound(VT,2)
write(*,720) ((VT(i,j),j=1,ubound(VT,2)),i=1,ubound(VT,1))
! stop
100 Format (' ',6(F10.4,1X))
200 Format (' ',3(F10.4,1X))
300 Format (' ',9(ES12.4,1X))
400 Format (' ',2(F10.4,1X))
500 Format (' ',3(F10.4,1X))
600 Format (' ',4(F10.4,1X))
700 Format (' ',5(F10.4,1X))
710 Format (' ',7(F10.4,1X))
720 Format (' ',8(F10.4,1X))
730 Format (' ',9(F10.4,1X))
contains
end program prog
++++++++
ABA_PARAM.INC
++++++++
implicit double precision (a-h,o-z)
parameter (j_sys_Dimension = 2)
parameter( n_vec_Length = 544 )
parameter( maxblk = n_vec_Length )
parameter(i_ipm_sta = -6)
character*5 j_ipm_Error
parameter(j_ipm_Error = "Error")
parameter(j_ipm_Aborted = 20)
++++++++
mydata.inc - only the relevant data
++++++++
DATA ((xmat(i,j),j=1,9),i=1,9) &
/1.5403 , -3.1706, 3.1706, -6.3411, -3.1706 , 3.1706 ,
&
3.3790, -1.6895, -1.6895, &
0.2701, 1.4679, -6.8261, 2.6690 , -4..1571,
-6.8261, &
-1.4222, -2.2152, 3.6374, &
-0.2701, 0.4679, 1.5507, 0.0000, -0.0000,
-0.0000, &
0.3885, 0.3885, 0.3885, &
0.5403, -0.0000 , -0.3027, 1.0288, -0..0000,
0.0000, &
-0.0000 , 0.0000, 0.0000, &
0.2701 , 0.4679 , -0.2678, 0.0255, 1..1376,
-0.0000, &
0.0000, 0.0000, 0.0000, &
-0.2701, 0.4679, 0.5507, -0.0523, 0.2340,
1.1950 , &
-0.0000 , 0.0000, 0.0000, &
-0.2879, 0.0000, 0.2889 , -0.6111, -0..8461,
0.6491, &
1.6417, -0.0000 , -0.0000, &
0.1439, 0.2493 , -0.2647 , -0.5584, 0.3023 ,
0.4464, &
0.3053, 1.9438, -0.0000, &
0.1439, -0.2493, -0.2768, -0.5573, 0.3450,
0.5842, &
-0.7036, -0.3306, -0.0000/
++++++++
svd_test.py
++++++++
import numpy as np
#% SVD test
A=np.array([[1.5403 , -3.1706, 3.1706, -6.3411,
-3.1706 , 3.1706 , 3.3790, -1.6895, -1.6895],
[0.2701, 1.4679, -6.8261, 2.6690 , -4.1571,
-6.8261, -1.4222, -2.2152, 3.6374],
[-0.2701, 0.4679, 1.5507, 0.0000, -0.0000,
-0.0000, 0.3885, 0.3885, 0.3885],
[ 0.5403, -0.0000 , -0.3027, 1.0288, -0.0000,
0.0000, -0.0000 , 0.0000, 0.0000],
[0.2701 , 0.4679 , -0.2678, 0.0255, 1.1376,
-0.0000, 0.0000, 0.0000, 0.0000],
[ -0.2701, 0.4679, 0.5507, -0.0523, 0.2340,
1.1950 , -0.0000 , 0.0000, 0.0000],
[ -0.2879, 0.0000, 0.2889 , -0.6111, -0.8461,
0.6491, 1.6417, -0.0000 , -0.0000],
[ 0.1439, 0.2493 , -0.2647 , -0.5584, 0.3023 ,
0.4464, 0.3053, 1.9438, -0.0000],
[ 0.1439, -0.2493, -0.2768, -0.5573, 0.3450,
0.5842, -0.7036, -0.3306, -0.0000]]);
[u,s,vh] = np.linalg.svd(A)
print u
print s
print vh
++++++++
It seems that you're using Linux. Is there any reason for not using
the blas/lapack packages provided by your distribution?
One is I could not find *.a files on my Ubuntu 9.04. There are *.so in
the python paths used by scipy e.g. flapack.so and clapack.so. I'm not
sure if it is okay to use those for fortran.
The second is I'm only developing and testing functionality on my
laptop but I intend to use the routines on departmental machines where
I can only work as a normal user. I therefore need to be able to
install and use these libraries with that level of permission - some
machines have them and some don't. Not sure what versions are
installed where they exist. This thus serves as a proof of concept
before I can go and fiddle with those machines.
It seems that your DGESVD subroutine has an extra argument. Try adding
the following interface block to your code (maybe before the second
include statement):
INTERFACE
SUBROUTINE DGESVD(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
WORK, LWORK, INFO )
INTEGER, PARAMETER :: WP = KIND(1.0D0)
CHARACTER(LEN=1), INTENT(IN) :: JOBU, JOBVT
INTEGER, INTENT(IN) :: M, N, LDA, LDU, LDVT, LWORK
INTEGER, INTENT(OUT) :: INFO
REAL(WP), INTENT(OUT) :: S(*)
REAL(WP), INTENT(INOUT) :: A(LDA,*)
REAL(WP), INTENT(OUT) :: U(LDU,*), VT(LDVT,*), WORK(*)
END SUBROUTINE DGESVD
END INTERFACE
Thank you very much! The interface is what was needed! It also forced
other problems to be exposed during compilation.
Regards,
Robert
.
- Follow-Ups:
- Re: LAPACK library usage in fortran 90
- From: Richard Maine
- Re: LAPACK library usage in fortran 90
- References:
- LAPACK library usage in fortran 90
- From: kiwanuka
- Re: LAPACK library usage in fortran 90
- From: jwm
- LAPACK library usage in fortran 90
- Prev by Date: Re: LAPACK library usage in fortran 90
- Next by Date: Re: LAPACK library usage in fortran 90
- Previous by thread: Re: LAPACK library usage in fortran 90
- Next by thread: Re: LAPACK library usage in fortran 90
- Index(es):
Relevant Pages
|
