![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SORBDB 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SORBDB + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorbdb.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorbdb.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorbdb.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SORBDB( 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 * REAL PHI( * ), THETA( * ) 00032 * REAL 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 *> SORBDB simultaneously bidiagonalizes the blocks of an M-by-M 00044 *> partitioned orthogonal matrix X: 00045 *> 00046 *> [ B11 | B12 0 0 ] 00047 *> [ X11 | X12 ] [ P1 | ] [ 0 | 0 -I 0 ] [ Q1 | ]**T 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 SORCSD 00055 *> for details.) 00056 *> 00057 *> The orthogonal 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 REAL array, dimension (LDX11,Q) 00108 *> On entry, the top-left block of the orthogonal 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 REAL array, dimension (LDX12,M-Q) 00128 *> On entry, the top-right block of the orthogonal 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 REAL array, dimension (LDX21,Q) 00148 *> On entry, the bottom-left block of the orthogonal 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 REAL array, dimension (LDX22,M-Q) 00166 *> On entry, the bottom-right block of the orthogonal 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 realOTHERcomputational 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 SORCSD for details. 00273 *> 00274 *> P1, P2, Q1, and Q2 are represented as products of elementary 00275 *> reflectors. See SORCSD for details on generating P1, P2, Q1, and Q2 00276 *> using SORGQR and SORGLQ. 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 SORBDB( 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 REAL PHI( * ), THETA( * ) 00302 REAL TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ), 00303 $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ), 00304 $ X21( LDX21, * ), X22( LDX22, * ) 00305 * .. 00306 * 00307 * ==================================================================== 00308 * 00309 * .. Parameters .. 00310 REAL REALONE 00311 PARAMETER ( REALONE = 1.0E0 ) 00312 REAL ONE 00313 PARAMETER ( ONE = 1.0E0 ) 00314 * .. 00315 * .. Local Scalars .. 00316 LOGICAL COLMAJOR, LQUERY 00317 INTEGER I, LWORKMIN, LWORKOPT 00318 REAL Z1, Z2, Z3, Z4 00319 * .. 00320 * .. External Subroutines .. 00321 EXTERNAL SAXPY, SLARF, SLARFGP, SSCAL, XERBLA 00322 * .. 00323 * .. External Functions .. 00324 REAL SNRM2 00325 LOGICAL LSAME 00326 EXTERNAL SNRM2, LSAME 00327 * .. 00328 * .. Intrinsic Functions 00329 INTRINSIC ATAN2, COS, MAX, SIN 00330 * .. 00331 * .. Executable Statements .. 00332 * 00333 * Test input arguments 00334 * 00335 INFO = 0 00336 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 00337 IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN 00338 Z1 = REALONE 00339 Z2 = REALONE 00340 Z3 = REALONE 00341 Z4 = REALONE 00342 ELSE 00343 Z1 = REALONE 00344 Z2 = -REALONE 00345 Z3 = REALONE 00346 Z4 = -REALONE 00347 END IF 00348 LQUERY = LWORK .EQ. -1 00349 * 00350 IF( M .LT. 0 ) THEN 00351 INFO = -3 00352 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 00353 INFO = -4 00354 ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR. 00355 $ Q .GT. M-Q ) THEN 00356 INFO = -5 00357 ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN 00358 INFO = -7 00359 ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN 00360 INFO = -7 00361 ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN 00362 INFO = -9 00363 ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN 00364 INFO = -9 00365 ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN 00366 INFO = -11 00367 ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN 00368 INFO = -11 00369 ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN 00370 INFO = -13 00371 ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN 00372 INFO = -13 00373 END IF 00374 * 00375 * Compute workspace 00376 * 00377 IF( INFO .EQ. 0 ) THEN 00378 LWORKOPT = M - Q 00379 LWORKMIN = M - Q 00380 WORK(1) = LWORKOPT 00381 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN 00382 INFO = -21 00383 END IF 00384 END IF 00385 IF( INFO .NE. 0 ) THEN 00386 CALL XERBLA( 'xORBDB', -INFO ) 00387 RETURN 00388 ELSE IF( LQUERY ) THEN 00389 RETURN 00390 END IF 00391 * 00392 * Handle column-major and row-major separately 00393 * 00394 IF( COLMAJOR ) THEN 00395 * 00396 * Reduce columns 1, ..., Q of X11, X12, X21, and X22 00397 * 00398 DO I = 1, Q 00399 * 00400 IF( I .EQ. 1 ) THEN 00401 CALL SSCAL( P-I+1, Z1, X11(I,I), 1 ) 00402 ELSE 00403 CALL SSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), 1 ) 00404 CALL SAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I,I-1), 00405 $ 1, X11(I,I), 1 ) 00406 END IF 00407 IF( I .EQ. 1 ) THEN 00408 CALL SSCAL( M-P-I+1, Z2, X21(I,I), 1 ) 00409 ELSE 00410 CALL SSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), 1 ) 00411 CALL SAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I,I-1), 00412 $ 1, X21(I,I), 1 ) 00413 END IF 00414 * 00415 THETA(I) = ATAN2( SNRM2( M-P-I+1, X21(I,I), 1 ), 00416 $ SNRM2( P-I+1, X11(I,I), 1 ) ) 00417 * 00418 CALL SLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) ) 00419 X11(I,I) = ONE 00420 CALL SLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) ) 00421 X21(I,I) = ONE 00422 * 00423 CALL SLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), 00424 $ X11(I,I+1), LDX11, WORK ) 00425 CALL SLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I), 00426 $ X12(I,I), LDX12, WORK ) 00427 CALL SLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I), 00428 $ X21(I,I+1), LDX21, WORK ) 00429 CALL SLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I), 00430 $ X22(I,I), LDX22, WORK ) 00431 * 00432 IF( I .LT. Q ) THEN 00433 CALL SSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I,I+1), 00434 $ LDX11 ) 00435 CALL SAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I,I+1), LDX21, 00436 $ X11(I,I+1), LDX11 ) 00437 END IF 00438 CALL SSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), LDX12 ) 00439 CALL SAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), LDX22, 00440 $ X12(I,I), LDX12 ) 00441 * 00442 IF( I .LT. Q ) 00443 $ PHI(I) = ATAN2( SNRM2( Q-I, X11(I,I+1), LDX11 ), 00444 $ SNRM2( M-Q-I+1, X12(I,I), LDX12 ) ) 00445 * 00446 IF( I .LT. Q ) THEN 00447 CALL SLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11, 00448 $ TAUQ1(I) ) 00449 X11(I,I+1) = ONE 00450 END IF 00451 CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, 00452 $ TAUQ2(I) ) 00453 X12(I,I) = ONE 00454 * 00455 IF( I .LT. Q ) THEN 00456 CALL SLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I), 00457 $ X11(I+1,I+1), LDX11, WORK ) 00458 CALL SLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I), 00459 $ X21(I+1,I+1), LDX21, WORK ) 00460 END IF 00461 CALL SLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 00462 $ X12(I+1,I), LDX12, WORK ) 00463 CALL SLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 00464 $ X22(I+1,I), LDX22, WORK ) 00465 * 00466 END DO 00467 * 00468 * Reduce columns Q + 1, ..., P of X12, X22 00469 * 00470 DO I = Q + 1, P 00471 * 00472 CALL SSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), LDX12 ) 00473 CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, 00474 $ TAUQ2(I) ) 00475 X12(I,I) = ONE 00476 * 00477 CALL SLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I), 00478 $ X12(I+1,I), LDX12, WORK ) 00479 IF( M-P-Q .GE. 1 ) 00480 $ CALL SLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12, 00481 $ TAUQ2(I), X22(Q+1,I), LDX22, WORK ) 00482 * 00483 END DO 00484 * 00485 * Reduce columns P + 1, ..., M - Q of X12, X22 00486 * 00487 DO I = 1, M - P - Q 00488 * 00489 CALL SSCAL( M-P-Q-I+1, Z2*Z4, X22(Q+I,P+I), LDX22 ) 00490 CALL SLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1), 00491 $ LDX22, TAUQ2(P+I) ) 00492 X22(Q+I,P+I) = ONE 00493 CALL SLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22, 00494 $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK ) 00495 * 00496 END DO 00497 * 00498 ELSE 00499 * 00500 * Reduce columns 1, ..., Q of X11, X12, X21, X22 00501 * 00502 DO I = 1, Q 00503 * 00504 IF( I .EQ. 1 ) THEN 00505 CALL SSCAL( P-I+1, Z1, X11(I,I), LDX11 ) 00506 ELSE 00507 CALL SSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), LDX11 ) 00508 CALL SAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I-1,I), 00509 $ LDX12, X11(I,I), LDX11 ) 00510 END IF 00511 IF( I .EQ. 1 ) THEN 00512 CALL SSCAL( M-P-I+1, Z2, X21(I,I), LDX21 ) 00513 ELSE 00514 CALL SSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), LDX21 ) 00515 CALL SAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I-1,I), 00516 $ LDX22, X21(I,I), LDX21 ) 00517 END IF 00518 * 00519 THETA(I) = ATAN2( SNRM2( M-P-I+1, X21(I,I), LDX21 ), 00520 $ SNRM2( P-I+1, X11(I,I), LDX11 ) ) 00521 * 00522 CALL SLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) ) 00523 X11(I,I) = ONE 00524 CALL SLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21, 00525 $ TAUP2(I) ) 00526 X21(I,I) = ONE 00527 * 00528 CALL SLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I), 00529 $ X11(I+1,I), LDX11, WORK ) 00530 CALL SLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I), 00531 $ X12(I,I), LDX12, WORK ) 00532 CALL SLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I), 00533 $ X21(I+1,I), LDX21, WORK ) 00534 CALL SLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21, 00535 $ TAUP2(I), X22(I,I), LDX22, WORK ) 00536 * 00537 IF( I .LT. Q ) THEN 00538 CALL SSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I+1,I), 1 ) 00539 CALL SAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I+1,I), 1, 00540 $ X11(I+1,I), 1 ) 00541 END IF 00542 CALL SSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), 1 ) 00543 CALL SAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), 1, 00544 $ X12(I,I), 1 ) 00545 * 00546 IF( I .LT. Q ) 00547 $ PHI(I) = ATAN2( SNRM2( Q-I, X11(I+1,I), 1 ), 00548 $ SNRM2( M-Q-I+1, X12(I,I), 1 ) ) 00549 * 00550 IF( I .LT. Q ) THEN 00551 CALL SLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) ) 00552 X11(I+1,I) = ONE 00553 END IF 00554 CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) ) 00555 X12(I,I) = ONE 00556 * 00557 IF( I .LT. Q ) THEN 00558 CALL SLARF( 'L', Q-I, P-I, X11(I+1,I), 1, TAUQ1(I), 00559 $ X11(I+1,I+1), LDX11, WORK ) 00560 CALL SLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, TAUQ1(I), 00561 $ X21(I+1,I+1), LDX21, WORK ) 00562 END IF 00563 CALL SLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I), 00564 $ X12(I,I+1), LDX12, WORK ) 00565 CALL SLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I), 00566 $ X22(I,I+1), LDX22, WORK ) 00567 * 00568 END DO 00569 * 00570 * Reduce columns Q + 1, ..., P of X12, X22 00571 * 00572 DO I = Q + 1, P 00573 * 00574 CALL SSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), 1 ) 00575 CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) ) 00576 X12(I,I) = ONE 00577 * 00578 CALL SLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I), 00579 $ X12(I,I+1), LDX12, WORK ) 00580 IF( M-P-Q .GE. 1 ) 00581 $ CALL SLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, TAUQ2(I), 00582 $ X22(I,Q+1), LDX22, WORK ) 00583 * 00584 END DO 00585 * 00586 * Reduce columns P + 1, ..., M - Q of X12, X22 00587 * 00588 DO I = 1, M - P - Q 00589 * 00590 CALL SSCAL( M-P-Q-I+1, Z2*Z4, X22(P+I,Q+I), 1 ) 00591 CALL SLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1, 00592 $ TAUQ2(P+I) ) 00593 X22(P+I,Q+I) = ONE 00594 * 00595 CALL SLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1, 00596 $ TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK ) 00597 * 00598 END DO 00599 * 00600 END IF 00601 * 00602 RETURN 00603 * 00604 * End of SORBDB 00605 * 00606 END 00607