6ff16a608b3c63d980b3d3e46ab4ff63c4841baa
[scilab.git] / scilab / modules / linear_algebra / src / fortran / DGELSY1.f
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