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