Annotation of rpl/lapack/lapack/zgelsx.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
! 2: $ WORK, RWORK, INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
! 11: DOUBLE PRECISION RCOND
! 12: * ..
! 13: * .. Array Arguments ..
! 14: INTEGER JPVT( * )
! 15: DOUBLE PRECISION RWORK( * )
! 16: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * This routine is deprecated and has been replaced by routine ZGELSY.
! 23: *
! 24: * ZGELSX computes the minimum-norm solution to a complex linear least
! 25: * squares problem:
! 26: * minimize || A * X - B ||
! 27: * using a complete orthogonal factorization of A. A is an M-by-N
! 28: * matrix which may be rank-deficient.
! 29: *
! 30: * Several right hand side vectors b and solution vectors x can be
! 31: * handled in a single call; they are stored as the columns of the
! 32: * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
! 33: * matrix X.
! 34: *
! 35: * The routine first computes a QR factorization with column pivoting:
! 36: * A * P = Q * [ R11 R12 ]
! 37: * [ 0 R22 ]
! 38: * with R11 defined as the largest leading submatrix whose estimated
! 39: * condition number is less than 1/RCOND. The order of R11, RANK,
! 40: * is the effective rank of A.
! 41: *
! 42: * Then, R22 is considered to be negligible, and R12 is annihilated
! 43: * by unitary transformations from the right, arriving at the
! 44: * complete orthogonal factorization:
! 45: * A * P = Q * [ T11 0 ] * Z
! 46: * [ 0 0 ]
! 47: * The minimum-norm solution is then
! 48: * X = P * Z' [ inv(T11)*Q1'*B ]
! 49: * [ 0 ]
! 50: * where Q1 consists of the first RANK columns of Q.
! 51: *
! 52: * Arguments
! 53: * =========
! 54: *
! 55: * M (input) INTEGER
! 56: * The number of rows of the matrix A. M >= 0.
! 57: *
! 58: * N (input) INTEGER
! 59: * The number of columns of the matrix A. N >= 0.
! 60: *
! 61: * NRHS (input) INTEGER
! 62: * The number of right hand sides, i.e., the number of
! 63: * columns of matrices B and X. NRHS >= 0.
! 64: *
! 65: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
! 66: * On entry, the M-by-N matrix A.
! 67: * On exit, A has been overwritten by details of its
! 68: * complete orthogonal factorization.
! 69: *
! 70: * LDA (input) INTEGER
! 71: * The leading dimension of the array A. LDA >= max(1,M).
! 72: *
! 73: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
! 74: * On entry, the M-by-NRHS right hand side matrix B.
! 75: * On exit, the N-by-NRHS solution matrix X.
! 76: * If m >= n and RANK = n, the residual sum-of-squares for
! 77: * the solution in the i-th column is given by the sum of
! 78: * squares of elements N+1:M in that column.
! 79: *
! 80: * LDB (input) INTEGER
! 81: * The leading dimension of the array B. LDB >= max(1,M,N).
! 82: *
! 83: * JPVT (input/output) INTEGER array, dimension (N)
! 84: * On entry, if JPVT(i) .ne. 0, the i-th column of A is an
! 85: * initial column, otherwise it is a free column. Before
! 86: * the QR factorization of A, all initial columns are
! 87: * permuted to the leading positions; only the remaining
! 88: * free columns are moved as a result of column pivoting
! 89: * during the factorization.
! 90: * On exit, if JPVT(i) = k, then the i-th column of A*P
! 91: * was the k-th column of A.
! 92: *
! 93: * RCOND (input) DOUBLE PRECISION
! 94: * RCOND is used to determine the effective rank of A, which
! 95: * is defined as the order of the largest leading triangular
! 96: * submatrix R11 in the QR factorization with pivoting of A,
! 97: * whose estimated condition number < 1/RCOND.
! 98: *
! 99: * RANK (output) INTEGER
! 100: * The effective rank of A, i.e., the order of the submatrix
! 101: * R11. This is the same as the order of the submatrix T11
! 102: * in the complete orthogonal factorization of A.
! 103: *
! 104: * WORK (workspace) COMPLEX*16 array, dimension
! 105: * (min(M,N) + max( N, 2*min(M,N)+NRHS )),
! 106: *
! 107: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
! 108: *
! 109: * INFO (output) INTEGER
! 110: * = 0: successful exit
! 111: * < 0: if INFO = -i, the i-th argument had an illegal value
! 112: *
! 113: * =====================================================================
! 114: *
! 115: * .. Parameters ..
! 116: INTEGER IMAX, IMIN
! 117: PARAMETER ( IMAX = 1, IMIN = 2 )
! 118: DOUBLE PRECISION ZERO, ONE, DONE, NTDONE
! 119: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, DONE = ZERO,
! 120: $ NTDONE = ONE )
! 121: COMPLEX*16 CZERO, CONE
! 122: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
! 123: $ CONE = ( 1.0D+0, 0.0D+0 ) )
! 124: * ..
! 125: * .. Local Scalars ..
! 126: INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
! 127: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
! 128: $ SMLNUM
! 129: COMPLEX*16 C1, C2, S1, S2, T1, T2
! 130: * ..
! 131: * .. External Subroutines ..
! 132: EXTERNAL XERBLA, ZGEQPF, ZLAIC1, ZLASCL, ZLASET, ZLATZM,
! 133: $ ZTRSM, ZTZRQF, ZUNM2R
! 134: * ..
! 135: * .. External Functions ..
! 136: DOUBLE PRECISION DLAMCH, ZLANGE
! 137: EXTERNAL DLAMCH, ZLANGE
! 138: * ..
! 139: * .. Intrinsic Functions ..
! 140: INTRINSIC ABS, DCONJG, MAX, MIN
! 141: * ..
! 142: * .. Executable Statements ..
! 143: *
! 144: MN = MIN( M, N )
! 145: ISMIN = MN + 1
! 146: ISMAX = 2*MN + 1
! 147: *
! 148: * Test the input arguments.
! 149: *
! 150: INFO = 0
! 151: IF( M.LT.0 ) THEN
! 152: INFO = -1
! 153: ELSE IF( N.LT.0 ) THEN
! 154: INFO = -2
! 155: ELSE IF( NRHS.LT.0 ) THEN
! 156: INFO = -3
! 157: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 158: INFO = -5
! 159: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
! 160: INFO = -7
! 161: END IF
! 162: *
! 163: IF( INFO.NE.0 ) THEN
! 164: CALL XERBLA( 'ZGELSX', -INFO )
! 165: RETURN
! 166: END IF
! 167: *
! 168: * Quick return if possible
! 169: *
! 170: IF( MIN( M, N, NRHS ).EQ.0 ) THEN
! 171: RANK = 0
! 172: RETURN
! 173: END IF
! 174: *
! 175: * Get machine parameters
! 176: *
! 177: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
! 178: BIGNUM = ONE / SMLNUM
! 179: CALL DLABAD( SMLNUM, BIGNUM )
! 180: *
! 181: * Scale A, B if max elements outside range [SMLNUM,BIGNUM]
! 182: *
! 183: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
! 184: IASCL = 0
! 185: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 186: *
! 187: * Scale matrix norm up to SMLNUM
! 188: *
! 189: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
! 190: IASCL = 1
! 191: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 192: *
! 193: * Scale matrix norm down to BIGNUM
! 194: *
! 195: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
! 196: IASCL = 2
! 197: ELSE IF( ANRM.EQ.ZERO ) THEN
! 198: *
! 199: * Matrix all zero. Return zero solution.
! 200: *
! 201: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
! 202: RANK = 0
! 203: GO TO 100
! 204: END IF
! 205: *
! 206: BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
! 207: IBSCL = 0
! 208: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 209: *
! 210: * Scale matrix norm up to SMLNUM
! 211: *
! 212: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
! 213: IBSCL = 1
! 214: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 215: *
! 216: * Scale matrix norm down to BIGNUM
! 217: *
! 218: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
! 219: IBSCL = 2
! 220: END IF
! 221: *
! 222: * Compute QR factorization with column pivoting of A:
! 223: * A * P = Q * R
! 224: *
! 225: CALL ZGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), RWORK,
! 226: $ INFO )
! 227: *
! 228: * complex workspace MN+N. Real workspace 2*N. Details of Householder
! 229: * rotations stored in WORK(1:MN).
! 230: *
! 231: * Determine RANK using incremental condition estimation
! 232: *
! 233: WORK( ISMIN ) = CONE
! 234: WORK( ISMAX ) = CONE
! 235: SMAX = ABS( A( 1, 1 ) )
! 236: SMIN = SMAX
! 237: IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
! 238: RANK = 0
! 239: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
! 240: GO TO 100
! 241: ELSE
! 242: RANK = 1
! 243: END IF
! 244: *
! 245: 10 CONTINUE
! 246: IF( RANK.LT.MN ) THEN
! 247: I = RANK + 1
! 248: CALL ZLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
! 249: $ A( I, I ), SMINPR, S1, C1 )
! 250: CALL ZLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
! 251: $ A( I, I ), SMAXPR, S2, C2 )
! 252: *
! 253: IF( SMAXPR*RCOND.LE.SMINPR ) THEN
! 254: DO 20 I = 1, RANK
! 255: WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
! 256: WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
! 257: 20 CONTINUE
! 258: WORK( ISMIN+RANK ) = C1
! 259: WORK( ISMAX+RANK ) = C2
! 260: SMIN = SMINPR
! 261: SMAX = SMAXPR
! 262: RANK = RANK + 1
! 263: GO TO 10
! 264: END IF
! 265: END IF
! 266: *
! 267: * Logically partition R = [ R11 R12 ]
! 268: * [ 0 R22 ]
! 269: * where R11 = R(1:RANK,1:RANK)
! 270: *
! 271: * [R11,R12] = [ T11, 0 ] * Y
! 272: *
! 273: IF( RANK.LT.N )
! 274: $ CALL ZTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO )
! 275: *
! 276: * Details of Householder rotations stored in WORK(MN+1:2*MN)
! 277: *
! 278: * B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
! 279: *
! 280: CALL ZUNM2R( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA,
! 281: $ WORK( 1 ), B, LDB, WORK( 2*MN+1 ), INFO )
! 282: *
! 283: * workspace NRHS
! 284: *
! 285: * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
! 286: *
! 287: CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
! 288: $ NRHS, CONE, A, LDA, B, LDB )
! 289: *
! 290: DO 40 I = RANK + 1, N
! 291: DO 30 J = 1, NRHS
! 292: B( I, J ) = CZERO
! 293: 30 CONTINUE
! 294: 40 CONTINUE
! 295: *
! 296: * B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
! 297: *
! 298: IF( RANK.LT.N ) THEN
! 299: DO 50 I = 1, RANK
! 300: CALL ZLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA,
! 301: $ DCONJG( WORK( MN+I ) ), B( I, 1 ),
! 302: $ B( RANK+1, 1 ), LDB, WORK( 2*MN+1 ) )
! 303: 50 CONTINUE
! 304: END IF
! 305: *
! 306: * workspace NRHS
! 307: *
! 308: * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
! 309: *
! 310: DO 90 J = 1, NRHS
! 311: DO 60 I = 1, N
! 312: WORK( 2*MN+I ) = NTDONE
! 313: 60 CONTINUE
! 314: DO 80 I = 1, N
! 315: IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN
! 316: IF( JPVT( I ).NE.I ) THEN
! 317: K = I
! 318: T1 = B( K, J )
! 319: T2 = B( JPVT( K ), J )
! 320: 70 CONTINUE
! 321: B( JPVT( K ), J ) = T1
! 322: WORK( 2*MN+K ) = DONE
! 323: T1 = T2
! 324: K = JPVT( K )
! 325: T2 = B( JPVT( K ), J )
! 326: IF( JPVT( K ).NE.I )
! 327: $ GO TO 70
! 328: B( I, J ) = T1
! 329: WORK( 2*MN+K ) = DONE
! 330: END IF
! 331: END IF
! 332: 80 CONTINUE
! 333: 90 CONTINUE
! 334: *
! 335: * Undo scaling
! 336: *
! 337: IF( IASCL.EQ.1 ) THEN
! 338: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
! 339: CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
! 340: $ INFO )
! 341: ELSE IF( IASCL.EQ.2 ) THEN
! 342: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
! 343: CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
! 344: $ INFO )
! 345: END IF
! 346: IF( IBSCL.EQ.1 ) THEN
! 347: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
! 348: ELSE IF( IBSCL.EQ.2 ) THEN
! 349: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
! 350: END IF
! 351: *
! 352: 100 CONTINUE
! 353: *
! 354: RETURN
! 355: *
! 356: * End of ZGELSX
! 357: *
! 358: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgelsx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack