 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 00002 $ WORK, LWORK, RWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00011 REAL RCOND 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER JPVT( * ) 00015 REAL RWORK( * ) 00016 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * CGELSY computes the minimum-norm solution to a complex linear least 00023 * squares problem: 00024 * minimize || A * X - B || 00025 * using a complete orthogonal factorization of A. A is an M-by-N 00026 * matrix which may be rank-deficient. 00027 * 00028 * Several right hand side vectors b and solution vectors x can be 00029 * handled in a single call; they are stored as the columns of the 00030 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution 00031 * matrix X. 00032 * 00033 * The routine first computes a QR factorization with column pivoting: 00034 * A * P = Q * [ R11 R12 ] 00035 * [ 0 R22 ] 00036 * with R11 defined as the largest leading submatrix whose estimated 00037 * condition number is less than 1/RCOND. The order of R11, RANK, 00038 * is the effective rank of A. 00039 * 00040 * Then, R22 is considered to be negligible, and R12 is annihilated 00041 * by unitary transformations from the right, arriving at the 00042 * complete orthogonal factorization: 00043 * A * P = Q * [ T11 0 ] * Z 00044 * [ 0 0 ] 00045 * The minimum-norm solution is then 00046 * X = P * Z**H [ inv(T11)*Q1**H*B ] 00047 * [ 0 ] 00048 * where Q1 consists of the first RANK columns of Q. 00049 * 00050 * This routine is basically identical to the original xGELSX except 00051 * three differences: 00052 * o The permutation of matrix B (the right hand side) is faster and 00053 * more simple. 00054 * o The call to the subroutine xGEQPF has been substituted by the 00055 * the call to the subroutine xGEQP3. This subroutine is a Blas-3 00056 * version of the QR factorization with column pivoting. 00057 * o Matrix B (the right hand side) is updated with Blas-3. 00058 * 00059 * Arguments 00060 * ========= 00061 * 00062 * M (input) INTEGER 00063 * The number of rows of the matrix A. M >= 0. 00064 * 00065 * N (input) INTEGER 00066 * The number of columns of the matrix A. N >= 0. 00067 * 00068 * NRHS (input) INTEGER 00069 * The number of right hand sides, i.e., the number of 00070 * columns of matrices B and X. NRHS >= 0. 00071 * 00072 * A (input/output) COMPLEX array, dimension (LDA,N) 00073 * On entry, the M-by-N matrix A. 00074 * On exit, A has been overwritten by details of its 00075 * complete orthogonal factorization. 00076 * 00077 * LDA (input) INTEGER 00078 * The leading dimension of the array A. LDA >= max(1,M). 00079 * 00080 * B (input/output) COMPLEX array, dimension (LDB,NRHS) 00081 * On entry, the M-by-NRHS right hand side matrix B. 00082 * On exit, the N-by-NRHS solution matrix X. 00083 * 00084 * LDB (input) INTEGER 00085 * The leading dimension of the array B. LDB >= max(1,M,N). 00086 * 00087 * JPVT (input/output) INTEGER array, dimension (N) 00088 * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted 00089 * to the front of AP, otherwise column i is a free column. 00090 * On exit, if JPVT(i) = k, then the i-th column of A*P 00091 * was the k-th column of A. 00092 * 00093 * RCOND (input) REAL 00094 * RCOND is used to determine the effective rank of A, which 00095 * is defined as the order of the largest leading triangular 00096 * submatrix R11 in the QR factorization with pivoting of A, 00097 * whose estimated condition number < 1/RCOND. 00098 * 00099 * RANK (output) INTEGER 00100 * The effective rank of A, i.e., the order of the submatrix 00101 * R11. This is the same as the order of the submatrix T11 00102 * in the complete orthogonal factorization of A. 00103 * 00104 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00105 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00106 * 00107 * LWORK (input) INTEGER 00108 * The dimension of the array WORK. 00109 * The unblocked strategy requires that: 00110 * LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS ) 00111 * where MN = min(M,N). 00112 * The block algorithm requires that: 00113 * LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS ) 00114 * where NB is an upper bound on the blocksize returned 00115 * by ILAENV for the routines CGEQP3, CTZRZF, CTZRQF, CUNMQR, 00116 * and CUNMRZ. 00117 * 00118 * If LWORK = -1, then a workspace query is assumed; the routine 00119 * only calculates the optimal size of the WORK array, returns 00120 * this value as the first entry of the WORK array, and no error 00121 * message related to LWORK is issued by XERBLA. 00122 * 00123 * RWORK (workspace) REAL array, dimension (2*N) 00124 * 00125 * INFO (output) INTEGER 00126 * = 0: successful exit 00127 * < 0: if INFO = -i, the i-th argument had an illegal value 00128 * 00129 * Further Details 00130 * =============== 00131 * 00132 * Based on contributions by 00133 * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA 00134 * E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 00135 * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 00136 * 00137 * ===================================================================== 00138 * 00139 * .. Parameters .. 00140 INTEGER IMAX, IMIN 00141 PARAMETER ( IMAX = 1, IMIN = 2 ) 00142 REAL ZERO, ONE 00143 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00144 COMPLEX CZERO, CONE 00145 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00146 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00147 * .. 00148 * .. Local Scalars .. 00149 LOGICAL LQUERY 00150 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN, 00151 $ NB, NB1, NB2, NB3, NB4 00152 REAL ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR, 00153 $ SMLNUM, WSIZE 00154 COMPLEX C1, C2, S1, S2 00155 * .. 00156 * .. External Subroutines .. 00157 EXTERNAL CCOPY, CGEQP3, CLAIC1, CLASCL, CLASET, CTRSM, 00158 $ CTZRZF, CUNMQR, CUNMRZ, SLABAD, XERBLA 00159 * .. 00160 * .. External Functions .. 00161 INTEGER ILAENV 00162 REAL CLANGE, SLAMCH 00163 EXTERNAL CLANGE, ILAENV, SLAMCH 00164 * .. 00165 * .. Intrinsic Functions .. 00166 INTRINSIC ABS, MAX, MIN, REAL, CMPLX 00167 * .. 00168 * .. Executable Statements .. 00169 * 00170 MN = MIN( M, N ) 00171 ISMIN = MN + 1 00172 ISMAX = 2*MN + 1 00173 * 00174 * Test the input arguments. 00175 * 00176 INFO = 0 00177 NB1 = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) 00178 NB2 = ILAENV( 1, 'CGERQF', ' ', M, N, -1, -1 ) 00179 NB3 = ILAENV( 1, 'CUNMQR', ' ', M, N, NRHS, -1 ) 00180 NB4 = ILAENV( 1, 'CUNMRQ', ' ', M, N, NRHS, -1 ) 00181 NB = MAX( NB1, NB2, NB3, NB4 ) 00182 LWKOPT = MAX( 1, MN+2*N+NB*(N+1), 2*MN+NB*NRHS ) 00183 WORK( 1 ) = CMPLX( LWKOPT ) 00184 LQUERY = ( LWORK.EQ.-1 ) 00185 IF( M.LT.0 ) THEN 00186 INFO = -1 00187 ELSE IF( N.LT.0 ) THEN 00188 INFO = -2 00189 ELSE IF( NRHS.LT.0 ) THEN 00190 INFO = -3 00191 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00192 INFO = -5 00193 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 00194 INFO = -7 00195 ELSE IF( LWORK.LT.( MN+MAX( 2*MN, N+1, MN+NRHS ) ) .AND. 00196 $ .NOT.LQUERY ) THEN 00197 INFO = -12 00198 END IF 00199 * 00200 IF( INFO.NE.0 ) THEN 00201 CALL XERBLA( 'CGELSY', -INFO ) 00202 RETURN 00203 ELSE IF( LQUERY ) THEN 00204 RETURN 00205 END IF 00206 * 00207 * Quick return if possible 00208 * 00209 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 00210 RANK = 0 00211 RETURN 00212 END IF 00213 * 00214 * Get machine parameters 00215 * 00216 SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) 00217 BIGNUM = ONE / SMLNUM 00218 CALL SLABAD( SMLNUM, BIGNUM ) 00219 * 00220 * Scale A, B if max entries outside range [SMLNUM,BIGNUM] 00221 * 00222 ANRM = CLANGE( 'M', M, N, A, LDA, RWORK ) 00223 IASCL = 0 00224 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00225 * 00226 * Scale matrix norm up to SMLNUM 00227 * 00228 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00229 IASCL = 1 00230 ELSE IF( ANRM.GT.BIGNUM ) THEN 00231 * 00232 * Scale matrix norm down to BIGNUM 00233 * 00234 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00235 IASCL = 2 00236 ELSE IF( ANRM.EQ.ZERO ) THEN 00237 * 00238 * Matrix all zero. Return zero solution. 00239 * 00240 CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00241 RANK = 0 00242 GO TO 70 00243 END IF 00244 * 00245 BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK ) 00246 IBSCL = 0 00247 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00248 * 00249 * Scale matrix norm up to SMLNUM 00250 * 00251 CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00252 IBSCL = 1 00253 ELSE IF( BNRM.GT.BIGNUM ) THEN 00254 * 00255 * Scale matrix norm down to BIGNUM 00256 * 00257 CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00258 IBSCL = 2 00259 END IF 00260 * 00261 * Compute QR factorization with column pivoting of A: 00262 * A * P = Q * R 00263 * 00264 CALL CGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), 00265 $ LWORK-MN, RWORK, INFO ) 00266 WSIZE = MN + REAL( WORK( MN+1 ) ) 00267 * 00268 * complex workspace: MN+NB*(N+1). real workspace 2*N. 00269 * Details of Householder rotations stored in WORK(1:MN). 00270 * 00271 * Determine RANK using incremental condition estimation 00272 * 00273 WORK( ISMIN ) = CONE 00274 WORK( ISMAX ) = CONE 00275 SMAX = ABS( A( 1, 1 ) ) 00276 SMIN = SMAX 00277 IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN 00278 RANK = 0 00279 CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00280 GO TO 70 00281 ELSE 00282 RANK = 1 00283 END IF 00284 * 00285 10 CONTINUE 00286 IF( RANK.LT.MN ) THEN 00287 I = RANK + 1 00288 CALL CLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), 00289 $ A( I, I ), SMINPR, S1, C1 ) 00290 CALL CLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), 00291 $ A( I, I ), SMAXPR, S2, C2 ) 00292 * 00293 IF( SMAXPR*RCOND.LE.SMINPR ) THEN 00294 DO 20 I = 1, RANK 00295 WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) 00296 WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 00297 20 CONTINUE 00298 WORK( ISMIN+RANK ) = C1 00299 WORK( ISMAX+RANK ) = C2 00300 SMIN = SMINPR 00301 SMAX = SMAXPR 00302 RANK = RANK + 1 00303 GO TO 10 00304 END IF 00305 END IF 00306 * 00307 * complex workspace: 3*MN. 00308 * 00309 * Logically partition R = [ R11 R12 ] 00310 * [ 0 R22 ] 00311 * where R11 = R(1:RANK,1:RANK) 00312 * 00313 * [R11,R12] = [ T11, 0 ] * Y 00314 * 00315 IF( RANK.LT.N ) 00316 $ CALL CTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ), 00317 $ LWORK-2*MN, INFO ) 00318 * 00319 * complex workspace: 2*MN. 00320 * Details of Householder rotations stored in WORK(MN+1:2*MN) 00321 * 00322 * B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS) 00323 * 00324 CALL CUNMQR( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA, 00325 $ WORK( 1 ), B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO ) 00326 WSIZE = MAX( WSIZE, 2*MN+REAL( WORK( 2*MN+1 ) ) ) 00327 * 00328 * complex workspace: 2*MN+NB*NRHS. 00329 * 00330 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) 00331 * 00332 CALL CTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, 00333 $ NRHS, CONE, A, LDA, B, LDB ) 00334 * 00335 DO 40 J = 1, NRHS 00336 DO 30 I = RANK + 1, N 00337 B( I, J ) = CZERO 00338 30 CONTINUE 00339 40 CONTINUE 00340 * 00341 * B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS) 00342 * 00343 IF( RANK.LT.N ) THEN 00344 CALL CUNMRZ( 'Left', 'Conjugate transpose', N, NRHS, RANK, 00345 $ N-RANK, A, LDA, WORK( MN+1 ), B, LDB, 00346 $ WORK( 2*MN+1 ), LWORK-2*MN, INFO ) 00347 END IF 00348 * 00349 * complex workspace: 2*MN+NRHS. 00350 * 00351 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS) 00352 * 00353 DO 60 J = 1, NRHS 00354 DO 50 I = 1, N 00355 WORK( JPVT( I ) ) = B( I, J ) 00356 50 CONTINUE 00357 CALL CCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 ) 00358 60 CONTINUE 00359 * 00360 * complex workspace: N. 00361 * 00362 * Undo scaling 00363 * 00364 IF( IASCL.EQ.1 ) THEN 00365 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00366 CALL CLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, 00367 $ INFO ) 00368 ELSE IF( IASCL.EQ.2 ) THEN 00369 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00370 CALL CLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, 00371 $ INFO ) 00372 END IF 00373 IF( IBSCL.EQ.1 ) THEN 00374 CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00375 ELSE IF( IBSCL.EQ.2 ) THEN 00376 CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00377 END IF 00378 * 00379 70 CONTINUE 00380 WORK( 1 ) = CMPLX( LWKOPT ) 00381 * 00382 RETURN 00383 * 00384 * End of CGELSY 00385 * 00386 END
1.7.3 
