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