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