![]() |
LAPACK
3.4.1
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 ONE 00313 PARAMETER ( ONE = (1.0D0,0.0D0) ) 00314 * .. 00315 * .. Local Scalars .. 00316 LOGICAL COLMAJOR, LQUERY 00317 INTEGER I, LWORKMIN, LWORKOPT 00318 DOUBLE PRECISION Z1, Z2, Z3, Z4 00319 * .. 00320 * .. External Subroutines .. 00321 EXTERNAL ZAXPY, ZLARF, ZLARFGP, ZSCAL, XERBLA 00322 EXTERNAL ZLACGV 00323 * 00324 * .. 00325 * .. External Functions .. 00326 DOUBLE PRECISION DZNRM2 00327 LOGICAL LSAME 00328 EXTERNAL DZNRM2, LSAME 00329 * .. 00330 * .. Intrinsic Functions 00331 INTRINSIC ATAN2, COS, MAX, MIN, SIN 00332 INTRINSIC DCMPLX, DCONJG 00333 * .. 00334 * .. Executable Statements .. 00335 * 00336 * Test input arguments 00337 * 00338 INFO = 0 00339 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 00340 IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN 00341 Z1 = REALONE 00342 Z2 = REALONE 00343 Z3 = REALONE 00344 Z4 = REALONE 00345 ELSE 00346 Z1 = REALONE 00347 Z2 = -REALONE 00348 Z3 = REALONE 00349 Z4 = -REALONE 00350 END IF 00351 LQUERY = LWORK .EQ. -1 00352 * 00353 IF( M .LT. 0 ) THEN 00354 INFO = -3 00355 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 00356 INFO = -4 00357 ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR. 00358 $ Q .GT. M-Q ) THEN 00359 INFO = -5 00360 ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN 00361 INFO = -7 00362 ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN 00363 INFO = -7 00364 ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN 00365 INFO = -9 00366 ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN 00367 INFO = -9 00368 ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN 00369 INFO = -11 00370 ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN 00371 INFO = -11 00372 ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN 00373 INFO = -13 00374 ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN 00375 INFO = -13 00376 END IF 00377 * 00378 * Compute workspace 00379 * 00380 IF( INFO .EQ. 0 ) THEN 00381 LWORKOPT = M - Q 00382 LWORKMIN = M - Q 00383 WORK(1) = LWORKOPT 00384 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN 00385 INFO = -21 00386 END IF 00387 END IF 00388 IF( INFO .NE. 0 ) THEN 00389 CALL XERBLA( 'xORBDB', -INFO ) 00390 RETURN 00391 ELSE IF( LQUERY ) THEN 00392 RETURN 00393 END IF 00394 * 00395 * Handle column-major and row-major separately 00396 * 00397 IF( COLMAJOR ) THEN 00398 * 00399 * Reduce columns 1, ..., Q of X11, X12, X21, and X22 00400 * 00401 DO I = 1, Q 00402 * 00403 IF( I .EQ. 1 ) THEN 00404 CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I), 1 ) 00405 ELSE 00406 CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ), 00407 $ X11(I,I), 1 ) 00408 CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)), 00409 $ 0.0D0 ), X12(I,I-1), 1, X11(I,I), 1 ) 00410 END IF 00411 IF( I .EQ. 1 ) THEN 00412 CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I), 1 ) 00413 ELSE 00414 CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ), 00415 $ X21(I,I), 1 ) 00416 CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)), 00417 $ 0.0D0 ), X22(I,I-1), 1, X21(I,I), 1 ) 00418 END IF 00419 * 00420 THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), 1 ), 00421 $ DZNRM2( P-I+1, X11(I,I), 1 ) ) 00422 * 00423 CALL ZLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) ) 00424 X11(I,I) = ONE 00425 CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) ) 00426 X21(I,I) = ONE 00427 * 00428 CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, DCONJG(TAUP1(I)), 00429 $ X11(I,I+1), LDX11, WORK ) 00430 CALL ZLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, 00431 $ DCONJG(TAUP1(I)), X12(I,I), LDX12, WORK ) 00432 CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, 00433 $ DCONJG(TAUP2(I)), X21(I,I+1), LDX21, WORK ) 00434 CALL ZLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, 00435 $ DCONJG(TAUP2(I)), X22(I,I), LDX22, WORK ) 00436 * 00437 IF( I .LT. Q ) THEN 00438 CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ), 00439 $ X11(I,I+1), LDX11 ) 00440 CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ), 00441 $ X21(I,I+1), LDX21, X11(I,I+1), LDX11 ) 00442 END IF 00443 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ), 00444 $ X12(I,I), LDX12 ) 00445 CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ), 00446 $ X22(I,I), LDX22, X12(I,I), LDX12 ) 00447 * 00448 IF( I .LT. Q ) 00449 $ PHI(I) = ATAN2( DZNRM2( Q-I, X11(I,I+1), LDX11 ), 00450 $ DZNRM2( M-Q-I+1, X12(I,I), LDX12 ) ) 00451 * 00452 IF( I .LT. Q ) THEN 00453 CALL ZLACGV( Q-I, X11(I,I+1), LDX11 ) 00454 CALL ZLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11, 00455 $ TAUQ1(I) ) 00456 X11(I,I+1) = ONE 00457 END IF 00458 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) 00459 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, 00460 $ TAUQ2(I) ) 00461 X12(I,I) = ONE 00462 * 00463 IF( I .LT. Q ) THEN 00464 CALL ZLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I), 00465 $ X11(I+1,I+1), LDX11, WORK ) 00466 CALL ZLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I), 00467 $ X21(I+1,I+1), LDX21, WORK ) 00468 END IF 00469 CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 00470 $ X12(I+1,I), LDX12, WORK ) 00471 CALL ZLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 00472 $ X22(I+1,I), LDX22, WORK ) 00473 * 00474 IF( I .LT. Q ) 00475 $ CALL ZLACGV( Q-I, X11(I,I+1), LDX11 ) 00476 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) 00477 * 00478 END DO 00479 * 00480 * Reduce columns Q + 1, ..., P of X12, X22 00481 * 00482 DO I = Q + 1, P 00483 * 00484 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I), 00485 $ LDX12 ) 00486 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) 00487 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, 00488 $ TAUQ2(I) ) 00489 X12(I,I) = ONE 00490 * 00491 CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 00492 $ X12(I+1,I), LDX12, WORK ) 00493 IF( M-P-Q .GE. 1 ) 00494 $ CALL ZLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12, 00495 $ TAUQ2(I), X22(Q+1,I), LDX22, WORK ) 00496 * 00497 CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) 00498 * 00499 END DO 00500 * 00501 * Reduce columns P + 1, ..., M - Q of X12, X22 00502 * 00503 DO I = 1, M - P - Q 00504 * 00505 CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ), 00506 $ X22(Q+I,P+I), LDX22 ) 00507 CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 ) 00508 CALL ZLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1), 00509 $ LDX22, TAUQ2(P+I) ) 00510 X22(Q+I,P+I) = ONE 00511 CALL ZLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22, 00512 $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK ) 00513 * 00514 CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 ) 00515 * 00516 END DO 00517 * 00518 ELSE 00519 * 00520 * Reduce columns 1, ..., Q of X11, X12, X21, X22 00521 * 00522 DO I = 1, Q 00523 * 00524 IF( I .EQ. 1 ) THEN 00525 CALL ZSCAL( P-I+1, DCMPLX( Z1, 0.0D0 ), X11(I,I), 00526 $ LDX11 ) 00527 ELSE 00528 CALL ZSCAL( P-I+1, DCMPLX( Z1*COS(PHI(I-1)), 0.0D0 ), 00529 $ X11(I,I), LDX11 ) 00530 CALL ZAXPY( P-I+1, DCMPLX( -Z1*Z3*Z4*SIN(PHI(I-1)), 00531 $ 0.0D0 ), X12(I-1,I), LDX12, X11(I,I), LDX11 ) 00532 END IF 00533 IF( I .EQ. 1 ) THEN 00534 CALL ZSCAL( M-P-I+1, DCMPLX( Z2, 0.0D0 ), X21(I,I), 00535 $ LDX21 ) 00536 ELSE 00537 CALL ZSCAL( M-P-I+1, DCMPLX( Z2*COS(PHI(I-1)), 0.0D0 ), 00538 $ X21(I,I), LDX21 ) 00539 CALL ZAXPY( M-P-I+1, DCMPLX( -Z2*Z3*Z4*SIN(PHI(I-1)), 00540 $ 0.0D0 ), X22(I-1,I), LDX22, X21(I,I), LDX21 ) 00541 END IF 00542 * 00543 THETA(I) = ATAN2( DZNRM2( M-P-I+1, X21(I,I), LDX21 ), 00544 $ DZNRM2( P-I+1, X11(I,I), LDX11 ) ) 00545 * 00546 CALL ZLACGV( P-I+1, X11(I,I), LDX11 ) 00547 CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 ) 00548 * 00549 CALL ZLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) ) 00550 X11(I,I) = ONE 00551 CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21, 00552 $ TAUP2(I) ) 00553 X21(I,I) = ONE 00554 * 00555 CALL ZLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I), 00556 $ X11(I+1,I), LDX11, WORK ) 00557 CALL ZLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I), 00558 $ X12(I,I), LDX12, WORK ) 00559 CALL ZLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I), 00560 $ X21(I+1,I), LDX21, WORK ) 00561 CALL ZLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21, 00562 $ TAUP2(I), X22(I,I), LDX22, WORK ) 00563 * 00564 CALL ZLACGV( P-I+1, X11(I,I), LDX11 ) 00565 CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 ) 00566 * 00567 IF( I .LT. Q ) THEN 00568 CALL ZSCAL( Q-I, DCMPLX( -Z1*Z3*SIN(THETA(I)), 0.0D0 ), 00569 $ X11(I+1,I), 1 ) 00570 CALL ZAXPY( Q-I, DCMPLX( Z2*Z3*COS(THETA(I)), 0.0D0 ), 00571 $ X21(I+1,I), 1, X11(I+1,I), 1 ) 00572 END IF 00573 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4*SIN(THETA(I)), 0.0D0 ), 00574 $ X12(I,I), 1 ) 00575 CALL ZAXPY( M-Q-I+1, DCMPLX( Z2*Z4*COS(THETA(I)), 0.0D0 ), 00576 $ X22(I,I), 1, X12(I,I), 1 ) 00577 * 00578 IF( I .LT. Q ) 00579 $ PHI(I) = ATAN2( DZNRM2( Q-I, X11(I+1,I), 1 ), 00580 $ DZNRM2( M-Q-I+1, X12(I,I), 1 ) ) 00581 * 00582 IF( I .LT. Q ) THEN 00583 CALL ZLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) ) 00584 X11(I+1,I) = ONE 00585 END IF 00586 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) ) 00587 X12(I,I) = ONE 00588 * 00589 IF( I .LT. Q ) THEN 00590 CALL ZLARF( 'L', Q-I, P-I, X11(I+1,I), 1, 00591 $ DCONJG(TAUQ1(I)), X11(I+1,I+1), LDX11, WORK ) 00592 CALL ZLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, 00593 $ DCONJG(TAUQ1(I)), X21(I+1,I+1), LDX21, WORK ) 00594 END IF 00595 CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, 00596 $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK ) 00597 CALL ZLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, 00598 $ DCONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK ) 00599 * 00600 END DO 00601 * 00602 * Reduce columns Q + 1, ..., P of X12, X22 00603 * 00604 DO I = Q + 1, P 00605 * 00606 CALL ZSCAL( M-Q-I+1, DCMPLX( -Z1*Z4, 0.0D0 ), X12(I,I), 1 ) 00607 CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) ) 00608 X12(I,I) = ONE 00609 * 00610 CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, 00611 $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK ) 00612 IF( M-P-Q .GE. 1 ) 00613 $ CALL ZLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, 00614 $ DCONJG(TAUQ2(I)), X22(I,Q+1), LDX22, WORK ) 00615 * 00616 END DO 00617 * 00618 * Reduce columns P + 1, ..., M - Q of X12, X22 00619 * 00620 DO I = 1, M - P - Q 00621 * 00622 CALL ZSCAL( M-P-Q-I+1, DCMPLX( Z2*Z4, 0.0D0 ), 00623 $ X22(P+I,Q+I), 1 ) 00624 CALL ZLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1, 00625 $ TAUQ2(P+I) ) 00626 X22(P+I,Q+I) = ONE 00627 * 00628 CALL ZLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1, 00629 $ DCONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22, 00630 $ WORK ) 00631 * 00632 END DO 00633 * 00634 END IF 00635 * 00636 RETURN 00637 * 00638 * End of ZUNBDB 00639 * 00640 END 00641