cgelsd problem at some sizes of the equation system

[simon]

Hello,

I use the driver CGELSD to find the solution of a complex equation system. The system has much more equations than unkwon.
My Problem is that for some number of equations/number of unkwon/nrhs the driver works well and I get the right solution. For some other combination I get one error message. The error depends on the actual sizes.

error messages:
- ** On entry to SGEMV parameter number 2 had an illegal value
- Speicherzugriffsfehler (segmentation fault)
- ** On entry to CLALSA parameter number 6 had an illegal value
- *** glibc detected *** error_search1: munmap_chunc(): invalid pointer: =xb7b80008 ***

In my test program "error_search1" I read the matrices A (src) and B (targ) from files, which I created before in matlab. Then I use "lwork=-1" to determine the sizes for the helper arrays. Then I use CGELSD to compute X (ws). At the end I write the solution to another file.

here is the program i use:

program error_search1

!! try to figure out the problem with cgelsd lapack driver
!! use sb_create_dummy_src_targ to create the files with the source and
!! the target points
!!
!! argument 1: filename to src file
!! argument 2: m of src
!! argument 3: n of src
!! argument 4: filename to targ file
!! argument 5: m of src
!! argument 6: n of src

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! declaration of variables
complex*8, dimension(:,:), allocatable :: src
complex*8, dimension(:,:), allocatable :: targ

complex*8, dimension(:,:), allocatable :: ws

integer :: msrc, nsrc, mtarg, ntarg

!! variables for cgelsd
integer :: m,n,nrhs,lda,ldb,rank,lwork,info
real, dimension(:), allocatable :: s, rwork
real :: rcond
complex*8, dimension(:), allocatable :: work
integer, dimension(:), allocatable :: iwork

!! helper
integer :: lrwork, liwork
character(100) :: arg
integer :: fileid = 20, stat

real :: time1, time2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! process command line parameters

!! allocate the memory for src and target
call getarg(2, arg)
read(unit=arg, fmt=*) msrc
print *, "msrc", msrc
call getarg(3, arg)
read(unit=arg, fmt=*) nsrc
print *, "nsrc", nsrc

allocate( src(msrc,nsrc) )

call getarg(5, arg)
read(unit=arg, fmt=*) mtarg
print *, "mtarg", mtarg
call getarg(6, arg)
read(unit=arg, fmt=*) ntarg
print *, "ntarg", ntarg

allocate( targ(mtarg, ntarg) )

allocate( ws(nsrc, ntarg) )

!! src matrix
call getarg(1, arg)
if (arg /= "") then
print *, "start loading src file ", arg
open(unit=fileid, file=arg, iostat=stat, status='old', access='direct', &
form='unformatted', action='read', recl=8*msrc*nsrc)
if (stat == 0) then
read (fileid, rec=1) src(1:,1:)
close (fileid)
else
print *, "Could not open file ", arg
print *, "STOP"
stop
end if
end if

!! targ matrix
call getarg(4, arg)
if (arg /= "") then
print *, "start loading targ file ", arg
open(unit=fileid, file=arg, iostat=stat, status='old', access='direct', &
form='unformatted', action='read', recl=8*mtarg*ntarg)
if (stat == 0) then
read (fileid, rec=1) targ(1:,1:)
close (fileid)
else
print *, "Could not open file ", arg
print *, "STOP"
stop
end if
end if

!! set first parameters
m = msrc
n = nsrc
nrhs = ntarg
lda = m
ldb = mtarg
allocate( s(min(m,n)) )
rcond = -1
lwork = -1
allocate( work(1), rwork(1), iwork(1) )

!! query for minimal dimensions for the helper arrays
call CGELSD(m,n,nrhs,src,lda,targ,ldb,s,rcond,rank,work,lwork,rwork,iwork,info)
lwork = int(work(1))
lrwork = int(rwork(1))
liwork = iwork(1)
deallocate( work, rwork, iwork )

call cpu_time(time1)
allocate( work(lwork), rwork(lrwork), iwork(liwork), stat=stat )
print *, "allocate status fuer hilfsarrays: ", stat

!! search the solution for the equation system
call CGELSD(m,n,nrhs,src,lda,targ,ldb,s,rcond,rank,work,lwork,rwork,iwork,info)

ws = targ(1:nsrc, 1:ntarg)
call cpu_time(time2)

print *, time2-time1, "s"

!! write ws to a file
open(unit=fileid, file="ws_res", iostat=stat, status='new', &
access='direct', &
form='unformatted', action='write', &
recl=8*size(ws) )
if (stat == 0) then
write (fileid, rec=1) ws(1:,1:)
close (fileid)
end if

deallocate(ws)
deallocate( work, rwork, iwork )
deallocate( s )
deallocate( src, targ )

end program error_search1

perhaps some of you have an idea, what could be the reason for such a strange behaviour or where I made a big mistake.

Thank you
simon
[/list]

[Julien Langou]

I looked quickly at your code and could not find any problem.
Can you give me an example of m, n, nrhs for which the code breaks?
Julien.

[simon]

these are two examples with the corresponding error message.

a)
m = 1000
n = 100
nrhs = 70
Speicherzugriffsfehler

b)
m = 1000
n = 100
nrhs = 100
** On entry to SGEMV parameter number 2 had an illegal value

[simon]

I get the feeling the helper arrays rwork and iwork are too small. If I allocate them with a bigger size e.g.

lrwork = 2*lwork
liwork = 2*lwork

the program works fine. The determination of lrwork and liwork is done by CGELSD with the option lwork = -1. Could it be that there is a bug, or did I make a mistake in my first call to CGELSD?

simon

[Julien Langou]

Here are my workspaces:

Code: Select all

M = 1000
N = 100
NRHS = 100
LWORK = 10200
LRWORK = 15776
LIWORK = 1700


Code: Select all

M = 1000
N = 100
NRHS = 70
LWORK = 7200
LRWORK = 13526
LIWORK = 1700


with my workspace query. Do you have the same?
If yes, we might have a problem in the workspace computation.
Julien.

[simon]

I get the same result from the workspace query.

If I compute lwork, lrwork and liwork by hand with the definitions for the minmal neccessary values given in the documentation I get the same result. So the query seems to work fine.

Could it be that the assumptions for the minmal neccessary helper array sizes are not general enough?

simon

[simon]

Now I use cgelss instead of cgelsd and everything works fine.