Return to dgeqpf.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

    1:       SUBROUTINE DGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO )
    2: *
    3: *  -- LAPACK deprecated computational routine (version 3.2.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     June 2010
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            INFO, LDA, M, N
   10: *     ..
   11: *     .. Array Arguments ..
   12:       INTEGER            JPVT( * )
   13:       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  This routine is deprecated and has been replaced by routine DGEQP3.
   20: *
   21: *  DGEQPF computes a QR factorization with column pivoting of a
   22: *  real M-by-N matrix A: A*P = Q*R.
   23: *
   24: *  Arguments
   25: *  =========
   26: *
   27: *  M       (input) INTEGER
   28: *          The number of rows of the matrix A. M >= 0.
   29: *
   30: *  N       (input) INTEGER
   31: *          The number of columns of the matrix A. N >= 0
   32: *
   33: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
   34: *          On entry, the M-by-N matrix A.
   35: *          On exit, the upper triangle of the array contains the
   36: *          min(M,N)-by-N upper triangular matrix R; the elements
   37: *          below the diagonal, together with the array TAU,
   38: *          represent the orthogonal matrix Q as a product of
   39: *          min(m,n) elementary reflectors.
   40: *
   41: *  LDA     (input) INTEGER
   42: *          The leading dimension of the array A. LDA >= max(1,M).
   43: *
   44: *  JPVT    (input/output) INTEGER array, dimension (N)
   45: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
   46: *          to the front of A*P (a leading column); if JPVT(i) = 0,
   47: *          the i-th column of A is a free column.
   48: *          On exit, if JPVT(i) = k, then the i-th column of A*P
   49: *          was the k-th column of A.
   50: *
   51: *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
   52: *          The scalar factors of the elementary reflectors.
   53: *
   54: *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
   55: *
   56: *  INFO    (output) INTEGER
   57: *          = 0:  successful exit
   58: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   59: *
   60: *  Further Details
   61: *  ===============
   62: *
   63: *  The matrix Q is represented as a product of elementary reflectors
   64: *
   65: *     Q = H(1) H(2) . . . H(n)
   66: *
   67: *  Each H(i) has the form
   68: *
   69: *     H = I - tau * v * v'
   70: *
   71: *  where tau is a real scalar, and v is a real vector with
   72: *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
   73: *
   74: *  The matrix P is represented in jpvt as follows: If
   75: *     jpvt(j) = i
   76: *  then the jth column of P is the ith canonical unit vector.
   77: *
   78: *  Partial column norm updating strategy modified by
   79: *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
   80: *    University of Zagreb, Croatia.
   81: *     June 2010
   82: *  For more details see LAPACK Working Note 176.
   83: *
   84: *  =====================================================================
   85: *
   86: *     .. Parameters ..
   87:       DOUBLE PRECISION   ZERO, ONE
   88:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
   89: *     ..
   90: *     .. Local Scalars ..
   91:       INTEGER            I, ITEMP, J, MA, MN, PVT
   92:       DOUBLE PRECISION   AII, TEMP, TEMP2, TOL3Z
   93: *     ..
   94: *     .. External Subroutines ..
   95:       EXTERNAL           DGEQR2, DLARF, DLARFG, DORM2R, DSWAP, XERBLA
   96: *     ..
   97: *     .. Intrinsic Functions ..
   98:       INTRINSIC          ABS, MAX, MIN, SQRT
   99: *     ..
  100: *     .. External Functions ..
  101:       INTEGER            IDAMAX
  102:       DOUBLE PRECISION   DLAMCH, DNRM2
  103:       EXTERNAL           IDAMAX, DLAMCH, DNRM2
  104: *     ..
  105: *     .. Executable Statements ..
  106: *
  107: *     Test the input arguments
  108: *
  109:       INFO = 0
  110:       IF( M.LT.0 ) THEN
  111:          INFO = -1
  112:       ELSE IF( N.LT.0 ) THEN
  113:          INFO = -2
  114:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  115:          INFO = -4
  116:       END IF
  117:       IF( INFO.NE.0 ) THEN
  118:          CALL XERBLA( 'DGEQPF', -INFO )
  119:          RETURN
  120:       END IF
  121: *
  122:       MN = MIN( M, N )
  123:       TOL3Z = SQRT(DLAMCH('Epsilon'))
  124: *
  125: *     Move initial columns up front
  126: *
  127:       ITEMP = 1
  128:       DO 10 I = 1, N
  129:          IF( JPVT( I ).NE.0 ) THEN
  130:             IF( I.NE.ITEMP ) THEN
  131:                CALL DSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
  132:                JPVT( I ) = JPVT( ITEMP )
  133:                JPVT( ITEMP ) = I
  134:             ELSE
  135:                JPVT( I ) = I
  136:             END IF
  137:             ITEMP = ITEMP + 1
  138:          ELSE
  139:             JPVT( I ) = I
  140:          END IF
  141:    10 CONTINUE
  142:       ITEMP = ITEMP - 1
  143: *
  144: *     Compute the QR factorization and update remaining columns
  145: *
  146:       IF( ITEMP.GT.0 ) THEN
  147:          MA = MIN( ITEMP, M )
  148:          CALL DGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
  149:          IF( MA.LT.N ) THEN
  150:             CALL DORM2R( 'Left', 'Transpose', M, N-MA, MA, A, LDA, TAU,
  151:      $                   A( 1, MA+1 ), LDA, WORK, INFO )
  152:          END IF
  153:       END IF
  154: *
  155:       IF( ITEMP.LT.MN ) THEN
  156: *
  157: *        Initialize partial column norms. The first n elements of
  158: *        work store the exact column norms.
  159: *
  160:          DO 20 I = ITEMP + 1, N
  161:             WORK( I ) = DNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
  162:             WORK( N+I ) = WORK( I )
  163:    20    CONTINUE
  164: *
  165: *        Compute factorization
  166: *
  167:          DO 40 I = ITEMP + 1, MN
  168: *
  169: *           Determine ith pivot column and swap if necessary
  170: *
  171:             PVT = ( I-1 ) + IDAMAX( N-I+1, WORK( I ), 1 )
  172: *
  173:             IF( PVT.NE.I ) THEN
  174:                CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
  175:                ITEMP = JPVT( PVT )
  176:                JPVT( PVT ) = JPVT( I )
  177:                JPVT( I ) = ITEMP
  178:                WORK( PVT ) = WORK( I )
  179:                WORK( N+PVT ) = WORK( N+I )
  180:             END IF
  181: *
  182: *           Generate elementary reflector H(i)
  183: *
  184:             IF( I.LT.M ) THEN
  185:                CALL DLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) )
  186:             ELSE
  187:                CALL DLARFG( 1, A( M, M ), A( M, M ), 1, TAU( M ) )
  188:             END IF
  189: *
  190:             IF( I.LT.N ) THEN
  191: *
  192: *              Apply H(i) to A(i:m,i+1:n) from the left
  193: *
  194:                AII = A( I, I )
  195:                A( I, I ) = ONE
  196:                CALL DLARF( 'LEFT', M-I+1, N-I, A( I, I ), 1, TAU( I ),
  197:      $                     A( I, I+1 ), LDA, WORK( 2*N+1 ) )
  198:                A( I, I ) = AII
  199:             END IF
  200: *
  201: *           Update partial column norms
  202: *
  203:             DO 30 J = I + 1, N
  204:                IF( WORK( J ).NE.ZERO ) THEN
  205: *
  206: *                 NOTE: The following 4 lines follow from the analysis in
  207: *                 Lapack Working Note 176.
  208: *                 
  209:                   TEMP = ABS( A( I, J ) ) / WORK( J )
  210:                   TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
  211:                   TEMP2 = TEMP*( WORK( J ) / WORK( N+J ) )**2
  212:                   IF( TEMP2 .LE. TOL3Z ) THEN 
  213:                      IF( M-I.GT.0 ) THEN
  214:                         WORK( J ) = DNRM2( M-I, A( I+1, J ), 1 )
  215:                         WORK( N+J ) = WORK( J )
  216:                      ELSE
  217:                         WORK( J ) = ZERO
  218:                         WORK( N+J ) = ZERO
  219:                      END IF
  220:                   ELSE
  221:                      WORK( J ) = WORK( J )*SQRT( TEMP )
  222:                   END IF
  223:                END IF
  224:    30       CONTINUE
  225: *
  226:    40    CONTINUE
  227:       END IF
  228:       RETURN
  229: *
  230: *     End of DGEQPF
  231: *
  232:       END