![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLATDF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZLATDF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlatdf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlatdf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlatdf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZLATDF( 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 * COMPLEX*16 RHS( * ), Z( LDZ, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> ZLATDF 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 ZGETC2. 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 ZGETC2 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 ZGECON, 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 DOUBLE PRECISION array, dimension (LDZ, N) 00074 *> On entry, the LU part of the factorization of the n-by-n 00075 *> matrix Z computed by ZGETC2: 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 DOUBLE PRECISION 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 DOUBLE PRECISION 00095 *> On entry, the sum of squares of computed contributions to 00096 *> the Dif-estimate under computation by ZTGSYL, 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 ZTGSY2 is called by CTGSYL. 00102 *> \endverbatim 00103 *> 00104 *> \param[in,out] RDSCAL 00105 *> \verbatim 00106 *> RDSCAL is DOUBLE PRECISION 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 ZTGSY2 is called by 00112 *> ZTGSYL. 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 complex16OTHERauxiliary 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 *>\n 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 ZLATDF( 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 DOUBLE PRECISION RDSCAL, RDSUM 00180 * .. 00181 * .. Array Arguments .. 00182 INTEGER IPIV( * ), JPIV( * ) 00183 COMPLEX*16 RHS( * ), Z( LDZ, * ) 00184 * .. 00185 * 00186 * ===================================================================== 00187 * 00188 * .. Parameters .. 00189 INTEGER MAXDIM 00190 PARAMETER ( MAXDIM = 2 ) 00191 DOUBLE PRECISION ZERO, ONE 00192 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00193 COMPLEX*16 CONE 00194 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 00195 * .. 00196 * .. Local Scalars .. 00197 INTEGER I, INFO, J, K 00198 DOUBLE PRECISION RTEMP, SCALE, SMINU, SPLUS 00199 COMPLEX*16 BM, BP, PMONE, TEMP 00200 * .. 00201 * .. Local Arrays .. 00202 DOUBLE PRECISION RWORK( MAXDIM ) 00203 COMPLEX*16 WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM ) 00204 * .. 00205 * .. External Subroutines .. 00206 EXTERNAL ZAXPY, ZCOPY, ZGECON, ZGESC2, ZLASSQ, ZLASWP, 00207 $ ZSCAL 00208 * .. 00209 * .. External Functions .. 00210 DOUBLE PRECISION DZASUM 00211 COMPLEX*16 ZDOTC 00212 EXTERNAL DZASUM, ZDOTC 00213 * .. 00214 * .. Intrinsic Functions .. 00215 INTRINSIC ABS, DBLE, SQRT 00216 * .. 00217 * .. Executable Statements .. 00218 * 00219 IF( IJOB.NE.2 ) THEN 00220 * 00221 * Apply permutations IPIV to RHS 00222 * 00223 CALL ZLASWP( 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 + DBLE( ZDOTC( N-J, Z( J+1, J ), 1, Z( J+1, 00237 $ J ), 1 ) ) 00238 SMINU = DBLE( ZDOTC( N-J, Z( J+1, J ), 1, RHS( J+1 ), 1 ) ) 00239 SPLUS = SPLUS*DBLE( RHS( J ) ) 00240 IF( SPLUS.GT.SMINU ) THEN 00241 RHS( J ) = BP 00242 ELSE IF( SMINU.GT.SPLUS ) THEN 00243 RHS( J ) = BM 00244 ELSE 00245 * 00246 * In this case the updating sums are equal and we can 00247 * choose RHS(J) +1 or -1. The first time this happens we 00248 * choose -1, thereafter +1. This is a simple way to get 00249 * good estimates of matrices like Byers well-known example 00250 * (see [1]). (Not done in BSOLVE.) 00251 * 00252 RHS( J ) = RHS( J ) + PMONE 00253 PMONE = CONE 00254 END IF 00255 * 00256 * Compute the remaining r.h.s. 00257 * 00258 TEMP = -RHS( J ) 00259 CALL ZAXPY( N-J, TEMP, Z( J+1, J ), 1, RHS( J+1 ), 1 ) 00260 10 CONTINUE 00261 * 00262 * Solve for U- part, lockahead for RHS(N) = +-1. This is not done 00263 * In BSOLVE and will hopefully give us a better estimate because 00264 * any ill-conditioning of the original matrix is transfered to U 00265 * and not to L. U(N, N) is an approximation to sigma_min(LU). 00266 * 00267 CALL ZCOPY( N-1, RHS, 1, WORK, 1 ) 00268 WORK( N ) = RHS( N ) + CONE 00269 RHS( N ) = RHS( N ) - CONE 00270 SPLUS = ZERO 00271 SMINU = ZERO 00272 DO 30 I = N, 1, -1 00273 TEMP = CONE / Z( I, I ) 00274 WORK( I ) = WORK( I )*TEMP 00275 RHS( I ) = RHS( I )*TEMP 00276 DO 20 K = I + 1, N 00277 WORK( I ) = WORK( I ) - WORK( K )*( Z( I, K )*TEMP ) 00278 RHS( I ) = RHS( I ) - RHS( K )*( Z( I, K )*TEMP ) 00279 20 CONTINUE 00280 SPLUS = SPLUS + ABS( WORK( I ) ) 00281 SMINU = SMINU + ABS( RHS( I ) ) 00282 30 CONTINUE 00283 IF( SPLUS.GT.SMINU ) 00284 $ CALL ZCOPY( N, WORK, 1, RHS, 1 ) 00285 * 00286 * Apply the permutations JPIV to the computed solution (RHS) 00287 * 00288 CALL ZLASWP( 1, RHS, LDZ, 1, N-1, JPIV, -1 ) 00289 * 00290 * Compute the sum of squares 00291 * 00292 CALL ZLASSQ( N, RHS, 1, RDSCAL, RDSUM ) 00293 RETURN 00294 END IF 00295 * 00296 * ENTRY IJOB = 2 00297 * 00298 * Compute approximate nullvector XM of Z 00299 * 00300 CALL ZGECON( 'I', N, Z, LDZ, ONE, RTEMP, WORK, RWORK, INFO ) 00301 CALL ZCOPY( N, WORK( N+1 ), 1, XM, 1 ) 00302 * 00303 * Compute RHS 00304 * 00305 CALL ZLASWP( 1, XM, LDZ, 1, N-1, IPIV, -1 ) 00306 TEMP = CONE / SQRT( ZDOTC( N, XM, 1, XM, 1 ) ) 00307 CALL ZSCAL( N, TEMP, XM, 1 ) 00308 CALL ZCOPY( N, XM, 1, XP, 1 ) 00309 CALL ZAXPY( N, CONE, RHS, 1, XP, 1 ) 00310 CALL ZAXPY( N, -CONE, XM, 1, RHS, 1 ) 00311 CALL ZGESC2( N, Z, LDZ, RHS, IPIV, JPIV, SCALE ) 00312 CALL ZGESC2( N, Z, LDZ, XP, IPIV, JPIV, SCALE ) 00313 IF( DZASUM( N, XP, 1 ).GT.DZASUM( N, RHS, 1 ) ) 00314 $ CALL ZCOPY( N, XP, 1, RHS, 1 ) 00315 * 00316 * Compute the sum of squares 00317 * 00318 CALL ZLASSQ( N, RHS, 1, RDSCAL, RDSUM ) 00319 RETURN 00320 * 00321 * End of ZLATDF 00322 * 00323 END