SRC/slatdf.f

   1 *> \brief \b SLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLATDF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slatdf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slatdf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slatdf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
  22 *                          JPIV )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            IJOB, LDZ, N
  26 *       REAL               RDSCAL, RDSUM
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IPIV( * ), JPIV( * )
  30 *       REAL               RHS( * ), Z( LDZ, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SLATDF uses the LU factorization of the n-by-n matrix Z computed by
  40 *> SGETC2 and computes a contribution to the reciprocal Dif-estimate
  41 *> by solving Z * x = b for x, and choosing the r.h.s. b such that
  42 *> the norm of x is as large as possible. On entry RHS = b holds the
  43 *> contribution from earlier solved sub-systems, and on return RHS = x.
  44 *>
  45 *> The factorization of Z returned by SGETC2 has the form Z = P*L*U*Q,
  46 *> where P and Q are permutation matrices. L is lower triangular with
  47 *> unit diagonal elements and U is upper triangular.
  48 *> \endverbatim
  49 *
  50 *  Arguments:
  51 *  ==========
  52 *
  53 *> \param[in] IJOB
  54 *> \verbatim
  55 *>          IJOB is INTEGER
  56 *>          IJOB = 2: First compute an approximative null-vector e
  57 *>              of Z using SGECON, e is normalized and solve for
  58 *>              Zx = +-e - f with the sign giving the greater value
  59 *>              of 2-norm(x). About 5 times as expensive as Default.
  60 *>          IJOB .ne. 2: Local look ahead strategy where all entries of
  61 *>              the r.h.s. b is chosen as either +1 or -1 (Default).
  62 *> \endverbatim
  63 *>
  64 *> \param[in] N
  65 *> \verbatim
  66 *>          N is INTEGER
  67 *>          The number of columns of the matrix Z.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] Z
  71 *> \verbatim
  72 *>          Z is REAL array, dimension (LDZ, N)
  73 *>          On entry, the LU part of the factorization of the n-by-n
  74 *>          matrix Z computed by SGETC2:  Z = P * L * U * Q
  75 *> \endverbatim
  76 *>
  77 *> \param[in] LDZ
  78 *> \verbatim
  79 *>          LDZ is INTEGER
  80 *>          The leading dimension of the array Z.  LDA >= max(1, N).
  81 *> \endverbatim
  82 *>
  83 *> \param[in,out] RHS
  84 *> \verbatim
  85 *>          RHS is REAL array, dimension N.
  86 *>          On entry, RHS contains contributions from other subsystems.
  87 *>          On exit, RHS contains the solution of the subsystem with
  88 *>          entries acoording to the value of IJOB (see above).
  89 *> \endverbatim
  90 *>
  91 *> \param[in,out] RDSUM
  92 *> \verbatim
  93 *>          RDSUM is REAL
  94 *>          On entry, the sum of squares of computed contributions to
  95 *>          the Dif-estimate under computation by STGSYL, where the
  96 *>          scaling factor RDSCAL (see below) has been factored out.
  97 *>          On exit, the corresponding sum of squares updated with the
  98 *>          contributions from the current sub-system.
  99 *>          If TRANS = 'T' RDSUM is not touched.
 100 *>          NOTE: RDSUM only makes sense when STGSY2 is called by STGSYL.
 101 *> \endverbatim
 102 *>
 103 *> \param[in,out] RDSCAL
 104 *> \verbatim
 105 *>          RDSCAL is REAL
 106 *>          On entry, scaling factor used to prevent overflow in RDSUM.
 107 *>          On exit, RDSCAL is updated w.r.t. the current contributions
 108 *>          in RDSUM.
 109 *>          If TRANS = 'T', RDSCAL is not touched.
 110 *>          NOTE: RDSCAL only makes sense when STGSY2 is called by
 111 *>                STGSYL.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] IPIV
 115 *> \verbatim
 116 *>          IPIV is INTEGER array, dimension (N).
 117 *>          The pivot indices; for 1 <= i <= N, row i of the
 118 *>          matrix has been interchanged with row IPIV(i).
 119 *> \endverbatim
 120 *>
 121 *> \param[in] JPIV
 122 *> \verbatim
 123 *>          JPIV is INTEGER array, dimension (N).
 124 *>          The pivot indices; for 1 <= j <= N, column j of the
 125 *>          matrix has been interchanged with column JPIV(j).
 126 *> \endverbatim
 127 *
 128 *  Authors:
 129 *  ========
 130 *
 131 *> \author Univ. of Tennessee
 132 *> \author Univ. of California Berkeley
 133 *> \author Univ. of Colorado Denver
 134 *> \author NAG Ltd.
 135 *
 136 *> \date June 2016
 137 *
 138 *> \ingroup realOTHERauxiliary
 139 *
 140 *> \par Further Details:
 141 *  =====================
 142 *>
 143 *>  This routine is a further developed implementation of algorithm
 144 *>  BSOLVE in [1] using complete pivoting in the LU factorization.
 145 *
 146 *> \par Contributors:
 147 *  ==================
 148 *>
 149 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 150 *>     Umea University, S-901 87 Umea, Sweden.
 151 *
 152 *> \par References:
 153 *  ================
 154 *>
 155 *> \verbatim
 156 *>
 157 *>
 158 *>  [1] Bo Kagstrom and Lars Westin,
 159 *>      Generalized Schur Methods with Condition Estimators for
 160 *>      Solving the Generalized Sylvester Equation, IEEE Transactions
 161 *>      on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.
 162 *>
 163 *>  [2] Peter Poromaa,
 164 *>      On Efficient and Robust Estimators for the Separation
 165 *>      between two Regular Matrix Pairs with Applications in
 166 *>      Condition Estimation. Report IMINF-95.05, Departement of
 167 *>      Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.
 168 *> \endverbatim
 169 *>
 170 *  =====================================================================
 171       SUBROUTINE SLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
 172      $                   JPIV )
 173 *
 174 *  -- LAPACK auxiliary routine (version 3.6.1) --
 175 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 176 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 177 *     June 2016
 178 *
 179 *     .. Scalar Arguments ..
 180       INTEGER            IJOB, LDZ, N
 181       REAL               RDSCAL, RDSUM
 182 *     ..
 183 *     .. Array Arguments ..
 184       INTEGER            IPIV( * ), JPIV( * )
 185       REAL               RHS( * ), Z( LDZ, * )
 186 *     ..
 187 *
 188 *  =====================================================================
 189 *
 190 *     .. Parameters ..
 191       INTEGER            MAXDIM
 192       PARAMETER          ( MAXDIM = 8 )
 193       REAL               ZERO, ONE
 194       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 195 *     ..
 196 *     .. Local Scalars ..
 197       INTEGER            I, INFO, J, K
 198       REAL               BM, BP, PMONE, SMINU, SPLUS, TEMP
 199 *     ..
 200 *     .. Local Arrays ..
 201       INTEGER            IWORK( MAXDIM )
 202       REAL               WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM )
 203 *     ..
 204 *     .. External Subroutines ..
 205       EXTERNAL           SAXPY, SCOPY, SGECON, SGESC2, SLASSQ, SLASWP,
 206      $                   SSCAL
 207 *     ..
 208 *     .. External Functions ..
 209       REAL               SASUM, SDOT
 210       EXTERNAL           SASUM, SDOT
 211 *     ..
 212 *     .. Intrinsic Functions ..
 213       INTRINSIC          ABS, SQRT
 214 *     ..
 215 *     .. Executable Statements ..
 216 *
 217       IF( IJOB.NE.2 ) THEN
 218 *
 219 *        Apply permutations IPIV to RHS
 220 *
 221          CALL SLASWP( 1, RHS, LDZ, 1, N-1, IPIV, 1 )
 222 *
 223 *        Solve for L-part choosing RHS either to +1 or -1.
 224 *
 225          PMONE = -ONE
 226 *
 227          DO 10 J = 1, N - 1
 228             BP = RHS( J ) + ONE
 229             BM = RHS( J ) - ONE
 230             SPLUS = ONE
 231 *
 232 *           Look-ahead for L-part RHS(1:N-1) = + or -1, SPLUS and
 233 *           SMIN computed more efficiently than in BSOLVE [1].
 234 *
 235             SPLUS = SPLUS + SDOT( N-J, Z( J+1, J ), 1, Z( J+1, J ), 1 )
 236             SMINU = SDOT( N-J, Z( J+1, J ), 1, RHS( J+1 ), 1 )
 237             SPLUS = SPLUS*RHS( J )
 238             IF( SPLUS.GT.SMINU ) THEN
 239                RHS( J ) = BP
 240             ELSE IF( SMINU.GT.SPLUS ) THEN
 241                RHS( J ) = BM
 242             ELSE
 243 *
 244 *              In this case the updating sums are equal and we can
 245 *              choose RHS(J) +1 or -1. The first time this happens
 246 *              we choose -1, thereafter +1. This is a simple way to
 247 *              get good estimates of matrices like Byers well-known
 248 *              example (see [1]). (Not done in BSOLVE.)
 249 *
 250                RHS( J ) = RHS( J ) + PMONE
 251                PMONE = ONE
 252             END IF
 253 *
 254 *           Compute the remaining r.h.s.
 255 *
 256             TEMP = -RHS( J )
 257             CALL SAXPY( N-J, TEMP, Z( J+1, J ), 1, RHS( J+1 ), 1 )
 258 *
 259    10    CONTINUE
 260 *
 261 *        Solve for U-part, look-ahead for RHS(N) = +-1. This is not done
 262 *        in BSOLVE and will hopefully give us a better estimate because
 263 *        any ill-conditioning of the original matrix is transfered to U
 264 *        and not to L. U(N, N) is an approximation to sigma_min(LU).
 265 *
 266          CALL SCOPY( N-1, RHS, 1, XP, 1 )
 267          XP( N ) = RHS( N ) + ONE
 268          RHS( N ) = RHS( N ) - ONE
 269          SPLUS = ZERO
 270          SMINU = ZERO
 271          DO 30 I = N, 1, -1
 272             TEMP = ONE / Z( I, I )
 273             XP( I ) = XP( I )*TEMP
 274             RHS( I ) = RHS( I )*TEMP
 275             DO 20 K = I + 1, N
 276                XP( I ) = XP( I ) - XP( K )*( Z( I, K )*TEMP )
 277                RHS( I ) = RHS( I ) - RHS( K )*( Z( I, K )*TEMP )
 278    20       CONTINUE
 279             SPLUS = SPLUS + ABS( XP( I ) )
 280             SMINU = SMINU + ABS( RHS( I ) )
 281    30    CONTINUE
 282          IF( SPLUS.GT.SMINU )
 283      $      CALL SCOPY( N, XP, 1, RHS, 1 )
 284 *
 285 *        Apply the permutations JPIV to the computed solution (RHS)
 286 *
 287          CALL SLASWP( 1, RHS, LDZ, 1, N-1, JPIV, -1 )
 288 *
 289 *        Compute the sum of squares
 290 *
 291          CALL SLASSQ( N, RHS, 1, RDSCAL, RDSUM )
 292 *
 293       ELSE
 294 *
 295 *        IJOB = 2, Compute approximate nullvector XM of Z
 296 *
 297          CALL SGECON( 'I', N, Z, LDZ, ONE, TEMP, WORK, IWORK, INFO )
 298          CALL SCOPY( N, WORK( N+1 ), 1, XM, 1 )
 299 *
 300 *        Compute RHS
 301 *
 302          CALL SLASWP( 1, XM, LDZ, 1, N-1, IPIV, -1 )
 303          TEMP = ONE / SQRT( SDOT( N, XM, 1, XM, 1 ) )
 304          CALL SSCAL( N, TEMP, XM, 1 )
 305          CALL SCOPY( N, XM, 1, XP, 1 )
 306          CALL SAXPY( N, ONE, RHS, 1, XP, 1 )
 307          CALL SAXPY( N, -ONE, XM, 1, RHS, 1 )
 308          CALL SGESC2( N, Z, LDZ, RHS, IPIV, JPIV, TEMP )
 309          CALL SGESC2( N, Z, LDZ, XP, IPIV, JPIV, TEMP )
 310          IF( SASUM( N, XP, 1 ).GT.SASUM( N, RHS, 1 ) )
 311      $      CALL SCOPY( N, XP, 1, RHS, 1 )
 312 *
 313 *        Compute the sum of squares
 314 *
 315          CALL SLASSQ( N, RHS, 1, RDSCAL, RDSUM )
 316 *
 317       END IF
 318 *
 319       RETURN
 320 *
 321 *     End of SLATDF
 322 *
 323       END