1*> \brief \b ZHETRS_3 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZHETRS_3 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrs_3.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrs_3.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrs_3.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZHETRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB, 22* INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER UPLO 26* INTEGER INFO, LDA, LDB, N, NRHS 27* .. 28* .. Array Arguments .. 29* INTEGER IPIV( * ) 30* COMPLEX*16 A( LDA, * ), B( LDB, * ), E( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> ZHETRS_3 solves a system of linear equations A * X = B with a complex 39*> Hermitian matrix A using the factorization computed 40*> by ZHETRF_RK or ZHETRF_BK: 41*> 42*> A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T), 43*> 44*> where U (or L) is unit upper (or lower) triangular matrix, 45*> U**H (or L**H) is the conjugate of U (or L), P is a permutation 46*> matrix, P**T is the transpose of P, and D is Hermitian and block 47*> diagonal with 1-by-1 and 2-by-2 diagonal blocks. 48*> 49*> This algorithm is using Level 3 BLAS. 50*> \endverbatim 51* 52* Arguments: 53* ========== 54* 55*> \param[in] UPLO 56*> \verbatim 57*> UPLO is CHARACTER*1 58*> Specifies whether the details of the factorization are 59*> stored as an upper or lower triangular matrix: 60*> = 'U': Upper triangular, form is A = P*U*D*(U**H)*(P**T); 61*> = 'L': Lower triangular, form is A = P*L*D*(L**H)*(P**T). 62*> \endverbatim 63*> 64*> \param[in] N 65*> \verbatim 66*> N is INTEGER 67*> The order of the matrix A. N >= 0. 68*> \endverbatim 69*> 70*> \param[in] NRHS 71*> \verbatim 72*> NRHS is INTEGER 73*> The number of right hand sides, i.e., the number of columns 74*> of the matrix B. NRHS >= 0. 75*> \endverbatim 76*> 77*> \param[in] A 78*> \verbatim 79*> A is COMPLEX*16 array, dimension (LDA,N) 80*> Diagonal of the block diagonal matrix D and factors U or L 81*> as computed by ZHETRF_RK and ZHETRF_BK: 82*> a) ONLY diagonal elements of the Hermitian block diagonal 83*> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k); 84*> (superdiagonal (or subdiagonal) elements of D 85*> should be provided on entry in array E), and 86*> b) If UPLO = 'U': factor U in the superdiagonal part of A. 87*> If UPLO = 'L': factor L in the subdiagonal part of A. 88*> \endverbatim 89*> 90*> \param[in] LDA 91*> \verbatim 92*> LDA is INTEGER 93*> The leading dimension of the array A. LDA >= max(1,N). 94*> \endverbatim 95*> 96*> \param[in] E 97*> \verbatim 98*> E is COMPLEX*16 array, dimension (N) 99*> On entry, contains the superdiagonal (or subdiagonal) 100*> elements of the Hermitian block diagonal matrix D 101*> with 1-by-1 or 2-by-2 diagonal blocks, where 102*> If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not referenced; 103*> If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced. 104*> 105*> NOTE: For 1-by-1 diagonal block D(k), where 106*> 1 <= k <= N, the element E(k) is not referenced in both 107*> UPLO = 'U' or UPLO = 'L' cases. 108*> \endverbatim 109*> 110*> \param[in] IPIV 111*> \verbatim 112*> IPIV is INTEGER array, dimension (N) 113*> Details of the interchanges and the block structure of D 114*> as determined by ZHETRF_RK or ZHETRF_BK. 115*> \endverbatim 116*> 117*> \param[in,out] B 118*> \verbatim 119*> B is COMPLEX*16 array, dimension (LDB,NRHS) 120*> On entry, the right hand side matrix B. 121*> On exit, the solution matrix X. 122*> \endverbatim 123*> 124*> \param[in] LDB 125*> \verbatim 126*> LDB is INTEGER 127*> The leading dimension of the array B. LDB >= max(1,N). 128*> \endverbatim 129*> 130*> \param[out] INFO 131*> \verbatim 132*> INFO is INTEGER 133*> = 0: successful exit 134*> < 0: if INFO = -i, the i-th argument had an illegal value 135*> \endverbatim 136* 137* Authors: 138* ======== 139* 140*> \author Univ. of Tennessee 141*> \author Univ. of California Berkeley 142*> \author Univ. of Colorado Denver 143*> \author NAG Ltd. 144* 145*> \ingroup complex16HEcomputational 146* 147*> \par Contributors: 148* ================== 149*> 150*> \verbatim 151*> 152*> June 2017, Igor Kozachenko, 153*> Computer Science Division, 154*> University of California, Berkeley 155*> 156*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, 157*> School of Mathematics, 158*> University of Manchester 159*> 160*> \endverbatim 161* 162* ===================================================================== 163 SUBROUTINE ZHETRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB, 164 $ INFO ) 165* 166* -- LAPACK computational routine -- 167* -- LAPACK is a software package provided by Univ. of Tennessee, -- 168* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 169* 170* .. Scalar Arguments .. 171 CHARACTER UPLO 172 INTEGER INFO, LDA, LDB, N, NRHS 173* .. 174* .. Array Arguments .. 175 INTEGER IPIV( * ) 176 COMPLEX*16 A( LDA, * ), B( LDB, * ), E( * ) 177* .. 178* 179* ===================================================================== 180* 181* .. Parameters .. 182 COMPLEX*16 ONE 183 PARAMETER ( ONE = ( 1.0D+0,0.0D+0 ) ) 184* .. 185* .. Local Scalars .. 186 LOGICAL UPPER 187 INTEGER I, J, K, KP 188 DOUBLE PRECISION S 189 COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM 190* .. 191* .. External Functions .. 192 LOGICAL LSAME 193 EXTERNAL LSAME 194* .. 195* .. External Subroutines .. 196 EXTERNAL ZDSCAL, ZSWAP, ZTRSM, XERBLA 197* .. 198* .. Intrinsic Functions .. 199 INTRINSIC ABS, DBLE, DCONJG, MAX 200* .. 201* .. Executable Statements .. 202* 203 INFO = 0 204 UPPER = LSAME( UPLO, 'U' ) 205 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 206 INFO = -1 207 ELSE IF( N.LT.0 ) THEN 208 INFO = -2 209 ELSE IF( NRHS.LT.0 ) THEN 210 INFO = -3 211 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 212 INFO = -5 213 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 214 INFO = -9 215 END IF 216 IF( INFO.NE.0 ) THEN 217 CALL XERBLA( 'ZHETRS_3', -INFO ) 218 RETURN 219 END IF 220* 221* Quick return if possible 222* 223 IF( N.EQ.0 .OR. NRHS.EQ.0 ) 224 $ RETURN 225* 226 IF( UPPER ) THEN 227* 228* Begin Upper 229* 230* Solve A*X = B, where A = U*D*U**H. 231* 232* P**T * B 233* 234* Interchange rows K and IPIV(K) of matrix B in the same order 235* that the formation order of IPIV(I) vector for Upper case. 236* 237* (We can do the simple loop over IPIV with decrement -1, 238* since the ABS value of IPIV(I) represents the row index 239* of the interchange with row i in both 1x1 and 2x2 pivot cases) 240* 241 DO K = N, 1, -1 242 KP = ABS( IPIV( K ) ) 243 IF( KP.NE.K ) THEN 244 CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 245 END IF 246 END DO 247* 248* Compute (U \P**T * B) -> B [ (U \P**T * B) ] 249* 250 CALL ZTRSM( 'L', 'U', 'N', 'U', N, NRHS, ONE, A, LDA, B, LDB ) 251* 252* Compute D \ B -> B [ D \ (U \P**T * B) ] 253* 254 I = N 255 DO WHILE ( I.GE.1 ) 256 IF( IPIV( I ).GT.0 ) THEN 257 S = DBLE( ONE ) / DBLE( A( I, I ) ) 258 CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB ) 259 ELSE IF ( I.GT.1 ) THEN 260 AKM1K = E( I ) 261 AKM1 = A( I-1, I-1 ) / AKM1K 262 AK = A( I, I ) / DCONJG( AKM1K ) 263 DENOM = AKM1*AK - ONE 264 DO J = 1, NRHS 265 BKM1 = B( I-1, J ) / AKM1K 266 BK = B( I, J ) / DCONJG( AKM1K ) 267 B( I-1, J ) = ( AK*BKM1-BK ) / DENOM 268 B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM 269 END DO 270 I = I - 1 271 END IF 272 I = I - 1 273 END DO 274* 275* Compute (U**H \ B) -> B [ U**H \ (D \ (U \P**T * B) ) ] 276* 277 CALL ZTRSM( 'L', 'U', 'C', 'U', N, NRHS, ONE, A, LDA, B, LDB ) 278* 279* P * B [ P * (U**H \ (D \ (U \P**T * B) )) ] 280* 281* Interchange rows K and IPIV(K) of matrix B in reverse order 282* from the formation order of IPIV(I) vector for Upper case. 283* 284* (We can do the simple loop over IPIV with increment 1, 285* since the ABS value of IPIV(I) represents the row index 286* of the interchange with row i in both 1x1 and 2x2 pivot cases) 287* 288 DO K = 1, N, 1 289 KP = ABS( IPIV( K ) ) 290 IF( KP.NE.K ) THEN 291 CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 292 END IF 293 END DO 294* 295 ELSE 296* 297* Begin Lower 298* 299* Solve A*X = B, where A = L*D*L**H. 300* 301* P**T * B 302* Interchange rows K and IPIV(K) of matrix B in the same order 303* that the formation order of IPIV(I) vector for Lower case. 304* 305* (We can do the simple loop over IPIV with increment 1, 306* since the ABS value of IPIV(I) represents the row index 307* of the interchange with row i in both 1x1 and 2x2 pivot cases) 308* 309 DO K = 1, N, 1 310 KP = ABS( IPIV( K ) ) 311 IF( KP.NE.K ) THEN 312 CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 313 END IF 314 END DO 315* 316* Compute (L \P**T * B) -> B [ (L \P**T * B) ] 317* 318 CALL ZTRSM( 'L', 'L', 'N', 'U', N, NRHS, ONE, A, LDA, B, LDB ) 319* 320* Compute D \ B -> B [ D \ (L \P**T * B) ] 321* 322 I = 1 323 DO WHILE ( I.LE.N ) 324 IF( IPIV( I ).GT.0 ) THEN 325 S = DBLE( ONE ) / DBLE( A( I, I ) ) 326 CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB ) 327 ELSE IF( I.LT.N ) THEN 328 AKM1K = E( I ) 329 AKM1 = A( I, I ) / DCONJG( AKM1K ) 330 AK = A( I+1, I+1 ) / AKM1K 331 DENOM = AKM1*AK - ONE 332 DO J = 1, NRHS 333 BKM1 = B( I, J ) / DCONJG( AKM1K ) 334 BK = B( I+1, J ) / AKM1K 335 B( I, J ) = ( AK*BKM1-BK ) / DENOM 336 B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 337 END DO 338 I = I + 1 339 END IF 340 I = I + 1 341 END DO 342* 343* Compute (L**H \ B) -> B [ L**H \ (D \ (L \P**T * B) ) ] 344* 345 CALL ZTRSM('L', 'L', 'C', 'U', N, NRHS, ONE, A, LDA, B, LDB ) 346* 347* P * B [ P * (L**H \ (D \ (L \P**T * B) )) ] 348* 349* Interchange rows K and IPIV(K) of matrix B in reverse order 350* from the formation order of IPIV(I) vector for Lower case. 351* 352* (We can do the simple loop over IPIV with decrement -1, 353* since the ABS value of IPIV(I) represents the row index 354* of the interchange with row i in both 1x1 and 2x2 pivot cases) 355* 356 DO K = N, 1, -1 357 KP = ABS( IPIV( K ) ) 358 IF( KP.NE.K ) THEN 359 CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 360 END IF 361 END DO 362* 363* END Lower 364* 365 END IF 366* 367 RETURN 368* 369* End of ZHETRS_3 370* 371 END 372