 |
LAPACK
3.4.0
LAPACK: Linear Algebra PACKage
|
|
LAPACK
3.4.0
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZUNBDB 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZUNBDB + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zunbdb.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zunbdb.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zunbdb.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 00022 * X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, 00023 * TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER SIGNS, TRANS 00027 * INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P, 00028 * $ Q 00029 * .. 00030 * .. Array Arguments .. 00031 * DOUBLE PRECISION PHI( * ), THETA( * ) 00032 * COMPLEX*16 TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ), 00033 * $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ), 00034 * $ X21( LDX21, * ), X22( LDX22, * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> ZUNBDB simultaneously bidiagonalizes the blocks of an M-by-M 00044 *> partitioned unitary matrix X: 00045 *> 00046 *> [ B11 | B12 0 0 ] 00047 *> [ X11 | X12 ] [ P1 | ] [ 0 | 0 -I 0 ] [ Q1 | ]**H 00048 *> X = [-----------] = [---------] [----------------] [---------] . 00049 *> [ X21 | X22 ] [ | P2 ] [ B21 | B22 0 0 ] [ | Q2 ] 00050 *> [ 0 | 0 0 I ] 00051 *> 00052 *> X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is 00053 *> not the case, then X must be transposed and/or permuted. This can be 00054 *> done in constant time using the TRANS and SIGNS options. See ZUNCSD 00055 *> for details.) 00056 *> 00057 *> The unitary matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by- 00058 *> (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are 00059 *> represented implicitly by Householder vectors. 00060 *> 00061 *> B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented 00062 *> implicitly by angles THETA, PHI. 00063 *> \endverbatim 00064 * 00065 * Arguments: 00066 * ========== 00067 * 00068 *> \param[in] TRANS 00069 *> \verbatim 00070 *> TRANS is CHARACTER 00071 *> = 'T': X, U1, U2, V1T, and V2T are stored in row-major 00072 *> order; 00073 *> otherwise: X, U1, U2, V1T, and V2T are stored in column- 00074 *> major order. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] SIGNS 00078 *> \verbatim 00079 *> SIGNS is CHARACTER 00080 *> = 'O': The lower-left block is made nonpositive (the 00081 *> "other" convention); 00082 *> otherwise: The upper-right block is made nonpositive (the 00083 *> "default" convention). 00084 *> \endverbatim 00085 *> 00086 *> \param[in] M 00087 *> \verbatim 00088 *> M is INTEGER 00089 *> The number of rows and columns in X. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] P 00093 *> \verbatim 00094 *> P is INTEGER 00095 *> The number of rows in X11 and X12. 0 <= P <= M. 00096 *> \endverbatim 00097 *> 00098 *> \param[in] Q 00099 *> \verbatim 00100 *> Q is INTEGER 00101 *> The number of columns in X11 and X21. 0 <= Q <= 00102 *> MIN(P,M-P,M-Q). 00103 *> \endverbatim 00104 *> 00105 *> \param[in,out] X11 00106 *> \verbatim 00107 *> X11 is COMPLEX*16 array, dimension (LDX11,Q) 00108 *> On entry, the top-left block of the unitary matrix to be 00109 *> reduced. On exit, the form depends on TRANS: 00110 *> If TRANS = 'N', then 00111 *> the columns of tril(X11) specify reflectors for P1, 00112 *> the rows of triu(X11,1) specify reflectors for Q1; 00113 *> else TRANS = 'T', and 00114 *> the rows of triu(X11) specify reflectors for P1, 00115 *> the columns of tril(X11,-1) specify reflectors for Q1. 00116 *> \endverbatim 00117 *> 00118 *> \param[in] LDX11 00119 *> \verbatim 00120 *> LDX11 is INTEGER 00121 *> The leading dimension of X11. If TRANS = 'N', then LDX11 >= 00122 *> P; else LDX11 >= Q. 00123 *> \endverbatim 00124 *> 00125 *> \param[in,out] X12 00126 *> \verbatim 00127 *> X12 is COMPLEX*16 array, dimension (LDX12,M-Q) 00128 *> On entry, the top-right block of the unitary matrix to 00129 *> be reduced. On exit, the form depends on TRANS: 00130 *> If TRANS = 'N', then 00131 *> the rows of triu(X12) specify the first P reflectors for 00132 *> Q2; 00133 *> else TRANS = 'T', and 00134 *> the columns of tril(X12) specify the first P reflectors 00135 *> for Q2. 00136 *> \endverbatim 00137 *> 00138 *> \param[in] LDX12 00139 *> \verbatim 00140 *> LDX12 is INTEGER 00141 *> The leading dimension of X12. If TRANS = 'N', then LDX12 >= 00142 *> P; else LDX11 >= M-Q. 00143 *> \endverbatim 00144 *> 00145 *> \param[in,out] X21 00146 *> \verbatim 00147 *> X21 is COMPLEX*16 array, dimension (LDX21,Q) 00148 *> On entry, the bottom-left block of the unitary matrix to 00149 *> be reduced. On exit, the form depends on TRANS: 00150 *> If TRANS = 'N', then 00151 *> the columns of tril(X21) specify reflectors for P2; 00152 *> else TRANS = 'T', and 00153 *> the rows of triu(X21) specify reflectors for P2. 00154 *> \endverbatim 00155 *> 00156 *> \param[in] LDX21 00157 *> \verbatim 00158 *> LDX21 is INTEGER 00159 *> The leading dimension of X21. If TRANS = 'N', then LDX21 >= 00160 *> M-P; else LDX21 >= Q. 00161 *> \endverbatim 00162 *> 00163 *> \param[in,out] X22 00164 *> \verbatim 00165 *> X22 is COMPLEX*16 array, dimension (LDX22,M-Q) 00166 *> On entry, the bottom-right block of the unitary matrix to 00167 *> be reduced. On exit, the form depends on TRANS: 00168 *> If TRANS = 'N', then 00169 *> the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last 00170 *> M-P-Q reflectors for Q2, 00171 *> else TRANS = 'T', and 00172 *> the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last 00173 *> M-P-Q reflectors for P2. 00174 *> \endverbatim 00175 *> 00176 *> \param[in] LDX22 00177 *> \verbatim 00178 *> LDX22 is INTEGER 00179 *> The leading dimension of X22. If TRANS = 'N', then LDX22 >= 00180 *> M-P; else LDX22 >= M-Q. 00181 *> \endverbatim 00182 *> 00183 *> \param[out] THETA 00184 *> \verbatim 00185 *> THETA is DOUBLE PRECISION array, dimension (Q) 00186 *> The entries of the bidiagonal blocks B11, B12, B21, B22 can 00187 *> be computed from the angles THETA and PHI. See Further 00188 *> Details. 00189 *> \endverbatim 00190 *> 00191 *> \param[out] PHI 00192 *> \verbatim 00193 *> PHI is DOUBLE PRECISION array, dimension (Q-1) 00194 *> The entries of the bidiagonal blocks B11, B12, B21, B22 can 00195 *> be computed from the angles THETA and PHI. See Further 00196 *> Details. 00197 *> \endverbatim 00198 *> 00199 *> \param[out] TAUP1 00200 *> \verbatim 00201 *> TAUP1 is COMPLEX*16 array, dimension (P) 00202 *> The scalar factors of the elementary reflectors that define 00203 *> P1. 00204 *> \endverbatim 00205 *> 00206 *> \param[out] TAUP2 00207 *> \verbatim 00208 *> TAUP2 is COMPLEX*16 array, dimension (M-P) 00209 *> The scalar factors of the elementary reflectors that define 00210 *> P2. 00211 *> \endverbatim 00212 *> 00213 *> \param[out] TAUQ1 00214 *> \verbatim 00215 *> TAUQ1 is COMPLEX*16 array, dimension (Q) 00216 *> The scalar factors of the elementary reflectors that define 00217 *> Q1. 00218 *> \endverbatim 00219 *> 00220 *> \param[out] TAUQ2 00221 *> \verbatim 00222 *> TAUQ2 is COMPLEX*16 array, dimension (M-Q) 00223 *> The scalar factors of the elementary reflectors that define 00224 *> Q2. 00225 *> \endverbatim 00226 *> 00227 *> \param[out] WORK 00228 *> \verbatim 00229 *> WORK is COMPLEX*16 array, dimension (LWORK) 00230 *> \endverbatim 00231 *> 00232 *> \param[in] LWORK 00233 *> \verbatim 00234 *> LWORK is INTEGER 00235 *> The dimension of the array WORK. LWORK >= M-Q. 00236 *> 00237 *> If LWORK = -1, then a workspace query is assumed; the routine 00238 *> only calculates the optimal size of the WORK array, returns 00239 *> this value as the first entry of the WORK array, and no error 00240 *> message related to LWORK is issued by XERBLA. 00241 *> \endverbatim 00242 *> 00243 *> \param[out] INFO 00244 *> \verbatim 00245 *> INFO is INTEGER 00246 *> = 0: successful exit. 00247 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00248 *> \endverbatim 00249 * 00250 * Authors: 00251 * ======== 00252 * 00253 *> \author Univ. of Tennessee 00254 *> \author Univ. of California Berkeley 00255 *> \author Univ. of Colorado Denver 00256 *> \author NAG Ltd. 00257 * 00258 *> \date November 2011 00259 * 00260 *> \ingroup complex16OTHERcomputational 00261 * 00262 *> \par Further Details: 00263 * ===================== 00264 *> 00265 *> \verbatim 00266 *> 00267 *> The bidiagonal blocks B11, B12, B21, and B22 are represented 00268 *> implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ..., 00269 *> PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are 00270 *> lower bidiagonal. Every entry in each bidiagonal band is a product 00271 *> of a sine or cosine of a THETA with a sine or cosine of a PHI. See 00272 *> [1] or ZUNCSD for details. 00273 *> 00274 *> P1, P2, Q1, and Q2 are represented as products of elementary 00275 *> reflectors. See ZUNCSD for details on generating P1, P2, Q1, and Q2 00276 *> using ZUNGQR and ZUNGLQ. 00277 *> \endverbatim 00278 * 00279 *> \par References: 00280 * ================ 00281 *> 00282 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. 00283 *> Algorithms, 50(1):33-65, 2009. 00284 *> 00285 * ===================================================================== 00286 SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 00287 $ X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, 00288 $ TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO ) 00289 * 00290 * -- LAPACK computational routine (version 3.4.0) -- 00291 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00292 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00293 * November 2011 00294 * 00295 * .. Scalar Arguments .. 00296 CHARACTER SIGNS, TRANS 00297 INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P, 00298 $ Q 00299 * .. 00300 * .. Array Arguments .. 00301 DOUBLE PRECISION PHI( * ), THETA( * ) 00302 COMPLEX*16 TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ), 00303 $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ), 00304 $ X21( LDX21, * ), X22( LDX22, * ) 00305 * .. 00306 * 00307 * ==================================================================== 00308 * 00309 * .. Parameters .. 00310 DOUBLE PRECISION REALONE 00311 PARAMETER ( REALONE = 1.0D0 ) 00312 COMPLEX*16 NEGONE, ONE 00313 PARAMETER ( NEGONE = (-1.0D0,0.0D0), 00314 $ ONE = (1.0D0,0.0D0) ) 00315 * .. 00316 * .. Local Scalars .. 00317 LOGICAL COLMAJOR, LQUERY 00318 INTEGER I, LWORKMIN, LWORKOPT 00319 DOUBLE PRECISION Z1, Z2, Z3, Z4 00320 * .. 00321 * .. External Subroutines .. 00322 EXTERNAL ZAXPY, ZLARF, ZLARFGP, ZSCAL, XERBLA 00323 EXTERNAL ZLACGV 00324 * 00325 * .. 00326 * .. External Functions .. 00327 DOUBLE PRECISION DZNRM2 00328 LOGICAL LSAME 00329 EXTERNAL DZNRM2, LSAME 00330 * .. 00331 * .. Intrinsic Functions 00332 INTRINSIC ATAN2, COS, MAX, MIN, SIN 00333 INTRINSIC DCMPLX, DCONJG 00334 * .. 00335 * .. Executable Statements .. 00336 * 00337 * Test input arguments 00338 * 00339 INFO = 0 00340 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 00341 IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN 00342 Z1 = REALONE 00343 Z2 = REALONE 00344 Z3 = REALONE 00345 Z4 = REALONE 00346 ELSE 00347 Z1 = REALONE 00348 Z2 = -REALONE 00349 Z3 = REALONE 00350 Z4 = -REALONE 00351 END IF 00352 LQUERY = LWORK .EQ. -1 00353 * 00354 IF( M .LT. 0 ) THEN 00355 INFO = -3 00356 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 00357 INFO = -4 00358 ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR. 00359 $ Q .GT. M-Q ) THEN 00360 INFO = -5 00361 ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN 00362 INFO = -7 00363 ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN 00364 INFO = -7 00365 ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN 00366 INFO = -9 00367 ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN 00368 INFO = -9 00369 ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN 00370 INFO = -11 00371 ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN 00372 INFO = -11 00373 ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN 00374 INFO = -13 00375 ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN 00376 INFO = -13 00377 END IF 00378 * 00379 * Compute workspace 00380 * 00381 IF( INFO .EQ. 0 ) THEN 00382 LWORKOPT = M - Q 00383 LWORKMIN = M - Q 00384 WORK(1) = LWORKOPT 00385 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN 00386 INFO = -21 00387 END IF 00388 END IF 00389 IF( INFO .NE. 0 ) THEN 00390 CALL XERBLA( 'xORBDB', -INFO ) 00391 RETURN 00392 ELSE IF( LQUERY ) THEN 00393 RETURN 00394 END IF 00395 * 00396 * Handle column-major and row-major separately 00397 * 00398 IF( COLMAJOR ) THEN 00399 * 00400 * Reduce columns 1, ..., Q of X11, X12, X21, and X22 00401 * 00402 DO I = 1, Q 00403 * 00404 IF( I .EQ. 1 ) THEN 00405 CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I), 1 ) 00406 ELSE 00407 CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ), 00408 $ X11(I,I), 1 ) 00409 CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)), 00410 $ 0.0D0 ), X12(I,I-1), 1, X11(I,I), 1 ) 00411 END IF 00412 IF( I .EQ. 1 ) THEN 00413 CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I), 1 ) 00414 ELSE 00415 CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ), 00416 $ X21(I,I), 1 ) 00417 CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)), 00418 $ 0.0D0 ), X22(I,I-1), 1, X21(I,I), 1 ) 00419 END IF 00420 * 00421 THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), 1 ), 00422 $ DZNRM2( P-I+1, X11(I,I), 1 ) ) 00423 * 00424 CALL ZLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) ) 00425 X11(I,I) = ONE 00426 CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) ) 00427 X21(I,I) = ONE 00428 * 00429 CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, DCONJG(TAUP1(I)), 00430 $ X11(I,I+1), LDX11, WORK ) 00431 CALL ZLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, 00432 $ DCONJG(TAUP1(I)), X12(I,I), LDX12, WORK ) 00433 CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, 00434 $ DCONJG(TAUP2(I)), X21(I,I+1), LDX21, WORK ) 00435 CALL ZLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, 00436 $ DCONJG(TAUP2(I)), X22(I,I), LDX22, WORK ) 00437 * 00438 IF( I .LT. Q ) THEN 00439 CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ), 00440 $ X11(I,I+1), LDX11 ) 00441 CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ), 00442 $ X21(I,I+1), LDX21, X11(I,I+1), LDX11 ) 00443 END IF 00444 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ), 00445 $ X12(I,I), LDX12 ) 00446 CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ), 00447 $ X22(I,I), LDX22, X12(I,I), LDX12 ) 00448 * 00449 IF( I .LT. Q ) 00450 $ PHI(I) = ATAN2( DZNRM2( Q-I, X11(I,I+1), LDX11 ), 00451 $ DZNRM2( M-Q-I+1, X12(I,I), LDX12 ) ) 00452 * 00453 IF( I .LT. Q ) THEN 00454 CALL ZLACGV( Q-I, X11(I,I+1), LDX11 ) 00455 CALL ZLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11, 00456 $ TAUQ1(I) ) 00457 X11(I,I+1) = ONE 00458 END IF 00459 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) 00460 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, 00461 $ TAUQ2(I) ) 00462 X12(I,I) = ONE 00463 * 00464 IF( I .LT. Q ) THEN 00465 CALL ZLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I), 00466 $ X11(I+1,I+1), LDX11, WORK ) 00467 CALL ZLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I), 00468 $ X21(I+1,I+1), LDX21, WORK ) 00469 END IF 00470 CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 00471 $ X12(I+1,I), LDX12, WORK ) 00472 CALL ZLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 00473 $ X22(I+1,I), LDX22, WORK ) 00474 * 00475 IF( I .LT. Q ) 00476 $ CALL ZLACGV( Q-I, X11(I,I+1), LDX11 ) 00477 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) 00478 * 00479 END DO 00480 * 00481 * Reduce columns Q + 1, ..., P of X12, X22 00482 * 00483 DO I = Q + 1, P 00484 * 00485 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I), 00486 $ LDX12 ) 00487 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) 00488 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, 00489 $ TAUQ2(I) ) 00490 X12(I,I) = ONE 00491 * 00492 CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 00493 $ X12(I+1,I), LDX12, WORK ) 00494 IF( M-P-Q .GE. 1 ) 00495 $ CALL ZLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12, 00496 $ TAUQ2(I), X22(Q+1,I), LDX22, WORK ) 00497 * 00498 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) 00499 * 00500 END DO 00501 * 00502 * Reduce columns P + 1, ..., M - Q of X12, X22 00503 * 00504 DO I = 1, M - P - Q 00505 * 00506 CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ), 00507 $ X22(Q+I,P+I), LDX22 ) 00508 CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 ) 00509 CALL ZLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1), 00510 $ LDX22, TAUQ2(P+I) ) 00511 X22(Q+I,P+I) = ONE 00512 CALL ZLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22, 00513 $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK ) 00514 * 00515 CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 ) 00516 * 00517 END DO 00518 * 00519 ELSE 00520 * 00521 * Reduce columns 1, ..., Q of X11, X12, X21, X22 00522 * 00523 DO I = 1, Q 00524 * 00525 IF( I .EQ. 1 ) THEN 00526 CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I), 00527 $ LDX11 ) 00528 ELSE 00529 CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ), 00530 $ X11(I,I), LDX11 ) 00531 CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)), 00532 $ 0.0D0 ), X12(I-1,I), LDX12, X11(I,I), LDX11 ) 00533 END IF 00534 IF( I .EQ. 1 ) THEN 00535 CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I), 00536 $ LDX21 ) 00537 ELSE 00538 CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ), 00539 $ X21(I,I), LDX21 ) 00540 CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)), 00541 $ 0.0D0 ), X22(I-1,I), LDX22, X21(I,I), LDX21 ) 00542 END IF 00543 * 00544 THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), LDX21 ), 00545 $ DZNRM2( P-I+1, X11(I,I), LDX11 ) ) 00546 * 00547 CALL ZLACGV( P-I+1, X11(I,I), LDX11 ) 00548 CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 ) 00549 * 00550 CALL ZLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) ) 00551 X11(I,I) = ONE 00552 CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21, 00553 $ TAUP2(I) ) 00554 X21(I,I) = ONE 00555 * 00556 CALL ZLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I), 00557 $ X11(I+1,I), LDX11, WORK ) 00558 CALL ZLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I), 00559 $ X12(I,I), LDX12, WORK ) 00560 CALL ZLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I), 00561 $ X21(I+1,I), LDX21, WORK ) 00562 CALL ZLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21, 00563 $ TAUP2(I), X22(I,I), LDX22, WORK ) 00564 * 00565 CALL ZLACGV( P-I+1, X11(I,I), LDX11 ) 00566 CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 ) 00567 * 00568 IF( I .LT. Q ) THEN 00569 CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ), 00570 $ X11(I+1,I), 1 ) 00571 CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ), 00572 $ X21(I+1,I), 1, X11(I+1,I), 1 ) 00573 END IF 00574 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ), 00575 $ X12(I,I), 1 ) 00576 CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ), 00577 $ X22(I,I), 1, X12(I,I), 1 ) 00578 * 00579 IF( I .LT. Q ) 00580 $ PHI(I) = ATAN2( DZNRM2( Q-I, X11(I+1,I), 1 ), 00581 $ DZNRM2( M-Q-I+1, X12(I,I), 1 ) ) 00582 * 00583 IF( I .LT. Q ) THEN 00584 CALL ZLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) ) 00585 X11(I+1,I) = ONE 00586 END IF 00587 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) ) 00588 X12(I,I) = ONE 00589 * 00590 IF( I .LT. Q ) THEN 00591 CALL ZLARF( 'L', Q-I, P-I, X11(I+1,I), 1, 00592 $ DCONJG(TAUQ1(I)), X11(I+1,I+1), LDX11, WORK ) 00593 CALL ZLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, 00594 $ DCONJG(TAUQ1(I)), X21(I+1,I+1), LDX21, WORK ) 00595 END IF 00596 CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, 00597 $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK ) 00598 CALL ZLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, 00599 $ DCONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK ) 00600 * 00601 END DO 00602 * 00603 * Reduce columns Q + 1, ..., P of X12, X22 00604 * 00605 DO I = Q + 1, P 00606 * 00607 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I), 1 ) 00608 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) ) 00609 X12(I,I) = ONE 00610 * 00611 CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, 00612 $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK ) 00613 IF( M-P-Q .GE. 1 ) 00614 $ CALL ZLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, 00615 $ DCONJG(TAUQ2(I)), X22(I,Q+1), LDX22, WORK ) 00616 * 00617 END DO 00618 * 00619 * Reduce columns P + 1, ..., M - Q of X12, X22 00620 * 00621 DO I = 1, M - P - Q 00622 * 00623 CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ), 00624 $ X22(P+I,Q+I), 1 ) 00625 CALL ZLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1, 00626 $ TAUQ2(P+I) ) 00627 X22(P+I,Q+I) = ONE 00628 * 00629 CALL ZLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1, 00630 $ DCONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22, 00631 $ WORK ) 00632 * 00633 END DO 00634 * 00635 END IF 00636 * 00637 RETURN 00638 * 00639 * End of ZUNBDB 00640 * 00641 END 00642
1.7.6.1 
