SUBROUTINE SLUBL( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine -- * Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, * Courant Institute, NAG Ltd., and Rice University * August 16, 1989 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) REAL A( LDA, * ) * .. * * Purpose * ======= * * SLUBL computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the left-looking Level 3 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) REAL array, dimension (LDA,N) * On entry, the m x n matrix to be factored. * On exit, the factors L and U; the unit diagonal elements of L * are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices. Row i of the matrix was interchanged with * row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations or to compute the inverse * of A. * * ===================================================================== * * .. Parameters .. REAL ONE PARAMETER ( ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, J, JB, NB * .. * .. External Subroutines .. EXTERNAL ENVIR, SGEMM, SLASWP, SLUL, STRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'SLUBL', -INFO ) RETURN END IF * * Quick return if possible. * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * * Determine the block size for this environment. * CALL ENVIR( 'Block', NB ) IF( NB.EQ.1 ) THEN CALL SLUL( M, N, A, LDA, IPIV, INFO ) ELSE * DO 20 J = 1, MIN( M, N ), NB JB = MIN( MIN( M, N )-J+1, NB ) * * Apply previous interchanges to current block. * CALL SLASWP( JB, A( 1, J ), LDA, 1, J-1, IPIV, 1 ) * * Compute superdiagonal block of U. * CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', J-1, $ JB, ONE, A, LDA, A( 1, J ), LDA ) * * Update diagonal and subdiagonal blocks. * CALL SGEMM( 'No transpose', 'No transpose', M-J+1, JB, J-1, $ -ONE, A( J, 1 ), LDA, A( 1, J ), LDA, ONE, $ A( J, J ), LDA ) * * Factorize diagonal and subdiagonal blocks and test for exact * singularity. * CALL SLUL( M-J+1, JB, A( J, J ), LDA, IPIV( J ), INFO ) DO 10 I = J, J + JB - 1 IPIV( I ) = J - 1 + IPIV( I ) 10 CONTINUE IF( INFO.EQ.0 ) THEN * * Apply interchanges in current block to previous blocks. * CALL SLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) ELSE * * If INFO is not zero, a zero pivot was found in SLUL. * Correct the index returned from SLUL and go on. * INFO = INFO + J - 1 END IF 20 CONTINUE * IF( N.GT.M ) THEN * * Apply interchanges from columns 1 to m to columns m+1 to n. * CALL SLASWP( N-M, A( 1, M+1 ), LDA, 1, M, IPIV, 1 ) * * Complete columns m+1 to n of U. * CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', M, N-M, $ ONE, A, LDA, A( 1, M+1 ), LDA ) END IF END IF RETURN * * End of SLUBL * END SUBROUTINE SLUL( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine -- * Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, * Courant Institute, NAG Ltd., and Rice University * August 16, 1989 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) REAL A( LDA, * ) * .. * * Purpose * ======= * * SLUL computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the left-looking Level 2 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) REAL array, dimension (LDA,N) * On entry, the m x n matrix to be factored. * On exit, the factors L and U; the unit diagonal elements of L * are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices. Row i of the matrix was interchanged with * row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations or to compute the inverse * of A. * * ===================================================================== * * .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Local Scalars .. INTEGER I, IP, J, JP, LEN REAL T * .. * .. External Functions .. INTEGER ISAMAX EXTERNAL ISAMAX * .. * .. External Subroutines .. EXTERNAL SGEMV, SSCAL, SSWAP, STRSV, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'SLUL', -INFO ) RETURN END IF * * Quick return if possible. * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * DO 20 J = 1, N * * Apply previous interchanges to J-th column. * LEN = MIN( J-1, M ) DO 10 I = 1, LEN IP = IPIV( I ) IF( IP.NE.I ) THEN T = A( I, J ) A( I, J ) = A( IP, J ) A( IP, J ) = T END IF 10 CONTINUE * * Compute elements 1:LEN of J-th column. * CALL STRSV( 'Lower', 'No transpose', 'Unit', LEN, A, LDA, $ A( 1, J ), 1 ) * * Update elements J:M of J-th column. * IF( J.LE.M ) THEN CALL SGEMV( 'No transpose', M-J+1, J-1, -ONE, A( J, 1 ), $ LDA, A( 1, J ), 1, ONE, A( J, J ), 1 ) * * Find pivot and test for singularity. * JP = J - 1 + ISAMAX( M-J+1, A( J, J ), 1 ) IPIV( J ) = JP IF( A( JP, J ).NE.ZERO ) THEN * * Apply interchange to columns 1:J. * IF( JP.NE.J ) $ CALL SSWAP( J, A( J, 1 ), LDA, A( JP, 1 ), LDA ) * * Compute elements J+1:M of J-th column. * IF( J.LT.M ) $ CALL SSCAL( M-J, ONE/A( J, J ), A( J+1, J ), 1 ) ELSE * * If A( JP, J ) is zero, set INFO to indicate that a zero * pivot has been found. INFO returns the index of the * last zero pivot to the calling routine if there is more * than one. * INFO = J END IF END IF 20 CONTINUE RETURN * * End of SLUL * END SUBROUTINE SLUBC( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine -- * Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, * Courant Institute, NAG Ltd., and Rice University * August 16, 1989 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) REAL A( LDA, * ) * .. * * Purpose * ======= * * SLUBC computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the Level 3 BLAS version of the Crout algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) REAL array, dimension (LDA,N) * On entry, the m x n matrix to be factored. * On exit, the factors L and U; the unit diagonal elements of L * are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices. Row i of the matrix was interchanged with * row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations or to compute the inverse * of A. * * ===================================================================== * * .. Parameters .. REAL ONE PARAMETER ( ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, J, JB, NB * .. * .. External Subroutines .. EXTERNAL ENVIR, SGEMM, SLASWP, SLUC, STRSM * .. * .. Intrinsic Functions .. INTRINSIC MIN * .. * .. Executable Statements .. * INFO = 0 * * Determine the block size for this environment. * CALL ENVIR( 'Block', NB ) IF( NB.EQ.1 ) THEN CALL SLUC( M, N, A, LDA, IPIV, INFO ) ELSE * DO 20 J = 1, MIN( M, N ), NB JB = MIN( MIN( M, N )-J+1, NB ) * * Update diagonal and subdiagonal blocks. * CALL SGEMM( 'No transpose', 'No transpose', M-J+1, JB, J-1, $ -ONE, A( J, 1 ), LDA, A( 1, J ), LDA, ONE, $ A( J, J ), LDA ) * * Factorize diagonal and subdiagonal blocks and test for * singularity. * CALL SLUC( M-J+1, JB, A( J, J ), LDA, IPIV( J ), INFO ) * * Update pivot indices and apply the interchanges to the * columns on either side of the current block. * DO 10 I = J, MIN( M, J+JB-1 ) IPIV( I ) = J - 1 + IPIV( I ) 10 CONTINUE CALL SLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) CALL SLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, IPIV, $ 1 ) * IF( INFO.NE.0 ) $ GO TO 30 * IF( J+JB.LE.N ) THEN * * Compute block row of U. * CALL SGEMM( 'No transpose', 'No transpose', JB, N-J-JB+1, $ J-1, -ONE, A( J, 1 ), LDA, A( 1, J+JB ), LDA, $ ONE, A( J, J+JB ), LDA ) CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), $ LDA ) END IF 20 CONTINUE END IF RETURN * 30 CONTINUE INFO = INFO + J - 1 RETURN * * End of SLUBC * END SUBROUTINE SLUC( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine -- * Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, * Courant Institute, NAG Ltd., and Rice University * August 16, 1989 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) REAL A( LDA, * ) * .. * * Purpose * ======= * * SLUC computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the Level 2 BLAS version of the Crout algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) REAL array, dimension (LDA,N) * On entry, the m x n matrix to be factored. * On exit, the factors L and U; the unit diagonal elements of L * are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices. Row i of the matrix was interchanged with * row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations or to compute the inverse * of A. * * ===================================================================== * * .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Local Scalars .. INTEGER J, JP * .. * .. External Functions .. INTEGER ISAMAX EXTERNAL ISAMAX * .. * .. External Subroutines .. EXTERNAL SGEMV, SSCAL, SSWAP * .. * .. Intrinsic Functions .. INTRINSIC MIN * .. * .. Executable Statements .. * INFO = 0 * DO 10 J = 1, MIN( M, N ) * * Update diagonal and subdiagonal elements in column J. * CALL SGEMV( 'No transpose', M-J+1, J-1, -ONE, A( J, 1 ), LDA, $ A( 1, J ), 1, ONE, A( J, J ), 1 ) * * Find pivot and test for singularity. * JP = J - 1 + ISAMAX( M-J+1, A( J, J ), 1 ) IPIV( J ) = JP IF( A( JP, J ).NE.ZERO ) THEN * * Apply interchange to columns 1:J. * IF( JP.NE.J ) $ CALL SSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA ) * * Compute elements J+1:M of J-th column. * IF( J.LT.M ) $ CALL SSCAL( M-J, ONE/A( J, J ), A( J+1, J ), 1 ) ELSE * * If A( JP, J ) is zero, set INFO to indicate that a zero * pivot has been found. INFO returns the index of the * last zero pivot to the calling routine if there is more * than one. * INFO = J END IF * IF( J+1.LE.N ) THEN * * Compute block row of U. * CALL SGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ), LDA, $ A( J, 1 ), LDA, ONE, A( J, J+1 ), LDA ) END IF 10 CONTINUE RETURN * * End of SLUC * END SUBROUTINE SLUBR( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine -- * Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, * Courant Institute, NAG Ltd., and Rice University * August 16, 1989 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) REAL A( LDA, * ) * .. * * Purpose * ======= * * SLUBR computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 3 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) REAL array, dimension (LDA,N) * On entry, the m x n matrix to be factored. * On exit, the factors L and U; the unit diagonal elements of L * are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices. Row i of the matrix was interchanged with * row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations or to compute the inverse * of A. * * ===================================================================== * * .. Parameters .. REAL ONE PARAMETER ( ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, J, JB, NB * .. * .. External Subroutines .. EXTERNAL ENVIR, SGEMM, SLASWP, SLUR, STRSM * .. * .. Intrinsic Functions .. INTRINSIC MIN * .. * .. Executable Statements .. * INFO = 0 * * Determine the block size for this environment. * CALL ENVIR( 'Block', NB ) IF( NB.EQ.1 ) THEN CALL SLUR( M, N, A, LDA, IPIV, INFO ) ELSE * DO 20 J = 1, MIN( M, N ), NB JB = MIN( MIN( M, N )-J+1, NB ) * * Factorize diagonal and subdiagonal blocks and test for * singularity. * CALL SLUR( M-J+1, JB, A( J, J ), LDA, IPIV( J ), INFO ) DO 10 I = J, J + JB - 1 IPIV( I ) = J - 1 + IPIV( I ) 10 CONTINUE IF( INFO.NE.0 ) $ GO TO 30 * * Apply interchanges to previous blocks. * CALL SLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) * * Apply interchanges to following blocks. * CALL SLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, IPIV, $ 1 ) IF( J+JB.LE.N ) THEN * * Compute block row of U. * CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), $ LDA ) IF( J+JB.LE.M ) THEN * * Update trailing submatrix. * CALL SGEMM( 'No transpose', 'No transpose', M-J-JB+1, $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA, $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ), $ LDA ) END IF END IF 20 CONTINUE END IF RETURN * 30 CONTINUE INFO = INFO + J - 1 RETURN * * End of SLUBR * END SUBROUTINE SLUR( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine -- * Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, * Courant Institute, NAG Ltd., and Rice University * August 16, 1989 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) REAL A( LDA, * ) * .. * * Purpose * ======= * * SLUR computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 2 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) REAL array, dimension (LDA,N) * On entry, the m x n matrix to be factored. * On exit, the factors L and U; the unit diagonal elements of L * are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER array, dimension (min(M,N)) * The pivot indices. Row i of the matrix was interchanged with * row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations or to compute the inverse * of A. * * ===================================================================== * * .. Parameters .. REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Local Scalars .. INTEGER J, JP * .. * .. External Functions .. INTEGER ISAMAX EXTERNAL ISAMAX * .. * .. External Subroutines .. EXTERNAL SGER, SSCAL, SSWAP * .. * .. Intrinsic Functions .. INTRINSIC MIN * .. * .. Executable Statements .. * INFO = 0 * DO 10 J = 1, MIN( M, N ) * * Find pivot and test for singularity. * JP = J - 1 + ISAMAX( M-J+1, A( J, J ), 1 ) IPIV( J ) = JP IF( A( JP, J ).NE.ZERO ) THEN * * Interchange rows J and JP. * IF( JP.NE.J ) $ CALL SSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA ) * * Compute elements J+1:M of J-th column. * IF( J.LT.M ) $ CALL SSCAL( M-J, ONE/A( J, J ), A( J+1, J ), 1 ) ELSE * * If A( JP, J ) is zero, set INFO to indicate that a zero * pivot has been found. INFO returns the index of the last * zero pivot to the calling routine if there is more than one. * INFO = J END IF * IF( J+1.LE.N ) THEN * * Update trailing submatrix. * CALL SGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA, $ A( J+1, J+1 ), LDA ) END IF 10 CONTINUE RETURN * * End of SLUR * END SUBROUTINE ENVIR( WHAT, NB ) * * -- LAPACK auxiliary routine -- * Argonne National Lab, Courant Institute, and N.A.G. Ltd. * April 1, 1989 * * .. Scalar Arguments .. CHARACTER WHAT INTEGER NB * .. * * Purpose * ======= * * ENVIR returns the block size to be used by the block algorithms. * * Arguments * ========= * * WHAT - CHARACTER * WHAT is a character code for the value to be returned. * * NB - INTEGER * The parameter NB returns an integer value from common * storage. * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Scalars in Common .. INTEGER NBLOCK, NPROC, NSHIFT * .. * .. Common blocks .. COMMON / CENVIR / NBLOCK, NPROC, NSHIFT * .. * .. Executable Statements .. IF( LSAME( WHAT, 'B' ) ) THEN NB = MAX( 1, NBLOCK ) ELSE IF( LSAME( WHAT, 'P' ) ) THEN NB = NPROC ELSE IF( LSAME( WHAT, 'S' ) ) THEN NB = NSHIFT END IF RETURN * * End of ENVIR * END