problem with DGELSx: variables A and B unchanged

Post here if you have a question about LAPACK or ScaLAPACK algorithm or data format

problem with DGELSx: variables A and B unchanged

Postby alexgc » Tue Dec 14, 2010 1:05 pm

hello!

I am trying to solve a linear least squares problem using LAPACK's specialized routines. Right after DGELS returns, the variables A and B are unchanged but info evaluates to 0. I am using PGI Fortran 7.2-5 with BLAS and LAPACK as released by PGI, on a Dell Precision T7400.

Has anyone else experienced this problem? Is there something I am not doing right? I even tried for a simple 2x2 matrix but the outcome was the same: A and B had the same value as before entering DGELS


Code: Select all
call DGELS(tr, 2, 2, 1, amat, 2, bvec, 2, ym, -1, infols)



Code: Select all
[alexgc@predeal SWE_NEW]$ ./sensobs
 RHS before
    1.000000000000000         1.000000000000000     
    1.000000000000000        -1.000000000000000     
 LHS before
    1.000000000000000     
    0.000000000000000     
 info (as returned by DGELS) =             0
 RHS after <=> factorization
    1.000000000000000         1.000000000000000     
    1.000000000000000        -1.000000000000000     
 LHS after <=> solution
    1.000000000000000     
    0.000000000000000     
FORTRAN STOP
alexgc
 
Posts: 1
Joined: Tue Dec 14, 2010 12:50 pm

Re: problem with DGELSx: variables A and B unchanged

Postby admin » Fri Dec 17, 2010 2:58 pm

Alex,
I am confused in the dimension of your matrices.
we need to know which matrix is which matrix

Code: Select all
*  NRHS    (input) INTEGER
*          The number of right hand sides, i.e., the number of
*          columns of the matrices B and X. NRHS >=0.


from your output your RHS has 2 colums, so the fourth arg NRHS should be two.

Also only B will be changed, A will stay unchanged.
Code: Select all
*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
*          On entry, the matrix B of right hand side vectors, stored
*          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
*          if TRANS = 'T'.
*          On exit, if INFO = 0, B is overwritten by the solution
*          vectors, stored columnwise


Below is a simple code that works, hope it helps
Julie

Code: Select all
      program main
      double precision A(2,2), B(2,1)
      integer m ,n, ldb, I, J

      integer info
      double precision WORK(66)
      lwork=-1
      A(1,1)=1
      A(2,1)=1
      A(1,2)=1
      A(2,2)=-1
         WRITE(*,*) "MATRIX A"
         WRITE(*,99997) ((A(I,J),J=1,2),I=1,2)
      B(1,1)=1
      B(2,1)=0
         WRITE(*,*) "MATRIX B"
         WRITE(*,99998) (B(I,1),I=1,2)
      lwork=66
      call dgels('T',2,2,1,A,2,B,2,work, lwork,info)
      write(*,*) "INFO for dgels", INFO
      IF (INFO .EQ.0) THEN
         write(*,*) "res vaut :"
         WRITE(*,99998) (B(I,1),I=1,2)
      ELSE
         write(*,*) "Problem in DGELS, INFO is not nul"
      END IF
99997 FORMAT ((3X,2(E16.8)))
99998 FORMAT ((3X,1(E16.8)))
      end
admin
Site Admin
 
Posts: 468
Joined: Wed Dec 08, 2004 7:07 pm

Re: problem with DGELSx: variables A and B unchanged

Postby akobotov » Tue Dec 21, 2010 9:42 pm

Hello Alex,

call DGELS(tr, 2, 2, 1, amat, 2, bvec, 2, ym, -1, infols)

You are supplying LWORK = -1. This means that you are doing workspace query. The function just writes recommended length of WORK in the first element of YM and immediately exits. For the workspace query case the result of unchanged matrices is expected. For your case the minimal length of YM and value of LWORK is 4.

From DGELS parameters description:
* LWORK (input) INTEGER
* The dimension of the array WORK.
* LWORK >= max( 1, MN + max( MN, NRHS ) ).
* For optimal performance,
* LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
* where MN = min(M,N) and NB is the optimum block size.
*
* If LWORK = -1, then a workspace query is assumed; the routine
* only calculates the optimal size of the WORK array, returns
* this value as the first entry of the WORK array, and no error
* message related to LWORK is issued by XERBLA.


W.B.R.,
Alexander
akobotov
 
Posts: 7
Joined: Wed Feb 03, 2010 7:38 am
Location: Intel Corp., Russia, Novosibirsk


Return to Algorithm / Data

Who is online

Users browsing this forum: No registered users and 0 guests

cron