1SUBROUTINE PZBSSOLVER1( N, M, IM, JM, DESCM, LAMBDA, X, IX, JX, & 2 DESCX, WORK, LWORK, IWORK, LIWORK, INFO ) 3! 4IMPLICIT NONE 5! 6! .. Scalar Arguments .. 7INTEGER N, IM, JM, IX, JX, LWORK, LIWORK, INFO 8! .. 9! .. Array Arguments .. 10INTEGER DESCM( * ), DESCX( * ), IWORK( * ) 11DOUBLE PRECISION M( * ), LAMBDA( * ), WORK( * ) 12COMPLEX*16 X( * ) 13! .. 14! 15! Purpose 16! ======= 17! 18! PZBSSOLVER1() computes all eigenvalues and (both right and 19! left) eigenvectors of 2n-by-2n complex matrix 20! 21! H = [ A, B; 22! -conj(B), -conj(A) ], 23! 24! where A is n-by-n Hermitian, and B is n-by-n symmetric. 25! 26! On entry, the information of H is provided in the lower triangular 27! part of 28! 29! M = [ real(A)+real(B), imag(A)-imag(B); 30! -imag(A)-imag(B), real(A)-real(B) ]. 31! 32! The matrix M is required to be positive definite. 33! 34! The structure of H leads to the following properties. 35! 36! 1) H is diagonalizable with n pairs of real eigenvalues 37! (lambda_i, -lambda_i). 38! 39! 2) The eigenvectors of H has the block structure 40! 41! X = [ X_1, conj(X_2); Y = [ X_1, -conj(X_2); 42! X_2, conj(X_1) ], -X_2, conj(X_1) ], 43! 44! and satisfy that 45! 46! X_1**T * X_2 = X_2**T * X_1, 47! X_1**H * X_1 - X_2**H * X_2 = I, 48! Y**H * X = I, 49! H * X = X * diag(lambda, -lambda), 50! Y**H * H = diag(lambda, -lambda) * Y**H. 51! 52! On exit, only the positive eigenvalues and the corresponding right 53! eigenvectors are returned. The eigenvalues are sorted in ascending 54! order. The eigenvectors are normalized (i.e., X = [ X_1; X_2 ] with 55! X_1**H * X_1 - X_2**H * X_2 = I). 56! 57! M is destroyed on exit. 58! 59! Notes 60! ===== 61! 62! Each global data object is described by an associated description 63! vector. This vector stores the information required to establish 64! the mapping between an object element and its corresponding process 65! and memory location. 66! 67! Let A be a generic term for any 2D block cyclicly distributed array. 68! Such a global array has an associated description vector DESCA. 69! In the following comments, the character _ should be read as 70! "of the global array". 71! 72! NOTATION STORED IN EXPLANATION 73! --------------- -------------- -------------------------------------- 74! DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 75! DTYPE_A = 1. 76! CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 77! the BLACS process grid A is distribu- 78! ted over. The context itself is glo- 79! bal, but the handle (the integer 80! value) may vary. 81! M_A (global) DESCA( M_ ) The number of rows in the global 82! array A. 83! N_A (global) DESCA( N_ ) The number of columns in the global 84! array A. 85! MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 86! the rows of the array. 87! NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 88! the columns of the array. 89! RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 90! row of the array A is distributed. 91! CSRC_A (global) DESCA( CSRC_ ) The process column over which the 92! first column of the array A is 93! distributed. 94! LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 95! array. LLD_A >= MAX(1,LOCr(M_A)). 96! 97! Let K be the number of rows or columns of a distributed matrix, 98! and assume that its process grid has dimension p x q. 99! LOCr( K ) denotes the number of elements of K that a process 100! would receive if K were distributed over the p processes of its 101! process column. 102! Similarly, LOCc( K ) denotes the number of elements of K that a 103! process would receive if K were distributed over the q processes of 104! its process row. 105! The values of LOCr() and LOCc() may be determined via a call to the 106! ScaLAPACK tool function, NUMROC: 107! LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 108! LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 109! An upper bound for these quantities may be computed by: 110! LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 111! LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 112! 113! Arguments 114! ========= 115! 116! N (global input) INTEGER 117! The number of rows and columns of A and B. 118! The size of M and X are N*N and (2N)*N, respectively. 119! N >= 0. 120! 121! M (local input and output) DOUBLE PRECISION pointer into the 122! local memory to an array of dimension 123! (LLD_M, LOCc(JM+2*N-1)). 124! On entry, the symmetric positive definite matrix M. Only the 125! lower triangular part of M is used to define the elements of 126! the symmetric matrix. 127! On exit, all entries of M are destroyed. 128! 129! IM (global input) INTEGER 130! The row index in the global array M indicating the first 131! row of sub( M ). 132! 133! JM (global input) INTEGER 134! The column index in the global array M indicating the 135! first column of sub( M ). 136! 137! DESCM (global and local input) INTEGER array of dimension DLEN_. 138! The array descriptor for the distributed matrix M. 139! 140! LAMBDA (global output) DOUBLE PRECISION array, dimension (N) 141! On normal exit LAMBDA contains the positive eigenvalues of H 142! in ascending order. 143! 144! X (local output) COMPLEX*16 array, 145! global dimension (2N, N), 146! local dimension ( LLD_X, LOCc(JX+N-1) ) 147! On normal exit X contains the normalized right eigenvectors 148! of H corresponding to the positive eigenvalues. 149! 150! IX (global input) INTEGER 151! X`s global row index, which points to the beginning of the 152! submatrix which is to be operated on. 153! In this version, only IX = JX = 1 is supported. 154! 155! JX (global input) INTEGER 156! X`s global column index, which points to the beginning of 157! the submatrix which is to be operated on. 158! In this version, only IX = JX = 1 is supported. 159! 160! DESCX (global and local input) INTEGER array of dimension DLEN_. 161! The array descriptor for the distributed matrix X. 162! DESCX( CTXT_ ) must equal DESCA( CTXT_ ) 163! 164! WORK (local workspace/output) DOUBLE PRECISION array, 165! dimension (LWORK) 166! On output, WORK( 1 ) returns the minimal amount of workspace 167! needed to guarantee completion. 168! If the input parameters are incorrect, WORK( 1 ) may also be 169! incorrect. 170! 171! LWORK (local input) INTEGER 172! The length of the workspace array WORK. 173! If LWORK = -1, the LWORK is global input and a workspace 174! query is assumed; the routine only calculates the minimum 175! size for the WORK/IWORK array. The required workspace is 176! returned as the first element of WORK/IWORK and no error 177! message is issued by PXERBLA. 178! 179! IWORK (local workspace/output) INTEGER array, 180! dimension (LIWORK) 181! On output, IWORK( 1 ) returns the minimal amount of workspace 182! needed to guarantee completion. 183! If the input parameters are incorrect, IWORK( 1 ) may also be 184! incorrect. 185! 186! LIWORK (local input) INTEGER 187! The length of the workspace array IWORK. 188! If LIWORK = -1, the LIWORK is global input and a workspace 189! query is assumed; the routine only calculates the minimum 190! size for the WORK/IWORK array. The required workspace is 191! returned as the first element of WORK/IWORK and no error 192! message is issued by PXERBLA. 193! 194! INFO (global output) INTEGER 195! = 0: successful exit 196! < 0: If the i-th argument is an array and the j-entry had 197! an illegal value, then INFO = -(i*100+j), if the i-th 198! argument is a scalar and had an illegal value, then 199! INFO = -i. 200! > 0: The eigensolver did not converge. 201! 202! Alignment requirements 203! ====================== 204! 205! This subroutine requires M and X to be distributed consistently in 206! the sense that DESCM( : ) = DESCX( : ) except for the fourth entry, 207! despite that X is complex and M is real. 208! In addition, square blocks ( i.e., MB = NB ) are required. 209! 210! ===================================================================== 211! 212! .. Parameters .. 213INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, & 214 LLD_, MB_, M_, NB_, N_, RSRC_ 215PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, & 216 CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, & 217 RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 218DOUBLE PRECISION ZERO, ONE, TWO 219PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 220! .. 221! .. Local Scalars .. 222LOGICAL LQUERY 223INTEGER ICTXT, NPROCS, NPROW, NPCOL, MYROW, MYCOL, NB, & 224 TWON, I, LWKOPT, LIWKOPT, LLWORK, LLIWORK, & 225 LOCALMAT, MROWS, MCOLS, LLDM, DIMV, NSPLIT, & 226 INDD, INDE, INDTAU, INDU, INDV, INDW, & 227 INDGAP, INDWORK, INDIBLOCK, INDISPLIT, & 228 INDIFAIL, INDICLUSTR, INDIWORK, ITMP(1) 229DOUBLE PRECISION ABSTOL, ORFAC, DTMP(1) 230DOUBLE PRECISION T_CHOL, T_FORMW, T_TRIDIAG1, T_TRIDIAG2, & 231 T_DIAG1, T_DIAG2, T_DIAG3, T_VEC1, T_VEC2, & 232 T_PREP 233! .. 234! .. Local Arrays .. 235INTEGER DESCU( DLEN_ ), DESCV( DLEN_ ), DESCW( DLEN_ ), & 236 ITMP2( 14 ) 237DOUBLE PRECISION DTMP2( 7 ) 238! .. 239! .. Intrinsic Functions .. 240INTRINSIC DBLE, DSQRT 241! .. 242! .. External Functions .. 243EXTERNAL NUMROC, PDLAMCH, MPI_WTIME 244INTEGER NUMROC 245DOUBLE PRECISION PDLAMCH, MPI_WTIME 246! .. 247! .. External Subroutines .. 248EXTERNAL PDSYRK, PDTRMM, PDTRSM, PDLACPY, PDLASET, & 249 PDPOTRF, PDSTEBZ, PDSTEIN, PXERBLA, & 250 BLACS_GRIDINFO, CHK1MAT, & 251 PDSORTEIG, PZBSAUX1, PDSSTRD, PDORMTR, & 252 PDGEADD, PZLASCL 253! .. 254! .. Executable Statements .. 255! 256T_PREP = MPI_WTIME() 257INFO = 0 258TWON = 2*N 259ICTXT = DESCM( CTXT_ ) 260CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 261NPROCS = NPROW*NPCOL 262IF ( NPROW .EQ. -1 ) THEN 263 INFO = -( 500+CTXT_ ) 264END IF 265! 266! Test the input arguments. 267! 268IF ( INFO .EQ. 0 .AND. N .LT. 0 ) THEN 269 INFO = -1 270END IF 271IF ( INFO .EQ. 0 ) & 272 CALL CHK1MAT( TWON, 1, TWON, 1, IM, JM, DESCM, 5, INFO ) 273IF ( INFO .EQ. 0 .AND. DESCM( MB_ ) .NE. DESCM( NB_ ) ) & 274 INFO = -( 500+MB_ ) 275IF ( INFO .EQ. 0 ) & 276 CALL CHK1MAT( TWON, 1, N, 1, IX, JX, DESCX, 10, INFO ) 277IF ( INFO .EQ. 0 .AND. DESCX( MB_ ) .NE. DESCX( NB_ ) ) & 278 INFO = -( 1000+MB_ ) 279! 280! Compute required workspace. 281! 282IF ( INFO .EQ. 0 ) THEN 283! 284! Set up local indices for the workspace. 285! 286 LQUERY = LWORK .EQ. -1 .OR. LIWORK .EQ. -1 287 ABSTOL = TWO*PDLAMCH( ICTXT, 'S' ) 288 ORFAC = PDLAMCH( ICTXT, 'P' ) 289 NB = DESCM( NB_ ) 290 LLDM = DESCM( LLD_ ) 291 MROWS = NUMROC( TWON, NB, MYROW, 0, NPROW ) 292 MCOLS = NUMROC( TWON, NB, MYCOL, 0, NPCOL ) 293 LOCALMAT = MAX( NB, LLDM )*MCOLS 294 DESCW( 1:DLEN_ ) = DESCM( 1:DLEN_ ) 295 DESCU( 1:DLEN_ ) = DESCM( 1:DLEN_ ) 296 DESCV( 1:DLEN_ ) = DESCM( 1:DLEN_ ) 297 INDE = 1 298 INDW = INDE + TWON 299 INDU = INDW + MAX( TWON+NPROCS, LOCALMAT ) 300 INDV = INDU + LOCALMAT 301 INDWORK = INDV + MAX( TWON, LOCALMAT ) 302 INDD = INDW 303 INDTAU = INDV 304 INDGAP = INDW + TWON 305 LLWORK = LWORK - INDWORK + 1 306 INDIBLOCK = 1 307 INDISPLIT = INDIBLOCK + TWON 308 INDIFAIL = INDISPLIT + TWON 309 INDICLUSTR = INDIFAIL + N 310 INDIWORK = INDICLUSTR + 2*NPROCS 311 LLIWORK = LIWORK - INDIWORK + 1 312! 313! Estimate the workspace required by external subroutines. 314! 315 CALL PDSSTRD( 'L', TWON, DTMP, 1, 1, DESCW, DTMP, DTMP, WORK, & 316 -1, ITMP ) 317 LWKOPT = INT( WORK( 1 ) ) 318 CALL PDORMTR( 'R', 'L', 'N', TWON, TWON, DTMP, 1, 1, DESCW, & 319 DTMP, DTMP, 1, 1, DESCU, WORK, -1, ITMP ) 320 LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) ) 321 CALL PDSTEBZ( ICTXT, 'I', 'B', TWON, ZERO, ZERO, N+1, TWON, & 322 ABSTOL, DTMP, DTMP, DIMV, NSPLIT, LAMBDA, ITMP, ITMP, & 323 DTMP2, -1, ITMP2, -1, ITMP ) 324 LWKOPT = MAX( LWKOPT, INT( DTMP2( 1 ) ) ) 325 LIWKOPT = ITMP2( 1 ) 326 CALL PDSTEIN( TWON, DTMP, DTMP, N, LAMBDA, ITMP, ITMP, & 327 ORFAC, DTMP, 1, 1, DESCU, WORK, -1, IWORK, -1, ITMP, & 328 ITMP, DTMP, ITMP(1) ) 329 LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) ) 330 LIWKOPT = MAX( LIWKOPT, IWORK( 1 ) ) 331! 332 LWKOPT = INDWORK - 1 + MAX( LWKOPT, LOCALMAT ) 333 LIWKOPT = INDIWORK - 1 + LIWKOPT 334 IF ( .NOT. LQUERY ) THEN 335 IF ( LWORK .LT. LWKOPT ) THEN 336 INFO = -12 337 ELSE IF ( LIWORK .LT. LIWKOPT ) THEN 338 INFO = -14 339 END IF 340 END IF 341END IF 342! 343IF ( INFO .NE. 0 ) THEN 344 CALL PXERBLA( ICTXT, 'PZBSSOLVER1', -INFO ) 345 RETURN 346END IF 347WORK( 1 ) = DBLE( LWKOPT ) 348IWORK( 1 ) = LIWKOPT 349IF ( LQUERY ) & 350 RETURN 351! 352! Quick return if possible. 353! 354IF ( N .EQ. 0 ) & 355 RETURN 356T_PREP = MPI_WTIME() - T_PREP 357! IF ( MYROW+MYCOL .EQ. 0 ) 358! $ WRITE( *, * ) 't_prep = ', T_PREP, ';' 359! 360! Compute the Cholesky factorization M = L * L**T. 361! In case of failure, a general solver is needed. 362! 363T_CHOL = MPI_WTIME() 364CALL PDPOTRF( 'L', TWON, M, 1, 1, DESCM, ITMP ) 365T_CHOL = MPI_WTIME() - T_CHOL 366! IF ( MYROW+MYCOL .EQ. 0 ) 367! $ WRITE( *, * ) 't_chol = ', T_CHOL, ';' 368! CALL PDLAPRNT( TWON, TWON, M, 1, 1, DESCM, 0, 0, 'L', 6, 369! $ WORK( INDWORK ) ) 370! IF ( MYROW+MYCOL .EQ. 0 ) 371! $ WRITE( *, * ) "L=tril(L);" 372IF ( ITMP(1) .NE. 0 ) THEN 373 INFO = -2 374 RETURN 375END IF 376! 377! Explicitly formulate W = L**T * J * L, 378! where 379! 380! J = [ 0, I_n; 381! -I_n, 0 ] 382! 383! is the standard symplectic matrix. 384! Only the lower triangular part of W is filled by the correct 385! values. 386! 387T_FORMW = MPI_WTIME() 388CALL PDLASET( 'A', TWON, TWON, ZERO, ZERO, WORK( INDW ), 1, 1, & 389 DESCW ) 390CALL PDGEADD( 'T', N, N, -ONE, M, N+1, 1, DESCM, ONE, & 391 WORK( INDW ), 1, 1, DESCW ) 392CALL PDTRADD( 'U', 'T', N, N, -ONE, M, N+1, N+1, DESCM, ONE, & 393 WORK( INDW ), N+1, 1, DESCW ) 394CALL PDTRMM( 'R', 'L', 'N', 'N', TWON, N, ONE, M, 1, 1, DESCM, & 395 WORK( INDW ), 1, 1, DESCW ) 396CALL PDLACPY( 'A', N, N, WORK( INDW ), 1, 1, DESCW, & 397 WORK( INDV ), 1, 1, DESCV ) 398CALL PDGEADD( 'T', N, N, -ONE, WORK( INDV ), 1, 1, DESCV, & 399 ONE, WORK( INDW ), 1, 1, DESCW ) 400T_FORMW = MPI_WTIME() - T_FORMW 401! IF ( MYROW+MYCOL .EQ. 0 ) 402! $ WRITE( *, * ) 't_formw = ', T_FORMW, ';' 403! CALL PDLAPRNT( TWON, TWON, WORK( INDW ), 1, 1, DESCW, 0, 0, 'W', 404! $ 6, WORK( INDWORK ) ) 405! 406! Tridiagonalization: U**T * W * U = T. 407! 408T_TRIDIAG1 = MPI_WTIME() 409CALL PDSSTRD( 'L', TWON, WORK( INDW ), 1, 1, DESCW, WORK( INDE ), & 410 WORK( INDTAU ), WORK( INDWORK ), LLWORK, ITMP ) 411! CALL PDLAPRNT( TWON, TWON, WORK( INDW ), 1, 1, DESCW, 0, 0, 'T', 412! $ 6, WORK( INDWORK ) ) 413DO I = 1, TWON-1 414 CALL PDELGET( 'A', ' ', WORK( INDE + I-1 ), WORK( INDW ), I+1, & 415 I, DESCW ) 416 WORK( INDE + I-1 ) = -WORK( INDE + I-1 ) 417! IF ( MYROW+MYCOL .EQ. 0 ) 418! $ WRITE( *, * ) 'E(', I, ', 1) =', WORK( INDE + I-1 ), ';' 419END DO 420T_TRIDIAG1 = MPI_WTIME() - T_TRIDIAG1 421! IF ( MYROW+MYCOL .EQ. 0 ) 422! $ WRITE( *, * ) 't_tridiag1 = ', T_TRIDIAG1, ';' 423! 424! Formulate sqrt(1/2) * L * U. 425! 426T_TRIDIAG2 = MPI_WTIME() 427CALL PDLASET( 'A', TWON, TWON, ZERO, ZERO, WORK( INDU ), 1, 1, & 428 DESCU ) 429CALL PDTRADD( 'L', 'N', TWON, TWON, DSQRT( ONE/TWO ), M, 1, 1, & 430 DESCM, ZERO, WORK( INDU ), 1, 1, DESCU ) 431CALL PDORMTR( 'R', 'L', 'N', TWON, TWON, WORK( INDW ), 1, 1, & 432 DESCW, WORK( INDTAU ), WORK( INDU ), 1, 1, DESCU, & 433 WORK( INDWORK ), LLWORK, ITMP ) 434T_TRIDIAG2 = MPI_WTIME() - T_TRIDIAG2 435! IF ( MYROW+MYCOL .EQ. 0 ) 436! $ WRITE( *, * ) 't_tridiag2 = ', T_TRIDIAG2, ';' 437! CALL PDLAPRNT( TWON, TWON, WORK( INDU ), 1, 1, DESCU, 0, 0, 'U', 438! $ 6, WORK( INDWORK ) ) 439! 440! Diagonalize the tridiagonal matrix 441! 442! D**H * (-i*T) * D = V * diag(lambda) * V**T, 443! 444! where D = diag{1,i,i^2,i^3,...}. 445! Only the positive half of the eigenpairs are computed. 446! 447DO I = 0, TWON-1 448 WORK( INDD + I ) = ZERO 449END DO 450! 451T_DIAG1 = MPI_WTIME() 452CALL PDSTEBZ( ICTXT, 'I', 'B', TWON, ZERO, ZERO, N+1, TWON, & 453 ABSTOL, WORK( INDD ), WORK( INDE ), DIMV, NSPLIT, LAMBDA, & 454 IWORK( INDIBLOCK ), IWORK( INDISPLIT ), WORK( INDWORK ), & 455 LLWORK, IWORK( INDIWORK ), LLIWORK, ITMP ) 456T_DIAG1 = MPI_WTIME() - T_DIAG1 457! IF ( MYROW+MYCOL .EQ. 0 ) 458! $ WRITE( *, * ) 't_diag1 = ', T_DIAG1, ';' 459IF ( ITMP(1) .NE. 0 ) THEN 460 INFO = ITMP(1) 461 WRITE( *, * ), '% PDSTEBZ fails with INFO =', INFO 462 RETURN 463END IF 464! IF ( MYROW+MYCOL .EQ. 0 ) THEN 465! DO I = 1, TWON 466! WRITE( *, * ) 'lambda0(', I, ', 1) =', LAMBDA( I ), ';' 467! END DO 468! END IF 469! 470T_DIAG2 = MPI_WTIME() 471CALL PDSTEIN( TWON, WORK( INDD ), WORK( INDE ), DIMV, LAMBDA, & 472 IWORK( INDIBLOCK ), IWORK( INDISPLIT ), ORFAC, & 473 WORK( INDV ), 1, 1, DESCV, WORK( INDWORK ), LLWORK, & 474 IWORK( INDIWORK ), LLIWORK, IWORK( INDIFAIL ), & 475 IWORK( INDICLUSTR ), WORK( INDGAP ), ITMP(1) ) 476IF ( ITMP(1) .NE. 0 ) THEN 477 INFO = ITMP(1) 478 WRITE( *, * ), '% PDSTEIN fails with INFO =', INFO 479 RETURN 480END IF 481T_DIAG2 = MPI_WTIME() - T_DIAG2 482! IF ( MYROW+MYCOL .EQ. 0 ) 483! $ WRITE( *, * ) 't_diag2 = ', T_DIAG2, ';' 484T_DIAG3 = MPI_WTIME() 485CALL PDSORTEIG( TWON, N, LAMBDA, WORK( INDV ), 1, 1, DESCV, 0 ) 486! 487! Orthogonalize the eigenvectors. 488! 489CALL PDSYRK( 'L', 'T', N, TWON, ONE, WORK( INDV ), 1, 1, DESCV, & 490 ZERO, WORK( INDW ), 1, 1, DESCW ) 491CALL PDPOTRF( 'L', N, WORK( INDW ), 1, 1, DESCW, ITMP ) 492CALL PDTRSM( 'R', 'L', 'T', 'N', TWON, N, ONE, WORK( INDW ), & 493 1, 1, DESCW, WORK( INDV ), 1, 1, DESCV ) 494T_DIAG3 = MPI_WTIME() - T_DIAG3 495! IF ( MYROW+MYCOL .EQ. 0 ) 496! $ WRITE( *, * ) 't_diag3 = ', T_DIAG3, ';' 497! CALL PDLAPRNT( TWON, N, WORK( INDV ), 1, 1, DESCV, 0, 0, 'V0', 498! $ 6, WORK( INDWORK ) ) 499! 500! Restore the eigenvectors of H: 501! 502! X = Q * L**{-T} * U * D * V * Lambda**{1/2}, 503! Y = Q * L * U * D * V * Lambda**{-1/2}, 504! 505! where 506! 507! Q = sqrt(1/2) * [ I_n, -i*I_n; 508! I_n, i*I_n ]. 509! 510T_VEC1 = MPI_WTIME() 511! 512! Scale V by Lambda**{-1/2}. 513! 514DO I = 1, N 515 CALL PDSCAL( TWON, ONE/DSQRT( LAMBDA( I ) ), & 516 WORK( INDV ), 1, I, DESCV, 1 ) 517END DO 518T_VEC1 = MPI_WTIME() - T_VEC1 519! IF ( MYROW+MYCOL .EQ. 0 ) 520! $ WRITE( *, * ) 't_vec1 = ', T_VEC1, ';' 521! CALL PDLAPRNT( TWON, TWON, WORK( INDU ), 1, 1, DESCU, 0, 0, 'Z', 522! $ 6, WORK( INDWORK ) ) 523! 524! Formulate 525! 526! Y = Q * L * U * D * V * Lambda**{-1/2}, 527! X = Q * L**{-T} * U * D * V * Lambda**{1/2}, 528! 529! Once Y is calculated, X is constructed from Y through 530! 531! Y = [ Y1; Y2 ], X = [ Y1; -Y2 ]. 532! 533! Use WORK and M as workspace for real and imaginary parts, 534! respectively. 535! 536T_VEC2 = MPI_WTIME() 537CALL PZBSAUX1( N, WORK( INDV ), 1, 1, DESCV, X, IX, JX, DESCX, & 538 WORK( INDU ), 1, 1, DESCU, WORK( INDWORK ), 1, 1, DESCW, & 539 M, IM, JM, DESCM ) 540CALL PZLASCL( 'G', -ONE, ONE, N, N, X, IX+N, JX, DESCX, ITMP ) 541T_VEC2 = MPI_WTIME() - T_VEC2 542! IF ( MYROW+MYCOL .EQ. 0 ) 543! $ WRITE( *, * ) 't_vec2 = ', T_VEC2, ';' 544! CALL PZLAPRNT( TWON, N, X, IX, JX, DESCX, 0, 0, 'X', 6, 545! $ WORK( INDWORK ) ) 546! 547WORK( 1 ) = DBLE( LWKOPT ) 548IWORK( 1 ) = LIWKOPT 549! 550RETURN 551! 552! End of PZBSSOLVER1(). 553! 554END 555