problem with DGELSx: variables A and B unchanged

[alexgc]

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


[admin]

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

[akobotov]

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