Commit 4a7e92be authored by Volker Blum's avatar Volker Blum Committed by Volker Blum
Browse files

Added more verbose step-by-step output to test programs. This is still not optimal,

but will do for now.
parent 1b97c9e8
...@@ -1163,7 +1163,7 @@ subroutine tridiag_band_real(na, nb, nblk, a, lda, d, e, mpi_comm_rows, mpi_comm ...@@ -1163,7 +1163,7 @@ subroutine tridiag_band_real(na, nb, nblk, a, lda, d, e, mpi_comm_rows, mpi_comm
deallocate(block_limits) deallocate(block_limits)
deallocate(global_id) deallocate(global_id)
end subroutine end subroutine tridiag_band_real
! -------------------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------------------
......
...@@ -88,12 +88,20 @@ program test_complex ...@@ -88,12 +88,20 @@ program test_complex
call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols ) call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols )
call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol ) call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol )
if (myid==0) then
print '(a)','| Past BLACS_Gridinfo.'
end if
! All ELPA routines need MPI communicators for communicating within ! All ELPA routines need MPI communicators for communicating within
! rows or columns of processes, these are set in get_elpa_row_col_comms. ! rows or columns of processes, these are set in get_elpa_row_col_comms.
call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, & call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, &
mpi_comm_rows, mpi_comm_cols) mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| Past split communicator setup for rows and columns.'
end if
! Determine the necessary size of the distributed matrices, ! Determine the necessary size of the distributed matrices,
! we use the Scalapack tools routine NUMROC for that. ! we use the Scalapack tools routine NUMROC for that.
...@@ -107,6 +115,10 @@ program test_complex ...@@ -107,6 +115,10 @@ program test_complex
call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info ) call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info )
if (myid==0) then
print '(a)','| Past scalapack descriptor setup.'
end if
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! Allocate matrices and set up a test matrix for the eigenvalue problem ! Allocate matrices and set up a test matrix for the eigenvalue problem
...@@ -133,8 +145,17 @@ program test_complex ...@@ -133,8 +145,17 @@ program test_complex
deallocate(xr) deallocate(xr)
a(:,:) = z(:,:) a(:,:) = z(:,:)
if (myid==0) then
print '(a)','| Random matrix block has been set up. (only processor 0 confirms this step)'
end if
call pztranc(na, na, CONE, z, 1, 1, sc_desc, CONE, a, 1, 1, sc_desc) ! A = A + Z**H call pztranc(na, na, CONE, z, 1, 1, sc_desc, CONE, a, 1, 1, sc_desc) ! A = A + Z**H
if (myid==0) then
print '(a)','| Random matrix has been symmetrized.'
end if
! Save original matrix A for later accuracy checks ! Save original matrix A for later accuracy checks
as = a as = a
...@@ -142,10 +163,20 @@ program test_complex ...@@ -142,10 +163,20 @@ program test_complex
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! Calculate eigenvalues/eigenvectors ! Calculate eigenvalues/eigenvectors
if (myid==0) then
print '(a)','| Entering one-step ELPA solver ... '
print *
end if
call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only
call solve_evp_complex(na, nev, a, na_rows, ev, z, na_rows, nblk, & call solve_evp_complex(na, nev, a, na_rows, ev, z, na_rows, nblk, &
mpi_comm_rows, mpi_comm_cols) mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| One-step ELPA solver complete.'
print *
end if
if(myid == 0) print *,'Time tridiag_complex :',time_evp_fwd if(myid == 0) print *,'Time tridiag_complex :',time_evp_fwd
if(myid == 0) print *,'Time solve_tridi :',time_evp_solve if(myid == 0) print *,'Time solve_tridi :',time_evp_solve
if(myid == 0) print *,'Time trans_ev_complex:',time_evp_back if(myid == 0) print *,'Time trans_ev_complex:',time_evp_back
......
...@@ -91,12 +91,20 @@ program test_complex_gen ...@@ -91,12 +91,20 @@ program test_complex_gen
call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols ) call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols )
call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol ) call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol )
if (myid==0) then
print '(a)','| Past BLACS_Gridinfo.'
end if
! All ELPA routines need MPI communicators for communicating within ! All ELPA routines need MPI communicators for communicating within
! rows or columns of processes, these are set in get_elpa_row_col_comms. ! rows or columns of processes, these are set in get_elpa_row_col_comms.
call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, & call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, &
mpi_comm_rows, mpi_comm_cols) mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| Past split communicator setup for rows and columns.'
end if
! Determine the necessary size of the distributed matrices, ! Determine the necessary size of the distributed matrices,
! we use the Scalapack tools routine NUMROC for that. ! we use the Scalapack tools routine NUMROC for that.
...@@ -110,6 +118,10 @@ program test_complex_gen ...@@ -110,6 +118,10 @@ program test_complex_gen
call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info ) call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info )
if (myid==0) then
print '(a)','| Past scalapack descriptor setup.'
end if
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! Allocate matrices and set up test matrices for the eigenvalue problem ! Allocate matrices and set up test matrices for the eigenvalue problem
...@@ -141,8 +153,17 @@ program test_complex_gen ...@@ -141,8 +153,17 @@ program test_complex_gen
deallocate(xr) deallocate(xr)
a(:,:) = z(:,:) a(:,:) = z(:,:)
if (myid==0) then
print '(a)','| Random matrix block has been set up. (only processor 0 confirms this step)'
end if
call pztranc(na, na, CONE, z, 1, 1, sc_desc, CONE, a, 1, 1, sc_desc) ! A = A + Z**H call pztranc(na, na, CONE, z, 1, 1, sc_desc, CONE, a, 1, 1, sc_desc) ! A = A + Z**H
if (myid==0) then
print '(a)','| Random matrix has been Hermite-icized.'
end if
! The matrix B in the generalized eigenvalue problem must be symmetric ! The matrix B in the generalized eigenvalue problem must be symmetric
! and positive definite - we use a simple diagonally dominant matrix ! and positive definite - we use a simple diagonally dominant matrix
...@@ -150,6 +171,10 @@ program test_complex_gen ...@@ -150,6 +171,10 @@ program test_complex_gen
call pzlaset('U', na, na, xc, CONE, b, 1, 1, sc_desc ) ! Upper part call pzlaset('U', na, na, xc, CONE, b, 1, 1, sc_desc ) ! Upper part
call pzlaset('L', na, na, conjg(xc), CONE, b, 1, 1, sc_desc ) ! Lower part call pzlaset('L', na, na, conjg(xc), CONE, b, 1, 1, sc_desc ) ! Lower part
if (myid==0) then
print '(a)','| Complex Hermitian diagonally dominant overlap matrix has been initialized.'
end if
! Save original matrices A and B for later accuracy checks ! Save original matrices A and B for later accuracy checks
as = a as = a
...@@ -161,13 +186,22 @@ program test_complex_gen ...@@ -161,13 +186,22 @@ program test_complex_gen
! 1. Calculate Cholesky factorization of Matrix B = U**T * U ! 1. Calculate Cholesky factorization of Matrix B = U**T * U
! and invert triangular matrix U ! and invert triangular matrix U
! !
! Please note: cholesky_real/invert_trm_real are not trimmed for speed. ! Please note: cholesky_complex/invert_trm_complex are not trimmed for speed.
! The only reason having them is that the Scalapack counterpart ! The only reason having them is that the Scalapack counterpart
! PDPOTRF very often fails on higher processor numbers for unknown reasons! ! PDPOTRF very often fails on higher processor numbers for unknown reasons!
call cholesky_complex(na, b, na_rows, nblk, mpi_comm_rows, mpi_comm_cols) call cholesky_complex(na, b, na_rows, nblk, mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| Cholesky factorization complete.'
end if
call invert_trm_complex(na, b, na_rows, nblk, mpi_comm_rows, mpi_comm_cols) call invert_trm_complex(na, b, na_rows, nblk, mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| Cholesky factor inverted.'
end if
ttt0 = MPI_Wtime() ttt0 = MPI_Wtime()
! 2. Calculate U**-T * A * U**-1 ! 2. Calculate U**-T * A * U**-1
...@@ -183,6 +217,11 @@ program test_complex_gen ...@@ -183,6 +217,11 @@ program test_complex_gen
call mult_ah_b_complex('U', 'U', na, na, b, na_rows, tmp2, na_rows, & call mult_ah_b_complex('U', 'U', na, na, b, na_rows, tmp2, na_rows, &
nblk, mpi_comm_rows, mpi_comm_cols, a, na_rows) nblk, mpi_comm_rows, mpi_comm_cols, a, na_rows)
ttt1 = MPI_Wtime() ttt1 = MPI_Wtime()
if (myid==0) then
print '(a)','| Matrix A transformed from generalized to orthogonal form using Cholesky factors.'
end if
if(myid == 0) print *,'Time U**-T*A*U**-1 :',ttt1-ttt0 if(myid == 0) print *,'Time U**-T*A*U**-1 :',ttt1-ttt0
! A is only set in the upper half, solve_evp_real needs a full matrix ! A is only set in the upper half, solve_evp_real needs a full matrix
...@@ -190,6 +229,10 @@ program test_complex_gen ...@@ -190,6 +229,10 @@ program test_complex_gen
call pztranc(na,na,CONE,a,1,1,sc_desc,CZERO,tmp1,1,1,sc_desc) call pztranc(na,na,CONE,a,1,1,sc_desc,CZERO,tmp1,1,1,sc_desc)
if (myid==0) then
print '(a)','| Lower half of A set by pztranc.'
end if
do i=1,na_cols do i=1,na_cols
! Get global column corresponding to i and number of local rows up to ! Get global column corresponding to i and number of local rows up to
! and including the diagonal, these are unchanged in A ! and including the diagonal, these are unchanged in A
...@@ -201,9 +244,19 @@ program test_complex_gen ...@@ -201,9 +244,19 @@ program test_complex_gen
! 3. Calculate eigenvalues/eigenvectors of U**-T * A * U**-1 ! 3. Calculate eigenvalues/eigenvectors of U**-T * A * U**-1
! Eigenvectors go to tmp1 ! Eigenvectors go to tmp1
if (myid==0) then
print '(a)','| Entering one-step ELPA solver ... '
print *
end if
call solve_evp_complex(na, nev, a, na_rows, ev, tmp1, na_rows, nblk, & call solve_evp_complex(na, nev, a, na_rows, ev, tmp1, na_rows, nblk, &
mpi_comm_rows, mpi_comm_cols) mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| One-step ELPA solver complete.'
print *
end if
if(myid == 0) print *,'Time tridiag_complex :',time_evp_fwd if(myid == 0) print *,'Time tridiag_complex :',time_evp_fwd
if(myid == 0) print *,'Time solve_tridi :',time_evp_solve if(myid == 0) print *,'Time solve_tridi :',time_evp_solve
if(myid == 0) print *,'Time trans_ev_complex:',time_evp_back if(myid == 0) print *,'Time trans_ev_complex:',time_evp_back
...@@ -217,6 +270,9 @@ program test_complex_gen ...@@ -217,6 +270,9 @@ program test_complex_gen
call mult_ah_b_complex('L', 'N', na, nev, tmp2, na_rows, tmp1, na_rows, & call mult_ah_b_complex('L', 'N', na, nev, tmp2, na_rows, tmp1, na_rows, &
nblk, mpi_comm_rows, mpi_comm_cols, z, na_rows) nblk, mpi_comm_rows, mpi_comm_cols, z, na_rows)
ttt1 = MPI_Wtime() ttt1 = MPI_Wtime()
if (myid==0) then
print '(a)','| Backtransform of eigenvectors to generalized form complete.'
end if
if(myid == 0) print *,'Time Back U**-1*Z :',ttt1-ttt0 if(myid == 0) print *,'Time Back U**-1*Z :',ttt1-ttt0
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
......
...@@ -85,12 +85,20 @@ program test_real ...@@ -85,12 +85,20 @@ program test_real
call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols ) call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols )
call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol ) call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol )
if (myid==0) then
print '(a)','| Past BLACS_Gridinfo.'
end if
! All ELPA routines need MPI communicators for communicating within ! All ELPA routines need MPI communicators for communicating within
! rows or columns of processes, these are set in get_elpa_row_col_comms. ! rows or columns of processes, these are set in get_elpa_row_col_comms.
call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, & call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, &
mpi_comm_rows, mpi_comm_cols) mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| Past split communicator setup for rows and columns.'
end if
! Determine the necessary size of the distributed matrices, ! Determine the necessary size of the distributed matrices,
! we use the Scalapack tools routine NUMROC for that. ! we use the Scalapack tools routine NUMROC for that.
...@@ -104,6 +112,10 @@ program test_real ...@@ -104,6 +112,10 @@ program test_real
call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info ) call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info )
if (myid==0) then
print '(a)','| Past scalapack descriptor setup.'
end if
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! Allocate matrices and set up a test matrix for the eigenvalue problem ! Allocate matrices and set up a test matrix for the eigenvalue problem
...@@ -125,8 +137,17 @@ program test_real ...@@ -125,8 +137,17 @@ program test_real
call RANDOM_NUMBER(z) call RANDOM_NUMBER(z)
a(:,:) = z(:,:) a(:,:) = z(:,:)
if (myid==0) then
print '(a)','| Random matrix block has been set up. (only processor 0 confirms this step)'
end if
call pdtran(na, na, 1.d0, z, 1, 1, sc_desc, 1.d0, a, 1, 1, sc_desc) ! A = A + Z**T call pdtran(na, na, 1.d0, z, 1, 1, sc_desc, 1.d0, a, 1, 1, sc_desc) ! A = A + Z**T
if (myid==0) then
print '(a)','| Random matrix has been symmetrized.'
end if
! Save original matrix A for later accuracy checks ! Save original matrix A for later accuracy checks
as = a as = a
...@@ -134,10 +155,20 @@ program test_real ...@@ -134,10 +155,20 @@ program test_real
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! Calculate eigenvalues/eigenvectors ! Calculate eigenvalues/eigenvectors
if (myid==0) then
print '(a)','| Entering one-step ELPA solver ... '
print *
end if
call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only
call solve_evp_real(na, nev, a, na_rows, ev, z, na_rows, nblk, & call solve_evp_real(na, nev, a, na_rows, ev, z, na_rows, nblk, &
mpi_comm_rows, mpi_comm_cols) mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| One-step ELPA solver complete.'
print *
end if
if(myid == 0) print *,'Time tridiag_real :',time_evp_fwd if(myid == 0) print *,'Time tridiag_real :',time_evp_fwd
if(myid == 0) print *,'Time solve_tridi :',time_evp_solve if(myid == 0) print *,'Time solve_tridi :',time_evp_solve
if(myid == 0) print *,'Time trans_ev_real:',time_evp_back if(myid == 0) print *,'Time trans_ev_real:',time_evp_back
......
...@@ -86,12 +86,20 @@ program test_real2 ...@@ -86,12 +86,20 @@ program test_real2
call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols ) call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols )
call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol ) call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol )
if (myid==0) then
print '(a)','| Past BLACS_Gridinfo.'
end if
! All ELPA routines need MPI communicators for communicating within ! All ELPA routines need MPI communicators for communicating within
! rows or columns of processes, these are set in get_elpa_row_col_comms. ! rows or columns of processes, these are set in get_elpa_row_col_comms.
call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, & call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, &
mpi_comm_rows, mpi_comm_cols) mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| Past split communicator setup for rows and columns.'
end if
! Determine the necessary size of the distributed matrices, ! Determine the necessary size of the distributed matrices,
! we use the Scalapack tools routine NUMROC for that. ! we use the Scalapack tools routine NUMROC for that.
...@@ -105,6 +113,10 @@ program test_real2 ...@@ -105,6 +113,10 @@ program test_real2
call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info ) call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info )
if (myid==0) then
print '(a)','| Past scalapack descriptor setup.'
end if
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! Allocate matrices and set up a test matrix for the eigenvalue problem ! Allocate matrices and set up a test matrix for the eigenvalue problem
...@@ -126,8 +138,17 @@ program test_real2 ...@@ -126,8 +138,17 @@ program test_real2
call RANDOM_NUMBER(z) call RANDOM_NUMBER(z)
a(:,:) = z(:,:) a(:,:) = z(:,:)
if (myid==0) then
print '(a)','| Random matrix block has been set up. (only processor 0 confirms this step)'
end if
call pdtran(na, na, 1.d0, z, 1, 1, sc_desc, 1.d0, a, 1, 1, sc_desc) ! A = A + Z**T call pdtran(na, na, 1.d0, z, 1, 1, sc_desc, 1.d0, a, 1, 1, sc_desc) ! A = A + Z**T
if (myid==0) then
print '(a)','| Random matrix has been symmetrized.'
end if
! Save original matrix A for later accuracy checks ! Save original matrix A for later accuracy checks
as = a as = a
...@@ -138,10 +159,20 @@ program test_real2 ...@@ -138,10 +159,20 @@ program test_real2
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! Calculate eigenvalues/eigenvectors ! Calculate eigenvalues/eigenvectors
if (myid==0) then
print '(a)','| Entering two-stage ELPA solver ... '
print *
end if
call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only call mpi_barrier(mpi_comm_world, mpierr) ! for correct timings only
call solve_evp_real_2stage(na, nev, a, na_rows, ev, z, na_rows, nblk, & call solve_evp_real_2stage(na, nev, a, na_rows, ev, z, na_rows, nblk, &
mpi_comm_rows, mpi_comm_cols, mpi_comm_world) mpi_comm_rows, mpi_comm_cols, mpi_comm_world)
if (myid==0) then
print '(a)','| Two-step ELPA solver complete.'
print *
end if
if(myid == 0) print *,'Time transform to tridi :',time_evp_fwd if(myid == 0) print *,'Time transform to tridi :',time_evp_fwd
if(myid == 0) print *,'Time solve tridi :',time_evp_solve if(myid == 0) print *,'Time solve tridi :',time_evp_solve
if(myid == 0) print *,'Time transform back EVs :',time_evp_back if(myid == 0) print *,'Time transform back EVs :',time_evp_back
......
...@@ -87,12 +87,20 @@ program test_real_gen ...@@ -87,12 +87,20 @@ program test_real_gen
call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols ) call BLACS_Gridinit( my_blacs_ctxt, 'C', np_rows, np_cols )
call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol ) call BLACS_Gridinfo( my_blacs_ctxt, nprow, npcol, my_prow, my_pcol )
if (myid==0) then
print '(a)','| Past BLACS_Gridinfo.'
end if
! All ELPA routines need MPI communicators for communicating within ! All ELPA routines need MPI communicators for communicating within
! rows or columns of processes, these are set in get_elpa_row_col_comms. ! rows or columns of processes, these are set in get_elpa_row_col_comms.
call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, & call get_elpa_row_col_comms(mpi_comm_world, my_prow, my_pcol, &
mpi_comm_rows, mpi_comm_cols) mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| Past split communicator setup for rows and columns.'
end if
! Determine the necessary size of the distributed matrices, ! Determine the necessary size of the distributed matrices,
! we use the Scalapack tools routine NUMROC for that. ! we use the Scalapack tools routine NUMROC for that.
...@@ -106,6 +114,10 @@ program test_real_gen ...@@ -106,6 +114,10 @@ program test_real_gen
call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info ) call descinit( sc_desc, na, na, nblk, nblk, 0, 0, my_blacs_ctxt, na_rows, info )
if (myid==0) then
print '(a)','| Past scalapack descriptor setup.'
end if
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! Allocate matrices and set up test matrices for the eigenvalue problem ! Allocate matrices and set up test matrices for the eigenvalue problem
...@@ -132,13 +144,26 @@ program test_real_gen ...@@ -132,13 +144,26 @@ program test_real_gen
call RANDOM_NUMBER(z) call RANDOM_NUMBER(z)
a(:,:) = z(:,:) a(:,:) = z(:,:)
if (myid==0) then
print '(a)','| Random matrix block has been set up. (only processor 0 confirms this step)'
end if
call pdtran(na, na, 1.d0, z, 1, 1, sc_desc, 1.d0, a, 1, 1, sc_desc) ! A = A + Z**T call pdtran(na, na, 1.d0, z, 1, 1, sc_desc, 1.d0, a, 1, 1, sc_desc) ! A = A + Z**T
if (myid==0) then
print '(a)','| Random matrix has been symmetrized.'
end if
! The matrix B in the generalized eigenvalue problem must be symmetric ! The matrix B in the generalized eigenvalue problem must be symmetric
! and positive definite - we use a simple diagonally dominant matrix ! and positive definite - we use a simple diagonally dominant matrix
call pdlaset('Full', na, na, 1.d0/na, 1.1d0, b, 1, 1, sc_desc ) call pdlaset('Full', na, na, 1.d0/na, 1.1d0, b, 1, 1, sc_desc )
if (myid==0) then
print '(a)','| Simple diagonally dominant overlap matrix has been initialized.'
end if
! Save original matrices A and B for later accuracy checks ! Save original matrices A and B for later accuracy checks
as = a as = a
...@@ -155,8 +180,17 @@ program test_real_gen ...@@ -155,8 +180,17 @@ program test_real_gen
! PDPOTRF very often fails on higher processor numbers for unknown reasons! ! PDPOTRF very often fails on higher processor numbers for unknown reasons!
call cholesky_real(na, b, na_rows, nblk, mpi_comm_rows, mpi_comm_cols) call cholesky_real(na, b, na_rows, nblk, mpi_comm_rows, mpi_comm_cols)
if (myid==0) then
print '(a)','| Cholesky factorization complete.'
end if
call invert_trm_real(na, b, na_rows, nblk, mpi_comm_rows, mpi_comm_cols) call invert_trm_real(na, b, na_rows, nblk, mpi_comm_rows, mpi_comm_cols)