001:       SUBROUTINE CORE_DTSTRF( M, N, IB, NB, U, LDU, A, LDA,
002:      $                       L, LDL, IPIV, WORK, LDWORK, INFO )
003: 
004: *********************************************************************
005: *     PLASMA core_blas routine (version 2.1.0)                      *
006: *     Author: Hatem Ltaief                                          *
007: *     Release Date: November, 15th 2009                             *
008: *     PLASMA is a software package provided by Univ. of Tennessee,  *
009: *     Univ. of California Berkeley and Univ. of Colorado Denver.    *
010: *********************************************************************
011: *
012: *     .. Scalar Arguments ..
013:       INTEGER            M, N, IB, NB, LDU, LDA, LDL, LDWORK, INFO
014: *     ..
015: *     .. Array Arguments ..
016:       DOUBLE PRECISION            U( LDU, * ), A( LDA, * )
017:       DOUBLE PRECISION            L( LDL, * ), WORK( LDWORK, * )
018:       INTEGER            IPIV( * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  CORE_DTSTRF computes an LU factorization of a real matrix formed
025: *  by an upper triangular NB-by-N tile U on top of a M-by-N tile A
026: *  using partial pivoting with row interchanges.
027: *
028: *  This is the right-looking Level 2.5 BLAS version of the algorithm.
029: *
030: *  Arguments
031: *  =========
032: *
033: *  M       (input) INTEGER
034: *          The number of rows of the tile A.  M >= 0.
035: *
036: *  N       (input) INTEGER
037: *          The number of columns of the tile A.  N >= 0.
038: *
039: *  IB      (input) INTEGER
040: *          The inner-blocking size.  IB >= 0.
041: *
042: *  U       (input/output) DOUBLE PRECISION array, dimension (LDU,N)
043: *          On entry, the NB-by-N upper triangular tile.
044: *          On exit, the new factor U from the factorization
045: *
046: *  LDU     (input) INTEGER
047: *          The leading dimension of the array U.  LDU >= max(1,NB).
048: *
049: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
050: *          On entry, the M-by-N tile to be factored.
051: *          On exit, the factor L from the factorization
052: *
053: *  LDA     (input) INTEGER
054: *          The leading dimension of the array A.  LDA >= max(1,M).
055: *
056: *  L       (input/output) DOUBLE PRECISION array, dimension (LDL,N)
057: *          On entry, the NB-by-N lower triangular tile.
058: *          On exit, the interchanged rows formthe tile A in case of pivoting.
059: *
060: *  LDL     (input) INTEGER
061: *          The leading dimension of the array L.  LDL >= max(1,NB).
062: *
063: *  IPIV    (output) INTEGER array, dimension (min(M,N))
064: *          The pivot indices; for 1 <= i <= min(M,N), row i of the
065: *          tile U was interchanged with row IPIV(i) of the tile A.
066: *
067: *  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,IB)
068: *
069: *  LDWORK  (input) INTEGER
070: *          The dimension of the array WORK.
071: *
072: *  INFO    (output) INTEGER
073: *          = 0: successful exit
074: *          < 0: if INFO = -k, the k-th argument had an illegal value
075: *          > 0: if INFO = k, U(k,k) is exactly zero. The factorization
076: *               has been completed, but the factor U is exactly
077: *               singular, and division by zero will occur if it is used
078: *               to solve a system of equations.
079: *
080: *  =====================================================================
081: *     ..
082: *     Internal variables ..
083:       INTEGER            I, II, J, SB, IM, IP, IDAMAX, IINFO
084:       EXTERNAL           IDAMAX
085: *     ..
086: *     .. Parameters ..
087:       DOUBLE PRECISION            MDONE, DZERO
088:       PARAMETER          ( MDONE = -1.0D+0 )
089:       PARAMETER          ( DZERO = 0.0D+0 )
090: *     ..
091: *     Test the input parameters.
092:       INFO = 0
093:       IF( M.LT.0 ) THEN
094:          INFO = -1
095:       ELSE IF( N.LT.0 ) THEN
096:          INFO = -2
097:       ELSE IF( IB.LT.0 ) THEN
098:          INFO = -3
099:       ELSE IF( LDU.LT.MAX( 1, M ) ) THEN
100:          INFO = -6
101:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
102:          INFO = -8
103:       ELSE IF( LDL.LT.MAX( 1, IB ) ) THEN
104:          INFO = -10
105:       END IF
106:       IF( INFO.NE.0 ) THEN
107:          CALL XERBLA( 'CORE_DTSTRF', -INFO )
108:          RETURN
109:       END IF
110: *
111: *     Quick return if possible.
112: *
113:       IF( M.EQ.0 .OR. N.EQ.0 .OR. IB.EQ.0 )
114:      $   RETURN
115: *
116:       IP = 1
117: *
118:       DO 20 II = 1, N, IB
119:          SB = MIN( N-II+1, IB )
120:          DO 10 I = 1, SB
121: *
122:             IM = IDAMAX( M, A( 1, II+I-1 ), 1 )
123:             IPIV( IP ) = II+I-1
124: *
125:             IF( ABS( A( IM, II+I-1 ) ) 
126:      $          .GT. ABS( U( II+I-1, II+I-1 ) ) ) THEN
127: *
128: *              Swap behind.
129: *
130:                CALL DSWAP( I-1, L( I, II ), LDL,
131:      $                    WORK( IM, 1 ), LDWORK )
132: *
133: *              Swap ahead.
134: *
135:                CALL DSWAP( SB-I+1, U( II+I-1, II+I-1 ), LDU,
136:      $                    A( IM, II+I-1 ), LDA )
137: *
138: *              Set IPIV.
139: *
140:                IPIV( IP ) = NB + IM
141: *
142:                DO 50 J = 1, I-1
143:                   A( IM, II+J-1 ) = DZERO
144:  50            CONTINUE
145: *
146:             END IF
147:             IF( INFO.EQ.0 .AND. ABS(A(IM,II+I-1)).EQ.DZERO 
148:      $          .AND. ABS(U(II+I-1, II+I-1)).EQ.DZERO ) THEN
149:                INFO = II+I-1
150:             END IF
151: *
152:             CALL DSCAL( M, 1/U( II+I-1, II+I-1 ), A( 1, II+I-1 ), 1 )
153:             CALL DCOPY( M, A( 1, II+I-1 ), 1, WORK( 1, I ), 1 )
154:             CALL DGER( M, SB-I, MDONE, A( 1, II+I-1 ), 1,
155:      $                U( II+I-1, II+I ), LDU,
156:      $                A( 1, II+I ), LDA )
157: *
158:             IP = IP+1
159:  10      CONTINUE
160: *
161: *        Apply the subpanel to the rest of the panel.
162: *
163:          IF( II+I-1.LE. N ) THEN
164:             DO 80 J = II, II+SB-1
165:                IF( IPIV( J ).LE.NB ) THEN
166:                   IPIV( J ) = IPIV( J )-II+1
167:                ENDIF
168:  80         CONTINUE
169: *
170:             CALL CORE_DSSSSM( NB, M, N-( II+SB-1 ), SB, SB,
171:      $                       U( II, II+SB ), LDU,
172:      $                       A( 1, II+SB ), LDA,
173:      $                       L( 1, II ), LDL,
174:      $                       WORK, LDWORK, IPIV( II ),
175:      $                       IINFO )
176: *
177:             DO 70 J = II, II+SB-1
178:                IF( IPIV( J ).LE.NB ) THEN
179:                   IPIV( J ) = IPIV( J )+II-1
180:                ENDIF
181:  70         CONTINUE
182:          END IF
183:  20   CONTINUE
184: *
185:       RETURN
186: *
187: *     End of CORE_DTSTRF.
188: *
189:       END
190: