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