![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> ZGELSX solves overdetermined or underdetermined systems for GE matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZGELSX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelsx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelsx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelsx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 00022 * WORK, RWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, LDB, M, N, NRHS, RANK 00026 * DOUBLE PRECISION RCOND 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER JPVT( * ) 00030 * DOUBLE PRECISION RWORK( * ) 00031 * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> This routine is deprecated and has been replaced by routine ZGELSY. 00041 *> 00042 *> ZGELSX computes the minimum-norm solution to a complex linear least 00043 *> squares problem: 00044 *> minimize || A * X - B || 00045 *> using a complete orthogonal factorization of A. A is an M-by-N 00046 *> matrix which may be rank-deficient. 00047 *> 00048 *> Several right hand side vectors b and solution vectors x can be 00049 *> handled in a single call; they are stored as the columns of the 00050 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution 00051 *> matrix X. 00052 *> 00053 *> The routine first computes a QR factorization with column pivoting: 00054 *> A * P = Q * [ R11 R12 ] 00055 *> [ 0 R22 ] 00056 *> with R11 defined as the largest leading submatrix whose estimated 00057 *> condition number is less than 1/RCOND. The order of R11, RANK, 00058 *> is the effective rank of A. 00059 *> 00060 *> Then, R22 is considered to be negligible, and R12 is annihilated 00061 *> by unitary transformations from the right, arriving at the 00062 *> complete orthogonal factorization: 00063 *> A * P = Q * [ T11 0 ] * Z 00064 *> [ 0 0 ] 00065 *> The minimum-norm solution is then 00066 *> X = P * Z**H [ inv(T11)*Q1**H*B ] 00067 *> [ 0 ] 00068 *> where Q1 consists of the first RANK columns of Q. 00069 *> \endverbatim 00070 * 00071 * Arguments: 00072 * ========== 00073 * 00074 *> \param[in] M 00075 *> \verbatim 00076 *> M is INTEGER 00077 *> The number of rows of the matrix A. M >= 0. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] N 00081 *> \verbatim 00082 *> N is INTEGER 00083 *> The number of columns of the matrix A. N >= 0. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] NRHS 00087 *> \verbatim 00088 *> NRHS is INTEGER 00089 *> The number of right hand sides, i.e., the number of 00090 *> columns of matrices B and X. NRHS >= 0. 00091 *> \endverbatim 00092 *> 00093 *> \param[in,out] A 00094 *> \verbatim 00095 *> A is COMPLEX*16 array, dimension (LDA,N) 00096 *> On entry, the M-by-N matrix A. 00097 *> On exit, A has been overwritten by details of its 00098 *> complete orthogonal factorization. 00099 *> \endverbatim 00100 *> 00101 *> \param[in] LDA 00102 *> \verbatim 00103 *> LDA is INTEGER 00104 *> The leading dimension of the array A. LDA >= max(1,M). 00105 *> \endverbatim 00106 *> 00107 *> \param[in,out] B 00108 *> \verbatim 00109 *> B is COMPLEX*16 array, dimension (LDB,NRHS) 00110 *> On entry, the M-by-NRHS right hand side matrix B. 00111 *> On exit, the N-by-NRHS solution matrix X. 00112 *> If m >= n and RANK = n, the residual sum-of-squares for 00113 *> the solution in the i-th column is given by the sum of 00114 *> squares of elements N+1:M in that column. 00115 *> \endverbatim 00116 *> 00117 *> \param[in] LDB 00118 *> \verbatim 00119 *> LDB is INTEGER 00120 *> The leading dimension of the array B. LDB >= max(1,M,N). 00121 *> \endverbatim 00122 *> 00123 *> \param[in,out] JPVT 00124 *> \verbatim 00125 *> JPVT is INTEGER array, dimension (N) 00126 *> On entry, if JPVT(i) .ne. 0, the i-th column of A is an 00127 *> initial column, otherwise it is a free column. Before 00128 *> the QR factorization of A, all initial columns are 00129 *> permuted to the leading positions; only the remaining 00130 *> free columns are moved as a result of column pivoting 00131 *> during the factorization. 00132 *> On exit, if JPVT(i) = k, then the i-th column of A*P 00133 *> was the k-th column of A. 00134 *> \endverbatim 00135 *> 00136 *> \param[in] RCOND 00137 *> \verbatim 00138 *> RCOND is DOUBLE PRECISION 00139 *> RCOND is used to determine the effective rank of A, which 00140 *> is defined as the order of the largest leading triangular 00141 *> submatrix R11 in the QR factorization with pivoting of A, 00142 *> whose estimated condition number < 1/RCOND. 00143 *> \endverbatim 00144 *> 00145 *> \param[out] RANK 00146 *> \verbatim 00147 *> RANK is INTEGER 00148 *> The effective rank of A, i.e., the order of the submatrix 00149 *> R11. This is the same as the order of the submatrix T11 00150 *> in the complete orthogonal factorization of A. 00151 *> \endverbatim 00152 *> 00153 *> \param[out] WORK 00154 *> \verbatim 00155 *> WORK is COMPLEX*16 array, dimension 00156 *> (min(M,N) + max( N, 2*min(M,N)+NRHS )), 00157 *> \endverbatim 00158 *> 00159 *> \param[out] RWORK 00160 *> \verbatim 00161 *> RWORK is DOUBLE PRECISION array, dimension (2*N) 00162 *> \endverbatim 00163 *> 00164 *> \param[out] INFO 00165 *> \verbatim 00166 *> INFO is INTEGER 00167 *> = 0: successful exit 00168 *> < 0: if INFO = -i, the i-th argument had an illegal value 00169 *> \endverbatim 00170 * 00171 * Authors: 00172 * ======== 00173 * 00174 *> \author Univ. of Tennessee 00175 *> \author Univ. of California Berkeley 00176 *> \author Univ. of Colorado Denver 00177 *> \author NAG Ltd. 00178 * 00179 *> \date November 2011 00180 * 00181 *> \ingroup complex16GEsolve 00182 * 00183 * ===================================================================== 00184 SUBROUTINE ZGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 00185 $ WORK, RWORK, INFO ) 00186 * 00187 * -- LAPACK driver routine (version 3.4.0) -- 00188 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00189 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00190 * November 2011 00191 * 00192 * .. Scalar Arguments .. 00193 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK 00194 DOUBLE PRECISION RCOND 00195 * .. 00196 * .. Array Arguments .. 00197 INTEGER JPVT( * ) 00198 DOUBLE PRECISION RWORK( * ) 00199 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 00200 * .. 00201 * 00202 * ===================================================================== 00203 * 00204 * .. Parameters .. 00205 INTEGER IMAX, IMIN 00206 PARAMETER ( IMAX = 1, IMIN = 2 ) 00207 DOUBLE PRECISION ZERO, ONE, DONE, NTDONE 00208 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, DONE = ZERO, 00209 $ NTDONE = ONE ) 00210 COMPLEX*16 CZERO, CONE 00211 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00212 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00213 * .. 00214 * .. Local Scalars .. 00215 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN 00216 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR, 00217 $ SMLNUM 00218 COMPLEX*16 C1, C2, S1, S2, T1, T2 00219 * .. 00220 * .. External Subroutines .. 00221 EXTERNAL XERBLA, ZGEQPF, ZLAIC1, ZLASCL, ZLASET, ZLATZM, 00222 $ ZTRSM, ZTZRQF, ZUNM2R 00223 * .. 00224 * .. External Functions .. 00225 DOUBLE PRECISION DLAMCH, ZLANGE 00226 EXTERNAL DLAMCH, ZLANGE 00227 * .. 00228 * .. Intrinsic Functions .. 00229 INTRINSIC ABS, DCONJG, MAX, MIN 00230 * .. 00231 * .. Executable Statements .. 00232 * 00233 MN = MIN( M, N ) 00234 ISMIN = MN + 1 00235 ISMAX = 2*MN + 1 00236 * 00237 * Test the input arguments. 00238 * 00239 INFO = 0 00240 IF( M.LT.0 ) THEN 00241 INFO = -1 00242 ELSE IF( N.LT.0 ) THEN 00243 INFO = -2 00244 ELSE IF( NRHS.LT.0 ) THEN 00245 INFO = -3 00246 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00247 INFO = -5 00248 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 00249 INFO = -7 00250 END IF 00251 * 00252 IF( INFO.NE.0 ) THEN 00253 CALL XERBLA( 'ZGELSX', -INFO ) 00254 RETURN 00255 END IF 00256 * 00257 * Quick return if possible 00258 * 00259 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 00260 RANK = 0 00261 RETURN 00262 END IF 00263 * 00264 * Get machine parameters 00265 * 00266 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 00267 BIGNUM = ONE / SMLNUM 00268 CALL DLABAD( SMLNUM, BIGNUM ) 00269 * 00270 * Scale A, B if max elements outside range [SMLNUM,BIGNUM] 00271 * 00272 ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK ) 00273 IASCL = 0 00274 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00275 * 00276 * Scale matrix norm up to SMLNUM 00277 * 00278 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00279 IASCL = 1 00280 ELSE IF( ANRM.GT.BIGNUM ) THEN 00281 * 00282 * Scale matrix norm down to BIGNUM 00283 * 00284 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00285 IASCL = 2 00286 ELSE IF( ANRM.EQ.ZERO ) THEN 00287 * 00288 * Matrix all zero. Return zero solution. 00289 * 00290 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00291 RANK = 0 00292 GO TO 100 00293 END IF 00294 * 00295 BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK ) 00296 IBSCL = 0 00297 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00298 * 00299 * Scale matrix norm up to SMLNUM 00300 * 00301 CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00302 IBSCL = 1 00303 ELSE IF( BNRM.GT.BIGNUM ) THEN 00304 * 00305 * Scale matrix norm down to BIGNUM 00306 * 00307 CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00308 IBSCL = 2 00309 END IF 00310 * 00311 * Compute QR factorization with column pivoting of A: 00312 * A * P = Q * R 00313 * 00314 CALL ZGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), RWORK, 00315 $ INFO ) 00316 * 00317 * complex workspace MN+N. Real workspace 2*N. Details of Householder 00318 * rotations stored in WORK(1:MN). 00319 * 00320 * Determine RANK using incremental condition estimation 00321 * 00322 WORK( ISMIN ) = CONE 00323 WORK( ISMAX ) = CONE 00324 SMAX = ABS( A( 1, 1 ) ) 00325 SMIN = SMAX 00326 IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN 00327 RANK = 0 00328 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00329 GO TO 100 00330 ELSE 00331 RANK = 1 00332 END IF 00333 * 00334 10 CONTINUE 00335 IF( RANK.LT.MN ) THEN 00336 I = RANK + 1 00337 CALL ZLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), 00338 $ A( I, I ), SMINPR, S1, C1 ) 00339 CALL ZLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), 00340 $ A( I, I ), SMAXPR, S2, C2 ) 00341 * 00342 IF( SMAXPR*RCOND.LE.SMINPR ) THEN 00343 DO 20 I = 1, RANK 00344 WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) 00345 WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 00346 20 CONTINUE 00347 WORK( ISMIN+RANK ) = C1 00348 WORK( ISMAX+RANK ) = C2 00349 SMIN = SMINPR 00350 SMAX = SMAXPR 00351 RANK = RANK + 1 00352 GO TO 10 00353 END IF 00354 END IF 00355 * 00356 * Logically partition R = [ R11 R12 ] 00357 * [ 0 R22 ] 00358 * where R11 = R(1:RANK,1:RANK) 00359 * 00360 * [R11,R12] = [ T11, 0 ] * Y 00361 * 00362 IF( RANK.LT.N ) 00363 $ CALL ZTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO ) 00364 * 00365 * Details of Householder rotations stored in WORK(MN+1:2*MN) 00366 * 00367 * B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS) 00368 * 00369 CALL ZUNM2R( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA, 00370 $ WORK( 1 ), B, LDB, WORK( 2*MN+1 ), INFO ) 00371 * 00372 * workspace NRHS 00373 * 00374 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) 00375 * 00376 CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, 00377 $ NRHS, CONE, A, LDA, B, LDB ) 00378 * 00379 DO 40 I = RANK + 1, N 00380 DO 30 J = 1, NRHS 00381 B( I, J ) = CZERO 00382 30 CONTINUE 00383 40 CONTINUE 00384 * 00385 * B(1:N,1:NRHS) := Y**H * B(1:N,1:NRHS) 00386 * 00387 IF( RANK.LT.N ) THEN 00388 DO 50 I = 1, RANK 00389 CALL ZLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA, 00390 $ DCONJG( WORK( MN+I ) ), B( I, 1 ), 00391 $ B( RANK+1, 1 ), LDB, WORK( 2*MN+1 ) ) 00392 50 CONTINUE 00393 END IF 00394 * 00395 * workspace NRHS 00396 * 00397 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS) 00398 * 00399 DO 90 J = 1, NRHS 00400 DO 60 I = 1, N 00401 WORK( 2*MN+I ) = NTDONE 00402 60 CONTINUE 00403 DO 80 I = 1, N 00404 IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN 00405 IF( JPVT( I ).NE.I ) THEN 00406 K = I 00407 T1 = B( K, J ) 00408 T2 = B( JPVT( K ), J ) 00409 70 CONTINUE 00410 B( JPVT( K ), J ) = T1 00411 WORK( 2*MN+K ) = DONE 00412 T1 = T2 00413 K = JPVT( K ) 00414 T2 = B( JPVT( K ), J ) 00415 IF( JPVT( K ).NE.I ) 00416 $ GO TO 70 00417 B( I, J ) = T1 00418 WORK( 2*MN+K ) = DONE 00419 END IF 00420 END IF 00421 80 CONTINUE 00422 90 CONTINUE 00423 * 00424 * Undo scaling 00425 * 00426 IF( IASCL.EQ.1 ) THEN 00427 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00428 CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, 00429 $ INFO ) 00430 ELSE IF( IASCL.EQ.2 ) THEN 00431 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00432 CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, 00433 $ INFO ) 00434 END IF 00435 IF( IBSCL.EQ.1 ) THEN 00436 CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00437 ELSE IF( IBSCL.EQ.2 ) THEN 00438 CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00439 END IF 00440 * 00441 100 CONTINUE 00442 * 00443 RETURN 00444 * 00445 * End of ZGELSX 00446 * 00447 END