LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgetri.f
Go to the documentation of this file.
00001 *> \brief \b SGETRI
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SGETRI + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgetri.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgetri.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgetri.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, LDA, LWORK, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       INTEGER            IPIV( * )
00028 *       REAL               A( LDA, * ), WORK( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> SGETRI computes the inverse of a matrix using the LU factorization
00038 *> computed by SGETRF.
00039 *>
00040 *> This method inverts U and then computes inv(A) by solving the system
00041 *> inv(A)*L = inv(U) for inv(A).
00042 *> \endverbatim
00043 *
00044 *  Arguments:
00045 *  ==========
00046 *
00047 *> \param[in] N
00048 *> \verbatim
00049 *>          N is INTEGER
00050 *>          The order of the matrix A.  N >= 0.
00051 *> \endverbatim
00052 *>
00053 *> \param[in,out] A
00054 *> \verbatim
00055 *>          A is REAL array, dimension (LDA,N)
00056 *>          On entry, the factors L and U from the factorization
00057 *>          A = P*L*U as computed by SGETRF.
00058 *>          On exit, if INFO = 0, the inverse of the original matrix A.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] LDA
00062 *> \verbatim
00063 *>          LDA is INTEGER
00064 *>          The leading dimension of the array A.  LDA >= max(1,N).
00065 *> \endverbatim
00066 *>
00067 *> \param[in] IPIV
00068 *> \verbatim
00069 *>          IPIV is INTEGER array, dimension (N)
00070 *>          The pivot indices from SGETRF; for 1<=i<=N, row i of the
00071 *>          matrix was interchanged with row IPIV(i).
00072 *> \endverbatim
00073 *>
00074 *> \param[out] WORK
00075 *> \verbatim
00076 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00077 *>          On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] LWORK
00081 *> \verbatim
00082 *>          LWORK is INTEGER
00083 *>          The dimension of the array WORK.  LWORK >= max(1,N).
00084 *>          For optimal performance LWORK >= N*NB, where NB is
00085 *>          the optimal blocksize returned by ILAENV.
00086 *>
00087 *>          If LWORK = -1, then a workspace query is assumed; the routine
00088 *>          only calculates the optimal size of the WORK array, returns
00089 *>          this value as the first entry of the WORK array, and no error
00090 *>          message related to LWORK is issued by XERBLA.
00091 *> \endverbatim
00092 *>
00093 *> \param[out] INFO
00094 *> \verbatim
00095 *>          INFO is INTEGER
00096 *>          = 0:  successful exit
00097 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00098 *>          > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is
00099 *>                singular and its inverse could not be computed.
00100 *> \endverbatim
00101 *
00102 *  Authors:
00103 *  ========
00104 *
00105 *> \author Univ. of Tennessee 
00106 *> \author Univ. of California Berkeley 
00107 *> \author Univ. of Colorado Denver 
00108 *> \author NAG Ltd. 
00109 *
00110 *> \date November 2011
00111 *
00112 *> \ingroup realGEcomputational
00113 *
00114 *  =====================================================================
00115       SUBROUTINE SGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
00116 *
00117 *  -- LAPACK computational routine (version 3.4.0) --
00118 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00119 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00120 *     November 2011
00121 *
00122 *     .. Scalar Arguments ..
00123       INTEGER            INFO, LDA, LWORK, N
00124 *     ..
00125 *     .. Array Arguments ..
00126       INTEGER            IPIV( * )
00127       REAL               A( LDA, * ), WORK( * )
00128 *     ..
00129 *
00130 *  =====================================================================
00131 *
00132 *     .. Parameters ..
00133       REAL               ZERO, ONE
00134       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00135 *     ..
00136 *     .. Local Scalars ..
00137       LOGICAL            LQUERY
00138       INTEGER            I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
00139      $                   NBMIN, NN
00140 *     ..
00141 *     .. External Functions ..
00142       INTEGER            ILAENV
00143       EXTERNAL           ILAENV
00144 *     ..
00145 *     .. External Subroutines ..
00146       EXTERNAL           SGEMM, SGEMV, SSWAP, STRSM, STRTRI, XERBLA
00147 *     ..
00148 *     .. Intrinsic Functions ..
00149       INTRINSIC          MAX, MIN
00150 *     ..
00151 *     .. Executable Statements ..
00152 *
00153 *     Test the input parameters.
00154 *
00155       INFO = 0
00156       NB = ILAENV( 1, 'SGETRI', ' ', N, -1, -1, -1 )
00157       LWKOPT = N*NB
00158       WORK( 1 ) = LWKOPT
00159       LQUERY = ( LWORK.EQ.-1 )
00160       IF( N.LT.0 ) THEN
00161          INFO = -1
00162       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00163          INFO = -3
00164       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
00165          INFO = -6
00166       END IF
00167       IF( INFO.NE.0 ) THEN
00168          CALL XERBLA( 'SGETRI', -INFO )
00169          RETURN
00170       ELSE IF( LQUERY ) THEN
00171          RETURN
00172       END IF
00173 *
00174 *     Quick return if possible
00175 *
00176       IF( N.EQ.0 )
00177      $   RETURN
00178 *
00179 *     Form inv(U).  If INFO > 0 from STRTRI, then U is singular,
00180 *     and the inverse is not computed.
00181 *
00182       CALL STRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
00183       IF( INFO.GT.0 )
00184      $   RETURN
00185 *
00186       NBMIN = 2
00187       LDWORK = N
00188       IF( NB.GT.1 .AND. NB.LT.N ) THEN
00189          IWS = MAX( LDWORK*NB, 1 )
00190          IF( LWORK.LT.IWS ) THEN
00191             NB = LWORK / LDWORK
00192             NBMIN = MAX( 2, ILAENV( 2, 'SGETRI', ' ', N, -1, -1, -1 ) )
00193          END IF
00194       ELSE
00195          IWS = N
00196       END IF
00197 *
00198 *     Solve the equation inv(A)*L = inv(U) for inv(A).
00199 *
00200       IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
00201 *
00202 *        Use unblocked code.
00203 *
00204          DO 20 J = N, 1, -1
00205 *
00206 *           Copy current column of L to WORK and replace with zeros.
00207 *
00208             DO 10 I = J + 1, N
00209                WORK( I ) = A( I, J )
00210                A( I, J ) = ZERO
00211    10       CONTINUE
00212 *
00213 *           Compute current column of inv(A).
00214 *
00215             IF( J.LT.N )
00216      $         CALL SGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
00217      $                     LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
00218    20    CONTINUE
00219       ELSE
00220 *
00221 *        Use blocked code.
00222 *
00223          NN = ( ( N-1 ) / NB )*NB + 1
00224          DO 50 J = NN, 1, -NB
00225             JB = MIN( NB, N-J+1 )
00226 *
00227 *           Copy current block column of L to WORK and replace with
00228 *           zeros.
00229 *
00230             DO 40 JJ = J, J + JB - 1
00231                DO 30 I = JJ + 1, N
00232                   WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
00233                   A( I, JJ ) = ZERO
00234    30          CONTINUE
00235    40       CONTINUE
00236 *
00237 *           Compute current block column of inv(A).
00238 *
00239             IF( J+JB.LE.N )
00240      $         CALL SGEMM( 'No transpose', 'No transpose', N, JB,
00241      $                     N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
00242      $                     WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
00243             CALL STRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
00244      $                  ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
00245    50    CONTINUE
00246       END IF
00247 *
00248 *     Apply column interchanges.
00249 *
00250       DO 60 J = N - 1, 1, -1
00251          JP = IPIV( J )
00252          IF( JP.NE.J )
00253      $      CALL SSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
00254    60 CONTINUE
00255 *
00256       WORK( 1 ) = IWS
00257       RETURN
00258 *
00259 *     End of SGETRI
00260 *
00261       END
 All Files Functions