LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlatdf.f
Go to the documentation of this file.
00001 *> \brief \b DLATDF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLATDF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlatdf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlatdf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlatdf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
00022 *                          JPIV )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            IJOB, LDZ, N
00026 *       DOUBLE PRECISION   RDSCAL, RDSUM
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IPIV( * ), JPIV( * )
00030 *       DOUBLE PRECISION   RHS( * ), Z( LDZ, * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DLATDF uses the LU factorization of the n-by-n matrix Z computed by
00040 *> DGETC2 and computes a contribution to the reciprocal Dif-estimate
00041 *> by solving Z * x = b for x, and choosing the r.h.s. b such that
00042 *> the norm of x is as large as possible. On entry RHS = b holds the
00043 *> contribution from earlier solved sub-systems, and on return RHS = x.
00044 *>
00045 *> The factorization of Z returned by DGETC2 has the form Z = P*L*U*Q,
00046 *> where P and Q are permutation matrices. L is lower triangular with
00047 *> unit diagonal elements and U is upper triangular.
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] IJOB
00054 *> \verbatim
00055 *>          IJOB is INTEGER
00056 *>          IJOB = 2: First compute an approximative null-vector e
00057 *>              of Z using DGECON, e is normalized and solve for
00058 *>              Zx = +-e - f with the sign giving the greater value
00059 *>              of 2-norm(x). About 5 times as expensive as Default.
00060 *>          IJOB .ne. 2: Local look ahead strategy where all entries of
00061 *>              the r.h.s. b is choosen as either +1 or -1 (Default).
00062 *> \endverbatim
00063 *>
00064 *> \param[in] N
00065 *> \verbatim
00066 *>          N is INTEGER
00067 *>          The number of columns of the matrix Z.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] Z
00071 *> \verbatim
00072 *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
00073 *>          On entry, the LU part of the factorization of the n-by-n
00074 *>          matrix Z computed by DGETC2:  Z = P * L * U * Q
00075 *> \endverbatim
00076 *>
00077 *> \param[in] LDZ
00078 *> \verbatim
00079 *>          LDZ is INTEGER
00080 *>          The leading dimension of the array Z.  LDA >= max(1, N).
00081 *> \endverbatim
00082 *>
00083 *> \param[in,out] RHS
00084 *> \verbatim
00085 *>          RHS is DOUBLE PRECISION array, dimension (N)
00086 *>          On entry, RHS contains contributions from other subsystems.
00087 *>          On exit, RHS contains the solution of the subsystem with
00088 *>          entries acoording to the value of IJOB (see above).
00089 *> \endverbatim
00090 *>
00091 *> \param[in,out] RDSUM
00092 *> \verbatim
00093 *>          RDSUM is DOUBLE PRECISION
00094 *>          On entry, the sum of squares of computed contributions to
00095 *>          the Dif-estimate under computation by DTGSYL, where the
00096 *>          scaling factor RDSCAL (see below) has been factored out.
00097 *>          On exit, the corresponding sum of squares updated with the
00098 *>          contributions from the current sub-system.
00099 *>          If TRANS = 'T' RDSUM is not touched.
00100 *>          NOTE: RDSUM only makes sense when DTGSY2 is called by STGSYL.
00101 *> \endverbatim
00102 *>
00103 *> \param[in,out] RDSCAL
00104 *> \verbatim
00105 *>          RDSCAL is DOUBLE PRECISION
00106 *>          On entry, scaling factor used to prevent overflow in RDSUM.
00107 *>          On exit, RDSCAL is updated w.r.t. the current contributions
00108 *>          in RDSUM.
00109 *>          If TRANS = 'T', RDSCAL is not touched.
00110 *>          NOTE: RDSCAL only makes sense when DTGSY2 is called by
00111 *>                DTGSYL.
00112 *> \endverbatim
00113 *>
00114 *> \param[in] IPIV
00115 *> \verbatim
00116 *>          IPIV is INTEGER array, dimension (N).
00117 *>          The pivot indices; for 1 <= i <= N, row i of the
00118 *>          matrix has been interchanged with row IPIV(i).
00119 *> \endverbatim
00120 *>
00121 *> \param[in] JPIV
00122 *> \verbatim
00123 *>          JPIV is INTEGER array, dimension (N).
00124 *>          The pivot indices; for 1 <= j <= N, column j of the
00125 *>          matrix has been interchanged with column JPIV(j).
00126 *> \endverbatim
00127 *
00128 *  Authors:
00129 *  ========
00130 *
00131 *> \author Univ. of Tennessee 
00132 *> \author Univ. of California Berkeley 
00133 *> \author Univ. of Colorado Denver 
00134 *> \author NAG Ltd. 
00135 *
00136 *> \date November 2011
00137 *
00138 *> \ingroup doubleOTHERauxiliary
00139 *
00140 *> \par Further Details:
00141 *  =====================
00142 *>
00143 *>  This routine is a further developed implementation of algorithm
00144 *>  BSOLVE in [1] using complete pivoting in the LU factorization.
00145 *
00146 *> \par Contributors:
00147 *  ==================
00148 *>
00149 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00150 *>     Umea University, S-901 87 Umea, Sweden.
00151 *
00152 *> \par References:
00153 *  ================
00154 *>
00155 *> \verbatim
00156 *>
00157 *>
00158 *>  [1] Bo Kagstrom and Lars Westin,
00159 *>      Generalized Schur Methods with Condition Estimators for
00160 *>      Solving the Generalized Sylvester Equation, IEEE Transactions
00161 *>      on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.
00162 *>
00163 *>  [2] Peter Poromaa,
00164 *>      On Efficient and Robust Estimators for the Separation
00165 *>      between two Regular Matrix Pairs with Applications in
00166 *>      Condition Estimation. Report IMINF-95.05, Departement of
00167 *>      Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.
00168 *> \endverbatim
00169 *>
00170 *  =====================================================================
00171       SUBROUTINE DLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
00172      $                   JPIV )
00173 *
00174 *  -- LAPACK auxiliary routine (version 3.4.0) --
00175 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00176 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00177 *     November 2011
00178 *
00179 *     .. Scalar Arguments ..
00180       INTEGER            IJOB, LDZ, N
00181       DOUBLE PRECISION   RDSCAL, RDSUM
00182 *     ..
00183 *     .. Array Arguments ..
00184       INTEGER            IPIV( * ), JPIV( * )
00185       DOUBLE PRECISION   RHS( * ), Z( LDZ, * )
00186 *     ..
00187 *
00188 *  =====================================================================
00189 *
00190 *     .. Parameters ..
00191       INTEGER            MAXDIM
00192       PARAMETER          ( MAXDIM = 8 )
00193       DOUBLE PRECISION   ZERO, ONE
00194       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00195 *     ..
00196 *     .. Local Scalars ..
00197       INTEGER            I, INFO, J, K
00198       DOUBLE PRECISION   BM, BP, PMONE, SMINU, SPLUS, TEMP
00199 *     ..
00200 *     .. Local Arrays ..
00201       INTEGER            IWORK( MAXDIM )
00202       DOUBLE PRECISION   WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM )
00203 *     ..
00204 *     .. External Subroutines ..
00205       EXTERNAL           DAXPY, DCOPY, DGECON, DGESC2, DLASSQ, DLASWP,
00206      $                   DSCAL
00207 *     ..
00208 *     .. External Functions ..
00209       DOUBLE PRECISION   DASUM, DDOT
00210       EXTERNAL           DASUM, DDOT
00211 *     ..
00212 *     .. Intrinsic Functions ..
00213       INTRINSIC          ABS, SQRT
00214 *     ..
00215 *     .. Executable Statements ..
00216 *
00217       IF( IJOB.NE.2 ) THEN
00218 *
00219 *        Apply permutations IPIV to RHS
00220 *
00221          CALL DLASWP( 1, RHS, LDZ, 1, N-1, IPIV, 1 )
00222 *
00223 *        Solve for L-part choosing RHS either to +1 or -1.
00224 *
00225          PMONE = -ONE
00226 *
00227          DO 10 J = 1, N - 1
00228             BP = RHS( J ) + ONE
00229             BM = RHS( J ) - ONE
00230             SPLUS = ONE
00231 *
00232 *           Look-ahead for L-part RHS(1:N-1) = + or -1, SPLUS and
00233 *           SMIN computed more efficiently than in BSOLVE [1].
00234 *
00235             SPLUS = SPLUS + DDOT( N-J, Z( J+1, J ), 1, Z( J+1, J ), 1 )
00236             SMINU = DDOT( N-J, Z( J+1, J ), 1, RHS( J+1 ), 1 )
00237             SPLUS = SPLUS*RHS( J )
00238             IF( SPLUS.GT.SMINU ) THEN
00239                RHS( J ) = BP
00240             ELSE IF( SMINU.GT.SPLUS ) THEN
00241                RHS( J ) = BM
00242             ELSE
00243 *
00244 *              In this case the updating sums are equal and we can
00245 *              choose RHS(J) +1 or -1. The first time this happens
00246 *              we choose -1, thereafter +1. This is a simple way to
00247 *              get good estimates of matrices like Byers well-known
00248 *              example (see [1]). (Not done in BSOLVE.)
00249 *
00250                RHS( J ) = RHS( J ) + PMONE
00251                PMONE = ONE
00252             END IF
00253 *
00254 *           Compute the remaining r.h.s.
00255 *
00256             TEMP = -RHS( J )
00257             CALL DAXPY( N-J, TEMP, Z( J+1, J ), 1, RHS( J+1 ), 1 )
00258 *
00259    10    CONTINUE
00260 *
00261 *        Solve for U-part, look-ahead for RHS(N) = +-1. This is not done
00262 *        in BSOLVE and will hopefully give us a better estimate because
00263 *        any ill-conditioning of the original matrix is transfered to U
00264 *        and not to L. U(N, N) is an approximation to sigma_min(LU).
00265 *
00266          CALL DCOPY( N-1, RHS, 1, XP, 1 )
00267          XP( N ) = RHS( N ) + ONE
00268          RHS( N ) = RHS( N ) - ONE
00269          SPLUS = ZERO
00270          SMINU = ZERO
00271          DO 30 I = N, 1, -1
00272             TEMP = ONE / Z( I, I )
00273             XP( I ) = XP( I )*TEMP
00274             RHS( I ) = RHS( I )*TEMP
00275             DO 20 K = I + 1, N
00276                XP( I ) = XP( I ) - XP( K )*( Z( I, K )*TEMP )
00277                RHS( I ) = RHS( I ) - RHS( K )*( Z( I, K )*TEMP )
00278    20       CONTINUE
00279             SPLUS = SPLUS + ABS( XP( I ) )
00280             SMINU = SMINU + ABS( RHS( I ) )
00281    30    CONTINUE
00282          IF( SPLUS.GT.SMINU )
00283      $      CALL DCOPY( N, XP, 1, RHS, 1 )
00284 *
00285 *        Apply the permutations JPIV to the computed solution (RHS)
00286 *
00287          CALL DLASWP( 1, RHS, LDZ, 1, N-1, JPIV, -1 )
00288 *
00289 *        Compute the sum of squares
00290 *
00291          CALL DLASSQ( N, RHS, 1, RDSCAL, RDSUM )
00292 *
00293       ELSE
00294 *
00295 *        IJOB = 2, Compute approximate nullvector XM of Z
00296 *
00297          CALL DGECON( 'I', N, Z, LDZ, ONE, TEMP, WORK, IWORK, INFO )
00298          CALL DCOPY( N, WORK( N+1 ), 1, XM, 1 )
00299 *
00300 *        Compute RHS
00301 *
00302          CALL DLASWP( 1, XM, LDZ, 1, N-1, IPIV, -1 )
00303          TEMP = ONE / SQRT( DDOT( N, XM, 1, XM, 1 ) )
00304          CALL DSCAL( N, TEMP, XM, 1 )
00305          CALL DCOPY( N, XM, 1, XP, 1 )
00306          CALL DAXPY( N, ONE, RHS, 1, XP, 1 )
00307          CALL DAXPY( N, -ONE, XM, 1, RHS, 1 )
00308          CALL DGESC2( N, Z, LDZ, RHS, IPIV, JPIV, TEMP )
00309          CALL DGESC2( N, Z, LDZ, XP, IPIV, JPIV, TEMP )
00310          IF( DASUM( N, XP, 1 ).GT.DASUM( N, RHS, 1 ) )
00311      $      CALL DCOPY( N, XP, 1, RHS, 1 )
00312 *
00313 *        Compute the sum of squares
00314 *
00315          CALL DLASSQ( N, RHS, 1, RDSCAL, RDSUM )
00316 *
00317       END IF
00318 *
00319       RETURN
00320 *
00321 *     End of DLATDF
00322 *
00323       END
 All Files Functions