1      SUBROUTINE DGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
2     $                   WORK, INFO )
3*
4*  -- LAPACK driver routine (instrumented to count ops, version 3.0) --
5*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6*     Courant Institute, Argonne National Lab, and Rice University
7*     September 30, 1994
8*
9*     .. Scalar Arguments ..
10      INTEGER            INFO, LDA, LDB, M, N, NRHS, RANK
11      DOUBLE PRECISION   RCOND
12*     ..
13*     .. Array Arguments ..
14      INTEGER            JPVT( * )
15      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
16*     ..
17*     Common blocks to return operation counts and timings
18*     .. Common blocks ..
19      COMMON             / LATIME / OPS, ITCNT
20      COMMON             / LSTIME / OPCNT, TIMNG
21*     ..
22*     .. Scalars in Common ..
23      DOUBLE PRECISION   ITCNT, OPS
24*     ..
25*     .. Arrays in Common ..
26      DOUBLE PRECISION   OPCNT( 6 ), TIMNG( 6 )
27*     ..
28*
29*  Purpose
30*  =======
31*
32*  DGELSX computes the minimum-norm solution to a real linear least
33*  squares problem:
34*      minimize || A * X - B ||
35*  using a complete orthogonal factorization of A.  A is an M-by-N
36*  matrix which may be rank-deficient.
37*
38*  Several right hand side vectors b and solution vectors x can be
39*  handled in a single call; they are stored as the columns of the
40*  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
41*  matrix X.
42*
43*  The routine first computes a QR factorization with column pivoting:
44*      A * P = Q * [ R11 R12 ]
45*                  [  0  R22 ]
46*  with R11 defined as the largest leading submatrix whose estimated
47*  condition number is less than 1/RCOND.  The order of R11, RANK,
48*  is the effective rank of A.
49*
50*  Then, R22 is considered to be negligible, and R12 is annihilated
51*  by orthogonal transformations from the right, arriving at the
52*  complete orthogonal factorization:
53*     A * P = Q * [ T11 0 ] * Z
54*                 [  0  0 ]
55*  The minimum-norm solution is then
56*     X = P * Z' [ inv(T11)*Q1'*B ]
57*                [        0       ]
58*  where Q1 consists of the first RANK columns of Q.
59*
60*  Arguments
61*  =========
62*
63*  M       (input) INTEGER
64*          The number of rows of the matrix A.  M >= 0.
65*
66*  N       (input) INTEGER
67*          The number of columns of the matrix A.  N >= 0.
68*
69*  NRHS    (input) INTEGER
70*          The number of right hand sides, i.e., the number of
71*          columns of matrices B and X. NRHS >= 0.
72*
73*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
74*          On entry, the M-by-N matrix A.
75*          On exit, A has been overwritten by details of its
76*          complete orthogonal factorization.
77*
78*  LDA     (input) INTEGER
79*          The leading dimension of the array A.  LDA >= max(1,M).
80*
81*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
82*          On entry, the M-by-NRHS right hand side matrix B.
83*          On exit, the N-by-NRHS solution matrix X.
84*          If m >= n and RANK = n, the residual sum-of-squares for
85*          the solution in the i-th column is given by the sum of
86*          squares of elements N+1:M in that column.
87*
88*  LDB     (input) INTEGER
89*          The leading dimension of the array B. LDB >= max(1,M,N).
90*
91*  JPVT    (input/output) INTEGER array, dimension (N)
92*          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
93*          initial column, otherwise it is a free column.  Before
94*          the QR factorization of A, all initial columns are
95*          permuted to the leading positions; only the remaining
96*          free columns are moved as a result of column pivoting
97*          during the factorization.
98*          On exit, if JPVT(i) = k, then the i-th column of A*P
99*          was the k-th column of A.
100*
101*  RCOND   (input) DOUBLE PRECISION
102*          RCOND is used to determine the effective rank of A, which
103*          is defined as the order of the largest leading triangular
104*          submatrix R11 in the QR factorization with pivoting of A,
105*          whose estimated condition number < 1/RCOND.
106*
107*  RANK    (output) INTEGER
108*          The effective rank of A, i.e., the order of the submatrix
109*          R11.  This is the same as the order of the submatrix T11
110*          in the complete orthogonal factorization of A.
111*
112*  WORK    (workspace) DOUBLE PRECISION array, dimension
113*                      (max( min(M,N)+3*N, 2*min(M,N)+NRHS )),
114*
115*  INFO    (output) INTEGER
116*          = 0:  successful exit
117*          < 0:  if INFO = -i, the i-th argument had an illegal value
118*
119*  =====================================================================
120*
121*     .. Parameters ..
122      INTEGER            IMAX, IMIN
123      PARAMETER          ( IMAX = 1, IMIN = 2 )
124      DOUBLE PRECISION   ZERO, ONE, DONE, NTDONE
125      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, DONE = ZERO,
126     $                   NTDONE = ONE )
127*     ..
128*     .. Local Scalars ..
129      INTEGER            GELSX, GEQPF, I, IASCL, IBSCL, ISMAX, ISMIN, J,
130     $                   K, LATZM, MN, ORM2R, TRSM, TZRQF
131      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
132     $                   SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2, TIM1,
133     $                   TIM2
134*     ..
135*     .. External Functions ..
136      DOUBLE PRECISION   DLAMCH, DLANGE, DOPBL3, DOPLA, DSECND
137      EXTERNAL           DLAMCH, DLANGE, DOPBL3, DOPLA, DSECND
138*     ..
139*     .. External Subroutines ..
140      EXTERNAL           DGEQPF, DLABAD, DLAIC1, DLASCL, DLASET, DLATZM,
141     $                   DORM2R, DTRSM, DTZRQF, XERBLA
142*     ..
143*     .. Intrinsic Functions ..
144      INTRINSIC          ABS, DBLE, MAX, MIN
145*     ..
146*     .. Data statements ..
147      DATA               GELSX / 1 / , GEQPF / 2 / , LATZM / 6 / ,
148     $                   ORM2R / 4 / , TRSM / 5 / , TZRQF / 3 /
149*     ..
150*     .. Executable Statements ..
151*
152      MN = MIN( M, N )
153      ISMIN = MN + 1
154      ISMAX = 2*MN + 1
155*
156*     Test the input arguments.
157*
158      INFO = 0
159      IF( M.LT.0 ) THEN
160         INFO = -1
161      ELSE IF( N.LT.0 ) THEN
162         INFO = -2
163      ELSE IF( NRHS.LT.0 ) THEN
164         INFO = -3
165      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
166         INFO = -5
167      ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
168         INFO = -7
169      END IF
170*
171      IF( INFO.NE.0 ) THEN
172         CALL XERBLA( 'DGELSX', -INFO )
173         RETURN
174      END IF
175*
176*     Quick return if possible
177*
178      IF( MIN( M, N, NRHS ).EQ.0 ) THEN
179         RANK = 0
180         RETURN
181      END IF
182*
183*     Get machine parameters
184*
185      OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( 2 )
186      SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
187      BIGNUM = ONE / SMLNUM
188      CALL DLABAD( SMLNUM, BIGNUM )
189*
190*     Scale A, B if max elements outside range [SMLNUM,BIGNUM]
191*
192      ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
193      IASCL = 0
194      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
195*
196*        Scale matrix norm up to SMLNUM
197*
198         OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( M*N )
199         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
200         IASCL = 1
201      ELSE IF( ANRM.GT.BIGNUM ) THEN
202*
203*        Scale matrix norm down to BIGNUM
204*
205         OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( M*N )
206         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
207         IASCL = 2
208      ELSE IF( ANRM.EQ.ZERO ) THEN
209*
210*        Matrix all zero. Return zero solution.
211*
212         CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
213         RANK = 0
214         GO TO 100
215      END IF
216*
217      BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
218      IBSCL = 0
219      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
220*
221*        Scale matrix norm up to SMLNUM
222*
223         OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( M*NRHS )
224         CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
225         IBSCL = 1
226      ELSE IF( BNRM.GT.BIGNUM ) THEN
227*
228*        Scale matrix norm down to BIGNUM
229*
230         OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( M*NRHS )
231         CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
232         IBSCL = 2
233      END IF
234*
235*     Compute QR factorization with column pivoting of A:
236*        A * P = Q * R
237*
238      OPCNT( GEQPF ) = OPCNT( GEQPF ) + DOPLA( 'DGEQPF', M, N, 0, 0, 0 )
239      TIM1 = DSECND( )
240      CALL DGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), INFO )
241      TIM2 = DSECND( )
242      TIMNG( GEQPF ) = TIMNG( GEQPF ) + ( TIM2-TIM1 )
243*
244*     workspace 3*N. Details of Householder rotations stored
245*     in WORK(1:MN).
246*
247*     Determine RANK using incremental condition estimation
248*
249      WORK( ISMIN ) = ONE
250      WORK( ISMAX ) = ONE
251      SMAX = ABS( A( 1, 1 ) )
252      SMIN = SMAX
253      IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
254         RANK = 0
255         CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
256         GO TO 100
257      ELSE
258         RANK = 1
259      END IF
260*
261   10 CONTINUE
262      IF( RANK.LT.MN ) THEN
263         I = RANK + 1
264         OPS = 0
265         CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
266     $                A( I, I ), SMINPR, S1, C1 )
267         CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
268     $                A( I, I ), SMAXPR, S2, C2 )
269         OPCNT( GELSX ) = OPCNT( GELSX ) + OPS + DBLE( 1 )
270*
271         IF( SMAXPR*RCOND.LE.SMINPR ) THEN
272            OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( RANK*2 )
273            DO 20 I = 1, RANK
274               WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
275               WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
276   20       CONTINUE
277            WORK( ISMIN+RANK ) = C1
278            WORK( ISMAX+RANK ) = C2
279            SMIN = SMINPR
280            SMAX = SMAXPR
281            RANK = RANK + 1
282            GO TO 10
283         END IF
284      END IF
285*
286*     Logically partition R = [ R11 R12 ]
287*                             [  0  R22 ]
288*     where R11 = R(1:RANK,1:RANK)
289*
290*     [R11,R12] = [ T11, 0 ] * Y
291*
292      IF( RANK.LT.N ) THEN
293         OPCNT( TZRQF ) = OPCNT( TZRQF ) +
294     $                    DOPLA( 'DTZRQF', RANK, N, 0, 0, 0 )
295         TIM1 = DSECND( )
296         CALL DTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO )
297         TIM2 = DSECND( )
298         TIMNG( TZRQF ) = TIMNG( TZRQF ) + ( TIM2-TIM1 )
299      END IF
300*
301*     Details of Householder rotations stored in WORK(MN+1:2*MN)
302*
303*     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
304*
305      OPCNT( ORM2R ) = OPCNT( ORM2R ) +
306     $                 DOPLA( 'DORMQR', M, NRHS, MN, 0, 0 )
307      TIM1 = DSECND( )
308      CALL DORM2R( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ),
309     $             B, LDB, WORK( 2*MN+1 ), INFO )
310      TIM2 = DSECND( )
311      TIMNG( ORM2R ) = TIMNG( ORM2R ) + ( TIM2-TIM1 )
312*
313*     workspace NRHS
314*
315*     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
316*
317      OPCNT( TRSM ) = OPCNT( TRSM ) + DOPBL3( 'DTRSM ', RANK, NRHS, 0 )
318      TIM1 = DSECND( )
319      CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
320     $            NRHS, ONE, A, LDA, B, LDB )
321      TIM2 = DSECND( )
322      TIMNG( TRSM ) = TIMNG( TRSM ) + ( TIM2-TIM1 )
323*
324      DO 40 I = RANK + 1, N
325         DO 30 J = 1, NRHS
326            B( I, J ) = ZERO
327   30    CONTINUE
328   40 CONTINUE
329*
330*     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
331*
332      IF( RANK.LT.N ) THEN
333         OPCNT( LATZM ) = OPCNT( LATZM ) +
334     $                    DBLE( 2*( ( N-RANK )*NRHS+NRHS+( N-RANK )*
335     $                    NRHS )*RANK )
336         TIM1 = DSECND( )
337         DO 50 I = 1, RANK
338            CALL DLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA,
339     $                   WORK( MN+I ), B( I, 1 ), B( RANK+1, 1 ), LDB,
340     $                   WORK( 2*MN+1 ) )
341   50    CONTINUE
342         TIM2 = DSECND( )
343         TIMNG( LATZM ) = TIMNG( LATZM ) + ( TIM2-TIM1 )
344      END IF
345*
346*     workspace NRHS
347*
348*     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
349*
350      DO 90 J = 1, NRHS
351         DO 60 I = 1, N
352            WORK( 2*MN+I ) = NTDONE
353   60    CONTINUE
354         DO 80 I = 1, N
355            IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN
356               IF( JPVT( I ).NE.I ) THEN
357                  K = I
358                  T1 = B( K, J )
359                  T2 = B( JPVT( K ), J )
360   70             CONTINUE
361                  B( JPVT( K ), J ) = T1
362                  WORK( 2*MN+K ) = DONE
363                  T1 = T2
364                  K = JPVT( K )
365                  T2 = B( JPVT( K ), J )
366                  IF( JPVT( K ).NE.I )
367     $               GO TO 70
368                  B( I, J ) = T1
369                  WORK( 2*MN+K ) = DONE
370               END IF
371            END IF
372   80    CONTINUE
373   90 CONTINUE
374*
375*     Undo scaling
376*
377      IF( IASCL.EQ.1 ) THEN
378         OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( N*NRHS+RANK*RANK )
379         CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
380         CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
381     $                INFO )
382      ELSE IF( IASCL.EQ.2 ) THEN
383         OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( N*NRHS+RANK*RANK )
384         CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
385         CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
386     $                INFO )
387      END IF
388      IF( IBSCL.EQ.1 ) THEN
389         OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( N*NRHS )
390         CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
391      ELSE IF( IBSCL.EQ.2 ) THEN
392         OPCNT( GELSX ) = OPCNT( GELSX ) + DBLE( N*NRHS )
393         CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
394      END IF
395*
396  100 CONTINUE
397*
398      RETURN
399*
400*     End of DGELSX
401*
402      END
403