6ff16a608b3c63d980b3d3e46ab4ff63c4841baa
1       SUBROUTINE DGELSY1( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
2      \$     WORK, LWORK, INFO )
3 *
4 *     -- LAPACK driver routine (version 3.0) --
5 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 *     Courant Institute, Argonne National Lab, and Rice University
7 *     June 30, 1999
8 *
9 *     .. Scalar Arguments ..
10       INTEGER            INFO, LDA, LDB, LWORK, 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 *
18 *     Purpose
19 *     =======
20 *
21 *     DGELSY1 computes a solution, with at least N-RANK zeros to a real
22 *     linear least squares problem:
23 *     minimize || A * X - B ||
24 *     using a complete orthogonal factorization of A.  A is an M-by-N
25 *     matrix which may be rank-deficient.
26 *
27 *     Several right hand side vectors b and solution vectors x can be
28 *     handled in a single call; they are stored as the columns of the
29 *     M-by-NRHS right hand side matrix B and the N-by-NRHS solution
30 *     matrix X.
31 *
32 *     The routine first computes a QR factorization with column
33 c     pivoting:
34 *     A * P = Q * [ R11 R12 ]
35 *                 [  0  R22 ]
36 *     with R11 defined as the largest leading submatrix whose estimated
37 *     condition number is less than 1/RCOND.  The order of R11, RANK,
38 *     is the effective rank of A.
39 *
40 *     Then, R22 is considered to be negligible,
41 *     The  solution return is then
42 *     X = P * [ inv(R11)*Q1'*B ]
43 *             [        0       ]
44 *     where Q1 consists of the first RANK columns of Q.
45 *
46 *     This routine is basically identical to the original xGELSX except
47 *     three differences:
48 *     o The call to the subroutine xGEQPF has been substituted by the
49 *     the call to the subroutine xGEQP3. This subroutine is a Blas-3
50 *     version of the QR factorization with column pivoting.
51 *     o Matrix B (the right hand side) is updated with Blas-3.
52 *     o The permutation of matrix B (the right hand side) is faster and
53 *     more simple.
54 *
55 *     Arguments
56 *     =========
57 *
58 *     M       (input) INTEGER
59 *     The number of rows of the matrix A.  M >= 0.
60 *
61 *     N       (input) INTEGER
62 *     The number of columns of the matrix A.  N >= 0.
63 *
64 *     NRHS    (input) INTEGER
65 *     The number of right hand sides, i.e., the number of
66 *     columns of matrices B and X. NRHS >= 0.
67 *
68 *     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
69 *     On entry, the M-by-N matrix A.
70 *     On exit, A has been overwritten by details of its
71 *     complete orthogonal factorization.
72 *
73 *     LDA     (input) INTEGER
74 *     The leading dimension of the array A.  LDA >= max(1,M).
75 *
76 *     B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS
77 c     )
78 *     On entry, the M-by-NRHS right hand side matrix B.
79 *     On exit, the N-by-NRHS solution matrix X.
80 *
81 *     LDB     (input) INTEGER
82 *     The leading dimension of the array B. LDB >= max(1,M,N).
83 *
84 *     JPVT    (input/output) INTEGER array, dimension (N)
85 *     On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
86 *     to the front of AP, otherwise column i is a free column.
87 *     On exit, if JPVT(i) = k, then the i-th column of AP
88 *     was the k-th column of A.
89 *
90 *     RCOND   (input) DOUBLE PRECISION
91 *     RCOND is used to determine the effective rank of A, which
92 *     is defined as the order of the largest leading triangular
93 *     submatrix R11 in the QR factorization with pivoting of A,
94 *     whose estimated condition number < 1/RCOND.
95 *
96 *     RANK    (output) INTEGER
97 *     The effective rank of A, i.e., the order of the submatrix
98 *     R11.  This is the same as the order of the submatrix T11
99 *     in the complete orthogonal factorization of A.
100 *
101 *     WORK    (workspace/output) DOUBLE PRECISION array, dimension
102 c     (LWORK)
103 *     On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
104 *
105 *     LWORK   (input) INTEGER
106 *     The dimension of the array WORK.
107 *     The unblocked strategy requires that:
108 *     LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
109 *     where MN = min( M, N ).
110 *     The block algorithm requires that:
111 *     LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
112 *     where NB is an upper bound on the blocksize returned
113 *     by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
114 *     and DORMRZ.
115 *
116 *     If LWORK = -1, then a workspace query is assumed; the routine
117 *     only calculates the optimal size of the WORK array, returns
118 *     this value as the first entry of the WORK array, and no error
119 *     message related to LWORK is issued by XERBLA.
120 *
121 *     INFO    (output) INTEGER
122 *     = 0: successful exit
123 *     < 0: If INFO = -i, the i-th argument had an illegal value.
124 *
125 *     Further Details
126 *     ===============
127 *
128 *     Based on contributions by
129 *     A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
130 *     E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I,
131 c     Spain
132 *     G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I,
133 c     Spain
134 *
135 *     ==================================================================
136 c     ===
137 *
138 *     .. Parameters ..
139       INTEGER            IMAX, IMIN
140       PARAMETER          ( IMAX = 1, IMIN = 2 )
141       DOUBLE PRECISION   ZERO, ONE
142       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
143 *     ..
144 *     .. Local Scalars ..
145       LOGICAL            LQUERY
146       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
147      \$     NB, NB1, NB2, NB3, NB4
148       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
149      \$     SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE
150 *     ..
151 *     .. External Functions ..
152       INTEGER            ILAENV
153       DOUBLE PRECISION   DLAMCH, DLANGE
154       EXTERNAL           ILAENV, DLAMCH, DLANGE
155 *     ..
156 *     .. External Subroutines ..
157       EXTERNAL           DCOPY, DGEQP3, DLABAD, DLAIC1, DLASCL, DLASET,
158      \$     DORMQR, DTRSM, XERBLA
159 *     ..
160 *     .. Intrinsic Functions ..
161       INTRINSIC          ABS, DBLE, MAX, MIN
162 *     ..
163 *     .. Executable Statements ..
164 *
165       MN = MIN( M, N )
166       ISMIN = MN + 1
167       ISMAX = 2*MN + 1
168 *
169 *     Test the input arguments.
170 *
171       INFO = 0
172       NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
173       NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
174       NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, NRHS, -1 )
175       NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, NRHS, -1 )
176       NB = MAX( NB1, NB2, NB3, NB4 )
177       LWKOPT = MAX( 1, MN+2*N+NB*( N+1 ), 2*MN+NB*NRHS )
178       WORK( 1 ) = DBLE( LWKOPT )
179       LQUERY = ( LWORK.EQ.-1 )
180       IF( M.LT.0 ) THEN
181          INFO = -1
182       ELSE IF( N.LT.0 ) THEN
183          INFO = -2
184       ELSE IF( NRHS.LT.0 ) THEN
185          INFO = -3
186       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
187          INFO = -5
188       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
189          INFO = -7
190       ELSE IF( LWORK.LT.MAX( 1, MN+3*N+1, 2*MN+NRHS ) .AND. .NOT.
191      \$        LQUERY ) THEN
192          INFO = -12
193       END IF
194 *
195       IF( INFO.NE.0 ) THEN
196          CALL XERBLA( 'DGELSY1', -INFO )
197          RETURN
198       ELSE IF( LQUERY ) THEN
199          RETURN
200       END IF
201 *
202 *     Quick return if possible
203 *
204       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
205          RANK = 0
206          RETURN
207       END IF
208 *
209 *     Get machine parameters
210 *
211       SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
212       BIGNUM = ONE / SMLNUM
213       CALL DLABAD( SMLNUM, BIGNUM )
214 *
215 *     Scale A, B if max entries outside range [SMLNUM,BIGNUM]
216 *
217       ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
218       IASCL = 0
219       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
220 *
221 *     Scale matrix norm up to SMLNUM
222 *
223          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
224          IASCL = 1
225       ELSE IF( ANRM.GT.BIGNUM ) THEN
226 *
227 *     Scale matrix norm down to BIGNUM
228 *
229          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
230          IASCL = 2
231       ELSE IF( ANRM.EQ.ZERO ) THEN
232 *
233 *     Matrix all zero. Return zero solution.
234 *
235          CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
236          RANK = 0
237          GO TO 70
238       END IF
239 *
240       BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
241       IBSCL = 0
242       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
243 *
244 *     Scale matrix norm up to SMLNUM
245 *
246          CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
247          IBSCL = 1
248       ELSE IF( BNRM.GT.BIGNUM ) THEN
249 *
250 *     Scale matrix norm down to BIGNUM
251 *
252          CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
253          IBSCL = 2
254       END IF
255 *
256 *     Compute QR factorization with column pivoting of A:
257 *     A * P = Q * R
258 *
259       CALL DGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ),
260      \$     LWORK-MN, INFO )
261       WSIZE = MN + WORK( MN+1 )
262 *
263 *     workspace: MN+2*N+NB*(N+1).
264 *     Details of Householder rotations stored in WORK(1:MN).
265 *
266 *     Determine RANK using incremental condition estimation
267 *
268       WORK( ISMIN ) = ONE
269       WORK( ISMAX ) = ONE
270       SMAX = ABS( A( 1, 1 ) )
271       SMIN = SMAX
272       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
273          RANK = 0
274          CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
275          GO TO 70
276       ELSE
277          RANK = 1
278       END IF
279 *
280  10   CONTINUE
281       IF( RANK.LT.MN ) THEN
282          I = RANK + 1
283          CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
284      \$        A( I, I ), SMINPR, S1, C1 )
285          CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
286      \$        A( I, I ), SMAXPR, S2, C2 )
287 *
288          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
289             DO 20 I = 1, RANK
290                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
291                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
292  20         CONTINUE
293             WORK( ISMIN+RANK ) = C1
294             WORK( ISMAX+RANK ) = C2
295             SMIN = SMINPR
296             SMAX = SMAXPR
297             RANK = RANK + 1
298             GO TO 10
299          END IF
300       END IF
301 *
302 *     workspace: 3*MN.
303 *
304 *     Logically partition R = [ R11 R12 ]
305 *     [  0  R22 ]
306 *     where R11 = R(1:RANK,1:RANK)
307 *
308 c     *     [R11,R12] = [ T11, 0 ] * Y
309 c     *
310 c     IF( RANK.LT.N )
311 c     \$   CALL DTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ),
312 c     \$                LWORK-2*MN, INFO )
313 c     *
314 c     *     workspace: 2*MN.
315 *     Details of Householder rotations stored in WORK(MN+1:2*MN)
316 *
317 *     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
318 *
319       CALL DORMQR( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ),
320      \$     B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
321       WSIZE = MAX( WSIZE, 2*MN+WORK( 2*MN+1 ) )
322 *
323 *     workspace: 2*MN+NB*NRHS.
324 *
325 *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
326 *
327       CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
328      \$     NRHS, ONE, A, LDA, B, LDB )
329 *
330       DO 40 J = 1, NRHS
331          DO 30 I = RANK + 1, N
332             B( I, J ) = ZERO
333  30      CONTINUE
334  40   CONTINUE
335 c     *
336 c     *     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
337 c     *
338 c     IF( RANK.LT.N ) THEN
339 c     CALL DORMRZ( 'Left', 'Transpose', N, NRHS, RANK, N-RANK, A,
340 c     \$                LDA, WORK( MN+1 ), B, LDB, WORK( 2*MN+1 ),
341 c     \$                LWORK-2*MN, INFO )
342 c     END IF
343 c     *
344 c     *     workspace: 2*MN+NRHS.
345 *
346 *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
347 *
348       DO 60 J = 1, NRHS
349          DO 50 I = 1, N
350             WORK( JPVT( I ) ) = B( I, J )
351  50      CONTINUE
352          CALL DCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
353  60   CONTINUE
354 *
355 *     workspace: N.
356 *
357 *     Undo scaling
358 *
359       IF( IASCL.EQ.1 ) THEN
360          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
361          CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
362      \$        INFO )
363       ELSE IF( IASCL.EQ.2 ) THEN
364          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
365          CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
366      \$        INFO )
367       END IF
368       IF( IBSCL.EQ.1 ) THEN
369          CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
370       ELSE IF( IBSCL.EQ.2 ) THEN
371          CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
372       END IF
373 *
374  70   CONTINUE
375       WORK( 1 ) = DBLE( LWKOPT )
376 *
377       RETURN
378 *
379 *     End of DGELSY1
380 *
381       END