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