00001 SUBROUTINE CLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, 00002 $ JPIV ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER IJOB, LDZ, N 00011 REAL RDSCAL, RDSUM 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IPIV( * ), JPIV( * ) 00015 COMPLEX RHS( * ), Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CLATDF computes the contribution to the reciprocal Dif-estimate 00022 * by solving for x in Z * x = b, where b is chosen such that the norm 00023 * of x is as large as possible. It is assumed that LU decomposition 00024 * of Z has been computed by CGETC2. On entry RHS = f holds the 00025 * contribution from earlier solved sub-systems, and on return RHS = x. 00026 * 00027 * The factorization of Z returned by CGETC2 has the form 00028 * Z = P * L * U * Q, where P and Q are permutation matrices. L is lower 00029 * triangular with unit diagonal elements and U is upper triangular. 00030 * 00031 * Arguments 00032 * ========= 00033 * 00034 * IJOB (input) INTEGER 00035 * IJOB = 2: First compute an approximative null-vector e 00036 * of Z using CGECON, e is normalized and solve for 00037 * Zx = +-e - f with the sign giving the greater value of 00038 * 2-norm(x). About 5 times as expensive as Default. 00039 * IJOB .ne. 2: Local look ahead strategy where 00040 * all entries of the r.h.s. b is choosen as either +1 or 00041 * -1. Default. 00042 * 00043 * N (input) INTEGER 00044 * The number of columns of the matrix Z. 00045 * 00046 * Z (input) REAL array, dimension (LDZ, N) 00047 * On entry, the LU part of the factorization of the n-by-n 00048 * matrix Z computed by CGETC2: Z = P * L * U * Q 00049 * 00050 * LDZ (input) INTEGER 00051 * The leading dimension of the array Z. LDA >= max(1, N). 00052 * 00053 * RHS (input/output) REAL array, dimension (N). 00054 * On entry, RHS contains contributions from other subsystems. 00055 * On exit, RHS contains the solution of the subsystem with 00056 * entries according to the value of IJOB (see above). 00057 * 00058 * RDSUM (input/output) REAL 00059 * On entry, the sum of squares of computed contributions to 00060 * the Dif-estimate under computation by CTGSYL, where the 00061 * scaling factor RDSCAL (see below) has been factored out. 00062 * On exit, the corresponding sum of squares updated with the 00063 * contributions from the current sub-system. 00064 * If TRANS = 'T' RDSUM is not touched. 00065 * NOTE: RDSUM only makes sense when CTGSY2 is called by CTGSYL. 00066 * 00067 * RDSCAL (input/output) REAL 00068 * On entry, scaling factor used to prevent overflow in RDSUM. 00069 * On exit, RDSCAL is updated w.r.t. the current contributions 00070 * in RDSUM. 00071 * If TRANS = 'T', RDSCAL is not touched. 00072 * NOTE: RDSCAL only makes sense when CTGSY2 is called by 00073 * CTGSYL. 00074 * 00075 * IPIV (input) INTEGER array, dimension (N). 00076 * The pivot indices; for 1 <= i <= N, row i of the 00077 * matrix has been interchanged with row IPIV(i). 00078 * 00079 * JPIV (input) INTEGER array, dimension (N). 00080 * The pivot indices; for 1 <= j <= N, column j of the 00081 * matrix has been interchanged with column JPIV(j). 00082 * 00083 * Further Details 00084 * =============== 00085 * 00086 * Based on contributions by 00087 * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00088 * Umea University, S-901 87 Umea, Sweden. 00089 * 00090 * This routine is a further developed implementation of algorithm 00091 * BSOLVE in [1] using complete pivoting in the LU factorization. 00092 * 00093 * [1] Bo Kagstrom and Lars Westin, 00094 * Generalized Schur Methods with Condition Estimators for 00095 * Solving the Generalized Sylvester Equation, IEEE Transactions 00096 * on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751. 00097 * 00098 * [2] Peter Poromaa, 00099 * On Efficient and Robust Estimators for the Separation 00100 * between two Regular Matrix Pairs with Applications in 00101 * Condition Estimation. Report UMINF-95.05, Department of 00102 * Computing Science, Umea University, S-901 87 Umea, Sweden, 00103 * 1995. 00104 * 00105 * ===================================================================== 00106 * 00107 * .. Parameters .. 00108 INTEGER MAXDIM 00109 PARAMETER ( MAXDIM = 2 ) 00110 REAL ZERO, ONE 00111 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00112 COMPLEX CONE 00113 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00114 * .. 00115 * .. Local Scalars .. 00116 INTEGER I, INFO, J, K 00117 REAL RTEMP, SCALE, SMINU, SPLUS 00118 COMPLEX BM, BP, PMONE, TEMP 00119 * .. 00120 * .. Local Arrays .. 00121 REAL RWORK( MAXDIM ) 00122 COMPLEX WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM ) 00123 * .. 00124 * .. External Subroutines .. 00125 EXTERNAL CAXPY, CCOPY, CGECON, CGESC2, CLASSQ, CLASWP, 00126 $ CSCAL 00127 * .. 00128 * .. External Functions .. 00129 REAL SCASUM 00130 COMPLEX CDOTC 00131 EXTERNAL SCASUM, CDOTC 00132 * .. 00133 * .. Intrinsic Functions .. 00134 INTRINSIC ABS, REAL, SQRT 00135 * .. 00136 * .. Executable Statements .. 00137 * 00138 IF( IJOB.NE.2 ) THEN 00139 * 00140 * Apply permutations IPIV to RHS 00141 * 00142 CALL CLASWP( 1, RHS, LDZ, 1, N-1, IPIV, 1 ) 00143 * 00144 * Solve for L-part choosing RHS either to +1 or -1. 00145 * 00146 PMONE = -CONE 00147 DO 10 J = 1, N - 1 00148 BP = RHS( J ) + CONE 00149 BM = RHS( J ) - CONE 00150 SPLUS = ONE 00151 * 00152 * Lockahead for L- part RHS(1:N-1) = +-1 00153 * SPLUS and SMIN computed more efficiently than in BSOLVE[1]. 00154 * 00155 SPLUS = SPLUS + REAL( CDOTC( N-J, Z( J+1, J ), 1, Z( J+1, 00156 $ J ), 1 ) ) 00157 SMINU = REAL( CDOTC( N-J, Z( J+1, J ), 1, RHS( J+1 ), 1 ) ) 00158 SPLUS = SPLUS*REAL( RHS( J ) ) 00159 IF( SPLUS.GT.SMINU ) THEN 00160 RHS( J ) = BP 00161 ELSE IF( SMINU.GT.SPLUS ) THEN 00162 RHS( J ) = BM 00163 ELSE 00164 * 00165 * In this case the updating sums are equal and we can 00166 * choose RHS(J) +1 or -1. The first time this happens we 00167 * choose -1, thereafter +1. This is a simple way to get 00168 * good estimates of matrices like Byers well-known example 00169 * (see [1]). (Not done in BSOLVE.) 00170 * 00171 RHS( J ) = RHS( J ) + PMONE 00172 PMONE = CONE 00173 END IF 00174 * 00175 * Compute the remaining r.h.s. 00176 * 00177 TEMP = -RHS( J ) 00178 CALL CAXPY( N-J, TEMP, Z( J+1, J ), 1, RHS( J+1 ), 1 ) 00179 10 CONTINUE 00180 * 00181 * Solve for U- part, lockahead for RHS(N) = +-1. This is not done 00182 * In BSOLVE and will hopefully give us a better estimate because 00183 * any ill-conditioning of the original matrix is transfered to U 00184 * and not to L. U(N, N) is an approximation to sigma_min(LU). 00185 * 00186 CALL CCOPY( N-1, RHS, 1, WORK, 1 ) 00187 WORK( N ) = RHS( N ) + CONE 00188 RHS( N ) = RHS( N ) - CONE 00189 SPLUS = ZERO 00190 SMINU = ZERO 00191 DO 30 I = N, 1, -1 00192 TEMP = CONE / Z( I, I ) 00193 WORK( I ) = WORK( I )*TEMP 00194 RHS( I ) = RHS( I )*TEMP 00195 DO 20 K = I + 1, N 00196 WORK( I ) = WORK( I ) - WORK( K )*( Z( I, K )*TEMP ) 00197 RHS( I ) = RHS( I ) - RHS( K )*( Z( I, K )*TEMP ) 00198 20 CONTINUE 00199 SPLUS = SPLUS + ABS( WORK( I ) ) 00200 SMINU = SMINU + ABS( RHS( I ) ) 00201 30 CONTINUE 00202 IF( SPLUS.GT.SMINU ) 00203 $ CALL CCOPY( N, WORK, 1, RHS, 1 ) 00204 * 00205 * Apply the permutations JPIV to the computed solution (RHS) 00206 * 00207 CALL CLASWP( 1, RHS, LDZ, 1, N-1, JPIV, -1 ) 00208 * 00209 * Compute the sum of squares 00210 * 00211 CALL CLASSQ( N, RHS, 1, RDSCAL, RDSUM ) 00212 RETURN 00213 END IF 00214 * 00215 * ENTRY IJOB = 2 00216 * 00217 * Compute approximate nullvector XM of Z 00218 * 00219 CALL CGECON( 'I', N, Z, LDZ, ONE, RTEMP, WORK, RWORK, INFO ) 00220 CALL CCOPY( N, WORK( N+1 ), 1, XM, 1 ) 00221 * 00222 * Compute RHS 00223 * 00224 CALL CLASWP( 1, XM, LDZ, 1, N-1, IPIV, -1 ) 00225 TEMP = CONE / SQRT( CDOTC( N, XM, 1, XM, 1 ) ) 00226 CALL CSCAL( N, TEMP, XM, 1 ) 00227 CALL CCOPY( N, XM, 1, XP, 1 ) 00228 CALL CAXPY( N, CONE, RHS, 1, XP, 1 ) 00229 CALL CAXPY( N, -CONE, XM, 1, RHS, 1 ) 00230 CALL CGESC2( N, Z, LDZ, RHS, IPIV, JPIV, SCALE ) 00231 CALL CGESC2( N, Z, LDZ, XP, IPIV, JPIV, SCALE ) 00232 IF( SCASUM( N, XP, 1 ).GT.SCASUM( N, RHS, 1 ) ) 00233 $ CALL CCOPY( N, XP, 1, RHS, 1 ) 00234 * 00235 * Compute the sum of squares 00236 * 00237 CALL CLASSQ( N, RHS, 1, RDSCAL, RDSUM ) 00238 RETURN 00239 * 00240 * End of CLATDF 00241 * 00242 END
1.7.2