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


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
 info (as returned by DGELS) =             0
 RHS after <=> factorization
    1.000000000000000         1.000000000000000     
    1.000000000000000        -1.000000000000000     
 LHS after <=> solution
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

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

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)
         WRITE(*,*) "MATRIX A"
         WRITE(*,99997) ((A(I,J),J=1,2),I=1,2)
         WRITE(*,*) "MATRIX B"
         WRITE(*,99998) (B(I,1),I=1,2)
      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)
         write(*,*) "Problem in DGELS, INFO is not nul"
      END IF
99997 FORMAT ((3X,2(E16.8)))
99998 FORMAT ((3X,1(E16.8)))
Site Admin
Posts: 546
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:
* 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.

Posts: 11
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 2 guests