Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / cgelsy.f
1 *> \brief <b> CGELSY solves overdetermined or underdetermined systems for GE matrices</b>
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 *            http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CGELSY + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgelsy.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgelsy.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgelsy.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE CGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
22 *                          WORK, LWORK, RWORK, INFO )
23 *
24 *       .. Scalar Arguments ..
25 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26 *       REAL               RCOND
27 *       ..
28 *       .. Array Arguments ..
29 *       INTEGER            JPVT( * )
30 *       REAL               RWORK( * )
31 *       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
32 *       ..
33 *
34 *
35 *> \par Purpose:
36 *  =============
37 *>
38 *> \verbatim
39 *>
40 *> CGELSY computes the minimum-norm solution to a complex linear least
41 *> squares problem:
42 *>     minimize || A * X - B ||
43 *> using a complete orthogonal factorization of A.  A is an M-by-N
44 *> matrix which may be rank-deficient.
45 *>
46 *> Several right hand side vectors b and solution vectors x can be
47 *> handled in a single call; they are stored as the columns of the
48 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
49 *> matrix X.
50 *>
51 *> The routine first computes a QR factorization with column pivoting:
52 *>     A * P = Q * [ R11 R12 ]
53 *>                 [  0  R22 ]
54 *> with R11 defined as the largest leading submatrix whose estimated
55 *> condition number is less than 1/RCOND.  The order of R11, RANK,
56 *> is the effective rank of A.
57 *>
58 *> Then, R22 is considered to be negligible, and R12 is annihilated
59 *> by unitary transformations from the right, arriving at the
60 *> complete orthogonal factorization:
61 *>    A * P = Q * [ T11 0 ] * Z
62 *>                [  0  0 ]
63 *> The minimum-norm solution is then
64 *>    X = P * Z**H [ inv(T11)*Q1**H*B ]
65 *>                 [        0         ]
66 *> where Q1 consists of the first RANK columns of Q.
67 *>
68 *> This routine is basically identical to the original xGELSX except
69 *> three differences:
70 *>   o The permutation of matrix B (the right hand side) is faster and
71 *>     more simple.
72 *>   o The call to the subroutine xGEQPF has been substituted by the
73 *>     the call to the subroutine xGEQP3. This subroutine is a Blas-3
74 *>     version of the QR factorization with column pivoting.
75 *>   o Matrix B (the right hand side) is updated with Blas-3.
76 *> \endverbatim
77 *
78 *  Arguments:
79 *  ==========
80 *
81 *> \param[in] M
82 *> \verbatim
83 *>          M is INTEGER
84 *>          The number of rows of the matrix A.  M >= 0.
85 *> \endverbatim
86 *>
87 *> \param[in] N
88 *> \verbatim
89 *>          N is INTEGER
90 *>          The number of columns of the matrix A.  N >= 0.
91 *> \endverbatim
92 *>
93 *> \param[in] NRHS
94 *> \verbatim
95 *>          NRHS is INTEGER
96 *>          The number of right hand sides, i.e., the number of
97 *>          columns of matrices B and X. NRHS >= 0.
98 *> \endverbatim
99 *>
100 *> \param[in,out] A
101 *> \verbatim
102 *>          A is COMPLEX array, dimension (LDA,N)
103 *>          On entry, the M-by-N matrix A.
104 *>          On exit, A has been overwritten by details of its
105 *>          complete orthogonal factorization.
106 *> \endverbatim
107 *>
108 *> \param[in] LDA
109 *> \verbatim
110 *>          LDA is INTEGER
111 *>          The leading dimension of the array A.  LDA >= max(1,M).
112 *> \endverbatim
113 *>
114 *> \param[in,out] B
115 *> \verbatim
116 *>          B is COMPLEX array, dimension (LDB,NRHS)
117 *>          On entry, the M-by-NRHS right hand side matrix B.
118 *>          On exit, the N-by-NRHS solution matrix X.
119 *> \endverbatim
120 *>
121 *> \param[in] LDB
122 *> \verbatim
123 *>          LDB is INTEGER
124 *>          The leading dimension of the array B. LDB >= max(1,M,N).
125 *> \endverbatim
126 *>
127 *> \param[in,out] JPVT
128 *> \verbatim
129 *>          JPVT is INTEGER array, dimension (N)
130 *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
131 *>          to the front of AP, otherwise column i is a free column.
132 *>          On exit, if JPVT(i) = k, then the i-th column of A*P
133 *>          was the k-th column of A.
134 *> \endverbatim
135 *>
136 *> \param[in] RCOND
137 *> \verbatim
138 *>          RCOND is REAL
139 *>          RCOND is used to determine the effective rank of A, which
140 *>          is defined as the order of the largest leading triangular
141 *>          submatrix R11 in the QR factorization with pivoting of A,
142 *>          whose estimated condition number < 1/RCOND.
143 *> \endverbatim
144 *>
145 *> \param[out] RANK
146 *> \verbatim
147 *>          RANK is INTEGER
148 *>          The effective rank of A, i.e., the order of the submatrix
149 *>          R11.  This is the same as the order of the submatrix T11
150 *>          in the complete orthogonal factorization of A.
151 *> \endverbatim
152 *>
153 *> \param[out] WORK
154 *> \verbatim
155 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
156 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
157 *> \endverbatim
158 *>
159 *> \param[in] LWORK
160 *> \verbatim
161 *>          LWORK is INTEGER
162 *>          The dimension of the array WORK.
163 *>          The unblocked strategy requires that:
164 *>            LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
165 *>          where MN = min(M,N).
166 *>          The block algorithm requires that:
167 *>            LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
168 *>          where NB is an upper bound on the blocksize returned
169 *>          by ILAENV for the routines CGEQP3, CTZRZF, CTZRQF, CUNMQR,
170 *>          and CUNMRZ.
171 *>
172 *>          If LWORK = -1, then a workspace query is assumed; the routine
173 *>          only calculates the optimal size of the WORK array, returns
174 *>          this value as the first entry of the WORK array, and no error
175 *>          message related to LWORK is issued by XERBLA.
176 *> \endverbatim
177 *>
178 *> \param[out] RWORK
179 *> \verbatim
180 *>          RWORK is REAL array, dimension (2*N)
181 *> \endverbatim
182 *>
183 *> \param[out] INFO
184 *> \verbatim
185 *>          INFO is INTEGER
186 *>          = 0: successful exit
187 *>          < 0: if INFO = -i, the i-th argument had an illegal value
188 *> \endverbatim
189 *
190 *  Authors:
191 *  ========
192 *
193 *> \author Univ. of Tennessee
194 *> \author Univ. of California Berkeley
195 *> \author Univ. of Colorado Denver
196 *> \author NAG Ltd.
197 *
198 *> \date November 2011
199 *
200 *> \ingroup complexGEsolve
201 *
202 *> \par Contributors:
203 *  ==================
204 *>
205 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
206 *>    E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
207 *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
208 *>
209 *  =====================================================================
210       SUBROUTINE CGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
211      $                   WORK, LWORK, RWORK, INFO )
212 *
213 *  -- LAPACK driver routine (version 3.4.0) --
214 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
215 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216 *     November 2011
217 *
218 *     .. Scalar Arguments ..
219       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
220       REAL               RCOND
221 *     ..
222 *     .. Array Arguments ..
223       INTEGER            JPVT( * )
224       REAL               RWORK( * )
225       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
226 *     ..
227 *
228 *  =====================================================================
229 *
230 *     .. Parameters ..
231       INTEGER            IMAX, IMIN
232       PARAMETER          ( IMAX = 1, IMIN = 2 )
233       REAL               ZERO, ONE
234       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
235       COMPLEX            CZERO, CONE
236       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
237      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
238 *     ..
239 *     .. Local Scalars ..
240       LOGICAL            LQUERY
241       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
242      $                   NB, NB1, NB2, NB3, NB4
243       REAL               ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
244      $                   SMLNUM, WSIZE
245       COMPLEX            C1, C2, S1, S2
246 *     ..
247 *     .. External Subroutines ..
248       EXTERNAL           CCOPY, CGEQP3, CLAIC1, CLASCL, CLASET, CTRSM,
249      $                   CTZRZF, CUNMQR, CUNMRZ, SLABAD, XERBLA
250 *     ..
251 *     .. External Functions ..
252       INTEGER            ILAENV
253       REAL               CLANGE, SLAMCH
254       EXTERNAL           CLANGE, ILAENV, SLAMCH
255 *     ..
256 *     .. Intrinsic Functions ..
257       INTRINSIC          ABS, MAX, MIN, REAL, CMPLX
258 *     ..
259 *     .. Executable Statements ..
260 *
261       MN = MIN( M, N )
262       ISMIN = MN + 1
263       ISMAX = 2*MN + 1
264 *
265 *     Test the input arguments.
266 *
267       INFO = 0
268       NB1 = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
269       NB2 = ILAENV( 1, 'CGERQF', ' ', M, N, -1, -1 )
270       NB3 = ILAENV( 1, 'CUNMQR', ' ', M, N, NRHS, -1 )
271       NB4 = ILAENV( 1, 'CUNMRQ', ' ', M, N, NRHS, -1 )
272       NB = MAX( NB1, NB2, NB3, NB4 )
273       LWKOPT = MAX( 1, MN+2*N+NB*(N+1), 2*MN+NB*NRHS )
274       WORK( 1 ) = CMPLX( LWKOPT )
275       LQUERY = ( LWORK.EQ.-1 )
276       IF( M.LT.0 ) THEN
277          INFO = -1
278       ELSE IF( N.LT.0 ) THEN
279          INFO = -2
280       ELSE IF( NRHS.LT.0 ) THEN
281          INFO = -3
282       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
283          INFO = -5
284       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
285          INFO = -7
286       ELSE IF( LWORK.LT.( MN+MAX( 2*MN, N+1, MN+NRHS ) ) .AND.
287      $   .NOT.LQUERY ) THEN
288          INFO = -12
289       END IF
290 *
291       IF( INFO.NE.0 ) THEN
292          CALL XERBLA( 'CGELSY', -INFO )
293          RETURN
294       ELSE IF( LQUERY ) THEN
295          RETURN
296       END IF
297 *
298 *     Quick return if possible
299 *
300       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
301          RANK = 0
302          RETURN
303       END IF
304 *
305 *     Get machine parameters
306 *
307       SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
308       BIGNUM = ONE / SMLNUM
309       CALL SLABAD( SMLNUM, BIGNUM )
310 *
311 *     Scale A, B if max entries outside range [SMLNUM,BIGNUM]
312 *
313       ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
314       IASCL = 0
315       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
316 *
317 *        Scale matrix norm up to SMLNUM
318 *
319          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
320          IASCL = 1
321       ELSE IF( ANRM.GT.BIGNUM ) THEN
322 *
323 *        Scale matrix norm down to BIGNUM
324 *
325          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
326          IASCL = 2
327       ELSE IF( ANRM.EQ.ZERO ) THEN
328 *
329 *        Matrix all zero. Return zero solution.
330 *
331          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
332          RANK = 0
333          GO TO 70
334       END IF
335 *
336       BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK )
337       IBSCL = 0
338       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
339 *
340 *        Scale matrix norm up to SMLNUM
341 *
342          CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
343          IBSCL = 1
344       ELSE IF( BNRM.GT.BIGNUM ) THEN
345 *
346 *        Scale matrix norm down to BIGNUM
347 *
348          CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
349          IBSCL = 2
350       END IF
351 *
352 *     Compute QR factorization with column pivoting of A:
353 *        A * P = Q * R
354 *
355       CALL CGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ),
356      $             LWORK-MN, RWORK, INFO )
357       WSIZE = MN + REAL( WORK( MN+1 ) )
358 *
359 *     complex workspace: MN+NB*(N+1). real workspace 2*N.
360 *     Details of Householder rotations stored in WORK(1:MN).
361 *
362 *     Determine RANK using incremental condition estimation
363 *
364       WORK( ISMIN ) = CONE
365       WORK( ISMAX ) = CONE
366       SMAX = ABS( A( 1, 1 ) )
367       SMIN = SMAX
368       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
369          RANK = 0
370          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
371          GO TO 70
372       ELSE
373          RANK = 1
374       END IF
375 *
376    10 CONTINUE
377       IF( RANK.LT.MN ) THEN
378          I = RANK + 1
379          CALL CLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
380      $                A( I, I ), SMINPR, S1, C1 )
381          CALL CLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
382      $                A( I, I ), SMAXPR, S2, C2 )
383 *
384          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
385             DO 20 I = 1, RANK
386                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
387                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
388    20       CONTINUE
389             WORK( ISMIN+RANK ) = C1
390             WORK( ISMAX+RANK ) = C2
391             SMIN = SMINPR
392             SMAX = SMAXPR
393             RANK = RANK + 1
394             GO TO 10
395          END IF
396       END IF
397 *
398 *     complex workspace: 3*MN.
399 *
400 *     Logically partition R = [ R11 R12 ]
401 *                             [  0  R22 ]
402 *     where R11 = R(1:RANK,1:RANK)
403 *
404 *     [R11,R12] = [ T11, 0 ] * Y
405 *
406       IF( RANK.LT.N )
407      $   CALL CTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ),
408      $                LWORK-2*MN, INFO )
409 *
410 *     complex workspace: 2*MN.
411 *     Details of Householder rotations stored in WORK(MN+1:2*MN)
412 *
413 *     B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
414 *
415       CALL CUNMQR( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA,
416      $             WORK( 1 ), B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
417       WSIZE = MAX( WSIZE, 2*MN+REAL( WORK( 2*MN+1 ) ) )
418 *
419 *     complex workspace: 2*MN+NB*NRHS.
420 *
421 *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
422 *
423       CALL CTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
424      $            NRHS, CONE, A, LDA, B, LDB )
425 *
426       DO 40 J = 1, NRHS
427          DO 30 I = RANK + 1, N
428             B( I, J ) = CZERO
429    30    CONTINUE
430    40 CONTINUE
431 *
432 *     B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS)
433 *
434       IF( RANK.LT.N ) THEN
435          CALL CUNMRZ( 'Left', 'Conjugate transpose', N, NRHS, RANK,
436      $                N-RANK, A, LDA, WORK( MN+1 ), B, LDB,
437      $                WORK( 2*MN+1 ), LWORK-2*MN, INFO )
438       END IF
439 *
440 *     complex workspace: 2*MN+NRHS.
441 *
442 *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
443 *
444       DO 60 J = 1, NRHS
445          DO 50 I = 1, N
446             WORK( JPVT( I ) ) = B( I, J )
447    50    CONTINUE
448          CALL CCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
449    60 CONTINUE
450 *
451 *     complex workspace: N.
452 *
453 *     Undo scaling
454 *
455       IF( IASCL.EQ.1 ) THEN
456          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
457          CALL CLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
458      $                INFO )
459       ELSE IF( IASCL.EQ.2 ) THEN
460          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
461          CALL CLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
462      $                INFO )
463       END IF
464       IF( IBSCL.EQ.1 ) THEN
465          CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
466       ELSE IF( IBSCL.EQ.2 ) THEN
467          CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
468       END IF
469 *
470    70 CONTINUE
471       WORK( 1 ) = CMPLX( LWKOPT )
472 *
473       RETURN
474 *
475 *     End of CGELSY
476 *
477       END