|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 00002 $ WORK, RWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK 00011 DOUBLE PRECISION RCOND 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER JPVT( * ) 00015 DOUBLE PRECISION RWORK( * ) 00016 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * This routine is deprecated and has been replaced by routine ZGELSY. 00023 * 00024 * ZGELSX computes the minimum-norm solution to a complex linear least 00025 * squares problem: 00026 * minimize || A * X - B || 00027 * using a complete orthogonal factorization of A. A is an M-by-N 00028 * matrix which may be rank-deficient. 00029 * 00030 * Several right hand side vectors b and solution vectors x can be 00031 * handled in a single call; they are stored as the columns of the 00032 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution 00033 * matrix X. 00034 * 00035 * The routine first computes a QR factorization with column pivoting: 00036 * A * P = Q * [ R11 R12 ] 00037 * [ 0 R22 ] 00038 * with R11 defined as the largest leading submatrix whose estimated 00039 * condition number is less than 1/RCOND. The order of R11, RANK, 00040 * is the effective rank of A. 00041 * 00042 * Then, R22 is considered to be negligible, and R12 is annihilated 00043 * by unitary transformations from the right, arriving at the 00044 * complete orthogonal factorization: 00045 * A * P = Q * [ T11 0 ] * Z 00046 * [ 0 0 ] 00047 * The minimum-norm solution is then 00048 * X = P * Z' [ inv(T11)*Q1'*B ] 00049 * [ 0 ] 00050 * where Q1 consists of the first RANK columns of Q. 00051 * 00052 * Arguments 00053 * ========= 00054 * 00055 * M (input) INTEGER 00056 * The number of rows of the matrix A. M >= 0. 00057 * 00058 * N (input) INTEGER 00059 * The number of columns of the matrix A. N >= 0. 00060 * 00061 * NRHS (input) INTEGER 00062 * The number of right hand sides, i.e., the number of 00063 * columns of matrices B and X. NRHS >= 0. 00064 * 00065 * A (input/output) COMPLEX*16 array, dimension (LDA,N) 00066 * On entry, the M-by-N matrix A. 00067 * On exit, A has been overwritten by details of its 00068 * complete orthogonal factorization. 00069 * 00070 * LDA (input) INTEGER 00071 * The leading dimension of the array A. LDA >= max(1,M). 00072 * 00073 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) 00074 * On entry, the M-by-NRHS right hand side matrix B. 00075 * On exit, the N-by-NRHS solution matrix X. 00076 * If m >= n and RANK = n, the residual sum-of-squares for 00077 * the solution in the i-th column is given by the sum of 00078 * squares of elements N+1:M in that column. 00079 * 00080 * LDB (input) INTEGER 00081 * The leading dimension of the array B. LDB >= max(1,M,N). 00082 * 00083 * JPVT (input/output) INTEGER array, dimension (N) 00084 * On entry, if JPVT(i) .ne. 0, the i-th column of A is an 00085 * initial column, otherwise it is a free column. Before 00086 * the QR factorization of A, all initial columns are 00087 * permuted to the leading positions; only the remaining 00088 * free columns are moved as a result of column pivoting 00089 * during the factorization. 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) DOUBLE PRECISION 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) COMPLEX*16 array, dimension 00105 * (min(M,N) + max( N, 2*min(M,N)+NRHS )), 00106 * 00107 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) 00108 * 00109 * INFO (output) INTEGER 00110 * = 0: successful exit 00111 * < 0: if INFO = -i, the i-th argument had an illegal value 00112 * 00113 * ===================================================================== 00114 * 00115 * .. Parameters .. 00116 INTEGER IMAX, IMIN 00117 PARAMETER ( IMAX = 1, IMIN = 2 ) 00118 DOUBLE PRECISION ZERO, ONE, DONE, NTDONE 00119 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, DONE = ZERO, 00120 $ NTDONE = ONE ) 00121 COMPLEX*16 CZERO, CONE 00122 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00123 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00124 * .. 00125 * .. Local Scalars .. 00126 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN 00127 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR, 00128 $ SMLNUM 00129 COMPLEX*16 C1, C2, S1, S2, T1, T2 00130 * .. 00131 * .. External Subroutines .. 00132 EXTERNAL XERBLA, ZGEQPF, ZLAIC1, ZLASCL, ZLASET, ZLATZM, 00133 $ ZTRSM, ZTZRQF, ZUNM2R 00134 * .. 00135 * .. External Functions .. 00136 DOUBLE PRECISION DLAMCH, ZLANGE 00137 EXTERNAL DLAMCH, ZLANGE 00138 * .. 00139 * .. Intrinsic Functions .. 00140 INTRINSIC ABS, DCONJG, MAX, MIN 00141 * .. 00142 * .. Executable Statements .. 00143 * 00144 MN = MIN( M, N ) 00145 ISMIN = MN + 1 00146 ISMAX = 2*MN + 1 00147 * 00148 * Test the input arguments. 00149 * 00150 INFO = 0 00151 IF( M.LT.0 ) THEN 00152 INFO = -1 00153 ELSE IF( N.LT.0 ) THEN 00154 INFO = -2 00155 ELSE IF( NRHS.LT.0 ) THEN 00156 INFO = -3 00157 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00158 INFO = -5 00159 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 00160 INFO = -7 00161 END IF 00162 * 00163 IF( INFO.NE.0 ) THEN 00164 CALL XERBLA( 'ZGELSX', -INFO ) 00165 RETURN 00166 END IF 00167 * 00168 * Quick return if possible 00169 * 00170 IF( MIN( M, N, NRHS ).EQ.0 ) THEN 00171 RANK = 0 00172 RETURN 00173 END IF 00174 * 00175 * Get machine parameters 00176 * 00177 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 00178 BIGNUM = ONE / SMLNUM 00179 CALL DLABAD( SMLNUM, BIGNUM ) 00180 * 00181 * Scale A, B if max elements outside range [SMLNUM,BIGNUM] 00182 * 00183 ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK ) 00184 IASCL = 0 00185 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00186 * 00187 * Scale matrix norm up to SMLNUM 00188 * 00189 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00190 IASCL = 1 00191 ELSE IF( ANRM.GT.BIGNUM ) THEN 00192 * 00193 * Scale matrix norm down to BIGNUM 00194 * 00195 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00196 IASCL = 2 00197 ELSE IF( ANRM.EQ.ZERO ) THEN 00198 * 00199 * Matrix all zero. Return zero solution. 00200 * 00201 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00202 RANK = 0 00203 GO TO 100 00204 END IF 00205 * 00206 BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK ) 00207 IBSCL = 0 00208 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00209 * 00210 * Scale matrix norm up to SMLNUM 00211 * 00212 CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00213 IBSCL = 1 00214 ELSE IF( BNRM.GT.BIGNUM ) THEN 00215 * 00216 * Scale matrix norm down to BIGNUM 00217 * 00218 CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00219 IBSCL = 2 00220 END IF 00221 * 00222 * Compute QR factorization with column pivoting of A: 00223 * A * P = Q * R 00224 * 00225 CALL ZGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), RWORK, 00226 $ INFO ) 00227 * 00228 * complex workspace MN+N. Real workspace 2*N. Details of Householder 00229 * rotations stored in WORK(1:MN). 00230 * 00231 * Determine RANK using incremental condition estimation 00232 * 00233 WORK( ISMIN ) = CONE 00234 WORK( ISMAX ) = CONE 00235 SMAX = ABS( A( 1, 1 ) ) 00236 SMIN = SMAX 00237 IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN 00238 RANK = 0 00239 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00240 GO TO 100 00241 ELSE 00242 RANK = 1 00243 END IF 00244 * 00245 10 CONTINUE 00246 IF( RANK.LT.MN ) THEN 00247 I = RANK + 1 00248 CALL ZLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), 00249 $ A( I, I ), SMINPR, S1, C1 ) 00250 CALL ZLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), 00251 $ A( I, I ), SMAXPR, S2, C2 ) 00252 * 00253 IF( SMAXPR*RCOND.LE.SMINPR ) THEN 00254 DO 20 I = 1, RANK 00255 WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) 00256 WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 00257 20 CONTINUE 00258 WORK( ISMIN+RANK ) = C1 00259 WORK( ISMAX+RANK ) = C2 00260 SMIN = SMINPR 00261 SMAX = SMAXPR 00262 RANK = RANK + 1 00263 GO TO 10 00264 END IF 00265 END IF 00266 * 00267 * Logically partition R = [ R11 R12 ] 00268 * [ 0 R22 ] 00269 * where R11 = R(1:RANK,1:RANK) 00270 * 00271 * [R11,R12] = [ T11, 0 ] * Y 00272 * 00273 IF( RANK.LT.N ) 00274 $ CALL ZTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO ) 00275 * 00276 * Details of Householder rotations stored in WORK(MN+1:2*MN) 00277 * 00278 * B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) 00279 * 00280 CALL ZUNM2R( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA, 00281 $ WORK( 1 ), B, LDB, WORK( 2*MN+1 ), INFO ) 00282 * 00283 * workspace NRHS 00284 * 00285 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) 00286 * 00287 CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, 00288 $ NRHS, CONE, A, LDA, B, LDB ) 00289 * 00290 DO 40 I = RANK + 1, N 00291 DO 30 J = 1, NRHS 00292 B( I, J ) = CZERO 00293 30 CONTINUE 00294 40 CONTINUE 00295 * 00296 * B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) 00297 * 00298 IF( RANK.LT.N ) THEN 00299 DO 50 I = 1, RANK 00300 CALL ZLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA, 00301 $ DCONJG( WORK( MN+I ) ), B( I, 1 ), 00302 $ B( RANK+1, 1 ), LDB, WORK( 2*MN+1 ) ) 00303 50 CONTINUE 00304 END IF 00305 * 00306 * workspace NRHS 00307 * 00308 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS) 00309 * 00310 DO 90 J = 1, NRHS 00311 DO 60 I = 1, N 00312 WORK( 2*MN+I ) = NTDONE 00313 60 CONTINUE 00314 DO 80 I = 1, N 00315 IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN 00316 IF( JPVT( I ).NE.I ) THEN 00317 K = I 00318 T1 = B( K, J ) 00319 T2 = B( JPVT( K ), J ) 00320 70 CONTINUE 00321 B( JPVT( K ), J ) = T1 00322 WORK( 2*MN+K ) = DONE 00323 T1 = T2 00324 K = JPVT( K ) 00325 T2 = B( JPVT( K ), J ) 00326 IF( JPVT( K ).NE.I ) 00327 $ GO TO 70 00328 B( I, J ) = T1 00329 WORK( 2*MN+K ) = DONE 00330 END IF 00331 END IF 00332 80 CONTINUE 00333 90 CONTINUE 00334 * 00335 * Undo scaling 00336 * 00337 IF( IASCL.EQ.1 ) THEN 00338 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00339 CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, 00340 $ INFO ) 00341 ELSE IF( IASCL.EQ.2 ) THEN 00342 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00343 CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, 00344 $ INFO ) 00345 END IF 00346 IF( IBSCL.EQ.1 ) THEN 00347 CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00348 ELSE IF( IBSCL.EQ.2 ) THEN 00349 CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00350 END IF 00351 * 00352 100 CONTINUE 00353 * 00354 RETURN 00355 * 00356 * End of ZGELSX 00357 * 00358 END
1.7.3 
