TABLE OF CONTENTS


ABINIT/m_slk [ Modules ]

[ Top ] [ Modules ]

NAME

 m_slk

FUNCTION

 This module contains the description of the variables used in the ScaLAPACK routines.

COPYRIGHT

 Copyright (C) 2004-2018 ABINIT group (CS,GZ,FB,MG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 module m_slk
24 
25  use defs_basis
26  use m_xmpi
27  use m_errors
28  use m_abicore
29 
30  use m_fstrings,      only : firstchar, toupper
31  use m_numeric_tools, only : print_arr
32 
33 #ifdef HAVE_LINALG_ELPA
34  use m_elpa
35 #endif
36 
37 #ifdef HAVE_MPI2
38  use mpi
39 #endif
40 
41  implicit none
42 
43 #ifdef HAVE_MPI1
44  include 'mpif.h'
45 #endif
46 
47  private
48 
49  ! parameters of the scaLAPACK array descriptor.
50  integer,private,parameter :: DLEN_ = 9    ! length
51  integer,private,parameter :: Dtype_ = 1   ! type
52  integer,private,parameter :: CTXT_ = 2    ! BLACS context
53  integer,private,parameter :: M_ = 3       ! nb global lines
54  integer,private,parameter :: N_ = 4       ! nb global columns
55  integer,private,parameter :: MB_ = 5      ! nb lines of a block
56  integer,private,parameter :: NB_ = 6      ! nb columns of a bloc
57  integer,private,parameter :: RSRC_ = 7    ! line of processors at the beginning
58  integer,private,parameter :: CSRC_ = 8    ! column of processors at the beginning
59  integer,private,parameter :: LLD_ = 9     ! local number of lines

m_slk/build_grid_scalapack [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  build_grid_scalapack

FUNCTION

  Set up the processor grid for ScaLAPACK as a function of the total number of processors attributed to the grid.

INPUTS

  nbprocs= total number of processors
  comm= MPI communicator

OUTPUT

  grid= the grid of processors used by Scalapack

PARENTS

      m_slk,m_xgScalapack

CHILDREN

SOURCE

245 subroutine build_grid_scalapack(grid,nbprocs,comm)
246 
247 
248 !This section has been created automatically by the script Abilint (TD).
249 !Do not modify the following lines by hand.
250 #undef ABI_FUNC
251 #define ABI_FUNC 'build_grid_scalapack'
252 !End of the abilint section
253 
254   implicit none
255 
256 !Arguments ------------------------------------
257  integer,intent(in) :: nbprocs,comm
258  type(grid_scalapack),intent(out) :: grid
259 !Local variables-------------------------------
260  integer  :: i
261 
262 ! *********************************************************************
263 
264  DBG_ENTER("COLL")
265 
266  grid%nbprocs=nbprocs
267 
268 !Search for a rectangular grid of processors
269  i=INT(SQRT(float(nbprocs)))
270  do while (MOD(nbprocs,i) /= 0)
271    i = i-1
272  end do
273  i=max(i,1)
274 
275  grid%dims(1) = i
276  grid%dims(2) = INT(nbprocs/i)
277 
278  grid%ictxt = comm
279 
280  call BLACS_GRIDINIT(grid%ictxt,'R',grid%dims(1),grid%dims(2))
281 
282  DBG_EXIT("COLL")
283 
284 end subroutine build_grid_scalapack

m_slk/build_processor_scalapack [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  build_processor_scalapack

FUNCTION

  Builds a processor descriptor for ScaLAPACK.
  Build of the data related to one processor in a grid

INPUTS

  grid= array representing the grid of processors
  myproc= selected processor
  comm= MPI communicator

OUTPUT

  processor= descriptor of a processor

PARENTS

      m_slk

CHILDREN

SOURCE

312 subroutine build_processor_scalapack(processor,grid,myproc,comm)
313 
314 
315 !This section has been created automatically by the script Abilint (TD).
316 !Do not modify the following lines by hand.
317 #undef ABI_FUNC
318 #define ABI_FUNC 'build_processor_scalapack'
319 !End of the abilint section
320 
321  implicit none
322 
323 !Arguments ------------------------------------
324  integer,intent(in) :: myproc,comm
325  type(processor_scalapack),intent(inout) :: processor
326  type(grid_scalapack),intent(in) :: grid
327 !Local variables-------------------------------
328 
329 ! *********************************************************************
330 
331  DBG_ENTER("COLL")
332 
333  processor%grid= grid
334 
335  processor%myproc = myproc
336 
337  processor%comm = comm
338 
339  call BLACS_GRIDINFO(grid%ictxt,processor%grid%dims(1), &
340 & processor%grid%dims(2),processor%coords(1), &
341 & processor%coords(2))
342 
343 !These values are the same as those computed by BLACS_GRIDINFO
344 !except in the case where the myproc argument is not the local proc
345  processor%coords(1) = INT((myproc) / grid%dims(2))
346  processor%coords(2) = MOD((myproc), grid%dims(2))
347 
348  DBG_EXIT("COLL")
349 
350 end subroutine build_processor_scalapack

m_slk/compute_eigen1 [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  compute_eigen1

FUNCTION

  Calculation of eigenvalues and eigenvectors.
  complex and real cases.

INPUTS

  comm= MPI communicator
  cplex=1 if matrix is real, 2 if complex
  nbli_global number of lines
  nbco_global number of columns
  matrix= the matrix to process
  vector= eigenvalues of the matrix
  istwf_k= option parameter that describes the storage of wfs

OUTPUT

  vector

SIDE EFFECTS

  results= ScaLAPACK matrix coming out of the operation
  eigen= eigenvalues of the matrix

PARENTS

      m_xgScalapack

CHILDREN

SOURCE

2689 subroutine compute_eigen1(comm,processor,cplex,nbli_global,nbco_global,matrix,vector,istwf_k)
2690 
2691 
2692 !This section has been created automatically by the script Abilint (TD).
2693 !Do not modify the following lines by hand.
2694 #undef ABI_FUNC
2695 #define ABI_FUNC 'compute_eigen1'
2696 !End of the abilint section
2697 
2698  implicit none
2699 
2700 !Arguments ------------------------------------
2701 !scalaras
2702  integer,intent(in) :: comm
2703  integer,intent(in) :: cplex,nbli_global,nbco_global
2704  integer,intent(in) :: istwf_k
2705  type(processor_scalapack),intent(in) :: processor
2706 !arrays
2707  real(dp),intent(inout) :: matrix(cplex*nbli_global,nbco_global)
2708  real(dp),intent(inout) :: vector(:)
2709 
2710 !Local variables-------------------------------
2711  integer :: i,j,ierr
2712  type(matrix_scalapack) :: sca_matrix1
2713  type(matrix_scalapack) :: sca_matrix2
2714  real(dp),allocatable :: r_tmp_evec(:,:)
2715  complex(dpc),allocatable :: z_tmp_evec(:,:)
2716 
2717 ! *************************************************************************
2718 
2719  ! ================================
2720  ! INITIALISATION SCALAPACK MATRIX
2721  ! ================================
2722  call init_matrix_scalapack(sca_matrix1,nbli_global,nbco_global,processor,istwf_k,10)
2723  call init_matrix_scalapack(sca_matrix2,nbli_global,nbco_global,processor,istwf_k,10)
2724 
2725  ! ==============================
2726  ! FILLING SCALAPACK MATRIX
2727  ! ==============================
2728  if ( istwf_k /= 2 ) then
2729    ABI_CHECK(cplex==2, "cplex != 2")
2730    ABI_MALLOC(z_tmp_evec,(nbli_global,nbco_global))
2731    z_tmp_evec=cmplx(0._DP,0._DP)
2732 #ifdef HAVE_LINALG_ELPA
2733    do j=1,nbco_global
2734       do i=j+1,nbli_global
2735          matrix(2*(i-1)+1,j) = matrix(2*(j-1)+1,i)
2736          matrix(2*(i-1)+2,j) = -matrix(2*(j-1)+2,i)
2737       end do
2738    end do
2739 #endif
2740    call matrix_from_complexmatrix(sca_matrix1,matrix,istwf_k)
2741  else
2742    ABI_CHECK(cplex==1, "cplex != 2")
2743    ABI_MALLOC(r_tmp_evec,(nbli_global,nbco_global))
2744    r_tmp_evec(:,:)=0._DP
2745 #ifdef HAVE_LINALG_ELPA
2746    do j=1,nbco_global
2747       do i=j+1,nbli_global
2748          matrix(i,j) = matrix(j,i)
2749       end do
2750    end do
2751 #endif
2752    call matrix_from_realmatrix(sca_matrix1,matrix,istwf_k)
2753  endif
2754 
2755  ! ================================
2756  ! COMPUTE EIGEN VALUES AND VECTORS : A * X = lambda  * X
2757  ! ================================
2758  call compute_eigen_problem(processor,sca_matrix1,&
2759 &                           sca_matrix2,vector,&
2760 &                           comm,istwf_k)
2761 
2762  ! ==============================
2763  ! CONCATENATE EIGEN VECTORS
2764  ! ==============================
2765  if ( istwf_k /= 2 ) then
2766    call matrix_to_complexmatrix(sca_matrix2,z_tmp_evec,istwf_k)
2767    call MPI_ALLREDUCE(z_tmp_evec, matrix, nbli_global*nbco_global, MPI_DOUBLE_complex,&
2768 &    MPI_SUM,comm,ierr)
2769  else
2770    call matrix_to_realmatrix(sca_matrix2,r_tmp_evec,istwf_k)
2771    call MPI_ALLREDUCE(r_tmp_evec, matrix, nbli_global*nbco_global, MPI_DOUBLE_PRECISION,&
2772 &    MPI_SUM,comm,ierr)
2773  endif
2774 
2775  ! ====================================
2776  ! DESTRUCTION SCALAPACK AND TMP MATRICES
2777  ! ====================================
2778  call destruction_matrix_scalapack(sca_matrix1)
2779  call destruction_matrix_scalapack(sca_matrix2)
2780 
2781  if ( istwf_k /= 2 ) then
2782    ABI_FREE(z_tmp_evec)
2783  else
2784    ABI_FREE(r_tmp_evec)
2785  endif
2786 
2787 end subroutine compute_eigen1

m_slk/compute_eigen2 [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  compute_eigen2

FUNCTION

  Calculation of eigenvalues and eigenvectors:
  A * X = lambda * B * X
  complex and real cases.

INPUTS

  comm= MPI communicator
  cplex=1 if matrix is real, 2 if complex
  nbli_global number of lines
  nbco_global number of columns
  matrix1= first ScaLAPACK matrix (matrix A)
  matrix2= second ScaLAPACK matrix (matrix B)
  vector=
  istwf_k= option parameter that describes the storage of wfs

OUTPUT

  None

SIDE EFFECTS

  results= ScaLAPACK matrix coming out of the operation
  eigen= eigenvalues of the matrix

PARENTS

      m_xgScalapack

CHILDREN

SOURCE

2825 subroutine compute_eigen2(comm,processor,cplex,nbli_global,nbco_global,matrix1,matrix2,vector,istwf_k)
2826 
2827 
2828 !This section has been created automatically by the script Abilint (TD).
2829 !Do not modify the following lines by hand.
2830 #undef ABI_FUNC
2831 #define ABI_FUNC 'compute_eigen2'
2832 !End of the abilint section
2833 
2834  implicit none
2835 
2836 !Arguments ------------------------------------
2837 !scalars
2838  integer,intent(in) :: cplex,nbli_global,nbco_global
2839  integer,intent(in) :: comm
2840  integer,intent(in) :: istwf_k
2841  type(processor_scalapack),intent(in) :: processor
2842 !arrays
2843  real(dp),intent(inout) :: matrix1(cplex*nbli_global,nbco_global)
2844  real(dp),intent(inout) :: matrix2(cplex*nbli_global,nbco_global)
2845  real(dp),intent(inout) :: vector(:)
2846 
2847 !Local variables-------------------------------
2848  integer :: i,j,ierr
2849  type(matrix_scalapack) :: sca_matrix1
2850  type(matrix_scalapack) :: sca_matrix2
2851  type(matrix_scalapack) :: sca_matrix3
2852  real(dp),allocatable :: r_tmp_evec(:,:)
2853  complex(dpc),allocatable :: z_tmp_evec(:,:)
2854 
2855 ! *************************************************************************
2856 
2857  ! ================================
2858  ! INITIALISATION SCALAPACK MATRIX
2859  ! ================================
2860  call init_matrix_scalapack(sca_matrix1,nbli_global,nbco_global,processor,istwf_k,10)
2861  call init_matrix_scalapack(sca_matrix2,nbli_global,nbco_global,processor,istwf_k,10)
2862  call init_matrix_scalapack(sca_matrix3,nbli_global,nbco_global,processor,istwf_k,10)
2863 
2864  ! ==============================
2865  ! FILLING SCALAPACK MATRIX
2866  ! ==============================
2867  if ( istwf_k /= 2 ) then
2868    ABI_CHECK(cplex==2, "cplex != 2")
2869    ABI_MALLOC(z_tmp_evec,(nbli_global,nbco_global))
2870    z_tmp_evec=cmplx(0._DP,0._DP)
2871 #ifdef HAVE_LINALG_ELPA
2872    do j=1,nbco_global
2873       do i=j+1,nbli_global
2874          matrix1(2*(i-1)+1,j) = matrix1(2*(j-1)+1,i)
2875          matrix1(2*(i-1)+2,j) = -matrix1(2*(j-1)+2,i)
2876          matrix2(2*(i-1)+1,j) = matrix2(2*(j-1)+1,i)
2877          matrix2(2*(i-1)+2,j) = -matrix2(2*(j-1)+2,i)
2878       end do
2879    end do
2880 #endif
2881    call matrix_from_complexmatrix(sca_matrix1,matrix1,istwf_k)
2882    call matrix_from_complexmatrix(sca_matrix2,matrix2,istwf_k)
2883  else
2884    ABI_CHECK(cplex==1, "cplex != 2")
2885    ABI_MALLOC(r_tmp_evec,(nbli_global,nbco_global))
2886    r_tmp_evec(:,:)=0._DP
2887 #ifdef HAVE_LINALG_ELPA
2888    do j=1,nbco_global
2889       do i=j+1,nbli_global
2890          matrix1(i,j) = matrix1(j,i)
2891          matrix2(i,j) = matrix2(j,i)
2892       end do
2893    end do
2894 #endif
2895    call matrix_from_realmatrix(sca_matrix1,matrix1,istwf_k)
2896    call matrix_from_realmatrix(sca_matrix2,matrix2,istwf_k)
2897  endif
2898 
2899  ! ================================
2900  ! COMPUTE EIGEN VALUES AND VECTORS : A * X = lambda * B * X
2901  ! ===============================
2902  call compute_generalized_eigen_problem(processor,sca_matrix1,sca_matrix2,&
2903 &                     sca_matrix3,vector,&
2904 &                     comm,istwf_k)
2905 
2906  ! ==============================
2907  ! CONCATENATE EIGEN VECTORS
2908  ! ==============================
2909  if ( istwf_k /= 2 ) then
2910    call matrix_to_complexmatrix(sca_matrix3,z_tmp_evec,istwf_k)
2911    call MPI_ALLREDUCE(z_tmp_evec, matrix1, nbli_global*nbco_global, MPI_DOUBLE_complex,&
2912 &    MPI_SUM,comm,ierr)
2913  else
2914    call matrix_to_realmatrix(sca_matrix3,r_tmp_evec,istwf_k)
2915    call MPI_ALLREDUCE(r_tmp_evec, matrix1, nbli_global*nbco_global, MPI_DOUBLE_PRECISION,&
2916 &    MPI_SUM,comm,ierr)
2917  endif
2918 
2919 
2920  ! ====================================
2921  ! DESTRUCTION SCALAPACK AND TMP MATRICES
2922  ! ====================================
2923  call destruction_matrix_scalapack(sca_matrix1)
2924  call destruction_matrix_scalapack(sca_matrix2)
2925  call destruction_matrix_scalapack(sca_matrix3)
2926 
2927 end subroutine compute_eigen2

m_slk/compute_eigen_problem [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  compute_eigen_problem

FUNCTION

  Calculation of eigenvalues and eigenvectors.
  A * X = lambda * X
  complex and real cases.

INPUTS

  processor= descriptor of a processor
  matrix= the matrix to process
  comm= MPI communicator
  istwf_k= option parameter that describes the storage of wfs

OUTPUT

  None

SIDE EFFECTS

  results= ScaLAPACK matrix coming out of the operation
  eigen= eigenvalues of the matrix

PARENTS

      m_slk

CHILDREN

SOURCE

2100 subroutine compute_eigen_problem(processor,matrix,results,eigen,comm,istwf_k)
2101 
2102 
2103 !This section has been created automatically by the script Abilint (TD).
2104 !Do not modify the following lines by hand.
2105 #undef ABI_FUNC
2106 #define ABI_FUNC 'compute_eigen_problem'
2107 !End of the abilint section
2108 
2109   implicit none
2110 
2111 #ifdef HAVE_LINALG_ELPA
2112   !Arguments ------------------------------------
2113   type(processor_scalapack),intent(in) :: processor
2114   type(matrix_scalapack),intent(inout) :: matrix
2115   type(matrix_scalapack),intent(inout) :: results
2116   DOUBLE PRECISION,intent(inout) :: eigen(:)
2117   integer,intent(in)  :: comm,istwf_k
2118   !Local variables ------------------------------
2119   type(elpa_hdl_t) :: elpa_hdl
2120 
2121 !************************************************************************
2122 
2123   call elpa_func_allocate(elpa_hdl,processor%comm,processor%coords(1),processor%coords(2))
2124   call elpa_func_set_matrix(elpa_hdl,matrix%sizeb_global(1),matrix%sizeb_blocs(1),&
2125 &                           matrix%sizeb_local(1),matrix%sizeb_local(2))
2126 
2127   if (istwf_k/=2) then
2128     call elpa_func_solve_evp_1stage(elpa_hdl,matrix%buffer_cplx,results%buffer_cplx,&
2129 &                                   eigen,matrix%sizeb_global(1))
2130   else
2131     call elpa_func_solve_evp_1stage(elpa_hdl,matrix%buffer_real,results%buffer_real,&
2132 &                                   eigen,matrix%sizeb_global(1))
2133   end if
2134 
2135   call elpa_func_deallocate(elpa_hdl)
2136 
2137 #else
2138   !Arguments ------------------------------------
2139   type(processor_scalapack),intent(in)       :: processor
2140   type(matrix_scalapack),intent(in)          :: matrix
2141   type(matrix_scalapack),intent(inout)       :: results
2142   DOUBLE PRECISION,intent(inout) :: eigen(:)
2143   integer,intent(in)  :: comm,istwf_k
2144   !Local variables-------------------------------
2145   integer            :: LRWORK,LIWORK,LCWORK,INFO
2146   character(len=500) :: msg
2147 
2148   integer         , dimension(1) :: IWORK_tmp
2149   DOUBLE PRECISION, dimension(1) :: RWORK_tmp
2150   complex(dpc)     , dimension(1) :: CWORK_tmp
2151 
2152   integer         , allocatable  :: IWORK(:)
2153   DOUBLE PRECISION, allocatable  :: RWORK(:)
2154   complex(dpc)     , allocatable  :: CWORK(:)
2155 
2156   integer,          allocatable :: ICLUSTR(:)
2157   integer,          allocatable :: IFAIL(:)
2158   DOUBLE PRECISION, allocatable :: GAP(:)
2159 
2160   DOUBLE PRECISION            :: ABSTOL,ORFAC
2161   integer,          parameter :: IZERO=0
2162 
2163   integer ::  M,NZ,IA,JA,IZ,JZ,ierr,TWORK_tmp(3),TWORK(3)
2164 
2165   DOUBLE PRECISION, external :: PDLAMCH
2166 
2167 ! *************************************************************************
2168 
2169   ! Initialisation
2170   INFO   = 0
2171   ABSTOL = zero
2172   ORFAC  = -1.D+0
2173 
2174   ! Allocation of the variables for the results of the calculations
2175   ABI_MALLOC(IFAIL,(matrix%sizeb_global(2)))
2176   ABI_MALLOC(ICLUSTR,(2*processor%grid%dims(1)*processor%grid%dims(2)))
2177   ABI_MALLOC(GAP,(processor%grid%dims(1)*processor%grid%dims(2)))
2178 
2179   ! Get the size of the work arrays
2180   if (istwf_k/=2) then
2181      call PZHEEVX('V','A','U',&
2182 &      matrix%sizeb_global(2),&
2183 &      matrix%buffer_cplx,1,1,matrix%descript%tab, &
2184 &      ZERO,ZERO,IZERO,IZERO,ABSTOL,&
2185 &      m,nz,eigen,ORFAC, &
2186 &      results%buffer_cplx,1,1,results%descript%tab, &
2187 &      CWORK_tmp,-1,RWORK_tmp,-1,IWORK_tmp,-1,&
2188 &      IFAIL,ICLUSTR,GAP,INFO)
2189   else
2190      call PDSYEVX('V','A','U',&
2191 &      matrix%sizeb_global(2),&
2192 &      matrix%buffer_real,1,1,matrix%descript%tab, &
2193 &      ZERO,ZERO,IZERO,IZERO,ABSTOL,&
2194 &      m,nz,eigen,ORFAC, &
2195 &      results%buffer_real,1,1,results%descript%tab, &
2196 &      RWORK_tmp,-1,IWORK_tmp,-1,&
2197 &      IFAIL,ICLUSTR,GAP,INFO)
2198   end if
2199 
2200   if (INFO/=0) then
2201      write(msg,'(A,I6)') "Problem to compute workspace to use ScaLAPACK, INFO=",INFO
2202      MSG_ERROR(msg)
2203   endif
2204 
2205   TWORK_tmp(1) = IWORK_tmp(1)
2206   TWORK_tmp(2) = INT(RWORK_tmp(1))
2207   TWORK_tmp(3) = INT(real(CWORK_tmp(1)))
2208 
2209  !! Get the maximum of the size of the work arrays processor%comm
2210   call MPI_ALLREDUCE(TWORK_tmp,TWORK,3,MPI_integer,MPI_MAX,comm,ierr)
2211 
2212   LIWORK = TWORK(1)
2213   LRWORK = TWORK(2) + matrix%sizeb_global(2) *(matrix%sizeb_global(2)-1)
2214   LCWORK = TWORK(3)
2215 
2216   ! Allocation of the work arrays
2217   if (LIWORK>0) then
2218     ABI_MALLOC(IWORK,(LIWORK))
2219     IWORK(:) = 0
2220   else
2221     ABI_MALLOC(IWORK,(1))
2222   end if
2223   if (LRWORK>0) then
2224     ABI_MALLOC(RWORK,(LRWORK))
2225     RWORK(:) = 0._dp
2226   else
2227     ABI_MALLOC(RWORK,(1))
2228   end if
2229   if (LCWORK>0) then
2230     ABI_MALLOC(CWORK,(LCWORK))
2231     CWORK(:) = (0._dp,0._dp)
2232   else
2233     ABI_MALLOC(CWORK,(1))
2234   end if
2235 
2236   ! Call the calculation routine
2237   if (istwf_k/=2) then
2238   !   write(std_out,*) 'I am using PZHEEVX'
2239      call PZHEEVX('V','A','U',&
2240 &      matrix%sizeb_global(2),&
2241 &      matrix%buffer_cplx,1,1,matrix%descript%tab, &
2242 &      ZERO,ZERO,IZERO,IZERO,ABSTOL,&
2243 &      m,nz,eigen,ORFAC, &
2244 &      results%buffer_cplx,1,1,results%descript%tab, &
2245 &      CWORK,LCWORK,RWORK,LRWORK,IWORK,LIWORK,&
2246 &      IFAIL,ICLUSTR,GAP,INFO)
2247   else
2248   !   write(std_out,*) ' I am using PDSYEVX'
2249      call PDSYEVX('V','A','U',&
2250 &      matrix%sizeb_global(2),&
2251 &      matrix%buffer_real,1,1,matrix%descript%tab, &
2252 &      ZERO,ZERO,IZERO,IZERO,ABSTOL,&
2253 &      m,nz,eigen,ORFAC, &
2254 &      results%buffer_real,1,1,results%descript%tab, &
2255 &      RWORK,LRWORK,IWORK,LIWORK,&
2256 &      IFAIL,ICLUSTR,GAP,INFO)
2257   endif
2258 
2259   if (INFO/=0) then
2260      write(msg,'(A,I6)') "Problem to compute eigenvalues and eigenvectors with ScaLAPACK, INFO=",INFO
2261      MSG_ERROR(msg)
2262   endif
2263 
2264   ABI_FREE(IFAIl)
2265   ABI_FREE(ICLUSTR)
2266   ABI_FREE(GAP)
2267   if (allocated(IWORK))  then
2268     ABI_FREE(IWORK)
2269   end if
2270   if (allocated(RWORK))  then
2271     ABI_FREE(RWORK)
2272   end if
2273   if (allocated(CWORK))  then
2274     ABI_FREE(CWORK)
2275   end if
2276 #endif
2277   return
2278 
2279 end subroutine compute_eigen_problem

m_slk/compute_generalized_eigen_problem [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  compute_generalized_eigen_problem

FUNCTION

  Calculation of eigenvalues and eigenvectors:
  A * X = lambda * B * X
  complex and real cases.

INPUTS

  processor= descriptor of a processor
  matrix1= first ScaLAPACK matrix (matrix A)
  matrix2= second ScaLAPACK matrix (matrix B)
  comm= MPI communicator
  istwf_k= option parameter that describes the storage of wfs

OUTPUT

  None

SIDE EFFECTS

  results= ScaLAPACK matrix coming out of the operation
  eigen= eigenvalues of the matrix

PARENTS

      m_slk,rayleigh_ritz

CHILDREN

SOURCE

2315 #ifdef HAVE_LINALG_ELPA
2316 
2317 subroutine solve_gevp_complex(na,nev,na_rows,na_cols,nblk,a,b,ev,z,tmp1,tmp2, &
2318                               my_prow,my_pcol,np_rows,np_cols,sc_desc,comm)
2319 
2320 
2321 !This section has been created automatically by the script Abilint (TD).
2322 !Do not modify the following lines by hand.
2323 #undef ABI_FUNC
2324 #define ABI_FUNC 'solve_gevp_complex'
2325 !End of the abilint section
2326 
2327   implicit none
2328 
2329   !-Arguments
2330   integer,intent(in) :: na
2331   integer,intent(in) :: nev
2332   integer,intent(in) :: na_rows,na_cols
2333   integer,intent(in) :: nblk
2334   integer,intent(in) :: my_pcol,my_prow
2335   integer,intent(in) :: np_cols,np_rows
2336   integer,intent(in) :: sc_desc(9)
2337   integer,intent(in) :: comm
2338   real*8 :: ev(na)
2339   complex*16 :: a(na_rows,na_cols),b(na_rows,na_cols),z(na_rows,na_cols)
2340   complex*16 :: tmp1(na_rows,na_cols),tmp2(na_rows,na_cols)
2341   !-Local variables
2342   integer :: i, n_col, n_row
2343   integer,external :: indxl2g,numroc
2344   complex*16, parameter :: CZERO = (0.d0,0.d0), CONE = (1.d0,0.d0)
2345   type(elpa_hdl_t) :: elpa_hdl
2346 
2347 ! *************************************************************************
2348 
2349   ! 0. Allocate ELPA handle
2350   call elpa_func_allocate(elpa_hdl,comm,my_prow,my_pcol)
2351   call elpa_func_set_matrix(elpa_hdl,na,nblk,na_rows,na_cols)
2352 
2353   ! 1. Calculate Cholesky factorization of Matrix B = U**T * U
2354   !    and invert triangular matrix U
2355   call elpa_func_cholesky(elpa_hdl,b)
2356   call elpa_func_invert_triangular(elpa_hdl,b)
2357   ! 2. Calculate U**-T * A * U**-1
2358   ! 2a. tmp1 = U**-T * A
2359   call elpa_func_hermitian_multiply(elpa_hdl,'U','L',na,b,a,na_rows,na_cols,tmp1,na_rows,na_cols)
2360   ! 2b. tmp2 = tmp1**T
2361   call pztranc(na,na,CONE,tmp1,1,1,sc_desc,CZERO,tmp2,1,1,sc_desc)
2362   ! 2c. A =  U**-T * tmp2 ( = U**-T * Aorig * U**-1 )
2363   call elpa_func_hermitian_multiply(elpa_hdl,'U','U',na,b,tmp2,na_rows,na_cols,a,na_rows,na_cols)
2364   ! A is only set in the upper half, solve_evp_real needs a full matrix
2365   ! Set lower half from upper half
2366   call pztranc(na,na,CONE,a,1,1,sc_desc,CZERO,tmp1,1,1,sc_desc)
2367   do i=1,na_cols
2368      ! Get global column corresponding to i and number of local rows up to
2369      ! and including the diagonal, these are unchanged in A
2370      n_col = indxl2g(i,     nblk, my_pcol, 0, np_cols)
2371      n_row = numroc (n_col, nblk, my_prow, 0, np_rows)
2372      a(n_row+1:na_rows,i) = tmp1(n_row+1:na_rows,i)
2373   enddo
2374   ! 3. Calculate eigenvalues/eigenvectors of U**-T * A * U**-1
2375   !    Eigenvectors go to tmp1
2376   call elpa_func_solve_evp_1stage(elpa_hdl,a,tmp1,ev,nev)
2377   ! 4. Backtransform eigenvectors: Z = U**-1 * tmp1
2378   ! hermitian_multiply needs the transpose of U**-1, thus tmp2 = (U**-1)**T
2379   call pztranc(na,na,CONE,b,1,1,sc_desc,CZERO,tmp2,1,1,sc_desc)
2380   call elpa_func_hermitian_multiply(elpa_hdl,'L','N',nev,tmp2,tmp1,na_rows,na_cols,z,na_rows,na_cols)
2381 
2382   call elpa_func_deallocate(elpa_hdl)
2383 
2384 end subroutine solve_gevp_complex
2385 
2386 !----------------------------------------------------------------------
2387 
2388 subroutine solve_gevp_real(na,nev,na_rows,na_cols,nblk,a,b,ev,z,tmp1,tmp2, &
2389                            my_prow,my_pcol,np_rows,np_cols,sc_desc,comm)
2390 
2391 
2392 !This section has been created automatically by the script Abilint (TD).
2393 !Do not modify the following lines by hand.
2394 #undef ABI_FUNC
2395 #define ABI_FUNC 'solve_gevp_real'
2396 !End of the abilint section
2397 
2398   implicit none
2399 
2400   !-Arguments
2401   integer,intent(in) :: na
2402   integer,intent(in) :: nev
2403   integer,intent(in) :: na_rows,na_cols
2404   integer,intent(in) :: nblk
2405   integer,intent(in) :: my_pcol,my_prow
2406   integer,intent(in) :: np_cols,np_rows
2407   integer,intent(in) :: sc_desc(9)
2408   integer,intent(in) :: comm
2409   real*8 :: ev(na)
2410   real*8 :: a(na_rows,na_cols),b(na_rows,na_cols),z(na_rows,na_cols)
2411   real*8::tmp1(na_rows,na_cols),tmp2(na_rows,na_cols)
2412   !-Local variables
2413   integer :: i, n_col, n_row
2414   integer,external :: indxl2g,numroc
2415   type(elpa_hdl_t) :: elpa_hdl
2416 
2417 ! *************************************************************************
2418 
2419   ! 0. Allocate ELPA handle
2420   call elpa_func_allocate(elpa_hdl,comm,my_prow,my_pcol)
2421   call elpa_func_set_matrix(elpa_hdl,na,nblk,na_rows,na_cols)
2422 
2423   ! 1. Calculate Cholesky factorization of Matrix B = U**T * U
2424   !    and invert triangular matrix U
2425   call elpa_func_cholesky(elpa_hdl,b)
2426   call elpa_func_invert_triangular(elpa_hdl,b)
2427   ! 2. Calculate U**-T * A * U**-1
2428   ! 2a. tmp1 = U**-T * A
2429   call elpa_func_hermitian_multiply(elpa_hdl,'U','L',na,b,a,na_rows,na_cols,tmp1,na_rows,na_cols)
2430   ! 2b. tmp2 = tmp1**T
2431   call pdtran(na,na,1.d0,tmp1,1,1,sc_desc,0.d0,tmp2,1,1,sc_desc)
2432   ! 2c. A =  U**-T * tmp2 ( = U**-T * Aorig * U**-1 )
2433   call elpa_func_hermitian_multiply(elpa_hdl,'U','U',na,b,tmp2,na_rows,na_cols,a,na_rows,na_cols)
2434   ! A is only set in the upper half, solve_evp_real needs a full matrix
2435   ! Set lower half from upper half
2436   call pdtran(na,na,1.d0,a,1,1,sc_desc,0.d0,tmp1,1,1,sc_desc)
2437   do i=1,na_cols
2438      ! Get global column corresponding to i and number of local rows up to
2439      ! and including the diagonal, these are unchanged in A
2440      n_col = indxl2g(i,     nblk, my_pcol, 0, np_cols)
2441      n_row = numroc (n_col, nblk, my_prow, 0, np_rows)
2442      a(n_row+1:na_rows,i) = tmp1(n_row+1:na_rows,i)
2443   enddo
2444   ! 3. Calculate eigenvalues/eigenvectors of U**-T * A * U**-1
2445   !    Eigenvectors go to tmp1
2446   call elpa_func_solve_evp_1stage(elpa_hdl,a,tmp1,ev,nev)
2447   ! 4. Backtransform eigenvectors: Z = U**-1 * tmp1
2448   !    hermitian_multiply needs the transpose of U**-1, thus tmp2 = (U**-1)**T
2449   call pdtran(na,na,1.d0,b,1,1,sc_desc,0.d0,tmp2,1,1,sc_desc)
2450   call elpa_func_hermitian_multiply(elpa_hdl,'L','N',nev,tmp2,tmp1,na_rows,na_cols,z,na_rows,na_cols)
2451 
2452   call elpa_func_deallocate(elpa_hdl)
2453 
2454  end subroutine solve_gevp_real
2455 #endif
2456 
2457 subroutine compute_generalized_eigen_problem(processor,matrix1,matrix2,results,eigen,comm,istwf_k)
2458 
2459 
2460 !This section has been created automatically by the script Abilint (TD).
2461 !Do not modify the following lines by hand.
2462 #undef ABI_FUNC
2463 #define ABI_FUNC 'compute_generalized_eigen_problem'
2464 !End of the abilint section
2465 
2466   implicit none
2467 
2468 #ifdef HAVE_LINALG_ELPA
2469 !Arguments ------------------------------------
2470   type(processor_scalapack),intent(in)       :: processor
2471   type(matrix_scalapack),intent(in)          :: matrix1,matrix2
2472   type(matrix_scalapack),intent(inout)       :: results
2473   DOUBLE PRECISION,intent(inout) :: eigen(:)
2474 
2475   integer,intent(in)  :: comm,istwf_k
2476 
2477 !Local
2478   type(matrix_scalapack)          :: tmp1,tmp2
2479   integer :: i,n_col, n_row
2480   integer,external :: indxl2g,numroc
2481 
2482   call  init_matrix_scalapack(tmp1,matrix1%sizeb_global(1),matrix1%sizeb_global(2),processor,istwf_k)
2483   call  init_matrix_scalapack(tmp2,matrix1%sizeb_global(1),matrix1%sizeb_global(2),processor,istwf_k)
2484   if (istwf_k/=2) then
2485      call solve_gevp_complex(matrix1%sizeb_global(1),matrix1%sizeb_global(2), &
2486 &          matrix1%sizeb_local(1),matrix1%sizeb_local(2),matrix1%sizeb_blocs(1), &
2487 &          matrix1%buffer_cplx,matrix2%buffer_cplx,eigen,results%buffer_cplx, &
2488 &          tmp1%buffer_cplx,tmp2%buffer_cplx, &
2489 &          processor%coords(1),processor%coords(2), &
2490 &          processor%grid%dims(1),processor%grid%dims(2), &
2491 &          matrix1%descript%tab,processor%comm)
2492   else
2493      call solve_gevp_real(matrix1%sizeb_global(1),matrix1%sizeb_global(2), &
2494 &          matrix1%sizeb_local(1),matrix1%sizeb_local(2),matrix1%sizeb_blocs(1), &
2495 &          matrix1%buffer_real,matrix2%buffer_real,eigen,results%buffer_real, &
2496 &          tmp1%buffer_real,tmp2%buffer_real, &
2497 &          processor%coords(1),processor%coords(2), &
2498 &          processor%grid%dims(1),processor%grid%dims(2), &
2499 &          matrix1%descript%tab,processor%comm)
2500   end if
2501   call destruction_matrix_scalapack(tmp1)
2502   call destruction_matrix_scalapack(tmp2)
2503 
2504 #else
2505 !Arguments ------------------------------------
2506   type(processor_scalapack),intent(in)       :: processor
2507   type(matrix_scalapack),intent(in)          :: matrix1,matrix2
2508   type(matrix_scalapack),intent(inout)       :: results
2509   DOUBLE PRECISION,intent(inout) :: eigen(:)
2510 
2511   integer,intent(in)  :: comm,istwf_k
2512 
2513 !Local variables-------------------------------
2514   integer            :: LRWORK,LIWORK,LCWORK,INFO
2515   character(len=500) :: msg
2516 
2517   integer         , dimension(1) :: IWORK_tmp
2518   DOUBLE PRECISION, dimension(1) :: RWORK_tmp
2519   complex(dpc)     , dimension(1) :: CWORK_tmp
2520 
2521   integer         , allocatable  :: IWORK(:)
2522   DOUBLE PRECISION, allocatable  :: RWORK(:)
2523   complex(dpc)     , allocatable  :: CWORK(:)
2524 
2525 
2526   integer,          allocatable :: ICLUSTR(:)
2527   integer,          allocatable :: IFAIL(:)
2528   DOUBLE PRECISION, allocatable :: GAP(:)
2529 
2530   DOUBLE PRECISION            :: ABSTOL,ORFAC
2531   integer         , parameter :: IZERO=0
2532 
2533   integer ::  M,NZ,IA,JA,IZ,JZ,ierr,TWORK_tmp(3),TWORK(3)
2534 
2535   DOUBLE PRECISION, external :: PDLAMCH
2536 
2537 ! *************************************************************************
2538 
2539   ! Initialisation
2540   INFO   = 0
2541   ABSTOL = zero
2542   ORFAC  = -1.D+0
2543 
2544   ! Allocate the arrays for the results of the calculation
2545   ABI_MALLOC(IFAIL  ,(matrix1%sizeb_global(2)))
2546   ABI_MALLOC(ICLUSTR,(2*processor%grid%dims(1)*processor%grid%dims(2)))
2547   ABI_MALLOC(GAP    ,(  processor%grid%dims(1)*processor%grid%dims(2)))
2548 
2549   ! Get the size of the work arrays
2550   if (istwf_k/=2) then
2551      call PZHEGVX(1,'V','A','U',&
2552 &      matrix1%sizeb_global(2),&
2553 &      matrix1%buffer_cplx,1,1,matrix1%descript%tab, &
2554 &      matrix2%buffer_cplx,1,1,matrix2%descript%tab, &
2555 &      ZERO,ZERO,IZERO,IZERO,ABSTOL,&
2556 &      m,nz,eigen,ORFAC, &
2557 &      results%buffer_cplx,1,1,results%descript%tab, &
2558 &      CWORK_tmp,-1,RWORK_tmp,-1,IWORK_tmp,-1,&
2559 &      IFAIL,ICLUSTR,GAP,INFO)
2560   else
2561      call PDSYGVX(1,'V','A','U',&
2562 &      matrix1%sizeb_global(2),&
2563 &      matrix1%buffer_real,1,1,matrix1%descript%tab, &
2564 &      matrix2%buffer_real,1,1,matrix2%descript%tab, &
2565 &      ZERO,ZERO,IZERO,IZERO,ABSTOL,&
2566 &      m,nz,eigen,ORFAC, &
2567 &      results%buffer_real,1,1,results%descript%tab, &
2568 &      RWORK_tmp,-1,IWORK_tmp,-1,&
2569 &      IFAIL,ICLUSTR,GAP,INFO)
2570   endif
2571 
2572   if (INFO/=0) then
2573      write(msg,'(A,I6)') "Problem to compute workspace to use ScaLAPACK, INFO=",INFO
2574      MSG_ERROR(msg)
2575   endif
2576 
2577   TWORK_tmp(1) = IWORK_tmp(1)
2578   TWORK_tmp(2) = INT(RWORK_tmp(1)) + matrix1%sizeb_global(2) *(matrix1%sizeb_global(2)-1)
2579   TWORK_tmp(3) = INT(real(CWORK_tmp(1)))
2580 
2581  ! Get the maximum of sizes of the work arrays processor%comm
2582   call MPI_ALLREDUCE(TWORK_tmp,TWORK,3,MPI_integer,MPI_MAX,comm,ierr)
2583 
2584   LIWORK = TWORK(1)
2585   LRWORK = TWORK(2)
2586   LCWORK = TWORK(3)
2587 
2588  ! Allocate the work arrays
2589   if (LIWORK>0) then
2590     ABI_MALLOC(IWORK,(LIWORK))
2591     IWORK(:) = 0
2592   else
2593     ABI_MALLOC(IWORK,(1))
2594   end if
2595   if (LRWORK>0) then
2596     ABI_MALLOC(RWORK,(LRWORK))
2597     RWORK(:) = 0._dp
2598   else
2599     ABI_MALLOC(RWORK,(1))
2600   end if
2601   if (LCWORK>0) then
2602     ABI_MALLOC(CWORK,(LCWORK))
2603     CWORK(:) = (0._dp,0._dp)
2604   else
2605     ABI_MALLOC(CWORK,(1))
2606   end if
2607 
2608   ! Call the calculation routine
2609   if (istwf_k/=2) then
2610   !   write(std_out,*) 'I am using PZHEGVX'
2611      call PZHEGVX(1,'V','A','U',&
2612 &      matrix1%sizeb_global(2),&
2613 &      matrix1%buffer_cplx,1,1,matrix1%descript%tab, &
2614 &      matrix2%buffer_cplx,1,1,matrix2%descript%tab, &
2615 &      ZERO,ZERO,IZERO,IZERO,ABSTOL,&
2616 &      m,nz,eigen,ORFAC, &
2617 &      results%buffer_cplx,1,1,results%descript%tab, &
2618 &      CWORK,LCWORK,RWORK,LRWORK,IWORK,LIWORK,&
2619 &      IFAIL,ICLUSTR,GAP,INFO)
2620   else
2621   !   write(std_out,*) 'I am using PDSYGVX'
2622      call PDSYGVX(1,'V','A','U',&
2623 &      matrix1%sizeb_global(2),&
2624 &      matrix1%buffer_real,1,1,matrix1%descript%tab, &
2625 &      matrix2%buffer_real,1,1,matrix2%descript%tab, &
2626 &      ZERO,ZERO,IZERO,IZERO,ABSTOL,&
2627 &      m,nz,eigen,ORFAC, &
2628 &      results%buffer_real,1,1,results%descript%tab, &
2629 &      RWORK,LRWORK,IWORK,LIWORK,&
2630 &      IFAIL,ICLUSTR,GAP,INFO)
2631   endif
2632 
2633   if (INFO/=0) then
2634      write(msg,'(A,I6)') "Problem to compute eigen problem with ScaLAPACK, INFO=",INFO
2635      MSG_ERROR(msg)
2636   endif
2637 
2638   ABI_FREE(IFAIl)
2639   ABI_FREE(ICLUSTR)
2640   ABI_FREE(GAP)
2641   if (allocated(IWORK))  then
2642     ABI_FREE(IWORK)
2643   end if
2644   if (allocated(RWORK))  then
2645     ABI_FREE(RWORK)
2646   end if
2647   if (allocated(CWORK))  then
2648     ABI_FREE(CWORK)
2649   end if
2650 #endif
2651   return
2652 
2653 end subroutine compute_generalized_eigen_problem

m_slk/descript_scalapack [ Types ]

[ Top ] [ m_slk ] [ Types ]

NAME

  descript_scalapack

FUNCTION

 Description of a ScaLAPACK matrix.

SOURCE

130  type,public :: descript_scalapack
131    integer :: tab(DLEN_)
132  end type descript_scalapack

m_slk/destruction_matrix_scalapack [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  destruction_matrix_scalapack

FUNCTION

  Destroys a matrix descriptor for ScaLAPACK.

PARENTS

      m_abilasi,m_exc_diago,m_slk

CHILDREN

SOURCE

585 subroutine destruction_matrix_scalapack(matrix)
586 
587 
588 !This section has been created automatically by the script Abilint (TD).
589 !Do not modify the following lines by hand.
590 #undef ABI_FUNC
591 #define ABI_FUNC 'destruction_matrix_scalapack'
592 !End of the abilint section
593 
594  implicit none
595 
596 !Arguments ------------------------------------
597  type(matrix_scalapack),intent(inout)    :: matrix
598 
599 ! *********************************************************************
600 
601 !MG This is dangerous as the pointers are not nullified when the object is initialized.
602 !Likely gfortran will bomb out here.
603 
604  NULLIFY(matrix%processor)
605 
606  matrix%sizeb_global = 0
607  if (allocated(matrix%buffer_cplx)) then
608    ABI_FREE(matrix%buffer_cplx)
609  end if
610  if (allocated(matrix%buffer_real)) then
611    ABI_FREE(matrix%buffer_real)
612  end if
613  if (allocated(matrix%ipiv)) then
614    ABI_FREE(matrix%ipiv)
615  end if
616 
617  matrix%sizeb_blocs = 0
618  matrix%sizeb_local = 0
619  matrix%descript%tab = 0
620 
621 end subroutine destruction_matrix_scalapack

m_slk/end_scalapack [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  end_scalapack

FUNCTION

  Removes a processor from the ScaLAPACK grid.

INPUTS

  None

OUTPUT

  None

SIDE EFFECTS

  processor= descriptor of a processor

PARENTS

      m_abi_linalg,m_abilasi,m_exc_diago

CHILDREN

SOURCE

435 subroutine end_scalapack(processor)
436 
437 
438 !This section has been created automatically by the script Abilint (TD).
439 !Do not modify the following lines by hand.
440 #undef ABI_FUNC
441 #define ABI_FUNC 'end_scalapack'
442 !End of the abilint section
443 
444  implicit none
445 
446 !Arguments ------------------------------------
447  type(processor_scalapack),intent(inout)    :: processor
448 
449 ! *********************************************************************
450 
451  call BLACS_GRIDEXIT(processor%grid%ictxt)
452  !call BLACS_EXIT(0)
453 
454 end subroutine end_scalapack

m_slk/glob_loc [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  glob_loc

FUNCTION

  Returns the global location of a matrix coefficient.

INPUTS

  matrix= the matrix to process
  idx= number of rows in the distributed matrix
  lico= block size index

PARENTS

SOURCE

873 function glob_loc(matrix,idx,lico)
874 
875 
876 !This section has been created automatically by the script Abilint (TD).
877 !Do not modify the following lines by hand.
878 #undef ABI_FUNC
879 #define ABI_FUNC 'glob_loc'
880 !End of the abilint section
881 
882  implicit none
883 
884 !Arguments ------------------------------------
885  integer, intent(in) :: idx, lico
886  type(matrix_scalapack),intent(in) :: matrix
887 
888 !Return value ---------------------------------
889  integer :: glob_loc
890 
891 !Local variables-------------------------------
892  integer,external :: NUMROC
893 
894 ! *********************************************************************
895 
896  glob_loc = NUMROC(idx,matrix%sizeb_blocs(lico), &
897 & matrix%processor%coords(lico),0, &
898 & matrix%processor%grid%dims(lico))
899 
900 end function glob_loc

m_slk/grid_scalapack [ Types ]

[ Top ] [ m_slk ] [ Types ]

NAME

  grid_scalapack

FUNCTION

  Grid of ScaLAPACK processors.

SOURCE

73  type,public :: grid_scalapack
74 
75    integer :: nbprocs
76    ! total number of processors
77 
78    integer   :: dims(2)
79 
80    integer :: ictxt
81    ! blacs context
82 
83  end type grid_scalapack
84 
85 #ifdef HAVE_LINALG_SCALAPACK
86  public :: build_grid_scalapack  ! Set up the processor grid for ScaLAPACK.
87 #endif

m_slk/idx_glob [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  idx_glob

FUNCTION

  Determination of the global indices of a term of the matrix with respect
  to its local indices.

INPUTS

  matrix= the matrix to process
  iloc= local row of the coefficient
  jloc= local column of the coefficient

OUTPUT

  i= row in the matrix
  j= column in the matrix

PARENTS

      m_slk,exc_diago

SOURCE

927 pure subroutine idx_glob(matrix,iloc,jloc,i,j)
928 
929 
930 !This section has been created automatically by the script Abilint (TD).
931 !Do not modify the following lines by hand.
932 #undef ABI_FUNC
933 #define ABI_FUNC 'idx_glob'
934 !End of the abilint section
935 
936  implicit none
937 
938 !Arguments ------------------------------------
939  integer, intent(in) :: iloc,jloc
940  integer, intent(out) :: i,j
941  type(matrix_scalapack),intent(in) :: matrix
942 
943 ! *********************************************************************
944 
945  i = loc_glob(matrix,matrix%processor,iloc,1)
946  j = loc_glob(matrix,matrix%processor,jloc,2)
947 
948 end subroutine idx_glob

m_slk/idx_loc [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  idx_loc

FUNCTION

  Determination of the local indices of a matrix coefficient with respect
  to its global indices, independently of the processor.

INPUTS

  matrix= the matrix to process
  i= row in the matrix
  j= column in the matrix

OUTPUT

  iloc= local row of the coefficient
  jloc= local column of the coefficient

PARENTS

      m_slk

CHILDREN

SOURCE

831 subroutine idx_loc(matrix,i,j,iloc,jloc)
832 
833 
834 !This section has been created automatically by the script Abilint (TD).
835 !Do not modify the following lines by hand.
836 #undef ABI_FUNC
837 #define ABI_FUNC 'idx_loc'
838 !End of the abilint section
839 
840  implicit none
841 
842 !Arguments ------------------------------------
843  integer, intent(in) :: i,j
844  integer, intent(out) :: iloc,jloc
845  type(matrix_scalapack),intent(in) :: matrix
846 
847 ! *********************************************************************
848 
849  iloc = glob_loc(matrix,i,1)
850  jloc = glob_loc(matrix,j,2)
851 
852 end subroutine idx_loc

m_slk/init_matrix_scalapack [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  init_matrix_scalapack

FUNCTION

  Initializes a matrix descriptor for ScaLAPACK.
  Initialisation of a SCALAPACK matrix (each proc initialize its own part of the matrix)

INPUTS

  processor= descriptor of a processor
  nbli_global= total number of lines
  nbco_global= total number of columns
  istwf_k= option parameter that describes the storage of wfs
  tbloc= custom block size

OUTPUT

  matrix= the matrix to process

PARENTS

      m_abilasi,m_exc_diago,m_slk,rayleigh_ritz

CHILDREN

SOURCE

484 subroutine init_matrix_scalapack(matrix,nbli_global,nbco_global,processor,istwf_k,tbloc)
485 
486 
487 !This section has been created automatically by the script Abilint (TD).
488 !Do not modify the following lines by hand.
489 #undef ABI_FUNC
490 #define ABI_FUNC 'init_matrix_scalapack'
491 !End of the abilint section
492 
493  implicit none
494 
495 !Arguments ------------------------------------
496  integer,intent(in) :: nbli_global,nbco_global,istwf_k
497  integer,intent(in),optional :: tbloc
498  type(matrix_scalapack),intent(inout) :: matrix
499  type(processor_scalapack),intent(in),target  :: processor
500 
501 !Local variables-------------------------------
502  INTEGER, PARAMETER :: SIZE_BLOCS = 24 ! As recommended by Intel MKL, a more sensible default than the previous value of 40
503  INTEGER            :: info,sizeb
504  integer,external :: NUMROC
505  character(len=500) :: msg
506 
507 ! *********************************************************************
508 
509  DBG_ENTER("COLL")
510 
511 #ifdef HAVE_LINALG_ELPA
512  sizeb  = 1
513 #else
514  sizeb = SIZE_BLOCS
515 #endif
516 
517 !Records of the matrix type :
518  matrix%processor => processor
519  matrix%sizeb_blocs(1) = MIN(sizeb,nbli_global)
520  matrix%sizeb_blocs(2) = MIN(sizeb,nbco_global)
521 
522 #ifdef HAVE_LINALG_ELPA
523  if(matrix%sizeb_blocs(1) .ne. matrix%sizeb_blocs(2)) then
524     matrix%sizeb_blocs(1) = MIN(matrix%sizeb_blocs(1),matrix%sizeb_blocs(2))
525     matrix%sizeb_blocs(2) = matrix%sizeb_blocs(1)
526  end if
527 #endif
528 
529  matrix%sizeb_global(1) = nbli_global
530  matrix%sizeb_global(2) = nbco_global
531 
532 !Size of the local buffer
533 !NUMROC computes the NUMber of Rows Or Columns of a distributed matrix owned by the process indicated by IPROC.
534  matrix%sizeb_local(1) = NUMROC(nbli_global,matrix%sizeb_blocs(1), &
535 & processor%coords(1),0, &
536 & processor%grid%dims(1))
537 
538  matrix%sizeb_local(2) = NUMROC(nbco_global,matrix%sizeb_blocs(2), &
539 & processor%coords(2),0, &
540 & processor%grid%dims(2))
541 
542  call idx_loc(matrix,matrix%sizeb_global(1),matrix%sizeb_global(2), &
543 & matrix%sizeb_local(1),matrix%sizeb_local(2))
544 
545 !Initialisation of the SCALAPACK description of the matrix
546  call DESCINIT(matrix%descript%tab, nbli_global, nbco_global, &
547 & matrix%sizeb_blocs(1), matrix%sizeb_blocs(2), 0,0 , &
548 & processor%grid%ictxt, MAX(1,matrix%sizeb_local(1)), &
549 & info)
550 
551  if (info /= 0) then
552    write(msg,'(2(a,i0))')" proc: ",processor%myproc,' error in the initialisation of the scalapack matrix: ',info
553    MSG_ERROR(msg)
554  end if
555 
556  if (istwf_k/=2) then
557    ABI_MALLOC(matrix%buffer_cplx,(matrix%sizeb_local(1),matrix%sizeb_local(2)))
558    matrix%buffer_cplx(:,:) = (0._DP,0._DP)
559  else
560    ABI_MALLOC(matrix%buffer_real,(matrix%sizeb_local(1),matrix%sizeb_local(2)))
561    matrix%buffer_real(:,:) = 0._DP
562  end if
563 
564  DBG_EXIT("COLL")
565 
566 end subroutine init_matrix_scalapack

m_slk/init_scalapack [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  init_scalapack

FUNCTION

  Initializes an instance of processor ScaLAPACK from an MPI communicator.

INPUTS

  comm= MPI communicator

OUTPUT

  processor= descriptor of a processor

PARENTS

      m_abi_linalg,m_abilasi,m_exc_diago

CHILDREN

SOURCE

375 subroutine init_scalapack(processor,comm)
376 
377 
378 !This section has been created automatically by the script Abilint (TD).
379 !Do not modify the following lines by hand.
380 #undef ABI_FUNC
381 #define ABI_FUNC 'init_scalapack'
382 !End of the abilint section
383 
384  implicit none
385 
386 !Arguments ------------------------------------
387  integer, intent(in) :: comm
388  type(processor_scalapack),intent(out) :: processor
389 
390 !Local variables-------------------------------
391  type(grid_scalapack) :: grid
392  integer :: nbproc,myproc,ierr
393 
394 ! *********************************************************************
395 
396  DBG_ENTER("COLL")
397 
398  call MPI_COMM_SIZE(comm, nbproc, ierr)
399  call MPI_COMM_RANK(comm, myproc, ierr)
400 
401  call build_grid_scalapack(grid, nbproc, comm)
402 
403  call build_processor_scalapack(processor, grid, myproc, comm)
404 
405  DBG_EXIT("COLL")
406 
407 end subroutine init_scalapack

m_slk/loc_glob [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  loc_glob

FUNCTION

  Determination of the global index from a local index (row or column)
  as a function of a given processor

INPUTS

  matrix= the matrix to process
  proc= descriptor of a processor
  idx= number of rows in the distributed matrix
  lico= block size index

PARENTS

SOURCE

971 pure function loc_glob(matrix,proc,idx,lico)
972 
973 
974 !This section has been created automatically by the script Abilint (TD).
975 !Do not modify the following lines by hand.
976 #undef ABI_FUNC
977 #define ABI_FUNC 'loc_glob'
978 !End of the abilint section
979 
980  implicit none
981 
982 !Arguments ------------------------------------
983  integer, intent(in) :: idx,lico
984  integer :: loc_glob
985  type(matrix_scalapack),intent(in) :: matrix
986  type(processor_scalapack),intent(in) :: proc
987 
988 !Local variables-------------------------------
989  integer :: nbcyc,reste,nblocs
990 
991 ! *********************************************************************
992 
993  nbcyc = INT((idx-1)/matrix%sizeb_blocs(lico))
994  reste = MOD(idx-1,matrix%sizeb_blocs(lico))
995  nblocs = nbcyc*proc%grid%dims(lico)+ proc%coords(lico)
996 
997  loc_glob = nblocs * matrix%sizeb_blocs(lico) + reste + 1
998 
999 end function loc_glob

m_slk/matrix_from_complexmatrix [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_from_complexmatrix

FUNCTION

  Routine to fill a SCALAPACK matrix from a global matrix.

INPUTS

  istwf_k= option parameter that describes the storage of wfs
  reference= a complex matrix

SIDE EFFECTS

  matrix= the matrix to process

PARENTS

      m_slk

CHILDREN

SOURCE

1267 subroutine matrix_from_complexmatrix(matrix,reference,istwf_k)
1268 
1269 
1270 !This section has been created automatically by the script Abilint (TD).
1271 !Do not modify the following lines by hand.
1272 #undef ABI_FUNC
1273 #define ABI_FUNC 'matrix_from_complexmatrix'
1274 !End of the abilint section
1275 
1276  implicit none
1277 
1278 !Arguments ------------------------------------
1279  integer,intent(in) :: istwf_k
1280  type(matrix_scalapack),intent(inout) :: matrix
1281 !arrays
1282  real(dp),intent(in) :: reference(:,:)
1283 
1284 !Local variables-------------------------------
1285  integer :: i,j,iglob,jglob,cptr
1286  real(dp) :: err
1287  complex(dpc) :: val
1288 
1289 ! *********************************************************************
1290 
1291 !err = 0._DP
1292 !cptr = 0
1293 
1294  do i=1,matrix%sizeb_local(1)
1295    do j=1,matrix%sizeb_local(2)
1296      call idx_glob(matrix,i,j,iglob,jglob)
1297 
1298      val = dcmplx(reference(2*iglob-1, jglob),reference(2*iglob, jglob))
1299      call matrix_set_local_cplx(matrix,i,j,val)
1300 
1301 !    cptr = cptr + 1
1302    end do
1303  end do
1304 
1305 !if (cptr /= 0) THEN
1306 !write(std_out,*) matrix%processor%myproc,"error Linf matrix scalapack", &
1307 !&  err,"on",cptr,"terms"
1308 !end if
1309 
1310 end subroutine matrix_from_complexmatrix

m_slk/matrix_from_global [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_from_global

FUNCTION

  Routine to fill a SCALAPACK matrix from a global matrix.

INPUTS

  istwf_k= option parameter that describes the storage of wfs
  reference= one-dimensional array

SIDE EFFECTS

  matrix= the matrix to process

PARENTS

CHILDREN

SOURCE

1024 subroutine matrix_from_global(matrix,reference,istwf_k)
1025 
1026 
1027 !This section has been created automatically by the script Abilint (TD).
1028 !Do not modify the following lines by hand.
1029 #undef ABI_FUNC
1030 #define ABI_FUNC 'matrix_from_global'
1031 !End of the abilint section
1032 
1033  implicit none
1034 
1035 !Arguments ------------------------------------
1036  integer,intent(in) :: istwf_k
1037  real(dp),intent(in) :: reference(*)
1038  type(matrix_scalapack),intent(inout) :: matrix
1039 
1040 !Local variables-------------------------------
1041  integer :: i,j,iglob,jglob,ind !cptr
1042  real(dp) :: err,val_real
1043  complex(dp)::val_cplx
1044  character(len=500) :: msg
1045 
1046 ! *********************************************************************
1047 
1048 !err = 0._DP
1049 !cptr = 0
1050 
1051  do i=1,matrix%sizeb_local(1)
1052    do j=1,matrix%sizeb_local(2)
1053      call idx_glob(matrix,i,j,iglob,jglob)
1054 
1055      if (istwf_k/=2) then
1056         ind = jglob*(jglob-1)+2*iglob-1
1057         val_cplx = dcmplx(reference(ind),reference(ind+1))
1058         call matrix_set_local_cplx(matrix,i,j,val_cplx)
1059     else
1060        ind = (jglob*(jglob-1))/2 + iglob
1061        val_real = reference(ind)
1062        call matrix_set_local_real(matrix,i,j,val_real)
1063 
1064 ! packed input now so this check has no sense anymore
1065 !        if(abs(reference(ind+1))>1.0d-10)then
1066 !          write(msg,'(2a,2i5,1es16.6,2a)')&
1067 ! &         '  For istwf_k=2, observed the following element of matrix :',ch10,&
1068 ! &         iglob,jglob,reference(ind+1),ch10,&
1069 ! &         '  with a non-negligible imaginary part.'
1070 !          MSG_BUG(msg)
1071 !        end if
1072 
1073      end if
1074 
1075 !    cptr = cptr + 1
1076    end do
1077  end do
1078 
1079 !if (cptr /= 0) then
1080 !write(std_out,*) matrix%processor%myproc,"error Linf matrix scalapack", &
1081 !&  err,"on",cptr,"terms"
1082 !endif
1083 
1084 end subroutine matrix_from_global

m_slk/matrix_from_global_sym [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_from_global_sym

FUNCTION

INPUTS

  istwf_k= option parameter that describes the storage of wfs
  reference= one-dimensional array

SIDE EFFECTS

  matrix= the matrix to process

PARENTS

CHILDREN

SOURCE

1108 subroutine matrix_from_global_sym(matrix,reference,istwf_k)
1109 
1110 
1111 !This section has been created automatically by the script Abilint (TD).
1112 !Do not modify the following lines by hand.
1113 #undef ABI_FUNC
1114 #define ABI_FUNC 'matrix_from_global_sym'
1115 !End of the abilint section
1116 
1117  implicit none
1118 
1119 !Arguments ------------------------------------
1120  integer,intent(in) :: istwf_k
1121  real(dp),intent(in) :: reference(:)
1122  type(matrix_scalapack),intent(inout)  :: matrix
1123 
1124 !Local variables-------------------------------
1125  integer :: i,j,iglob,jglob,ind !cptr
1126  real(dp) :: err
1127  complex(dp)::val_cplx
1128  real(dp)   ::val_real
1129  character(len=500) :: msg
1130 
1131 ! *********************************************************************
1132 
1133 !err = 0._DP
1134 !cptr = 0
1135 
1136  do i=1,matrix%sizeb_local(1)
1137    do j=1,matrix%sizeb_local(2)
1138      call idx_glob(matrix,i,j,iglob,jglob)
1139      if(jglob < iglob) then
1140         ind = iglob*(iglob-1)+2*jglob-1
1141      else
1142         ind = jglob*(jglob-1)+2*iglob-1
1143      end if
1144      if (istwf_k/=2) then
1145         val_cplx = dcmplx(reference(ind),reference(ind+1))
1146         if(jglob < iglob) then
1147            call matrix_set_local_cplx(matrix,i,j,conjg(val_cplx))
1148         else
1149            call matrix_set_local_cplx(matrix,i,j,val_cplx)
1150         end if
1151      else
1152         ind = (ind + 1) / 2
1153         val_real = reference(ind)
1154         call matrix_set_local_real(matrix,i,j,val_real)
1155         ! if(abs(reference(ind+1))>1.0d-10)then
1156         !    write(msg,'(2a,2i5,1es16.6,2a)')&
1157         !         &         '  For istwf_k=2, observed the following element of matrix :',ch10,&
1158         !         &         iglob,jglob,reference(ind+1),ch10,&
1159         !         &         '  with a non-negligible imaginary part.'
1160         !    MSG_BUG(msg)
1161         ! end if
1162      end if
1163 
1164      ! cptr = cptr + 1
1165    end do
1166  end do
1167 
1168 !if (cptr /= 0) then
1169 !write(std_out,*) matrix%processor%myproc,"error Linf matrix scalapack", &
1170 !&  err,"on",cptr,"terms"
1171 !endif
1172 
1173 end subroutine matrix_from_global_sym

m_slk/matrix_from_realmatrix [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_from_realmatrix

FUNCTION

  Routine to fill a SCALAPACK matrix from a real global matrix.

INPUTS

  istwf_k= option parameter that describes the storage of wfs
  reference= a real matrix

SIDE EFFECTS

  matrix= the matrix to process

PARENTS

      m_slk

CHILDREN

SOURCE

1199 subroutine matrix_from_realmatrix(matrix,reference,istwf_k)
1200 
1201 
1202 !This section has been created automatically by the script Abilint (TD).
1203 !Do not modify the following lines by hand.
1204 #undef ABI_FUNC
1205 #define ABI_FUNC 'matrix_from_realmatrix'
1206 !End of the abilint section
1207 
1208  implicit none
1209 
1210 !Arguments ------------------------------------
1211  integer,intent(in) :: istwf_k
1212  type(matrix_scalapack),intent(inout) :: matrix
1213 !arrays
1214  real(dp),intent(in) :: reference(:,:)
1215 
1216 !Local variables-------------------------------
1217  integer :: i,j,iglob,jglob,cptr
1218  real(dp) :: err,val
1219 
1220 ! *********************************************************************
1221 
1222 !err = 0._DP
1223 !cptr = 0
1224 
1225  do i=1,matrix%sizeb_local(1)
1226    do j=1,matrix%sizeb_local(2)
1227      call idx_glob(matrix,i,j,iglob,jglob)
1228 
1229      val = reference(iglob, jglob)
1230      call matrix_set_local_real(matrix,i,j,val)
1231 
1232 !    cptr = cptr + 1
1233    end do
1234  end do
1235 
1236 !if (cptr /= 0) THEN
1237 !write(std_out,*) matrix%processor%myproc,"error Linf matrix scalapack", &
1238 !&  err,"on",cptr,"terms"
1239 !end if
1240 
1241 end subroutine matrix_from_realmatrix

m_slk/matrix_get_local_cplx [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_get_local_cplx

FUNCTION

  Returns a local matrix coefficient of complex type.
 Access to a component thanks to its local indices

INPUTS

  matrix= the matrix to process
  i= row in the matrix
  j= column in the matrix

OUTPUT

  The value of the local matrix.

PARENTS

SOURCE

648 pure function matrix_get_local_cplx(matrix,i,j)
649 
650 
651 !This section has been created automatically by the script Abilint (TD).
652 !Do not modify the following lines by hand.
653 #undef ABI_FUNC
654 #define ABI_FUNC 'matrix_get_local_cplx'
655 !End of the abilint section
656 
657  implicit none
658 
659 !Arguments ------------------------------------
660  integer, intent(in) :: i,j
661  complex(dp) :: matrix_get_local_cplx
662  type(matrix_scalapack),intent(in) :: matrix
663 
664 ! *********************************************************************
665 
666  matrix_get_local_cplx = matrix%buffer_cplx(i,j)
667 
668 end function matrix_get_local_cplx

m_slk/matrix_get_local_real [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_get_local_real

FUNCTION

  Returns a local matrix coefficient of double precision type.

INPUTS

  matrix= the matrix to process
  i= row in the matrix
  j= column in the matrix

PARENTS

SOURCE

689 pure function matrix_get_local_real(matrix,i,j)
690 
691 
692 !This section has been created automatically by the script Abilint (TD).
693 !Do not modify the following lines by hand.
694 #undef ABI_FUNC
695 #define ABI_FUNC 'matrix_get_local_real'
696 !End of the abilint section
697 
698  implicit none
699 
700 !Arguments ------------------------------------
701  integer, intent(in) :: i,j
702  type(matrix_scalapack),intent(in) :: matrix
703  real(dp) :: matrix_get_local_real
704 
705 ! *********************************************************************
706 
707  matrix_get_local_real = matrix%buffer_real(i,j)
708 
709 end function matrix_get_local_real

m_slk/matrix_scalapack [ Types ]

[ Top ] [ m_slk ] [ Types ]

NAME

FUNCTION

 The local part of a ScaLAPACK matrix

SOURCE

145  type,public :: matrix_scalapack
146 
147    integer :: sizeb_local(2)
148      ! dimensions of the buffer
149 
150    integer :: sizeb_global(2)
151 
152    integer :: sizeb_blocs(2)
153      ! size of the block of consecutive data
154 
155    integer,allocatable :: ipiv(:)
156 
157    real(dp),allocatable  :: buffer_real(:,:)
158      ! buffer of the local part of the matrix
159 
160    complex(dpc),allocatable :: buffer_cplx(:,:)
161      ! buffer of the local part of the matrix
162 
163    type(processor_scalapack),pointer :: processor
164 
165    type(descript_scalapack) :: descript
166 
167  end type matrix_scalapack
168 
169 
170 #ifdef HAVE_LINALG_SCALAPACK
171  public :: init_matrix_scalapack           ! Initialisation of a SCALAPACK matrix (each proc initialize its own part of the matrix)
172  public :: destruction_matrix_scalapack    ! Destroys a matrix descriptor for ScaLAPACK.
173  public :: matrix_get_local_cplx           ! Returns a local matrix coefficient of complex type.
174  public :: matrix_get_local_real           ! Returns a local matrix coefficient of double precision type.
175  public :: matrix_set_local_cplx           ! Sets a local matrix coefficient of complex type.
176  public :: matrix_set_local_real           ! Sets a local matrix coefficient of double precision type.
177  public :: idx_loc                         ! Local indices of a matrix coefficient with respect to its global indices, independently of the processor.
178  public :: glob_loc                        ! Returns the global location of a matrix coefficient.
179  public :: idx_glob                        ! Determination of the global indices of a term of the matrix with respect to its local indices.
180  public :: loc_glob                        ! Determination of the global index from a local index (row or column) as a function of a given processor
181  public :: matrix_from_global              ! Fills a SCALAPACK matrix with respect to a full matrix.
182  public :: matrix_from_global_sym          ! Fills a SCALAPACK matrix with respect to a full matrix.
183  public :: matrix_from_realmatrix          ! Fills a SCALAPACK matrix with respect to a full matrix.
184  public :: matrix_from_complexmatrix       ! Fills a SCALAPACK matrix with respect to a full matrix.
185  public :: matrix_to_global                ! Inserts a ScaLAPACK matrix into a global one.
186  public :: matrix_to_realmatrix            ! Inserts a ScaLAPACK matrix into a real matrix.
187  public :: matrix_to_complexmatrix         ! Inserts a ScaLAPACK matrix into a complex matrix.
188  public :: matrix_to_reference             ! Fill a full matrix with respect to a SCALAPACK matrix.
189  public :: slk_matrix_from_global_dpc_2D   ! Fill a complex SCALAPACK matrix with respect to a global matrix.
190  public :: slk_matrix_from_global_dpc_1Dp  ! Fill a complex SCALAPACK matrix with respect to a global matrix.
191                                            !   target: double precision complex matrix in packed form.
192  public :: slk_matrix_to_global_dpc_2D     ! Fill a global matrix with respect to a SCALAPACK matrix.
193                                            !   target: Two-dimensional Double precision complex matrix.
194  !public :: my_locr
195  !public :: my_locc
196 
197  public :: slk_pzgemm                        ! C := alpha*A*B - beta*C
198  public :: compute_eigen_problem             ! Calculation of eigenvalues and eigenvectors. A * X = lambda * X complex and real cases.
199  public :: compute_generalized_eigen_problem ! compute_generalized_eigen_problem
200  public :: compute_eigen1                    ! Calculation of eigenvalues and eigenvectors.  complex and real cases.
201  public :: compute_eigen2                    ! Calculation of eigenvalues and eigenvectors: A * X = lambda * B * X complex and real cases.
202  public :: slk_pzheev                        ! Eigenvalues and, optionally, eigenvectors of an Hermitian matrix A.  A * X = lambda * X
203  public :: slk_pzheevx                       ! Eigenvalues and, optionally, eigenvectors of a complex hermitian matrix A. A * X = lambda *  X
204  public :: slk_pzhegvx                       ! genvalues and, optionally, eigenvectors of a complex generalized  Hermitian-definite eigenproblem, of the form
205                                              ! sub( A )*x=(lambda)*sub( B )*x,  sub( A )*sub( B )x=(lambda)*x,  or sub( B )*sub( A )*x=(lambda)*x.
206  public :: slk_zinvert                       ! Compute the inverse of a complex matrix in double precision
207  public :: slk_zdhp_invert                   ! Computes the inverse of a Hermitian positive definite matrix.
208  public :: slk_write                         ! Writes a square scaLAPACK distributed matrix on an external file using MPI-IO.
209  public :: slk_read                          ! Read a square scaLAPACK distributed matrix from an external file using MPI-IO.
210  public :: slk_single_fview_read_mask        ! Returns an MPI datatype that can be used to read a scaLAPACK distributed matrix from
211                                              !   a binary file using MPI-IO. The view is created using the user-defined mask function
212  public :: slk_symmetrize                    ! Symmetrizes a square scaLAPACK-distributed matrix.
213  public :: slk_single_fview_read             ! Returns an MPI datatype to read a scaLAPACK distributed matrix from a binary file using MPI-IO.
214  public :: slk_single_fview_write            ! Returns an MPI datatype to write a scaLAPACK distributed matrix to a binary file using MPI-IO.
215  public :: slk_bsize_and_type                ! Returns the byte size and the MPI datatype associated to the matrix elements
216                                              !  that are stored in the ScaLAPACK_matrix
217 #endif
218 
219 CONTAINS  !==============================================================================

m_slk/matrix_set_local_cplx [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_set_local_cplx

FUNCTION

  Sets a local matrix coefficient of complex type.
 -------------------------------------------------------
  Positioning of a component of a matrix thanks to its local indices
 -------------------------------------------------------

INPUTS

  i= row in the matrix
  j= column in the matrix
  value= the value to set

SIDE EFFECTS

  matrix%buffer_cplx(i,j) filled with value

PARENTS

      m_slk

SOURCE

737 pure subroutine matrix_set_local_cplx(matrix,i,j,value)
738 
739 
740 !This section has been created automatically by the script Abilint (TD).
741 !Do not modify the following lines by hand.
742 #undef ABI_FUNC
743 #define ABI_FUNC 'matrix_set_local_cplx'
744 !End of the abilint section
745 
746  implicit none
747 
748 !Arguments ------------------------------------
749  integer, intent(in) :: i,j
750  complex(dp), intent(in) :: value
751  type(matrix_scalapack),intent(inout) :: matrix
752 
753 ! *********************************************************************
754 
755  matrix%buffer_cplx(i,j) = value
756 
757 end subroutine matrix_set_local_cplx

m_slk/matrix_set_local_real [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_set_local_real

FUNCTION

  Sets a local matrix coefficient of double precision type.

INPUTS

  i= row in the matrix
  j= column in the matrix
  value= the value to set

SIDE EFFECTS

  matrix%buffer_real(i,j) set to value

PARENTS

      m_slk

SOURCE

782 pure subroutine matrix_set_local_real(matrix,i,j,value)
783 
784 
785 !This section has been created automatically by the script Abilint (TD).
786 !Do not modify the following lines by hand.
787 #undef ABI_FUNC
788 #define ABI_FUNC 'matrix_set_local_real'
789 !End of the abilint section
790 
791  implicit none
792 
793 !Arguments ------------------------------------
794  integer, intent(in) :: i,j
795  real(dp), intent(in) :: value
796  type(matrix_scalapack),intent(inout) :: matrix
797 
798 ! *********************************************************************
799 
800  matrix%buffer_real(i,j) = value
801 
802 end subroutine matrix_set_local_real

m_slk/matrix_to_complexmatrix [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_to_complexmatrix

FUNCTION

  Inserts a ScaLAPACK matrix into a complex matrix.

INPUTS

  matrix= the matrix to process
  istwf_k= option parameter that describes the storage of wfs

SIDE EFFECTS

  reference= the matrix to fill

PARENTS

      m_slk

CHILDREN

SOURCE

1477 subroutine matrix_to_complexmatrix(matrix,reference,istwf_k)
1478 
1479 
1480 !This section has been created automatically by the script Abilint (TD).
1481 !Do not modify the following lines by hand.
1482 #undef ABI_FUNC
1483 #define ABI_FUNC 'matrix_to_complexmatrix'
1484 !End of the abilint section
1485 
1486  implicit none
1487 
1488 !Arguments ------------------------------------
1489  integer,intent(in) :: istwf_k
1490  type(matrix_scalapack),intent(in) :: matrix
1491 !arrays
1492  complex(dpc),intent(inout) :: reference(:,:)
1493 
1494 !Local variables-------------------------------
1495  integer  :: i,j,iglob,jglob,cptr
1496  real(dp) :: err
1497 
1498 ! *********************************************************************
1499 
1500 !err = 0._DP
1501 !cptr = 0
1502 
1503  do i=1,matrix%sizeb_local(1)
1504    do j=1,matrix%sizeb_local(2)
1505      call idx_glob(matrix,i,j,iglob,jglob)
1506 
1507      reference(iglob,jglob) = matrix_get_local_cplx(matrix,i,j)
1508 !    cptr = cptr + 1
1509    end do
1510  end do
1511 
1512 !if (cptr /= 0) THEN
1513 !write(std_out,*) matrix%processor%myproc,"erreur Linf matrix scalapack", &
1514 !&  err,"on",cptr,"terms"
1515 !end if
1516 
1517 end subroutine matrix_to_complexmatrix

m_slk/matrix_to_global [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_to_global

FUNCTION

  Inserts a ScaLAPACK matrix into a global one.

INPUTS

  matrix= the matrix to process
  istwf_k= option parameter that describes the storage of wfs
  nband_k= number of bands at this k point for that spin polarization

SIDE EFFECTS

  reference= one-dimensional array

PARENTS

CHILDREN

SOURCE

1336 subroutine matrix_to_global(matrix,reference,istwf_k)
1337 
1338 
1339 !This section has been created automatically by the script Abilint (TD).
1340 !Do not modify the following lines by hand.
1341 #undef ABI_FUNC
1342 #define ABI_FUNC 'matrix_to_global'
1343 !End of the abilint section
1344 
1345  implicit none
1346 
1347 !Arguments ------------------------------------
1348  integer,intent(in) :: istwf_k          !,nband_k
1349  real(dp),intent(inout) :: reference(*) !(nband_k*(nband_k+1))
1350  type(matrix_scalapack),intent(in) :: matrix
1351 
1352 !Local variables-------------------------------
1353  integer  :: i,j,iglob,jglob,ind,cptr
1354  real(dp) :: err
1355 
1356 ! *********************************************************************
1357 
1358 !err = 0._DP
1359 !cptr = 0
1360 
1361  do i=1,matrix%sizeb_local(1)
1362    do j=1,matrix%sizeb_local(2)
1363      call idx_glob(matrix,i,j,iglob,jglob)
1364 
1365      ind = jglob*(jglob-1)+2*iglob-1
1366      if (ind <= matrix%sizeb_global(2)*(matrix%sizeb_global(2)+1)) then
1367         if (istwf_k/=2) then
1368          reference(ind)   = real(matrix_get_local_cplx(matrix,i,j))
1369          reference(ind+1) = IMAG(matrix_get_local_cplx(matrix,i,j))
1370        else
1371           ind=(ind+1)/2 !real packed storage
1372           reference(ind) = matrix_get_local_real(matrix,i,j)
1373        end if
1374      end if
1375 !    cptr = cptr + 1
1376    end do
1377  end do
1378 
1379 !if (cptr /= 0) then
1380 !write(std_out,*) matrix%processor%myproc,"erreur Linf matrix scalapack", &
1381 !&  err,"on",cptr,"terms"
1382 !endif
1383 
1384 end subroutine matrix_to_global

m_slk/matrix_to_realmatrix [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_to_realmatrix

FUNCTION

  Inserts a ScaLAPACK matrix into a real matrix.

INPUTS

  matrix= the matrix to process
  istwf_k= option parameter that describes the storage of wfs

SIDE EFFECTS

  reference= the matrix to fill

PARENTS

      m_slk

CHILDREN

SOURCE

1410 subroutine matrix_to_realmatrix(matrix,reference,istwf_k)
1411 
1412 
1413 !This section has been created automatically by the script Abilint (TD).
1414 !Do not modify the following lines by hand.
1415 #undef ABI_FUNC
1416 #define ABI_FUNC 'matrix_to_realmatrix'
1417 !End of the abilint section
1418 
1419  implicit none
1420 
1421 !Arguments ------------------------------------
1422  integer,intent(in) :: istwf_k
1423  type(matrix_scalapack),intent(in) :: matrix
1424 !arrays
1425  real(dp),intent(inout) :: reference(:,:)
1426 
1427 !Local variables-------------------------------
1428  integer      :: i,j,iglob,jglob,cptr
1429  real(dp)     :: err
1430  complex(dpc) :: zvar
1431 
1432 ! *********************************************************************
1433 
1434 !err = 0._DP
1435 !cptr = 0
1436 
1437  do i=1,matrix%sizeb_local(1)
1438    do j=1,matrix%sizeb_local(2)
1439      call idx_glob(matrix,i,j,iglob,jglob)
1440 
1441      reference(iglob,jglob) = matrix_get_local_real(matrix,i,j)
1442 !    cptr = cptr + 1
1443    end do
1444  end do
1445 
1446 !if (cptr /= 0) THEN
1447 !write(std_out,*) matrix%processor%myproc,"erreur Linf matrix scalapack", &
1448 !&  err,"on",cptr,"terms"
1449 !end if
1450 
1451 end subroutine matrix_to_realmatrix

m_slk/matrix_to_reference [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  matrix_to_reference

FUNCTION

  Routine to fill a full matrix with respect to a SCALAPACK matrix.

INPUTS

  matrix= the matrix to process
  istwf_k= option parameter that describes the storage of wfs

SIDE EFFECTS

  reference= one-dimensional array

PARENTS

CHILDREN

SOURCE

1542 subroutine matrix_to_reference(matrix,reference,istwf_k)
1543 
1544 
1545 !This section has been created automatically by the script Abilint (TD).
1546 !Do not modify the following lines by hand.
1547 #undef ABI_FUNC
1548 #define ABI_FUNC 'matrix_to_reference'
1549 !End of the abilint section
1550 
1551  implicit none
1552 
1553 !Arguments ------------------------------------
1554  integer,intent(in) :: istwf_k
1555  type(matrix_scalapack),intent(in) :: matrix
1556 !arrays
1557  real(dp),intent(inout) :: reference(:,:)
1558 
1559 !Local variables-------------------------------
1560  integer  :: i,j,iglob,jglob,ind,cptr
1561  real(dp) :: err
1562 
1563 ! *********************************************************************
1564 
1565 !err = 0._DP
1566 !cptr = 0
1567 
1568  do i=1,matrix%sizeb_local(1)
1569    do j=1,matrix%sizeb_local(2)
1570      call idx_glob(matrix,i,j,iglob,jglob)
1571 
1572      if (istwf_k/=2) then
1573        ind=(iglob-1)*2+1
1574        reference(ind,jglob)   = real(matrix_get_local_cplx(matrix,i,j))
1575        reference(ind+1,jglob) = IMAG(matrix_get_local_cplx(matrix,i,j))
1576      else
1577         ind=iglob
1578         reference(ind,jglob)   = matrix_get_local_real(matrix,i,j)
1579         !reference(ind+1,jglob) = 0._dp
1580      end if
1581 
1582 !    cptr = cptr + 1
1583    end do
1584  end do
1585 
1586 !if (cptr /= 0) then
1587 !write(std_out,*) matrix%processor%myproc,"error Linf matrix scalapack", &
1588 !&  err,"on",cptr,"terms"
1589 !endif
1590 
1591 end subroutine matrix_to_reference

m_slk/my_locc [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

 my_locc

FUNCTION

  Method of matrix_scalapack wrapping the scaLAPACK tool LOCc.

INPUTS

  Slk_mat<matrix_scalapack>

OUTPUT

  my_locc= For the meaning see NOTES below.

NOTES

  Let K be the number of rows or columns of a distributed matrix, and assume that its process grid has dimension p x q.
  LOCr( K ) denotes the number of elements of K that a process would receive if K were distributed over the p
  processes of its process column.
  Similarly,  LOCc(  K  ) denotes the number of elements of K that a process would receive if K were distributed over
  the q processes of its process row.
  The values of LOCr() and LOCc() may be determined via a call to the ScaLAPACK tool function, NUMROC:
          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).  An upper bound for these quantities may  be  computed
  by:
          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A

PARENTS

SOURCE

1973 function my_locc(Slk_mat)
1974 
1975 
1976 !This section has been created automatically by the script Abilint (TD).
1977 !Do not modify the following lines by hand.
1978 #undef ABI_FUNC
1979 #define ABI_FUNC 'my_locc'
1980 !End of the abilint section
1981 
1982  implicit none
1983 
1984 !Arguments ------------------------------------
1985 !scalars
1986  integer :: my_locc
1987  type(matrix_scalapack),intent(in) :: Slk_mat
1988 
1989 !Local variables-------------------------------
1990  integer :: N, NB_A, MYCOL, CSRC_A, NPCOL
1991  integer,external :: NUMROC
1992 
1993 ! *************************************************************************
1994 
1995  N      = Slk_mat%descript%tab(N_ )      ! The number of columns in the global matrix.
1996  NB_A   = Slk_mat%descript%tab(NB_)      ! The number of columns in a block.
1997  MYCOL  = Slk_mat%processor%coords(2)    ! The column index of my processor
1998  CSRC_A = Slk_mat%descript%tab(CSRC_)    ! The column of the processors at the beginning.
1999  NPCOL  = Slk_mat%processor%grid%dims(2) ! The number of processors per column in the Scalapack grid.
2000 
2001  my_locc = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL )
2002 
2003 end function my_locc

m_slk/my_locr [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

 my_locr

FUNCTION

  Method of matrix_scalapack wrapping the scaLAPACK tool LOCr.

INPUTS

  Slk_mat<matrix_scalapack>

OUTPUT

  my_locr= For the meaning see NOTES below.

NOTES

  Let K be the number of rows or columns of a distributed matrix, and assume that its process grid has dimension p x q.
  LOCr( K ) denotes the number of elements of K that a process would receive if K were distributed over the p
  processes of its process column.
  Similarly,  LOCc(  K  ) denotes the number of elements of K that a process would receive if K were distributed over
  the q processes of its process row.
  The values of LOCr() and LOCc() may be determined via a call to the ScaLAPACK tool function, NUMROC:
          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).  An upper bound for these quantities may  be  computed
  by:
          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A

PARENTS

SOURCE

1908 function my_locr(Slk_mat)
1909 
1910 
1911 !This section has been created automatically by the script Abilint (TD).
1912 !Do not modify the following lines by hand.
1913 #undef ABI_FUNC
1914 #define ABI_FUNC 'my_locr'
1915 !End of the abilint section
1916 
1917  implicit none
1918 
1919 !Arguments ------------------------------------
1920 !scalars
1921  integer :: my_locr
1922  type(matrix_scalapack),intent(in) :: Slk_mat
1923 
1924 !Local variables-------------------------------
1925  integer :: M, MB_A, MYROW, RSRC_A, NPROW
1926  integer,external :: NUMROC
1927 
1928 ! *************************************************************************
1929 
1930  M      = Slk_mat%descript%tab(M_ )      ! The number of rows in the global matrix.
1931  MB_A   = Slk_mat%descript%tab(MB_)      ! The number of rows in a block.
1932  MYROW  = Slk_mat%processor%coords(1)    ! The row index of my processor
1933  RSRC_A = Slk_mat%descript%tab(RSRC_)    ! The row of the processors at the beginning.
1934  NPROW  = Slk_mat%processor%grid%dims(1) ! The number of processors per row in the Scalapack grid.
1935 
1936  my_locr = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW )
1937 
1938 end function my_locr

m_slk/no_scalapack [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  no_scalapack

FUNCTION

   Empty placeholder.

PARENTS

CHILDREN

SOURCE

5170 #else
5171 
5172 subroutine no_scalapack()
5173 
5174 
5175 !This section has been created automatically by the script Abilint (TD).
5176 !Do not modify the following lines by hand.
5177 #undef ABI_FUNC
5178 #define ABI_FUNC 'no_scalapack'
5179 !End of the abilint section
5180 
5181  implicit none
5182 
5183 ! *************************************************************************
5184 
5185 end subroutine no_scalapack

m_slk/processor_scalapack [ Types ]

[ Top ] [ m_slk ] [ Types ]

NAME

  processor_scalapack

FUNCTION

  One processor in the grid.

SOURCE

101  type,public :: processor_scalapack
102 
103    integer :: myproc
104    ! number of the processor
105    integer :: comm
106    ! MPI communicator MPI underlying the grid BLACS
107    integer :: coords(2)
108 
109    type(grid_scalapack) :: grid   ! the grid to which the processor is associated
110  end type processor_scalapack
111 
112 #ifdef HAVE_LINALG_SCALAPACK
113  public :: build_processor_scalapack     ! Builds a processor descriptor for ScaLAPACK.
114  public :: init_scalapack                ! Initializes an instance of processor ScaLAPACK from an MPI communicator.
115  public :: end_scalapack                 ! Removes a processor from the ScaLAPACK grid.
116 #endif

m_slk/slk_bsize_and_type [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_bsize_and_type

FUNCTION

  Returns the byte size and the MPI datatype associated to the matrix elements
  that are stored in the ScaLAPACK_matrix

INPUTS

  Slk_mat<matrix_scalapack>=Structured datatype defining the scaLAPACK distribution with the local buffer

OUTPUT

  bsize_elm=Byte size of the matrix element.
  mpi_type_elm=MPI datatype of the matrix element.

PARENTS

      m_slk

CHILDREN

SOURCE

4985 subroutine slk_bsize_and_type(Slk_mat,bsize_elm,mpi_type_elm)
4986 
4987 
4988 !This section has been created automatically by the script Abilint (TD).
4989 !Do not modify the following lines by hand.
4990 #undef ABI_FUNC
4991 #define ABI_FUNC 'slk_bsize_and_type'
4992 !End of the abilint section
4993 
4994  implicit none
4995 
4996 !Arguments ------------------------------------
4997 !scalars
4998  type(matrix_scalapack),intent(in) :: Slk_mat
4999  integer,intent(out) :: bsize_elm,mpi_type_elm
5000 
5001 !Local variables ------------------------------
5002 !scalars
5003  integer :: ierr
5004  character(len=500) :: msg
5005 
5006 ! ************************************************************************
5007 
5008 !@matrix_scalapack
5009  ierr=0
5010  if (allocated(Slk_mat%buffer_cplx)) then
5011    ierr=ierr+1
5012    mpi_type_elm = MPI_DOUBLE_complex
5013    bsize_elm    = xmpi_bsize_dpc
5014  end if
5015 
5016  if (allocated(Slk_mat%buffer_real)) then
5017    ierr=ierr+1
5018    mpi_type_elm = MPI_DOUBLE_PRECISION
5019    bsize_elm    = xmpi_bsize_dp
5020  end if
5021 !
5022 !One and only one buffer should be allocated.
5023  if (ierr/=1) then
5024    write(msg,'(a,i0)')" ScaLAPACK buffers are not allocated correctly, ierr= ",ierr
5025    MSG_ERROR(msg)
5026  end if
5027 
5028 end subroutine slk_bsize_and_type

m_slk/slk_matrix_from_global_dpc_1Dp [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_matrix_from_global_dpc_1Dp

FUNCTION

  Routine to fill a complex SCALAPACK matrix with respect to a global matrix.
  target: double precision complex matrix in packed form.

INPUTS

  glob_pmat(n*(n+1)/2)=One-dimensional array containing the global matrix A packed columnwise in a linear array.
    The j-th column of A is stored in the array glob_pmat as follows:
      if uplo = "U", glob_pmat(i + (j-1)*j/2)       = A(i,j) for 1<=i<=j;
      if uplo = "L", glob_pmat(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
      where n is the number of rows or columns in the global matrix.
  uplo=String specifying whether only the upper or lower triangular part of the global matrix is used:
    = "U":  Upper triangular
    = "L":  Lower triangular

SIDE EFFECTS

  Slk_mat<matrix_scalapack>=The distributed matrix.
    %buffer_cplx=Local buffer containg the value this node is dealing with.

PARENTS

      m_abilasi

CHILDREN

SOURCE

1719 subroutine slk_matrix_from_global_dpc_1Dp(Slk_mat,uplo,glob_pmat)
1720 
1721 
1722 !This section has been created automatically by the script Abilint (TD).
1723 !Do not modify the following lines by hand.
1724 #undef ABI_FUNC
1725 #define ABI_FUNC 'slk_matrix_from_global_dpc_1Dp'
1726 !End of the abilint section
1727 
1728  implicit none
1729 
1730 !Arguments ------------------------------------
1731 !scalars
1732  character(len=*),intent(in) :: uplo
1733  type(matrix_scalapack),intent(inout)  :: Slk_mat
1734 !array
1735  complex(dpc),intent(in) :: glob_pmat(:)
1736 
1737 !Local variables-------------------------------
1738  integer :: ii,jj,iglob,jglob,ind,n
1739  real(dp) :: szm
1740 
1741 !************************************************************************
1742 
1743  ABI_CHECK(allocated(Slk_mat%buffer_cplx),"%buffer_cplx not allocated")
1744 
1745  szm = SIZE(glob_pmat)
1746  n = NINT( (-1 + SQRT(one+8*szm) )*half )
1747  if (n*(n+1)/2 /= SIZE(glob_pmat)) then
1748    MSG_ERROR("Buggy compiler")
1749  end if
1750 
1751  select case (uplo(1:1))
1752 
1753    case ("U","u") ! Only the upper triangle of the global matrix is used.
1754      do jj=1,Slk_mat%sizeb_local(2)
1755        do ii=1,Slk_mat%sizeb_local(1)
1756          call idx_glob(Slk_mat,ii,jj,iglob,jglob)
1757 
1758          if (jglob>=iglob) then
1759            ind = iglob + jglob*(jglob-1)/2
1760            Slk_mat%buffer_cplx(ii,jj) =        glob_pmat(ind)
1761          else
1762            ind = jglob + iglob*(iglob-1)/2
1763            Slk_mat%buffer_cplx(ii,jj) = DCONJG( glob_pmat(ind) )
1764          end if
1765 
1766        end do
1767      end do
1768 
1769    case ("L","l") ! Only the lower triangle of the global matrix is used.
1770      do jj=1,Slk_mat%sizeb_local(2)
1771        do ii=1,Slk_mat%sizeb_local(1)
1772          call idx_glob(Slk_mat,ii,jj,iglob,jglob)
1773 
1774          if (jglob<=iglob) then
1775            ind = iglob + (jglob-1)*(2*n-jglob)/2
1776            Slk_mat%buffer_cplx(ii,jj) =        glob_pmat(ind)
1777          else
1778            ind = jglob + (iglob-1)*(2*n-iglob)/2
1779            Slk_mat%buffer_cplx(ii,jj) = DCONJG( glob_pmat(ind) )
1780          end if
1781 
1782        end do
1783      end do
1784 
1785      case default
1786      MSG_BUG(" Wrong uplo: "//TRIM(uplo))
1787  end select
1788 
1789 end subroutine slk_matrix_from_global_dpc_1Dp

m_slk/slk_matrix_from_global_dpc_2D [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_matrix_from_global_dpc_2D

FUNCTION

  Routine to fill a complex SCALAPACK matrix with respect to a global matrix.
  target: Two-dimensional double precision complex matrix

INPUTS

  glob_mat=Two-dimensional array containing the global matrix.
  uplo=String specifying whether only the upper or lower triangular part of the global matrix is used:
    = "U":  Upper triangular
    = "L":  Lower triangular
    = "A":  Full matrix (used for general complex matrices)

SIDE EFFECTS

  Slk_mat<matrix_scalapack>=The distributed matrix.
    %buffer_cplx=Local buffer containg the value this node is dealing with.

PARENTS

      m_abilasi

CHILDREN

SOURCE

1622 subroutine slk_matrix_from_global_dpc_2D(Slk_mat,uplo,glob_mat)
1623 
1624 
1625 !This section has been created automatically by the script Abilint (TD).
1626 !Do not modify the following lines by hand.
1627 #undef ABI_FUNC
1628 #define ABI_FUNC 'slk_matrix_from_global_dpc_2D'
1629 !End of the abilint section
1630 
1631  implicit none
1632 
1633 !Arguments ------------------------------------
1634 !scalars
1635  character(len=*),intent(in) :: uplo
1636  type(matrix_scalapack),intent(inout)  :: Slk_mat
1637 !array
1638  complex(dpc),intent(in) :: glob_mat(:,:)
1639 
1640 !Local variables-------------------------------
1641  integer :: ii,jj,iglob,jglob
1642 
1643 !************************************************************************
1644 
1645  ABI_CHECK(allocated(Slk_mat%buffer_cplx),"%buffer_cplx not allocated")
1646 
1647  select case (uplo(1:1))
1648 
1649    case ("A","a") ! Full global matrix is used.
1650      do jj=1,Slk_mat%sizeb_local(2)
1651        do ii=1,Slk_mat%sizeb_local(1)
1652          call idx_glob(Slk_mat,ii,jj,iglob,jglob)
1653          Slk_mat%buffer_cplx(ii,jj) = glob_mat(iglob,jglob)
1654        end do
1655      end do
1656 
1657    case ("U","u") ! Only the upper triangle of the global matrix is used.
1658      do jj=1,Slk_mat%sizeb_local(2)
1659        do ii=1,Slk_mat%sizeb_local(1)
1660          call idx_glob(Slk_mat,ii,jj,iglob,jglob)
1661          if (jglob>=iglob) then
1662            Slk_mat%buffer_cplx(ii,jj) =        glob_mat(iglob,jglob)
1663          else
1664            Slk_mat%buffer_cplx(ii,jj) = DCONJG( glob_mat(jglob,iglob) )
1665          end if
1666        end do
1667      end do
1668 
1669    case ("L","l") ! Only the lower triangle of the global matrix is used.
1670      do jj=1,Slk_mat%sizeb_local(2)
1671        do ii=1,Slk_mat%sizeb_local(1)
1672          call idx_glob(Slk_mat,ii,jj,iglob,jglob)
1673          if (jglob<=iglob) then
1674            Slk_mat%buffer_cplx(ii,jj) =        glob_mat(iglob,jglob)
1675          else
1676            Slk_mat%buffer_cplx(ii,jj) = DCONJG( glob_mat(jglob,iglob) )
1677          end if
1678        end do
1679      end do
1680 
1681      case default
1682      MSG_BUG(" Wrong uplo: "//TRIM(uplo))
1683  end select
1684 
1685 end subroutine slk_matrix_from_global_dpc_2D

m_slk/slk_matrix_to_global_dpc_2D [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_matrix_to_global_dpc_2D

FUNCTION

  Routine to fill a global matrix with respect to a SCALAPACK matrix.
  target: Two-dimensional Double precision complex matrix.

INPUTS

  Slk_mat<matrix_scalapack>=The distributed matrix.
  uplo=String specifying whether the upper or lower triangular part of the global matrix has to be filled:
    = "U":  Upper triangular
    = "L":  Lower triangular
    = "A":  Full matrix is filled (used for general complex matrices)

SIDE EFFECTS

  glob_mat=The global matrix where the entries owned by this processors have been overwritten.
  Note that the remaing entries not treated by this node are not changed.

PARENTS

      m_abilasi

CHILDREN

SOURCE

1820 subroutine slk_matrix_to_global_dpc_2D(Slk_mat,uplo,glob_mat)
1821 
1822 
1823 !This section has been created automatically by the script Abilint (TD).
1824 !Do not modify the following lines by hand.
1825 #undef ABI_FUNC
1826 #define ABI_FUNC 'slk_matrix_to_global_dpc_2D'
1827 !End of the abilint section
1828 
1829  implicit none
1830 
1831 !Arguments ------------------------------------
1832 !scalaras
1833  character(len=*),intent(in) :: uplo
1834  type(matrix_scalapack),intent(in) :: Slk_mat
1835 !arrays
1836  complex(dpc),intent(inout) :: glob_mat(:,:)
1837 
1838 !Local variables-------------------------------
1839  integer :: ii,jj,iglob,jglob
1840 
1841 !************************************************************************
1842 
1843  select case (uplo(1:1))
1844 
1845    case ("A","a")  ! Full global matrix has to be filled.
1846      do jj=1,Slk_mat%sizeb_local(2)
1847        do ii=1,Slk_mat%sizeb_local(1)
1848          call idx_glob(Slk_mat,ii,jj,iglob,jglob)
1849          glob_mat(iglob,jglob) = Slk_mat%buffer_cplx(ii,jj)
1850        end do
1851      end do
1852 
1853    case ("U","u") ! Only the upper triangle of the global matrix is filled.
1854      do jj=1,Slk_mat%sizeb_local(2)
1855        do ii=1,Slk_mat%sizeb_local(1)
1856          call idx_glob(Slk_mat,ii,jj,iglob,jglob)
1857          if (jglob>=iglob) glob_mat(iglob,jglob) = Slk_mat%buffer_cplx(ii,jj)
1858        end do
1859      end do
1860 
1861    case ("L","l") ! Only the lower triangle of the global matrix is filled.
1862      do jj=1,Slk_mat%sizeb_local(2)
1863        do ii=1,Slk_mat%sizeb_local(1)
1864          call idx_glob(Slk_mat,ii,jj,iglob,jglob)
1865          if (jglob<=iglob) glob_mat(iglob,jglob) = Slk_mat%buffer_cplx(ii,jj)
1866        end do
1867      end do
1868 
1869      case default
1870      MSG_BUG(" Wrong uplo: "//TRIM(uplo))
1871  end select
1872 
1873 end subroutine slk_matrix_to_global_dpc_2D

m_slk/slk_my_rclist [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_my_rclist

FUNCTION

  Returns a list with the (row|column) indices of the global matrix that are treated by this node.

INPUTS

  Slk_mat<matrix_scalapack>=Structured datatype defining the scaLAPACK distribution.
  rc_str= "C" if the list of columns is wanted. "R" for rows.

OUTPUT

  how_many=Number of (rows|columns) treated by this node.

SIDE EFFECTS

  rc_list
    input: pointer to null
    output: rc_list(ii=1,how_many) gives the global indices treated by this node.

TODO

  Likely there's a much faster way to retrieve the list of indices using scaLAPACK primitives.

PARENTS

CHILDREN

SOURCE

5061 subroutine slk_my_rclist(Slk_mat,rc_str,how_many,rc_list)
5062 
5063 
5064 !This section has been created automatically by the script Abilint (TD).
5065 !Do not modify the following lines by hand.
5066 #undef ABI_FUNC
5067 #define ABI_FUNC 'slk_my_rclist'
5068 !End of the abilint section
5069 
5070  implicit none
5071 
5072 !Arguments ------------------------------------
5073 !scalars
5074  integer,intent(out) :: how_many
5075  character(len=*),intent(in) :: rc_str
5076  type(matrix_scalapack),intent(in) :: Slk_mat
5077 !scalars
5078  integer,pointer :: rc_list(:)
5079 
5080 !Local variables ------------------------------
5081 !scalars
5082  integer :: col_loc,row_loc,row_glob,col_glob
5083  integer :: nseen,nrow_glob,ncol_glob
5084  !character(len=500) :: msg
5085 ! arrays
5086  integer,allocatable :: seen(:)
5087 
5088 !************************************************************************
5089 
5090 !@matrix_scalapack
5091  how_many=0
5092 
5093  nrow_glob = Slk_mat%sizeb_global(1)
5094  ncol_glob = Slk_mat%sizeb_global(2)
5095 
5096  select case (toupper(rc_str(1:1)))
5097 
5098  case ("C")
5099    ABI_MALLOC(seen,(ncol_glob))
5100    nseen=0
5101    !
5102    do col_loc=1,Slk_mat%sizeb_local(2)
5103      do row_loc=1,Slk_mat%sizeb_local(1)
5104        call idx_glob(Slk_mat,row_loc,col_loc,row_glob,col_glob)
5105        !
5106        if (col_glob==1.and.row_loc==1) then
5107          nseen = nseen+1
5108          seen(nseen) = col_glob
5109        else
5110         if ( ALL(col_glob /= seen(1:nseen)) ) then
5111          nseen = nseen+1
5112          seen(nseen) = col_glob
5113         end if
5114        end if
5115      end do
5116    end do
5117 
5118    how_many = nseen
5119    ABI_MALLOC(rc_list,(nseen))
5120    rc_list = seen(1:nseen)
5121    ABI_FREE(seen)
5122 
5123  case ("R")
5124    ABI_MALLOC(seen,(nrow_glob))
5125    nseen=0
5126    !
5127    do col_loc=1,Slk_mat%sizeb_local(2)
5128      do row_loc=1,Slk_mat%sizeb_local(1)
5129        call idx_glob(Slk_mat,row_loc,col_loc,row_glob,col_glob)
5130        !
5131        if (col_glob==1.and.row_loc==1) then
5132          nseen = nseen+1
5133          seen(nseen) = row_glob
5134        else
5135         if ( ALL(row_glob /= seen(1:nseen)) ) then
5136          nseen = nseen+1
5137          seen(nseen) = row_glob
5138         end if
5139        end if
5140      end do
5141    end do
5142 
5143    how_many = nseen
5144    ABI_MALLOC(rc_list,(nseen))
5145    rc_list = seen(1:nseen)
5146    ABI_FREE(seen)
5147 
5148  case default
5149    MSG_ERROR(" Wrong rc_str: "//TRIM(rc_str))
5150  end select
5151 
5152 end subroutine slk_my_rclist

m_slk/slk_pzgemm [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_pzgemm

FUNCTION

  Extended matrix*matrix product
  C := alpha*A*B - beta*C

  For a simple matrix vector product, one can simply pass
  alpha = (1.,0.) and beta (0.,0.)

INPUTS

  matrix1= first ScaLAPACK matrix (matrix A)
  matrix2= second ScaLAPACK matrix (matrix B)
  alpha= scalar multiplicator for the A*B product
  beta= scalar multiplicator for the C matrix

OUTPUT

  None

NOTES

 The matrices matrix1 and matrix2 must have no common elements; otherwise, results are unpredictable.

SIDE EFFECTS

  results= ScaLAPACK matrix coming out of the operation

PARENTS

      m_exc_diago

CHILDREN

SOURCE

2041 subroutine slk_pzgemm(transa,transb,matrix1,alpha,matrix2,beta,results)
2042 
2043 
2044 !This section has been created automatically by the script Abilint (TD).
2045 !Do not modify the following lines by hand.
2046 #undef ABI_FUNC
2047 #define ABI_FUNC 'slk_pzgemm'
2048 !End of the abilint section
2049 
2050  implicit none
2051 
2052 !Arguments ------------------------------------
2053  character(len=*),intent(in) :: transa,transb
2054  type(matrix_scalapack),intent(in) :: matrix1,matrix2
2055  type(matrix_scalapack),intent(inout) :: results
2056  complex(dpc),intent(in) :: alpha, beta
2057 
2058 !************************************************************************
2059 
2060  call PZGEMM(transa,transb,matrix1%sizeb_global(1),matrix2%sizeb_global(2),&
2061 &  matrix1%sizeb_global(2),alpha,matrix1%buffer_cplx,1,1,&
2062 &  matrix1%descript%tab,matrix2%buffer_cplx,1,1,         &
2063 &  matrix2%descript%tab,beta,results%buffer_cplx,1,1,    &
2064 &  results%descript%tab)
2065 
2066 end subroutine slk_pzgemm

m_slk/slk_pzheev [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

 slk_pzheev

FUNCTION

  slk_pzheev provides an object-oriented interface to the ScaLAPACK routine PHZEEV which computes selected
  eigenvalues and, optionally, eigenvectors of an Hermitian matrix A.
   A * X = lambda * X

INPUTS

  JOBZ    (global input) CHARACTER*1
          Specifies whether or not to compute the eigenvectors:
          = "N":  Compute eigenvalues only.
          = "V":  Compute eigenvalues and eigenvectors.
  UPLO    (global input) CHARACTER*1
          Specifies whether the upper or lower triangular part of the symmetric matrix A is stored:
          = "U":  Upper triangular
          = "L":  Lower triangular

  Slk_mat<type(matrix_scalapack)>=The object storing the local buffer, the array descriptor, the context
    and other quantities needed to call ScaLAPACK routines.
  Slk_vec<matrix_scalapack>=The distributed eigenvectors. Not referenced if JOBZ="N"

OUTPUT

  W       (global output) DOUBLE PRECISION array, dimension (N) where N is the rank of the global matrix.
          On normal exit, the first M entries contain the selected eigenvalues in ascending order.

SIDE EFFECTS

  If JOBZ="V", the local buffer Slk_vec%buffer_cplx will contain part of the distributed eigenvectors.

PARENTS

      m_abilasi,m_exc_diago

CHILDREN

SOURCE

2970 subroutine slk_pzheev(jobz,uplo,Slk_mat,Slk_vec,w)
2971 
2972 
2973 !This section has been created automatically by the script Abilint (TD).
2974 !Do not modify the following lines by hand.
2975 #undef ABI_FUNC
2976 #define ABI_FUNC 'slk_pzheev'
2977 !End of the abilint section
2978 
2979  implicit none
2980 
2981 !Arguments ------------------------------------
2982 !scalars
2983  character(len=*),intent(in) :: jobz,uplo
2984  type(matrix_scalapack),intent(inout) :: Slk_mat
2985  type(matrix_scalapack),intent(inout) :: Slk_vec
2986 !arrays
2987  real(dp),intent(out) :: w(:)
2988 
2989 !Local variables ------------------------------
2990 !scalars
2991  integer :: lwork,lrwork,info,nn
2992  character(len=500) :: msg
2993 !arrays
2994  real(dp),allocatable :: rwork(:)
2995  complex(dpc),allocatable :: work(:)
2996 
2997 !************************************************************************
2998 
2999  ABI_CHECK(allocated(Slk_mat%buffer_cplx),"buffer_cplx not allocated")
3000 
3001 !Get optimal size of workspace.
3002  lwork=-1; lrwork=-1
3003  ABI_MALLOC(work,(1))
3004  ABI_MALLOC(rwork,(1))
3005 
3006  call PZHEEV(jobz,uplo,Slk_mat%sizeb_global(2),Slk_mat%buffer_cplx,1,1,Slk_mat%descript%tab, &
3007 & w,Slk_vec%buffer_cplx,1,1,Slk_vec%descript%tab,work,lwork,rwork,lrwork,info)
3008 
3009  ABI_CHECK(info==0,"Error during the calculation of the workspace size")
3010 
3011  lwork = NINT(real(work(1))); lrwork= NINT(rwork(1)) !*2
3012  ABI_FREE(work)
3013  ABI_FREE(rwork)
3014  !
3015  ! MG: Nov 23 2011. On my mac with the official scalapack package, rwork(1) is not large enough and causes a SIGFAULT.
3016  nn = Slk_mat%sizeb_global(1)
3017  if (firstchar(jobz,(/'V'/))) then
3018   if (lrwork < 2*nn + 2*nn-2) lrwork = 2*nn + 2*nn-2
3019  else if (firstchar(jobz,(/'N'/))) then
3020   if (lrwork < 2*nn) lrwork = 2*nn
3021  end if
3022  !write(std_out,*)lwork,lrwork
3023 
3024  !Solve the problem.
3025  ABI_MALLOC(work,(lwork))
3026  ABI_MALLOC(rwork,(lrwork))
3027 
3028  call PZHEEV(jobz,uplo,Slk_mat%sizeb_global(2),Slk_mat%buffer_cplx,1,1,Slk_mat%descript%tab, &
3029 & w,Slk_vec%buffer_cplx,1,1,Slk_vec%descript%tab,work,lwork,rwork,lrwork,info)
3030 
3031  if (info/=0) then
3032    write(msg,'(a,i0)')"PZHEEV returned info= ",info
3033    MSG_ERROR(msg)
3034  end if
3035 
3036  ABI_FREE(work)
3037  ABI_FREE(rwork)
3038 
3039 end subroutine slk_pzheev

m_slk/slk_pzheevx [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_pzheevx

FUNCTION

  slk_pzheevx provides an object-oriented interface to the ScaLAPACK routine PZHEEVX that
  computes selected eigenvalues and, optionally, eigenvectors of a complex hermitian matrix A.
   A * X = lambda *  X

INPUTS

  Slk_mat<matrix_scalapack>=ScaLAPACK matrix (matrix A)

  Slk_vec<matrix_scalapack>=The distributed eigenvectors X. Not referenced if JOBZ="N"

  JOBZ    (global input) CHARACTER*1
          Specifies whether or not to compute the eigenvectors:
          = "N":  Compute eigenvalues only.
          = "V":  Compute eigenvalues and eigenvectors.

  RANGE   (global input) CHARACTER*1
          = "A": all eigenvalues will be found.
          = "V": all eigenvalues in the interval [VL,VU] will be found.
          = "I": the IL-th through IU-th eigenvalues will be found.

  UPLO    (global input) CHARACTER*1
          Specifies whether the upper or lower triangular part of the Hermitian matrix A is stored:
          = "U":  Upper triangular
          = "L":  Lower triangular

  VL      (global input) DOUBLE PRECISION
          If  RANGE="V",the lower bound of the interval to be searched for eigenvalues. Not referenced if RANGE =
          "A" or "I"

  VU      (global input) DOUBLE PRECISION
          If RANGE="V", the upper bound of the interval to be searched for eigenvalues.  Not referenced  if  RANGE  =
          "A" or "I".

  IL     (global input) integer
         If  RANGE="I",  the  index  (from smallest to largest) of the smallest eigenvalue to be returned.  IL >= 1.
         Not referenced if RANGE = "A" or "V".

  IU     (global input) integer
         If RANGE="I", the index (from smallest to largest) of the largest eigenvalue to be returned.  min(IL,N)  <=
         IU <= N.  Not referenced if RANGE = "A" or "V"

  ABSTOL  (global input) DOUBLE PRECISION
          If JOBZ="V", setting ABSTOL to PDLAMCH( CONTEXT, "U") yields the most orthogonal eigenvectors.
          The  absolute error tolerance for the eigenvalues.  An approximate eigenvalue is accepted as converged when
          it is determined to lie in an interval [a,b] of width less than or equal to

           ABSTOL + EPS *   max( |a|,|b| ) ,

          where EPS is the machine precision.  If ABSTOL is less than or equal to zero, then EPS*norm(T) will be used
          in  its  place, where norm(T) is the 1-norm of the tridiagonal matrix obtained by reducing A to tridiagonal form.
          Eigenvalues will be computed  most  accurately  when  ABSTOL  is  set  to  twice  the  underflow  threshold
          2*PDLAMCH("S")  not  zero.   If  this routine returns with ((MOD(INFO,2).NE.0) .OR.  (MOD(INFO/8,2).NE.0)),
          indicating that some eigenvalues or eigenvectors did not converge, try setting ABSTOL to 2*PDLAMCH("S").

OUTPUT

  mene_found= (global output) Total number of eigenvalues found.  0 <= mene_found <= N.

  eigen(N)= (global output) Eigenvalues of A where N is the dimension of M
            On normal exit, the first mene_found entries contain the selected eigenvalues in ascending order.

SIDE EFFECTS

  If JOBZ="V", the local buffer Slk_vec%buffer_cplx will contain part of the distributed eigenvectors.

  Slk_mat<ScaLAPACK_matrix>=
    %buffer_cplx is destroyed when the routine returns

PARENTS

      m_abilasi,m_exc_diago

CHILDREN

SOURCE

3121 subroutine slk_pzheevx(jobz,range,uplo,Slk_mat,vl,vu,il,iu,abstol,Slk_vec,mene_found,eigen)
3122 
3123 
3124 !This section has been created automatically by the script Abilint (TD).
3125 !Do not modify the following lines by hand.
3126 #undef ABI_FUNC
3127 #define ABI_FUNC 'slk_pzheevx'
3128 !End of the abilint section
3129 
3130  implicit none
3131 
3132 !Arguments ------------------------------------
3133  integer,intent(in) :: il,iu
3134  integer,intent(out) :: mene_found
3135  real(dp),intent(in) :: abstol,vl,vu
3136  character(len=*),intent(in) :: jobz,range,uplo
3137  type(matrix_scalapack),intent(inout) :: Slk_mat
3138  type(matrix_scalapack),intent(inout) :: Slk_vec
3139 !arrays
3140  real(dp),intent(out) :: eigen(*)
3141 
3142 !Local variables-------------------------------
3143 !scalars
3144  integer  :: lwork,lrwork,liwork,info
3145  integer ::  nvec_calc,ierr
3146  real(dp) :: orfac
3147  character(len=500) :: msg
3148 !arrays
3149  integer :: ibuff(3),max_ibuff(3)
3150  integer,allocatable  :: iwork(:),iclustr(:),ifail(:)
3151  real(dp),allocatable  :: rwork(:),gap(:)
3152  complex(dpc),allocatable :: work(:)
3153 
3154 !************************************************************************
3155 
3156 ! abstol = PDLAMCH(Slk_vecprocessor%grid%ictxt,'U')
3157 
3158   orfac  = -one ! Only for eigenvectors: use default value 10d-3.
3159 ! Vectors within orfac*norm(A) will be reorthogonalized.
3160 
3161 ! Allocate the arrays for the results of the calculation
3162   ABI_MALLOC(gap,( Slk_mat%processor%grid%dims(1) * Slk_mat%processor%grid%dims(2)))
3163 
3164   if (firstchar(jobz,(/"V","v"/))) then
3165     ABI_MALLOC(ifail,(Slk_mat%sizeb_global(2)))
3166     ABI_MALLOC(iclustr,( 2*Slk_mat%processor%grid%dims(1) * Slk_mat%processor%grid%dims(2)))
3167   end if
3168 !
3169 ! Get the optimal size of the work arrays.
3170   lwork=-1; lrwork=-1; liwork=-1
3171   ABI_MALLOC(work,(1))
3172   ABI_MALLOC(iwork,(1))
3173   ABI_MALLOC(rwork,(3))
3174 ! This is clearly seen in the source in which rwork(1:3) is accessed
3175 ! during the calcuation of the workspace size.
3176 
3177   call PZHEEVX(jobz,range,uplo, Slk_mat%sizeb_global(2),Slk_mat%buffer_cplx,1,1,Slk_mat%descript%tab,&
3178 &  vl,vu,il,iu,abstol,mene_found,nvec_calc,eigen,orfac,&
3179 &  Slk_vec%buffer_cplx,1,1,Slk_vec%descript%tab,&
3180 &  work,lwork,rwork,lrwork,iwork,liwork,ifail,iclustr,gap,info)
3181 
3182   write(msg,'(a,i7)')" Problem to compute workspace to use ScaLAPACK, info: ",info
3183   ABI_CHECK(info==0,msg)
3184 
3185   lwork  = NINT(real(work(1)),kind=dp)
3186   lrwork = NINT(rwork(1))
3187   liwork = iwork(1)
3188 
3189   ABI_FREE(work)
3190   ABI_FREE(rwork)
3191   ABI_FREE(iwork)
3192 !
3193 ! FROM THE SCALAPACK MAN PAGE:
3194 ! The computed eigenvectors may not be orthogonal if the minimal workspace is supplied and ORFAC is too
3195 ! small. If you  want to guarantee orthogonality (at the cost of potentially poor performance) you should
3196 ! add the following to LRWORK: (CLUSTERSIZE-1)*N where CLUSTERSIZE is  the  number  of  eigenvalues  in  the
3197 ! largest cluster, where a cluster is defined as a set of close eigenvalues: { W(K),...,W(K+CLUSTERSIZE-1) |
3198 ! W(J+1) <= W(J) + ORFAC*2*norm(A) }.
3199 
3200   if ( firstchar(jobz,(/"V","v"/)) ) then
3201     lrwork = INT( lrwork + Slk_mat%sizeb_global(2) *(Slk_mat%sizeb_global(2)-1) )
3202   end if
3203 
3204 ! ibuff(1) = lwork
3205 ! ibuff(2) = lrwork !INT(lrwork + Slk_mat%sizeb_global(2) *(Slk_mat%sizeb_global(2)-1)
3206 ! ibuff(3) = liwork
3207 
3208 ! Get the maximum of sizes of the work arrays processor%comm
3209 ! call MPI_ALLREDUCE(ibuff,max_ibuff,3,MPI_integer,MPI_MAX,comm,ierr)
3210 
3211 ! lwork  = max_ibuff(1)
3212 ! lrwork = max_ibuff(2)
3213 ! liwork = max_ibuff(3)
3214 
3215   ABI_MALLOC(work ,(lwork ))
3216   ABI_MALLOC(rwork,(lrwork))
3217   ABI_MALLOC(iwork,(liwork))
3218 !
3219 ! Call the scaLAPACK routine.
3220 
3221 ! write(std_out,*) 'I am using PZHEEVX'
3222   call PZHEEVX(jobz,range,uplo, Slk_mat%sizeb_global(2),Slk_mat%buffer_cplx,1,1,Slk_mat%descript%tab,&
3223 &  vl,vu,il,iu,abstol,mene_found,nvec_calc, eigen,orfac,&
3224 &  Slk_vec%buffer_cplx,1,1,Slk_vec%descript%tab,&
3225 &  work,lwork,rwork,lrwork,iwork,liwork,ifail,iclustr,gap,info)
3226 
3227 ! Handle the possible error.
3228   if (info < 0) then
3229     write(msg,'(a,i7,a)')" The ",-info,"-th argument of PZHEEVX had an illegal value."
3230     if (info==-25) msg = " LRWORK is too small to compute all the eigenvectors requested, no computation is performed"
3231     MSG_ERROR(msg)
3232   end if
3233 
3234   if (info > 0) then
3235     write(msg,'(a,i7)') " PZHEEVX returned info: ",info
3236     call wrtout(std_out,msg,"PERS")
3237     if (MOD(info,2)/=0)then
3238       write(msg,'(3a)')&
3239 &      " One or more eigenvectors failed to converge. ",ch10,&
3240 &      " Their indices are stored in IFAIL. Ensure ABSTOL=2.0*PDLAMCH('U')"
3241       call wrtout(std_out,msg,"PERS")
3242     end if
3243     if (MOD(info/2,2)/=0) then
3244       write(msg,'(5a)')&
3245 &      " Eigenvectors corresponding to one or more clusters of eigenvalues ",ch10,&
3246 &      " could not be reorthogonalized because of insufficient workspace. ",ch10,&
3247 &      " The indices of the clusters are stored in the array ICLUSTR."
3248       call wrtout(std_out,msg,"PERS")
3249     end if
3250     if (MOD(info/4,2)/=0) then
3251       write(msg,'(3a)')" Space limit prevented PZHEEVX from computing all of the eigenvectors between VL and VU. ",ch10,&
3252 &      " The number of eigenvectors  computed  is returned in NZ."
3253       call wrtout(std_out,msg,"PERS")
3254     end if
3255     if (MOD(info/8,2)/=0) then
3256       msg = " PZSTEBZ  failed to compute eigenvalues. Ensure ABSTOL=2.0*PDLAMCH('U')"
3257       call wrtout(std_out,msg,"PERS")
3258     end if
3259     MSG_ERROR("Cannot continue")
3260   end if
3261 
3262 ! Check the number of eigenvalues found wrt to the number of vectors calculated.
3263   if ( firstchar(jobz,(/'V','v'/)) .and. mene_found/=nvec_calc) then
3264     write(msg,'(5a)')&
3265 &    " The user supplied insufficient space and PZHEEVX is not able to detect this before beginning computation. ",ch10,&
3266 &    " To get all the  eigenvectors requested, the user must supply both sufficient space to hold the ",ch10,&
3267 &    " eigenvectors in Z (M .LE. DESCZ(N_)) and sufficient workspace to compute them. "
3268     MSG_ERROR(msg)
3269   end if
3270 
3271   ABI_FREE(work)
3272   ABI_FREE(rwork)
3273   ABI_FREE(iwork)
3274   ABI_FREE(gap)
3275 
3276   if ( firstchar(jobz,(/"V","v"/)) ) then
3277     ABI_FREE(ifail)
3278     ABI_FREE(iclustr)
3279   end if
3280 
3281 end subroutine slk_pzheevx

m_slk/slk_pzhegvx [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_pzhegvx

FUNCTION

  slk_pzhegvx provides an object-oriented interface to the ScaLAPACK routine PZHEGVX that
  computes selected eigenvalues and, optionally, eigenvectors of a complex generalized
  Hermitian-definite eigenproblem, of the form
  sub( A )*x=(lambda)*sub( B )*x,  sub( A )*sub( B )x=(lambda)*x,  or sub( B )*sub( A )*x=(lambda)*x.
  Here sub( A ) denoting A( IA:IA+N-1, JA:JA+N-1 ) is assumed to be
  Hermitian, and sub( B ) denoting B( IB:IB+N-1, JB:JB+N-1 ) is assumed
  to be Hermitian positive definite.

INPUTS

  Slk_matA<matrix_scalapack>=ScaLAPACK matrix (matrix A)
  Slk_matB<matrix_scalapack>=ScaLAPACK matrix (matrix B)
  Slk_vec<matrix_scalapack>=The distributed eigenvectors X. Not referenced if JOBZ="N"

  IBtype   (global input) integer
          Specifies the problem type to be solved:
          = 1:  sub( A )*x = (lambda)*sub( B )*x
          = 2:  sub( A )*sub( B )*x = (lambda)*x
          = 3:  sub( B )*sub( A )*x = (lambda)*x

  JOBZ    (global input) CHARACTER*1
          Specifies whether or not to compute the eigenvectors:
          = "N":  Compute eigenvalues only.
          = "V":  Compute eigenvalues and eigenvectors.

  RANGE   (global input) CHARACTER*1
          = "A": all eigenvalues will be found.
          = "V": all eigenvalues in the interval [VL,VU] will be found.
          = "I": the IL-th through IU-th eigenvalues will be found.

  UPLO    (global input) CHARACTER*1
          Specifies whether the upper or lower triangular part of the Hermitian matrix sub(A) and sub(B) is stored:
          = "U":  Upper triangular
          = "L":  Lower triangular

  VL      (global input) DOUBLE PRECISION
          If  RANGE="V",the lower bound of the interval to be searched for eigenvalues. Not referenced if RANGE =
          "A" or "I"

  VU      (global input) DOUBLE PRECISION
          If RANGE="V", the upper bound of the interval to be searched for eigenvalues.  Not referenced  if  RANGE  =
          "A" or "I".

  IL     (global input) integer
         If  RANGE="I",  the  index  (from smallest to largest) of the smallest eigenvalue to be returned.  IL >= 1.
         Not referenced if RANGE = "A" or "V".

  IU     (global input) integer
         If RANGE="I", the index (from smallest to largest) of the largest eigenvalue to be returned.  min(IL,N)  <=
         IU <= N.  Not referenced if RANGE = "A" or "V"

  ABSTOL  (global input) DOUBLE PRECISION
          If JOBZ="V", setting ABSTOL to PDLAMCH( CONTEXT, "U") yields the most orthogonal eigenvectors.
          The  absolute error tolerance for the eigenvalues.  An approximate eigenvalue is accepted as converged when
          it is determined to lie in an interval [a,b] of width less than or equal to

           ABSTOL + EPS *   max( |a|,|b| ) ,

          where EPS is the machine precision.  If ABSTOL is less than or equal to zero, then EPS*norm(T) will be used
          in  its  place, where norm(T) is the 1-norm of the tridiagonal matrix obtained by reducing A to tridiagonal form.
          Eigenvalues will be computed  most  accurately  when  ABSTOL  is  set  to  twice  the  underflow  threshold
          2*PDLAMCH("S")  not  zero.   If  this routine returns with ((MOD(INFO,2).NE.0) .OR.  (MOD(INFO/8,2).NE.0)),
          indicating that some eigenvalues or eigenvectors did not converge, try setting ABSTOL to 2*PDLAMCH("S").

OUTPUT

  mene_found= (global output) Total number of eigenvalues found.  0 <= mene_found <= N.

  eigen(N)= (global output) Eigenvalues of A where N is the dimension of M
            On normal exit, the first mene_found entries contain the selected eigenvalues in ascending order.

SIDE EFFECTS

  Slk_vec<matrix_scalapack>:
   %buffer_cplx local output (global dimension (N,N)
     If JOBZ = 'V', then on normal exit the first M columns of Z
     contain the orthonormal eigenvectors of the matrix
     corresponding to the selected eigenvalues.
     If JOBZ = 'N', then Z is not referenced.

  Slk_matA<matrix_scalapack>:
    %buffer_cplx
      (local input/local output) complex(DPC) pointer into the
      local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
      On entry, this array contains the local pieces of the
      N-by-N Hermitian distributed matrix sub( A ). If UPLO = 'U',
      the leading N-by-N upper triangular part of sub( A ) contains
      the upper triangular part of the matrix.  If UPLO = 'L', the
      leading N-by-N lower triangular part of sub( A ) contains
      the lower triangular part of the matrix.

      On exit, if JOBZ = 'V', then if INFO = 0, sub( A ) contains
      the distributed matrix Z of eigenvectors.  The eigenvectors
      are normalized as follows:
      if IBtype = 1 or 2, Z**H*sub( B )*Z = I;
      if IBtype = 3, Z**H*inv( sub( B ) )*Z = I.
      If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
      or the lower triangle (if UPLO='L') of sub( A ), including
      the diagonal, is destroyed.

  Slk_matB=
    %buffer_cplx
      (local input/local output) complex*(DPC) pointer into the
      local memory to an array of dimension (LLD_B, LOCc(JB+N-1)).
      On entry, this array contains the local pieces of the
      N-by-N Hermitian distributed matrix sub( B ). If UPLO = 'U',
      the leading N-by-N upper triangular part of sub( B ) contains
      the upper triangular part of the matrix.  If UPLO = 'L', the
      leading N-by-N lower triangular part of sub( B ) contains
      the lower triangular part of the matrix.

      On exit, if INFO <= N, the part of sub( B ) containing the
      matrix is overwritten by the triangular factor U or L from
      the Cholesky factorization sub( B ) = U**H*U or
      sub( B ) = L*L**H.

PARENTS

      m_exc_diago

CHILDREN

SOURCE

3411 subroutine slk_pzhegvx(ibtype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,mene_found,eigen)
3412 
3413 
3414 !This section has been created automatically by the script Abilint (TD).
3415 !Do not modify the following lines by hand.
3416 #undef ABI_FUNC
3417 #define ABI_FUNC 'slk_pzhegvx'
3418 !End of the abilint section
3419 
3420  implicit none
3421 
3422 !Arguments ------------------------------------
3423  integer,intent(in) :: il,iu,ibtype
3424  integer,intent(out) :: mene_found
3425  real(dp),intent(in) :: abstol,vl,vu
3426  character(len=*),intent(in) :: jobz,range,uplo
3427  type(matrix_scalapack),intent(inout) :: Slk_matA
3428  type(matrix_scalapack),intent(inout) :: Slk_matB
3429  type(matrix_scalapack),intent(inout) :: Slk_vec
3430 !arrays
3431  real(dp),intent(out) :: eigen(*)
3432 
3433 !Local variables-------------------------------
3434 !scalars
3435  integer  :: lwork,lrwork,liwork,info,nvec_calc,ierr
3436  real(dp) :: orfac
3437  logical :: ltest
3438  character(len=500) :: msg
3439 !arrays
3440  integer :: ibuff(3),max_ibuff(3)
3441  integer :: desca(DLEN_),descb(DLEN_),descz(DLEN_)
3442  integer,allocatable  :: iwork(:),iclustr(:),ifail(:)
3443  real(dp),allocatable  :: rwork(:),gap(:)
3444  complex(dpc),allocatable :: work(:)
3445 
3446 !************************************************************************
3447 
3448 ! abstol = PDLAMCH(Slk_vecprocessor%grid%ictxt,'U')
3449 
3450  orfac  = -one ! Only for eigenvectors: use default value 10d-3.
3451 ! Vectors within orfac*norm(A) will be reorthogonalized.
3452 
3453 ! ======================
3454 ! Alignment requirements
3455 ! ======================
3456 ! The distributed submatrices A(IA:*, JA:*), C(IC:IC+M-1,JC:JC+N-1),
3457 ! and B( IB:IB+N-1, JB:JB+N-1 ) must verify some alignment properties,
3458 
3459  desca = Slk_matA%descript%tab
3460  descb = Slk_matB%descript%tab
3461  if (firstchar(jobz,(/"V","v"/))) then
3462    descz = Slk_vec%descript%tab
3463  else
3464    descz = Slk_matA%descript%tab
3465  end if
3466 
3467  ltest = .TRUE.
3468  ltest = ltest .and. (DESCA(MB_) == DESCA(NB_))
3469 !IA = IB = IZ
3470 !JA = IB = JZ
3471  ltest = ltest .and.ALL(DESCA(M_   )==(/DESCB(M_   ),DESCZ(M_   )/))
3472  ltest = ltest .and.ALL(DESCA(N_   )==(/DESCB(N_   ),DESCZ(N_   )/))
3473  ltest = ltest .and.ALL(DESCA(MB_  )==(/DESCB(MB_  ),DESCZ(MB_  )/))
3474  ltest = ltest .and.ALL(DESCA(NB_  )==(/DESCB(NB_  ),DESCZ(NB_  )/))
3475  ltest = ltest .and.ALL(DESCA(RSRC_)==(/DESCB(RSRC_),DESCZ(RSRC_)/))
3476  ltest = ltest .and.ALL(DESCA(CSRC_)==(/DESCB(CSRC_),DESCZ(CSRC_)/))
3477 !MOD( IA-1, DESCA( MB_ ) ) = 0
3478 !MOD( JA-1, DESCA( NB_ ) ) = 0
3479 !MOD( IB-1, DESCB( MB_ ) ) = 0
3480 !MOD( JB-1, DESCB( NB_ ) ) = 0
3481 
3482  if (.not.ltest) then
3483    msg = " Alignment requirements not satisfied, check the caller"
3484    MSG_ERROR(msg)
3485  end if
3486 
3487 !Allocate the arrays for the results of the calculation
3488  ABI_MALLOC(gap,( Slk_matA%processor%grid%dims(1) * Slk_matA%processor%grid%dims(2)))
3489 
3490  if (firstchar(jobz,(/"V","v"/))) then
3491    ABI_MALLOC(ifail,(Slk_matA%sizeb_global(2)))
3492    ABI_MALLOC(iclustr,( 2*Slk_matA%processor%grid%dims(1) * Slk_matA%processor%grid%dims(2)))
3493  else
3494    ABI_MALLOC(ifail,(1))
3495  end if
3496 !
3497 !Get the optimal size of the work arrays.
3498  lwork=-1; lrwork=-1; liwork=-1
3499  ABI_MALLOC(work,(1))
3500  ABI_MALLOC(iwork,(1))
3501  ABI_MALLOC(rwork,(3))
3502 !This is clearly seen in the source in which rwork(1:3) is accessed
3503 !during the calcuation of the workspace size.
3504 
3505  call pzhegvx(ibtype,jobz,range,uplo, Slk_matA%sizeb_global(2),Slk_matA%buffer_cplx,1,1,Slk_matA%descript%tab,&
3506 & Slk_matB%buffer_cplx,1,1,Slk_matB%descript%tab,&
3507 & vl,vu,il,iu,abstol,mene_found,nvec_calc,eigen,orfac,&
3508 & Slk_vec%buffer_cplx,1,1,Slk_vec%descript%tab,&
3509 & work,lwork,rwork,lrwork,iwork,liwork,ifail,iclustr,gap,info)
3510 
3511  write(msg,'(a,i0)')" Problem to compute workspace to use ScaLAPACK, info: ",info
3512  ABI_CHECK(info==0,msg)
3513 
3514  lwork  = NINT(real(work(1)),kind=dp)
3515  lrwork = NINT(rwork(1))
3516  liwork = iwork(1)
3517 
3518  ABI_FREE(work)
3519  ABI_FREE(rwork)
3520  ABI_FREE(iwork)
3521 !
3522 !FROM THE SCALAPACK MAN PAGE:
3523 !The computed eigenvectors may not be orthogonal if the minimal workspace is supplied and ORFAC is too
3524 !small. If you  want to guarantee orthogonality (at the cost of potentially poor performance) you should
3525 !add the following to LRWORK: (CLUSTERSIZE-1)*N where CLUSTERSIZE is  the  number  of  eigenvalues  in  the
3526 !largest cluster, where a cluster is defined as a set of close eigenvalues: { W(K),...,W(K+CLUSTERSIZE-1) |
3527 !W(J+1) <= W(J) + ORFAC*2*norm(A) }.
3528 
3529  if ( firstchar(jobz,(/"V","v"/)) ) then
3530    lrwork = INT( lrwork + Slk_matA%sizeb_global(2) *(Slk_matA%sizeb_global(2)-1) )
3531  end if
3532 
3533 !ibuff(1) = lwork
3534 !ibuff(2) = lrwork !INT(lrwork + Slk_matA%sizeb_global(2) *(Slk_matA%sizeb_global(2)-1)
3535 !ibuff(3) = liwork
3536 
3537 !Get the maximum of sizes of the work arrays processor%comm
3538 !call MPI_ALLREDUCE(ibuff,max_ibuff,3,MPI_integer,MPI_MAX,comm,ierr)
3539 
3540 !lwork  = max_ibuff(1)
3541 !lrwork = max_ibuff(2)
3542 !liwork = max_ibuff(3)
3543 
3544  ABI_MALLOC(work ,(lwork ))
3545  ABI_MALLOC(rwork,(lrwork))
3546  ABI_MALLOC(iwork,(liwork))
3547 !
3548 !Call the scaLAPACK routine.
3549 
3550 !write(std_out,*) 'I am using PZHEGVX'
3551  call pzhegvx(ibtype,jobz,range,uplo, Slk_matA%sizeb_global(2),Slk_matA%buffer_cplx,1,1,Slk_matA%descript%tab,&
3552 & Slk_matB%buffer_cplx,1,1,Slk_matB%descript%tab,&
3553 & vl,vu,il,iu,abstol,mene_found,nvec_calc, eigen,orfac,&
3554 & Slk_vec%buffer_cplx,1,1,Slk_vec%descript%tab,&
3555 & work,lwork,rwork,lrwork,iwork,liwork,ifail,iclustr,gap,info)
3556 
3557 !Handle the possible error.
3558  if (info < 0) then
3559    write(msg,'(a,i7,a)')" The ",-info,"-th argument of PZHEGVX had an illegal value."
3560    if (info==-25) msg = " LRWORK is too small to compute all the eigenvectors requested, no computation is performed"
3561    MSG_ERROR(msg)
3562  end if
3563 
3564  if (info > 0) then
3565    write(msg,'(a,i7)') " PZHEGVX returned info: ",info
3566    call wrtout(std_out,msg,"PERS")
3567    if (MOD(info,2)/=0)then
3568      write(msg,'(3a)')&
3569 &     " One or more eigenvectors failed to converge. ",ch10,&
3570 &     " Their indices are stored in IFAIL. Ensure ABSTOL=2.0*PDLAMCH('U')"
3571      call wrtout(std_out,msg,"PERS")
3572    end if
3573    if (MOD(info/2,2)/=0) then
3574      write(msg,'(5a)')&
3575 &     " Eigenvectors corresponding to one or more clusters of eigenvalues ",ch10,&
3576 &     " could not be reorthogonalized because of insufficient workspace. ",ch10,&
3577 &     " The indices of the clusters are stored in the array ICLUSTR."
3578      call wrtout(std_out,msg,"PERS")
3579    end if
3580    if (MOD(info/4,2)/=0) then
3581      write(msg,'(3a)')&
3582 &     " Space limit prevented PZHEGVX from computing all of the eigenvectors between VL and VU. ",ch10,&
3583 &     " The number of eigenvectors  computed  is returned in NZ."
3584      call wrtout(std_out,msg,"PERS")
3585    end if
3586    if (MOD(info/8,2)/=0) then
3587      msg = " PZSTEBZ  failed to compute eigenvalues. Ensure ABSTOL=2.0*PDLAMCH('U')"
3588      call wrtout(std_out,msg,"PERS")
3589    end if
3590    if (MOD(info/16,2)/=0) then
3591      write(msg,'(3a)')&
3592 &     " B was not positive definite.",ch10,&
3593 &     " IFAIL(1) indicates the order of the smallest minor which is not positive definite."
3594      call wrtout(std_out,msg,"PERS")
3595    end if
3596    MSG_ERROR("Cannot continue")
3597  end if
3598 
3599 !Check the number of eigenvalues found wrt to the number of vectors calculated.
3600  if ( firstchar(jobz,(/'V','v'/)) .and. mene_found/=nvec_calc) then
3601    write(msg,'(5a)')&
3602 &   " The user supplied insufficient space and PZHEGVX is not able to detect this before beginning computation. ",ch10,&
3603 &   " To get all the  eigenvectors requested, the user must supply both sufficient space to hold the ",ch10,&
3604 &   " eigenvectors in Z (M .LE. DESCZ(N_)) and sufficient workspace to compute them. "
3605    MSG_ERROR(msg)
3606  end if
3607 
3608  ABI_FREE(work)
3609  ABI_FREE(rwork)
3610  ABI_FREE(iwork)
3611  ABI_FREE(gap)
3612 
3613  if ( firstchar(jobz,(/"V","v"/)) )  then
3614    ABI_FREE(iclustr)
3615  end if
3616  ABI_FREE(ifail)
3617 
3618 end subroutine slk_pzhegvx

m_slk/slk_read [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_read

FUNCTION

  Routine to read a square scaLAPACK distributed matrix from an external file using MPI-IO.

INPUTS

  uplo=String specifying whether only the upper or lower triangular part of the global matrix is stored on disk:
    = "U":  Upper triangular is stored
    = "L":  Lower triangular is stored
    = "A":  Full matrix (used for general complex matrices)
  symtype=Symmetry type of the matrix stored on disk (used only if uplo = "L" or "A").
    = "H" for Hermitian matrix
    = "S" for symmetric matrix.
    = "N" if matrix has no symmetry (not compatible with uplo="L" or uplo="U".
  is_fortran_file=.FALSE. is C stream is used. .TRUE. for writing Fortran binary files.
  [fname]= Mutually exclusive with mpi_fh. The name of the external file from which the matrix will be read.
           The file is open and closed inside the routine with MPI flags specified by flags.
  [mpi_fh]=File handler associated to the file (already open in the caller). Not compatible with fname.
  [flags]=MPI-IO flags used to open the file in MPI_FILE_OPEN. Default is MPI_MODE_RDONLY. Referenced only when fname is used.

SIDE EFFECTS

  Slk_mat<matrix_scalapack>=Structured datatype defining the scaLAPACK distribution with the local buffer
    supposed to be allocated.
    %buffer_cplx=Local buffer containg the distributed matrix stored on the external file.
  If fname is present then the file is opened and closed inside the routine. Any exception is fatal.
  [offset]=
    input:  Offset used to access the content of the file. Default is zero.
    output: New offset incremented with the byte size of the matrix that has been read (Fortran
            markers are included if is_fortran_file=.TRUE.)

TODO

  * Generalize the implementation adding the reading of the real buffer.

PARENTS

      m_exc_diago

CHILDREN

SOURCE

4068 subroutine slk_read(Slk_mat,uplo,symtype,is_fortran_file,fname,mpi_fh,offset,flags)
4069 
4070 
4071 !This section has been created automatically by the script Abilint (TD).
4072 !Do not modify the following lines by hand.
4073 #undef ABI_FUNC
4074 #define ABI_FUNC 'slk_read'
4075 !End of the abilint section
4076 
4077  implicit none
4078 
4079 !Arguments ------------------------------------
4080 !scalars
4081  integer,optional,intent(in) :: flags,mpi_fh
4082  integer(XMPI_OFFSET_KIND),optional,intent(inout) :: offset
4083  character(len=*),optional,intent(in) :: fname
4084  character(len=*),intent(in) :: uplo,symtype
4085  logical,intent(in) :: is_fortran_file
4086  type(matrix_scalapack),intent(inout) :: Slk_mat
4087 
4088 !Local variables ------------------------------
4089 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
4090 !scalars
4091  integer :: nrows_glob,offset_err,slk_type,etype
4092  integer(XMPI_OFFSET_KIND) :: my_offset
4093  logical :: do_open
4094  integer :: comm,my_flags,my_fh,buffer_size,ierr,col_glob
4095  integer :: nfrec,bsize_elm,mpi_type_elm
4096  complex(dpc) :: ctest
4097  logical,parameter :: check_frm=.TRUE.
4098  integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecord(:)
4099 !arrays
4100  character(len=500) :: msg
4101 #endif
4102 
4103 !************************************************************************
4104 
4105 !@matrix_scalapack
4106 
4107 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
4108 
4109  do_open = PRESENT(fname)
4110  if (do_open) then
4111    ABI_CHECK(.not.PRESENT(fname),"fname should not be present")
4112  else
4113    ABI_CHECK(PRESENT(mpi_fh),"mpi_fh should be present")
4114  end if
4115 
4116  my_offset=0; if (PRESENT(offset)) my_offset=offset
4117 
4118  ABI_CHECK(allocated(Slk_mat%buffer_cplx),"%buffer_cplx not allocated")
4119  if (firstchar(uplo,(/"U","L"/)) .and. Slk_mat%sizeb_global(1) /= Slk_mat%sizeb_global(2) ) then
4120    MSG_ERROR("rectangular matrices are not compatible with the specified uplo")
4121  end if
4122 
4123  nrows_glob = Slk_mat%sizeb_global(1)
4124 
4125  buffer_size= PRODUCT(Slk_mat%sizeb_local(1:2))
4126 
4127  call wrtout(std_out,"slk_read: Using MPI-IO","PERS")
4128 
4129  comm = Slk_mat%processor%comm
4130 
4131  if (do_open) then ! Open the file.
4132    my_flags=MPI_MODE_RDONLY; if (PRESENT(flags)) my_flags = flags
4133    call MPI_FILE_OPEN(comm, fname, my_flags, MPI_INFO_NULL, my_fh, ierr)
4134    ABI_CHECK_MPI(ierr,"FILE_OPEN "//TRIM(fname))
4135  else
4136    my_fh = mpi_fh
4137  end if
4138 
4139  call slk_single_fview_read(Slk_mat,uplo,etype,slk_type,offset_err,is_fortran_file=is_fortran_file)
4140 
4141  if (offset_err/=0) then
4142    write(msg,"(3a)")&
4143 &   "Global position index cannot be stored in standard Fortran integer ",ch10,&
4144 &   "scaLAPACK matrix cannot be read with a single MPI-IO call."
4145    MSG_ERROR(msg)
4146  end if
4147 
4148  call MPI_FILE_SET_VIEW(my_fh, my_offset, etype, slk_type, 'native', MPI_INFO_NULL, ierr)
4149  ABI_CHECK_MPI(ierr,"SET_VIEW")
4150 
4151  call MPI_FILE_READ_ALL(my_fh, Slk_mat%buffer_cplx, buffer_size, MPI_DOUBLE_complex, MPI_STATUS_IGNORE, ierr)
4152  ABI_CHECK_MPI(ierr,"READ_ALL")
4153  !
4154  ! Symmetrize local buffer if uplo /= "All"
4155  call slk_symmetrize(Slk_mat,uplo,symtype)
4156 
4157 !BEGINDEBUG
4158 !call MPI_FILE_READ_AT(mpi_fh,my_offset+xmpio_bsize_frm,ctest,1,MPI_DOUBLE_complex,MPI_STATUS_IGNORE,ierr)
4159 !write(std_out,*)"ctest",ctest
4160 !call MPI_FILE_READ_AT(mpi_fh,my_offset+2*xmpio_bsize_frm,ctest,1,MPI_DOUBLE_complex,MPI_STATUS_IGNORE,ierr)
4161 !write(std_out,*)"ctest",ctest
4162 !ENDDEBUG
4163 
4164  !call print_arr(Slk_mat%buffer_cplx,max_r=10,max_c=10,unit=std_out,mode_paral="PERS")
4165  !
4166  ! Close the file and release the MPI filetype.
4167  call MPI_type_FREE(slk_type,ierr)
4168  ABI_CHECK_MPI(ierr,"MPI_type_FREE")
4169 
4170  call slk_bsize_and_type(Slk_mat,bsize_elm,mpi_type_elm)
4171 
4172 !It seems that personal call makes the code stuck
4173 !if (is_fortran_file .and. check_frm .and. Slk_mat%Processor%myproc==0) then ! Master checks the Fortran markers.
4174  if (is_fortran_file .and. check_frm) then ! Master checks the Fortran markers.
4175    call wrtout(std_out,"Checking Fortran record markers...","PERS", do_flush=.True.)
4176    nfrec = Slk_mat%sizeb_global(2)
4177    ABI_MALLOC(bsize_frecord,(nfrec))
4178    if (firstchar(uplo,(/"A"/))) then
4179      bsize_frecord = Slk_mat%sizeb_global(1) * bsize_elm
4180    else if (firstchar(uplo,(/"U"/))) then
4181      bsize_frecord = (/(col_glob * bsize_elm, col_glob=1,nfrec)/)
4182    else if (firstchar(uplo,(/"L"/))) then
4183      bsize_frecord = (/(col_glob * bsize_elm, col_glob=nfrec,1,-1)/)
4184    else
4185      MSG_ERROR("Wrong uplo")
4186    end if
4187    call xmpio_check_frmarkers(my_fh,my_offset,xmpio_collective,nfrec,bsize_frecord,ierr)
4188    ABI_CHECK(ierr==0,"Wrong Fortran record markers")
4189    ABI_FREE(bsize_frecord)
4190  end if
4191 
4192  if (do_open) then ! Close the file.
4193    call MPI_FILE_CLOSE(my_fh, ierr)
4194    ABI_CHECK_MPI(ierr,"FILE_CLOSE")
4195  end if
4196 !
4197 !Increment the offset
4198  if (PRESENT(offset)) then
4199    if (firstchar(uplo,(/"A"/))) then
4200      offset = offset + PRODUCT(Slk_mat%sizeb_global(1:2)) * bsize_elm
4201      if (is_fortran_file) offset = offset + Slk_mat%sizeb_global(2)*2*xmpio_bsize_frm
4202    else if (firstchar(uplo,(/"U","L"/))) then
4203      offset = offset + ( (Slk_mat%sizeb_global(2) * (Slk_mat%sizeb_global(2))+1)/2 ) * bsize_elm
4204      if (is_fortran_file) offset = offset + Slk_mat%sizeb_global(2)*2*xmpio_bsize_frm
4205    else
4206      MSG_ERROR("Wrong uplo")
4207    end if
4208  end if
4209 
4210  call xmpi_barrier(comm)
4211  RETURN
4212 
4213 #else
4214  MSG_ERROR("MPI-IO support not enabled")
4215 #endif
4216 
4217 end subroutine slk_read

m_slk/slk_single_fview_read [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_single_fview_read

FUNCTION

  Return an MPI datatype that can be used to read a scaLAPACK distributed matrix from
  a binary file using MPI-IO.

INPUTS

  Slk_mat<matrix_scalapack>=Structured datatype defining the scaLAPACK distribution with the local buffer.
  uplo=String specifying whether only the upper or lower triangular part of the global matrix is stored on disk:
    = "U":  Upper triangular is stored
    = "L":  Lower triangular is stored
    = "A":  Full matrix (used for general complex matrices)
  [is_fortran_file]=.FALSE. is C stream is used. .TRUE. for writing Fortran binary files
    with record markers. In this case etype is set xmpio_mpi_type_frm provided that
    the mpi_type of the matrix element is commensurate with xmpio_mpi_type_frm. Defaults to .TRUE.

OUTPUT

  etype=Elementary data type (handle) defining the elementary unit used to access the file.
  slk_type=New MPI type that can be used to instantiate the MPI-IO view for the Fortran file.
     Note that the view assumes that the file pointer points to the FIRST Fortran record marker.
  offset_err=Error code. A non-zero value signals that the global matrix is too large
    for a single MPI-IO access (see notes below).

NOTES

  With (signed) Fortran integers, the maximum size of the file that
  that can be read in one-shot is around 2Gb when etype is set to byte.
  Using a larger etype might create portability problems (real data on machines using
  integer*16 for the marker) since etype must be a multiple of the Fortran record marker
  Due to the above reason, block_displ is given in bytes and must be stored in a integer
  of kind XMPI_ADDRESS_KIND. If the displacement is too large, the routine returns
  offset_err=1 so that the caller will know that several MPI-IO reads are needed to
  read the local buffer.

PARENTS

      m_slk

CHILDREN

SOURCE

4575 subroutine slk_single_fview_read(Slk_mat,uplo,etype,slk_type,offset_err,is_fortran_file)
4576 
4577 
4578 !This section has been created automatically by the script Abilint (TD).
4579 !Do not modify the following lines by hand.
4580 #undef ABI_FUNC
4581 #define ABI_FUNC 'slk_single_fview_read'
4582 !End of the abilint section
4583 
4584  implicit none
4585 
4586 !Arguments ------------------------------------
4587 !scalars
4588  integer,intent(out) :: offset_err,slk_type,etype
4589  character(len=*),intent(in) :: uplo
4590  logical,optional,intent(in) :: is_fortran_file
4591  type(matrix_scalapack),intent(in) :: Slk_mat
4592 
4593 !Local variables ------------------------------
4594 !scalars
4595  integer :: jloc,iloc,iglob,jglob,nrows_glob,ncols_glob,mpi_err,nel
4596  integer :: bsize_frm,mpi_type_elm,ij_loc,bsize_etype,bsize_elm,cpad_bsize
4597  integer(XMPI_OFFSET_KIND) :: ijp_glob,my_offset,cpad_frm
4598 !arrays
4599  character(len=500) :: msg
4600  integer,allocatable :: block_length(:),block_type(:)
4601  integer(XMPI_ADDRESS_KIND),allocatable :: block_displ(:)
4602 
4603 !************************************************************************
4604 
4605 #ifdef HAVE_MPI_IO
4606 !@matrix_scalapack
4607  bsize_frm = xmpio_bsize_frm    ! Byte size of the Fortran record marker.
4608  if (PRESENT(is_fortran_file)) then
4609    if (.not.is_fortran_file) bsize_frm = 0
4610  end if
4611 
4612  call slk_bsize_and_type(Slk_mat,bsize_elm,mpi_type_elm)
4613 !
4614 !Global dimensions.
4615  nrows_glob=Slk_mat%sizeb_global(1)
4616  ncols_glob=Slk_mat%sizeb_global(2)
4617 !
4618 !Number of matrix elements treated by this node.
4619  nel = PRODUCT(Slk_mat%sizeb_local(1:2))
4620 !
4621 !Cannot use MPI_type_CREATE_INDEXED_BLOCK since it is not correctly implemented in several MPI libraries.
4622 !etype has to be set to MPI_BYTE, since the displacement in MPI structures is always in byte.
4623 ! MSG_WARNING("Using MPI_type_STRUCT for the MPI-IO file view")
4624 
4625  etype = MPI_BYTE
4626  call MPI_type_SIZE(etype,bsize_etype,mpi_err)
4627 !
4628 !Define the mapping between scaLAPACK buffer and the storage on file.
4629  ABI_MALLOC(block_length,(nel+2))
4630  ABI_MALLOC(block_displ,(nel+2))
4631  ABI_MALLOC(block_type,(nel+2))
4632  block_length(1)=1
4633  block_displ (1)=0
4634  block_type  (1)=MPI_LB
4635 !
4636 !Note that the view assumes that the file pointer points to the first Fortran record marker.
4637  offset_err=0
4638  select case (uplo(1:1))
4639  case ("A","a") ! The entire global matrix is stored on disk. TODO can use contigous vectors for better access.
4640    ij_loc=0
4641    do jloc=1,Slk_mat%sizeb_local(2)
4642      do iloc=1,Slk_mat%sizeb_local(1)
4643        call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4644        ij_loc  = ij_loc+1
4645        my_offset = 2*(jglob-1)*bsize_frm + bsize_frm + (jglob-1)*nrows_glob*bsize_elm + (iglob-1) * bsize_elm
4646        my_offset = my_offset / bsize_etype
4647        if (xmpio_max_address(my_offset)) offset_err=1   ! Test for possible wraparounds
4648        block_displ (ij_loc+1) = my_offset
4649        block_type  (ij_loc+1) = mpi_type_elm
4650        block_length(ij_loc+1) = 1
4651      end do
4652    end do
4653 
4654  case ("U","u") ! Only the upper triangle of the global matrix is stored on disk.
4655    ij_loc=0
4656    do jloc=1,Slk_mat%sizeb_local(2)
4657      do iloc=1,Slk_mat%sizeb_local(1)
4658        call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4659        if (jglob>=iglob) then
4660          ijp_glob = iglob + jglob*(jglob-1)/2  ! Index for packed form
4661          cpad_frm = 2*(jglob-1)*bsize_frm
4662        else
4663          ijp_glob = jglob + iglob*(iglob-1)/2  ! Index for packed form
4664          cpad_frm = 2*(iglob-1)*bsize_frm
4665        end if
4666        ij_loc = ij_loc+1
4667        my_offset = cpad_frm + bsize_frm + (ijp_glob-1) * bsize_elm
4668        my_offset = my_offset / bsize_etype
4669        if (xmpio_max_address(my_offset)) offset_err=1  ! Test for possible wraparounds
4670        block_displ (ij_loc+1) = my_offset
4671        block_type  (ij_loc+1) = mpi_type_elm
4672        block_length(ij_loc+1) = 1
4673      end do
4674    end do
4675 
4676  case ("L","l") ! Only the lower triangle of the global matrix is stored on disk.
4677    ij_loc=0
4678    do jloc=1,Slk_mat%sizeb_local(2)
4679      do iloc=1,Slk_mat%sizeb_local(1)
4680        call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4681        if (jglob<=iglob) then
4682          ijp_glob = iglob + (jglob-1)*(2*nrows_glob-jglob)/2 ! Index for packed form
4683          cpad_frm = 2*(jglob-1)*bsize_frm
4684        else
4685          ijp_glob = jglob + (iglob-1)*(2*nrows_glob-iglob)/2 ! Index for packed form
4686          cpad_frm = 2*(iglob-1)*bsize_frm
4687        end if
4688        ij_loc = ij_loc+1
4689        my_offset = cpad_frm + bsize_frm + (ijp_glob-1) * bsize_elm
4690        my_offset = my_offset / bsize_etype
4691        if (xmpio_max_address(my_offset)) offset_err=1   ! block_displ is usually integer*4. Test for possible wraparounds
4692        block_displ  (ij_loc+1) = my_offset
4693        block_type  (ij_loc+1) = mpi_type_elm
4694        block_length(ij_loc+1) = 1
4695      end do
4696    end do
4697 
4698    if (offset_err/=0) then  ! just warn, let the caller handle the exception.
4699      write(msg,"(3a)")&
4700 &     " Global position index cannot be stored in standard Fortran integer ",ch10,&
4701 &     " scaLAPACK matrix cannot be read with a single MPI-IO call ."
4702      MSG_WARNING(msg)
4703    end if
4704 
4705  case default
4706    MSG_BUG(" Wrong uplo: "//TRIM(uplo))
4707  end select
4708 
4709  block_length(nel+2)= 1
4710  block_displ (nel+2)= ncols_glob * (nrows_glob*bsize_elm + 2*bsize_frm) / bsize_etype
4711  block_type  (nel+2)= MPI_UB
4712 
4713  call xmpio_type_struct(nel+2,block_length,block_displ,block_type,slk_type,mpi_err)
4714  ABI_CHECK_MPI(mpi_err,"MPI_type_STRUCT")
4715 
4716  ABI_FREE(block_length)
4717  ABI_FREE(block_displ)
4718  ABI_FREE(block_type)
4719 
4720  call MPI_type_COMMIT(slk_type,mpi_err)
4721  ABI_CHECK_MPI(mpi_err,"MPI_type_COMMIT")
4722 
4723 #else
4724  MSG_ERROR("MPI-IO is mandatatory in slk_single_fview_read")
4725 #endif
4726 
4727 end subroutine slk_single_fview_read

m_slk/slk_single_fview_read_mask [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_single_fview_read_mask

FUNCTION

  Return an MPI datatype that can be used to read a scaLAPACK distributed matrix from
  a binary file using MPI-IO. The view is created using the user-defined mask function
  mask_of_glob. The storage of the data on file is described via the user-defined function offset_of_glob.

INPUTS

  Slk_mat<matrix_scalapack>=Structured datatype defining the scaLAPACK matrix.
  mask_of_glob(row_glob,col_glob,size_glob) is an integer function that accepts in input
     the global indeces of the matrix size_glob(1:2) are the global dimensions.
     Return 0 if (row_glob,col_glob) should not be read.
  offset_of_glob(row_glob,col_glob,size_glob,nsblocks,sub_block,bsize_elm,bsize_frm)
  nsblocks=Number of sub-blocks (will be passed to offset_of_glob)
  sub_block(2,2,nsblocks)=Global coordinates of the extremal points delimiting the sub-blocs
   e.g. sub_block(:,1,1) gives the coordinates of the left upper corner of the first block.
        sub_block(:,2,1) gives the coordinates of the right lower corner of the first block.
  [is_fortran_file]=.FALSE. is C stream is used. Set to .TRUE. for writing Fortran binary files
    with record marker.

OUTPUT

  my_nel=Number of elements that will be read by this node.
  etype=Elementary data type (handle) defining the elementary unit used to access the file.
    This is the elementary type that must be used to creae the view (MPI_BYTE is used).
  slk_type=New MPI type that can be used to instantiate the MPI-IO view for the Fortran file.
     Note that the view assumes that the file pointer points to the FIRST Fortran record marker.
  offset_err=Error code. A returned non-zero value signals that the global matrix is too large
    for a single MPI-IO access. See notes in other slk_single_fview_* routines.

SIDE EFFECTS

  myel2loc(:,:)
    input: pointer to NULL
    output: myel2loc(2,my_nel):  myel2loc(:,el) gives (iloc,jloc) for el=1,my_nel.

PARENTS

      m_exc_diago

CHILDREN

SOURCE

4265 subroutine slk_single_fview_read_mask(Slk_mat,mask_of_glob,offset_of_glob,nsblocks,sub_block,&
4266 & my_nel,myel2loc,etype,slk_type,offset_err,is_fortran_file)
4267 
4268 
4269 !This section has been created automatically by the script Abilint (TD).
4270 !Do not modify the following lines by hand.
4271 #undef ABI_FUNC
4272 #define ABI_FUNC 'slk_single_fview_read_mask'
4273 !End of the abilint section
4274 
4275  implicit none
4276 
4277 !Arguments ------------------------------------
4278 !scalars
4279  integer,intent(in) :: nsblocks
4280  integer,intent(out) :: my_nel,offset_err,slk_type,etype
4281  logical,optional,intent(in) :: is_fortran_file
4282  type(matrix_scalapack),intent(in) :: Slk_mat
4283 !arrays
4284  integer,intent(in) :: sub_block(2,2,nsblocks)
4285  integer,pointer :: myel2loc(:,:)
4286 
4287  interface
4288    function mask_of_glob(row_glob,col_glob,size_glob)
4289      use defs_basis
4290      integer :: mask_of_glob
4291      integer,intent(in) :: row_glob,col_glob
4292      integer,intent(in) :: size_glob(2)
4293    end function mask_of_glob
4294  end interface
4295 
4296  interface
4297    function offset_of_glob(row_glob,col_glob,size_glob,nsblocks,sub_block,bsize_elm,bsize_frm)
4298      use defs_basis
4299      use m_xmpi
4300      integer(XMPI_OFFSET_KIND) :: offset_of_glob
4301      integer,intent(in) :: row_glob,col_glob,bsize_elm,bsize_frm,nsblocks
4302      integer,intent(in) :: size_glob(2),sub_block(2,2,nsblocks)
4303    end function offset_of_glob
4304  end interface
4305 
4306 !Local variables ------------------------------
4307 !scalars
4308  integer :: el,jloc,iloc,iglob,jglob,mpi_err,sweep
4309  integer :: bsize_frm,mpi_type_elm,bsize_elm
4310  integer(XMPI_OFFSET_KIND) :: tmp_off,max_displ
4311 !arrays
4312  character(len=500) :: msg
4313  integer,allocatable :: block_length(:),block_type(:)
4314  integer(XMPI_ADDRESS_KIND),allocatable :: block_displ(:)
4315 
4316 !************************************************************************
4317 
4318 #ifdef HAVE_MPI_IO
4319 !@matrix_scalapack
4320  bsize_frm = xmpio_bsize_frm  ! Byte size of the Fortran record marker.
4321  if (PRESENT(is_fortran_file)) then
4322    if (.not.is_fortran_file) bsize_frm = 0
4323  end if
4324 
4325  ! Byte size of the matrix element.
4326  call slk_bsize_and_type(Slk_mat,bsize_elm,mpi_type_elm)
4327  !
4328  ! Find the number of local matrix elements to be read, then create the table myel2loc.
4329  do sweep=1,2
4330    if (sweep==2) then
4331       ABI_MALLOC(myel2loc,(2,my_nel))
4332    end if
4333    my_nel=0
4334 
4335    do jloc=1,Slk_mat%sizeb_local(2)
4336      do iloc=1,Slk_mat%sizeb_local(1)
4337        call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4338        if ( mask_of_glob(iglob,jglob,Slk_mat%sizeb_global)/= 0) then ! Will fill this entry.
4339          my_nel  = my_nel+1
4340          if (sweep==2) myel2loc(:,my_nel) = (/iloc,jloc/)
4341        end if
4342      end do
4343    end do
4344  end do
4345  !
4346  etype = MPI_BYTE
4347  !
4348  ! Define the mapping between scaLAPACK buffer and the storage on file.
4349  ! Note that the view assumes that the file pointer points to the first Fortran record marker.
4350  ABI_MALLOC(block_length,(my_nel+2))
4351  ABI_MALLOC(block_displ,(my_nel+2))
4352  ABI_MALLOC(block_type,(my_nel+2))
4353  block_length(1)=1
4354  block_displ (1)=0
4355  block_type  (1)=MPI_LB
4356 
4357  offset_err=0; max_displ=0
4358  do el=1,my_nel
4359    iloc = myel2loc(1,el)
4360    jloc = myel2loc(2,el)
4361    call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4362    tmp_off = offset_of_glob(iglob,jglob,Slk_mat%sizeb_global,nsblocks,sub_block,bsize_elm,bsize_frm)
4363    if (xmpio_max_address(tmp_off)) offset_err=1   ! Test for possible wraparounds.
4364    max_displ = MAX(max_displ,tmp_off)
4365    block_displ (el+1) = tmp_off
4366    block_type  (el+1) = mpi_type_elm
4367    block_length(el+1) = 1
4368    !write(std_out,*)" iglob, jglob, tmp_off ",iglob, jglob, tmp_off
4369  end do
4370  !write(std_out,*)" MAX displ is ",MAXVAL(block_displ)
4371 
4372  if (offset_err/=0) then  ! just warn, let the caller handle the exception.
4373    write(msg,"(3a)")&
4374 &   " Global position index cannot be stored in standard Fortran integer ",ch10,&
4375 &   " scaLAPACK matrix cannot be read with a single MPI-IO call ."
4376    MSG_WARNING(msg)
4377  end if
4378 
4379  block_length(my_nel+2) = 1
4380  block_displ (my_nel+2) = max_displ
4381  block_type  (my_nel+2) = MPI_UB
4382 
4383  call xmpio_type_struct(my_nel+2,block_length,block_displ,block_type,slk_type,mpi_err)
4384  ABI_CHECK_MPI(mpi_err,"MPI_type_STRUCT")
4385 
4386  ABI_FREE(block_length)
4387  ABI_FREE(block_displ)
4388  ABI_FREE(block_type)
4389 
4390  call MPI_type_COMMIT(slk_type,mpi_err)
4391  ABI_CHECK_MPI(mpi_err,"MPI_type_COMMIT")
4392 
4393 #else
4394  MSG_ERROR("MPI-IO is mandatatory in slk_single_fview_read_mask")
4395 #endif
4396 
4397 end subroutine slk_single_fview_read_mask

m_slk/slk_single_fview_write [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_single_fview_write

FUNCTION

  Returns an MPI datatype that can be used to write a scaLAPACK distributed matrix to
  a binary file using MPI-IO.

INPUTS

  Slk_mat<matrix_scalapack>=Structured datatype defining the scaLAPACK distribution with the local buffer.
  uplo=String specifying whether only the upper or lower triangular part of the global matrix is stored on disk:
    = "U":  Upper triangular is stored
    = "L":  Lower triangular is stored
    = "A":  Full matrix (used for general complex matrices)
  [is_fortran_file]=.FALSE. is C stream is used. .TRUE. for writing Fortran binary files
    with record marker. In this case etype is set xmpio_mpi_type_frm provided that
    the mpi_type of the matrix element is commensurate with xmpio_mpi_type_frm. Defaults to .TRUE.
  glob_subarray(2,2) = Used to select the subarray of the global matrix. Used only when uplo="All"
     glob_subarray(:,1)=starting global coordinates of the subarray in each dimension (array of nonnegative integers >=1, <=array_of_sizes)
     glob_subarray(:,2)=Number of elements in each dimension of the subarray (array of positive integers)

OUTPUT

  nelw=Number of elements to be written.
  etype=Elementary data type (handle) defining the elementary unit used to access the file.
  slk_type=New MPI type that can be used to instantiate the MPI-IO view for the Fortran file.
     Note that the view assumes that the file pointer points to the FIRST Fortran record marker.
  offset_err=Error code. A non-zero value signals that the global matrix is too large
    for a single MPI-IO access (see notes below).

SIDE EFFECTS

  elw2slk(:,:) =
    input:  pointer to null().
    output: elw2slk(2,nelw) contains the local coordinates of the matrix elements to be written.
     (useful only if the upper or lower triangle of the global matrix has to be written or when
      uplo="all" but a global subarray is written.

NOTES

  With (signed) Fortran integers, the maximum size of the file that
  that can be read in one-shot is around 2Gb when etype is set to byte.
  Using a larger etype might create portability problems (real data on machines using
  integer*16 for the marker) since etype must be a multiple of the Fortran record marker
  Due to the above reason, block_displ is given in bytes and must be stored in a Fortran
  integer of kind XMPI_ADDRESS_KIND. If the displacement is too large, the routine returns
  offset_err=1 so that the caller will know that several MPI-IO reads are needed to
  write the local buffer.

PARENTS

      m_slk

CHILDREN

SOURCE

4785 subroutine slk_single_fview_write(Slk_mat,uplo,nelw,elw2slk,etype,slk_type,offset_err,is_fortran_file,glob_subarray)
4786 
4787 
4788 !This section has been created automatically by the script Abilint (TD).
4789 !Do not modify the following lines by hand.
4790 #undef ABI_FUNC
4791 #define ABI_FUNC 'slk_single_fview_write'
4792 !End of the abilint section
4793 
4794  implicit none
4795 
4796 !Arguments ------------------------------------
4797 !scalars
4798  integer,intent(out) :: offset_err,slk_type,etype,nelw
4799  character(len=*),intent(in) :: uplo
4800  logical,optional,intent(in) :: is_fortran_file
4801  type(matrix_scalapack),intent(in) :: Slk_mat
4802 !arrays
4803  integer,pointer :: elw2slk(:,:)
4804  integer,optional,intent(in) :: glob_subarray(2,2)
4805 
4806 !Local variables ------------------------------
4807 !scalars
4808  integer :: jloc,iloc,iglob,jglob,nrows_glob,ncols_glob,mpi_err,nel_max
4809  integer :: grow_min,grow_max,gcol_min,gcol_max
4810  integer :: bsize_frm,mpi_type_elm,ij_loc,bsize_elm,cpad_bsize
4811  integer(XMPI_OFFSET_KIND) :: ijp_glob,my_offset,cpad_frm
4812 !arrays
4813  character(len=500) :: msg
4814  integer :: starts(2),sub_sizes(2)
4815  integer,allocatable :: block_length(:),block_type(:)
4816  integer(XMPI_ADDRESS_KIND),allocatable :: block_displ(:)
4817 
4818 !************************************************************************
4819 
4820 #ifdef HAVE_MPI_IO
4821 !@matrix_scalapack
4822  bsize_frm = xmpio_bsize_frm    ! Byte size of the Fortran record marker.
4823  if (PRESENT(is_fortran_file)) then
4824    if (.not.is_fortran_file) bsize_frm = 0
4825  end if
4826 
4827  if (PRESENT(glob_subarray).and..not.firstchar(uplo,(/"A"/))) then
4828    MSG_ERROR("glob_subarray should not be used when uplo/=All")
4829  end if
4830 
4831  call slk_bsize_and_type(Slk_mat,bsize_elm,mpi_type_elm)
4832  !
4833  ! Global dimensions.
4834  nrows_glob=Slk_mat%sizeb_global(1)
4835  ncols_glob=Slk_mat%sizeb_global(2)
4836  !
4837  ! Number of matrix elements treated by this node.
4838  nel_max = PRODUCT(Slk_mat%sizeb_local(1:2))
4839 
4840  ABI_MALLOC(elw2slk,(2,nel_max))
4841  elw2slk=0
4842 !
4843 !Cannot use MPI_type_CREATE_INDEXED_BLOCK since it is not correctly implemented in several MPI libraries.
4844 !etype has to be set to MPI_BYTE, since the displacement in MPI structures is always in byte.
4845  etype = MPI_BYTE
4846 !
4847 !Define the mapping between scaLAPACK buffer and the storage on file.
4848  ABI_MALLOC(block_length,(nel_max+2))
4849  ABI_MALLOC(block_displ,(nel_max+2))
4850  ABI_MALLOC(block_type,(nel_max+2))
4851  block_length(1)=1
4852  block_displ (1)=0
4853  block_type  (1)=MPI_LB
4854 
4855 !
4856 !Note that the view assumes that the file pointer points to the first Fortran record marker.
4857  offset_err=0
4858  select case (uplo(1:1))
4859 
4860  case ("A","a") ! The entire global matrix is written on disk. TODO can use contigous vectors for better access.
4861 
4862    grow_min=1; grow_max=nrows_glob
4863    gcol_min=1; gcol_max=ncols_glob
4864    if (PRESENT(glob_subarray)) then ! subarray access.
4865      grow_min = glob_subarray(1,1)
4866      gcol_min = glob_subarray(2,1)
4867      grow_max = grow_min + glob_subarray(1,2) -1
4868      gcol_max = gcol_min + glob_subarray(2,2) -1
4869    end if
4870 
4871    ij_loc=0
4872    do jloc=1,Slk_mat%sizeb_local(2)
4873      do iloc=1,Slk_mat%sizeb_local(1)
4874        call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4875        if (iglob>=grow_min.and.iglob<=grow_max .and. &  ! glob_subarray element.
4876            jglob>=gcol_min.and.jglob<=gcol_max) then
4877          ij_loc  = ij_loc+1
4878          my_offset = 2*(jglob-1)*bsize_frm + bsize_frm + (jglob-1)*nrows_glob*bsize_elm + (iglob-1) * bsize_elm
4879          if (xmpio_max_address(my_offset)) offset_err=1   ! Test for possible wraparounds
4880          block_displ (ij_loc+1) = my_offset
4881          block_type  (ij_loc+1) = mpi_type_elm
4882          block_length(ij_loc+1) = 1
4883          elw2slk(:,ij_loc) = (/iloc,jloc/) ! useless when subarray are not used but oh well!
4884        end if
4885      end do
4886    end do
4887 
4888  case ("U","u") ! Only the upper triangle of the global matrix is stored on disk.
4889    ij_loc=0
4890    do jloc=1,Slk_mat%sizeb_local(2)
4891      do iloc=1,Slk_mat%sizeb_local(1)
4892        call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4893        if (jglob>=iglob) then
4894          ijp_glob = iglob + jglob*(jglob-1)/2  ! Index for packed form
4895          cpad_frm = 2*(jglob-1)*bsize_frm
4896          ij_loc = ij_loc+1
4897          my_offset = cpad_frm + bsize_frm + (ijp_glob-1) * bsize_elm
4898          if (xmpio_max_address(my_offset)) offset_err=1   ! Test for possible wraparounds
4899          block_displ (ij_loc+1) = my_offset
4900          block_type  (ij_loc+1) = mpi_type_elm
4901          block_length(ij_loc+1) = 1
4902          elw2slk(:,ij_loc) = (/iloc,jloc/)
4903        end if
4904      end do
4905    end do
4906 
4907  case ("L","l") ! Only the lower triangle of the global matrix is stored on disk.
4908    ij_loc=0
4909    do jloc=1,Slk_mat%sizeb_local(2)
4910      do iloc=1,Slk_mat%sizeb_local(1)
4911        call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4912        if (jglob<=iglob) then
4913          ijp_glob = iglob + (jglob-1)*(2*nrows_glob-jglob)/2 ! Index for packed form
4914          cpad_frm = 2*(jglob-1)*bsize_frm
4915          ij_loc = ij_loc+1
4916          my_offset = cpad_frm + bsize_frm + (ijp_glob-1) * bsize_elm
4917          if (xmpio_max_address(my_offset)) offset_err=1   ! block_displ is usually integer*4. Test for possible wraparounds
4918          block_displ (ij_loc+1) = my_offset
4919          block_type  (ij_loc+1) = mpi_type_elm
4920          block_length(ij_loc+1) = 1
4921          elw2slk(:,ij_loc) = (/iloc,jloc/)
4922        end if
4923      end do
4924    end do
4925 
4926  case default
4927    MSG_BUG(" Wrong uplo: "//TRIM(uplo))
4928  end select
4929 
4930  if (offset_err/=0) then  ! just warn, let the caller handle the exception.
4931    write(msg,"(3a)")&
4932 &   "Global position index cannot be stored in standard Fortran integer ",ch10,&
4933 &   "scaLAPACK matrix cannot be read with a single MPI-IO call ."
4934    MSG_WARNING(msg)
4935  end if
4936  !
4937  ! Final number of matrix elements that will be written by this node.
4938  nelw = ij_loc
4939 
4940  block_length(nelw+2)= 1
4941  block_displ (nelw+2)= ncols_glob * (nrows_glob*bsize_elm + 2*bsize_frm)
4942  block_type  (nelw+2)= MPI_UB
4943 
4944  call xmpio_type_struct(nelw+2,block_length,block_displ,block_type,slk_type,mpi_err)
4945  ABI_CHECK_MPI(mpi_err,"MPI_type_STRUCT")
4946 
4947  ABI_FREE(block_length)
4948  ABI_FREE(block_displ)
4949  ABI_FREE(block_type)
4950 
4951  call MPI_type_COMMIT(slk_type,mpi_err)
4952  ABI_CHECK_MPI(mpi_err,"MPI_type_COMMIT")
4953 
4954 #else
4955  MSG_ERROR("MPI-IO is mandatatory in slk_single_fview_read")
4956 #endif
4957 
4958 end subroutine slk_single_fview_write

m_slk/slk_symmetrize [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_symmetrize

FUNCTION

  This routine symmetrize a square scaLAPACK-distributed matrix.

INPUTS

  uplo=String specifying whether only the upper or lower triangular part of the global matrix has been read
    = "U":  Upper triangular has been read.
    = "L":  Lower triangular has been read.
    = "A":  Full matrix (used for general complex matrices)
  symtype=Symmetry type of the matrix (used only if uplo = "L" or "A").
    = "H" for Hermitian matrix
    = "S" for symmetric matrix.
    = "N" if matrix has no symmetry (not compatible with uplo="L" or uplo="U".

SIDE EFFECTS

  Slk_mat<matrix_scalapack>=Structured datatype defining the scaLAPACK distribution with the local buffer
    supposed to be allocated.
    %buffer_cplx=Local buffer containg the distributed matrix stored on the external file.

PARENTS

      m_slk

CHILDREN

SOURCE

4431 subroutine slk_symmetrize(Slk_mat,uplo,symtype)
4432 
4433 
4434 !This section has been created automatically by the script Abilint (TD).
4435 !Do not modify the following lines by hand.
4436 #undef ABI_FUNC
4437 #define ABI_FUNC 'slk_symmetrize'
4438 !End of the abilint section
4439 
4440  implicit none
4441 
4442 !Arguments ------------------------------------
4443 !scalars
4444  character(len=*),intent(in) :: symtype
4445  character(len=*),intent(in) :: uplo
4446  type(matrix_scalapack),intent(inout) :: Slk_mat
4447 
4448 !Local variables ------------------------------
4449 !scalars
4450  integer :: jloc,iloc,iglob,jglob,ij_loc
4451  logical :: is_hermitian,is_real,is_cplx,is_symmetric
4452  character(len=500) :: msg
4453 
4454 !************************************************************************
4455 
4456 !@matrix_scalapack
4457  is_cplx = (allocated(Slk_mat%buffer_cplx))
4458  is_real = (allocated(Slk_mat%buffer_real))
4459 !
4460 !One and only one buffer should be allocated.
4461  if (is_real.and.is_cplx) then
4462    write(msg,'(a,2l1)')" ScaLAPACK buffers are not allocated correctly, is_real=, is_cplx ",is_real,is_cplx
4463    MSG_ERROR(msg)
4464  end if
4465 
4466  if (is_real) RETURN
4467 
4468  is_hermitian=.FALSE.; is_symmetric=.FALSE.
4469  select case (symtype(1:1))
4470    case ("H","h")
4471      is_hermitian=.TRUE.
4472    case ("S","s")
4473      is_symmetric=.TRUE.
4474    case("N","n")
4475      if ( ALL(uplo(1:1)/=(/"A","a"/)) ) then
4476        msg = " Found symtype= "//TRIM(symtype)//", but uplo= "//TRIM(uplo)
4477        MSG_ERROR(msg)
4478      end if
4479      RETURN  ! Nothing to do.
4480    case default
4481      MSG_ERROR("Wrong symtype "//TRIM(symtype))
4482  end select
4483 
4484  !write(std_out,*)"is_cplx",is_cplx
4485  !write(std_out,*)"is_hermitian",is_hermitian
4486 
4487  select case (uplo(1:1))
4488 
4489  case ("A","a") ! Full global matrix has been read, nothing to do.
4490    RETURN
4491 
4492  case ("U","u") ! Only the upper triangle of the global matrix was read.
4493 
4494    if (is_cplx .and. is_hermitian) then
4495      ij_loc=0
4496      do jloc=1,Slk_mat%sizeb_local(2)
4497        do iloc=1,Slk_mat%sizeb_local(1)
4498          call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4499          ij_loc = ij_loc+1
4500          if (jglob<iglob) then ! Diagonal elements are not forced to be real.
4501            Slk_mat%buffer_cplx(iloc,jloc) =  DCONJG(Slk_mat%buffer_cplx(iloc,jloc))
4502          end if
4503 !        if (iglob==jglob) Slk_mat%buffer_cplx(iloc,jloc) =  real(Slk_mat%buffer_cplx(iloc,jloc))
4504        end do
4505      end do
4506    end if
4507 
4508  case ("L","l") ! Only the lower triangle of the global matrix was read.
4509 
4510    if (is_cplx .and. is_hermitian) then
4511      ij_loc=0
4512      do jloc=1,Slk_mat%sizeb_local(2)
4513        do iloc=1,Slk_mat%sizeb_local(1)
4514          call idx_glob(Slk_mat,iloc,jloc,iglob,jglob)
4515          ij_loc = ij_loc+1
4516          if (jglob>iglob) then ! diagonal elements are not forced to be real.
4517            Slk_mat%buffer_cplx(iloc,jloc) =  DCONJG(Slk_mat%buffer_cplx(iloc,jloc))
4518          end if
4519 !        if (iglob==jglob) Slk_mat%buffer_cplx(iloc,jloc) =  real(Slk_mat%buffer_cplx(iloc,jloc))
4520        end do
4521      end do
4522    end if
4523 
4524  case default
4525    MSG_BUG(" Wrong uplo: "//TRIM(uplo))
4526  end select
4527 
4528 end subroutine slk_symmetrize

m_slk/slk_write [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

  slk_write

FUNCTION

  Routine to write a squared scaLAPACK distributed matrix on an external file using MPI-IO.

INPUTS

  Slk_mat<matrix_scalapack>=Structured datatype defining the scaLAPACK distribution with the local buffer
    containing the distributed matrix.
  uplo=String specifying whether only the upper or lower triangular part of the global matrix is used:
    = "U":  Upper triangular
    = "L":  Lower triangular
    = "A":  Full matrix (used for general complex matrices)
  is_fortran_file=.FALSE. is C stream is used. .TRUE. for writing Fortran binary files.
  [fname]= Mutually exclusive with mpi_fh. The name of the external file on which the matrix will be written.
           The file is open and closed inside the routine with MPI flags specified by flags.
  [mpi_fh]=File handler associated to the file (already open in the caller). Not compatible with fname.
  [flags]=MPI-IO flags used to open the file in MPI_FILE_OPEN.
    Default is MPI_MODE_CREATE + MPI_MODE_WRONLY + MPI_MODE_EXCL.
  [glob_subarray(2,2)] = Used to select the subarray of the global matrix. Used only when uplo="All"
     NOTE that each node should call the routine with the same value.
     glob_subarray(:,1)=starting global coordinates of the subarray in each dimension (array of nonnegative integers >=1, <=array_of_sizes)
     glob_subarray(:,2)=Number of elements in each dimension of the subarray (array of positive integers)

OUTPUT

  Only writing. The global scaLAPACK matrix is written on file fname.
  If fname is present then the file is open and closed inside the routine. Any exception is fatal.

SIDE EFFECTS

  [offset]=
    input:  Offset used to access the content of the file. Default is zero.
    output: New offset incremented with the byte size of the matrix that has been read (Fortran
            markers are included if is_fortran_file=.TRUE.)

TODO

  * Generalize the implementation adding the writing the real buffer.

PARENTS

      m_exc_diago

CHILDREN

SOURCE

3847 subroutine slk_write(Slk_mat,uplo,is_fortran_file,fname,mpi_fh,offset,flags,glob_subarray)
3848 
3849 
3850 !This section has been created automatically by the script Abilint (TD).
3851 !Do not modify the following lines by hand.
3852 #undef ABI_FUNC
3853 #define ABI_FUNC 'slk_write'
3854 !End of the abilint section
3855 
3856  implicit none
3857 
3858 !Arguments ------------------------------------
3859 !scalars
3860  integer,optional,intent(in) :: flags
3861  integer,optional,intent(inout) :: mpi_fh
3862  integer(XMPI_OFFSET_KIND),optional,intent(inout) :: offset
3863  logical,intent(in) :: is_fortran_file
3864  character(len=*),optional,intent(in) :: fname
3865  character(len=*),intent(in) :: uplo
3866  type(matrix_scalapack),intent(in) :: Slk_mat
3867 !array
3868  integer,optional,intent(in) :: glob_subarray(2,2)
3869 
3870 !Local variables ------------------------------
3871 !scalars
3872  integer :: jloc,iloc,iglob,jglob,nrows_glob,ncols_glob,elw,nrows_w,ncols_w
3873  integer :: slk_type,offset_err,etype,nfrec,bsize_elm,mpi_type_elm
3874  integer(XMPI_OFFSET_KIND) :: my_offset
3875  logical :: do_open
3876 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
3877  integer :: comm,my_flags,my_fh,buffer_size
3878  integer ::ierr,ij_loc,nelw,col_glob
3879 !arrays
3880  integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecord(:)
3881  integer,pointer :: elw2slk(:,:)
3882  complex(dpc),allocatable :: buffer1_cplx(:)
3883 #endif
3884  character(len=500) :: msg
3885 
3886 !************************************************************************
3887 
3888 !@matrix_scalapack
3889  ABI_CHECK(allocated(Slk_mat%buffer_cplx),"%buffer_cplx not allocated")
3890 
3891  if (firstchar(uplo,(/"U","L"/)) .and. Slk_mat%sizeb_global(1) /= Slk_mat%sizeb_global(2) ) then
3892    MSG_ERROR("rectangular matrices are not compatible with the specified uplo")
3893  end if
3894 
3895  if (PRESENT(glob_subarray).and..not.firstchar(uplo,(/"A"/))) then
3896    MSG_ERROR("glob_subarray should not be used when uplo/=All")
3897  end if
3898 
3899  do_open = PRESENT(fname)
3900  if (do_open) then
3901    ABI_CHECK(.not.PRESENT(fname),"fname should not be present")
3902  else
3903    ABI_CHECK(PRESENT(mpi_fh),"mpi_fh should be present")
3904  end if
3905 
3906  my_offset=0; if (PRESENT(offset)) my_offset=offset
3907 
3908 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
3909  comm = Slk_mat%processor%comm
3910 
3911  nrows_glob=Slk_mat%sizeb_global(1)
3912  ncols_glob=Slk_mat%sizeb_global(1)
3913  buffer_size= PRODUCT(Slk_mat%sizeb_local(1:2))
3914 
3915  call slk_bsize_and_type(Slk_mat,bsize_elm,mpi_type_elm)
3916 
3917  if (do_open) then !Open the file.
3918    my_flags=MPI_MODE_CREATE + MPI_MODE_WRONLY + MPI_MODE_APPEND
3919    if (PRESENT(flags)) my_flags = flags
3920 
3921    call MPI_FILE_OPEN(comm, fname, my_flags, MPI_INFO_NULL, my_fh, ierr)
3922    msg = " MPI_FILE_OPEN "//TRIM(fname)
3923    ABI_CHECK_MPI(ierr,msg)
3924  else
3925    my_fh = mpi_fh
3926  end if
3927 
3928  if (PRESENT(glob_subarray)) then
3929    call slk_single_fview_write(Slk_mat,uplo,nelw,elw2slk,etype,slk_type,offset_err,&
3930 &    is_fortran_file=is_fortran_file,glob_subarray=glob_subarray)
3931  else
3932    call slk_single_fview_write(Slk_mat,uplo,nelw,elw2slk,etype,slk_type,offset_err,&
3933 &    is_fortran_file=is_fortran_file)
3934  end if
3935 
3936  if (offset_err/=0) then
3937    write(msg,"(3a)")&
3938 &   " Global position index cannot be stored in standard Fortran integer ",ch10,&
3939 &   " scaLAPACK matrix cannot be read with a single MPI-IO call."
3940    MSG_ERROR(msg)
3941  end if
3942 
3943  call MPI_FILE_SET_VIEW(my_fh, my_offset, etype, slk_type, 'native', MPI_INFO_NULL, ierr)
3944  ABI_CHECK_MPI(ierr,"SET_VIEW")
3945 
3946  call MPI_type_FREE(slk_type,ierr)
3947  ABI_CHECK_MPI(ierr,"MPI_type_FREE")
3948 
3949  if (nelw==buffer_size) then ! Dump Slk_mat% immediately.
3950    call MPI_FILE_WRITE_ALL(my_fh, Slk_mat%buffer_cplx, buffer_size, MPI_DOUBLE_complex, MPI_STATUS_IGNORE, ierr)
3951    ABI_CHECK_MPI(ierr,"WRITE_ALL")
3952  else ! Have to extract the data to be written.
3953    ABI_MALLOC(buffer1_cplx,(nelw))
3954    do elw=1,nelw
3955      iloc = elw2slk(1,elw)
3956      jloc = elw2slk(2,elw)
3957      buffer1_cplx(elw) = Slk_mat%buffer_cplx(iloc,jloc)
3958    end do
3959    call MPI_FILE_WRITE_ALL(my_fh, buffer1_cplx, nelw, MPI_DOUBLE_complex, MPI_STATUS_IGNORE, ierr)
3960    ABI_CHECK_MPI(ierr,"WRITE_ALL")
3961    ABI_FREE(buffer1_cplx)
3962  end if
3963 
3964  ABI_FREE(elw2slk)
3965  !
3966  ! Number of columns and rows that have been written.
3967  ! Used to write the Fortran markers and to increment the offset.
3968  nrows_w = nrows_glob
3969  ncols_w = ncols_glob
3970  if (PRESENT(glob_subarray)) then
3971    nrows_w = glob_subarray(1,2) - glob_subarray(1,1) + 1
3972    ncols_w = glob_subarray(2,2) - glob_subarray(2,1) + 1
3973    if (.not.firstchar(uplo,(/"A"/))) then
3974      MSG_ERROR("glob_subarray should not be used when uplo/=All")
3975    end if
3976  end if
3977 
3978  !TODO check whether slk_single_fview_write can report an offset to reduce the extent.
3979  if (is_fortran_file) then ! Collective writing of the Fortran markers.
3980    nfrec = ncols_w
3981    ABI_MALLOC(bsize_frecord,(nfrec))
3982    if (firstchar(uplo,(/"A"/))) then
3983      bsize_frecord = nrows_w * bsize_elm
3984    else if (firstchar(uplo,(/"U"/))) then
3985      bsize_frecord = (/(col_glob * bsize_elm, col_glob=1,nfrec)/)
3986    else if (firstchar(uplo,(/"L"/))) then
3987      bsize_frecord = (/(col_glob * bsize_elm, col_glob=nfrec,1,-1)/)
3988    else
3989      MSG_ERROR("Wrong uplo")
3990    end if
3991    call xmpio_write_frmarkers(mpi_fh,my_offset,xmpio_collective,nfrec,bsize_frecord,ierr)
3992    ABI_CHECK(ierr==0,"Error while writing Fortran markers")
3993    ABI_FREE(bsize_frecord)
3994  end if
3995  !
3996  if (do_open) then ! Close the file.
3997    call MPI_FILE_CLOSE(my_fh, ierr)
3998    ABI_CHECK_MPI(ierr,"FILE_CLOSE")
3999  end if
4000  !
4001  ! Increment the offset
4002  if (PRESENT(offset)) then
4003    if (firstchar(uplo,(/"A"/))) then
4004      offset = offset + nrows_w*ncols_w*bsize_elm
4005      if (is_fortran_file) offset = offset + ncols_w*2*xmpio_bsize_frm
4006    else if (firstchar(uplo,(/"U","L"/))) then
4007      offset = offset + ( (Slk_mat%sizeb_global(2) * (Slk_mat%sizeb_global(2))+1)/2 ) * bsize_elm
4008      if (is_fortran_file) offset = offset + Slk_mat%sizeb_global(2)*2*xmpio_bsize_frm
4009    else
4010      MSG_ERROR("Wrong uplo")
4011    end if
4012  end if
4013 
4014  call xmpi_barrier(comm)
4015  RETURN
4016 
4017 #else
4018   MSG_ERROR("MPI-IO support not activated")
4019 #endif
4020 
4021 end subroutine slk_write

m_slk/slk_zdhp_invert [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

 slk_zdhp_invert

FUNCTION

  slk_zdhp_invert provides an object-oriented interface to the ScaLAPACK set of routines used to compute
  the inverse of a Hermitian positive definite matrix.

INPUTS

  uplo(global input)
    = 'U':  Upper triangle of sub( A ) is stored;
    = 'L':  Lower triangle of sub( A ) is stored.

SIDE EFFECTS

  Slk_mat<type(matrix_scalapack)>=The object storing the local buffer, the array descriptor, the context
    and other quantities needed to call ScaLAPACK routines.
    On entry, this array contains the local pieces of the N-by-N Hermitian distributed matrix sub( A ) to be factored.
    If UPLO = 'U', the leading N-by-N upper triangular part of sub( A ) contains the upper triangular part of the matrix,
    and its strictly lower triangular part is not referenced.
    If UPLO = 'L', the leading N-by-N lower triangular part of sub( A ) contains the lower triangular part of the distribu-
    ted matrix, and its strictly upper triangular part is not referenced.
    On exit, the local pieces of the upper or lower triangle of the (Hermitian) inverse of sub( A )

PARENTS

      m_abilasi

CHILDREN

SOURCE

3753 subroutine slk_zdhp_invert(Slk_mat,uplo)
3754 
3755 
3756 !This section has been created automatically by the script Abilint (TD).
3757 !Do not modify the following lines by hand.
3758 #undef ABI_FUNC
3759 #define ABI_FUNC 'slk_zdhp_invert'
3760 !End of the abilint section
3761 
3762  implicit none
3763 
3764 !Arguments ------------------------------------
3765 !scalars
3766  character(len=*),intent(in) :: uplo
3767  type(matrix_scalapack),intent(inout) :: Slk_mat
3768 
3769 !Local variables ------------------------------
3770 !scalars
3771  integer :: info
3772  character(len=500) :: msg
3773 
3774 !************************************************************************
3775 
3776  ABI_CHECK(allocated(Slk_mat%buffer_cplx),"buffer_cplx not allocated")
3777 
3778  ! *  ZPOTRF computes the Cholesky factorization of a complex Hermitian positive definite.
3779  ! *     A = U**H * U,   if UPLO = 'U', or
3780  ! *     A = L  * L**H,  if UPLO = 'L',
3781  call PZPOTRF(uplo,Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx,1,1,Slk_mat%descript%tab,info)
3782 
3783  if (info/=0) then
3784    write(msg,'(a,i0)')" PZPOTRF returned info= ",info
3785    MSG_ERROR(msg)
3786  end if
3787  !
3788  ! PZPOTRI computes the inverse of a complex Hermitian positive definite
3789  ! distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1) using the
3790  ! Cholesky factorization sub( A ) = U**H*U or L*L**H computed by PZPOTRF.
3791  call PZPOTRI(uplo,Slk_mat%sizeb_global,Slk_mat%buffer_cplx,1,1,Slk_mat%descript%tab,info)
3792 
3793  if (info/=0) then
3794    write(msg,'(a,i0)')" PZPOTRI returned info= ",info
3795    MSG_ERROR(msg)
3796  end if
3797 
3798 end subroutine slk_zdhp_invert

m_slk/slk_zinvert [ Functions ]

[ Top ] [ m_slk ] [ Functions ]

NAME

 slk_zinvert

FUNCTION

  slk_zinvert provides an object-oriented interface to the ScaLAPACK set of routines used to compute
  the inverse of a complex matrix in double precision

SIDE EFFECTS

  Slk_mat<type(matrix_scalapack)>=The object storing the local buffer, the array descriptor, the context
    and other quantities needed to call ScaLAPACK routines. In input, the matrix to invert, in output
    the matrix inverted and distributed among the nodes.

PARENTS

      m_abilasi,m_exc_diago

CHILDREN

SOURCE

3643 subroutine slk_zinvert(Slk_mat)
3644 
3645 
3646 !This section has been created automatically by the script Abilint (TD).
3647 !Do not modify the following lines by hand.
3648 #undef ABI_FUNC
3649 #define ABI_FUNC 'slk_zinvert'
3650 !End of the abilint section
3651 
3652  implicit none
3653 
3654 !Arguments ------------------------------------
3655 !scalars
3656  type(matrix_scalapack),intent(inout) :: Slk_mat
3657 
3658 !Local variables ------------------------------
3659 !scalars
3660  integer :: lwork,info
3661  integer :: ipiv_size,liwork
3662  character(len=500) :: msg
3663 !array
3664  integer,allocatable :: ipiv(:),iwork(:)
3665  complex(dpc),allocatable :: work(:)
3666 
3667 !************************************************************************
3668 
3669  ABI_CHECK(allocated(Slk_mat%buffer_cplx),"buffer_cplx not allocated")
3670 
3671 !IMPORTANT NOTE: PZGETRF requires square block decomposition i.e.,  MB_A = NB_A.
3672  if ( Slk_mat%descript%tab(MB_)/=Slk_mat%descript%tab(NB_) ) then
3673    msg =" PZGETRF requires square block decomposition i.e.,  MB_A = NB_A."
3674    MSG_ERROR(msg)
3675  end if
3676 
3677  ipiv_size = my_locr(Slk_mat) + Slk_mat%descript%tab(MB_)
3678  ABI_MALLOC(ipiv,(ipiv_size))
3679 
3680  call PZGETRF(Slk_mat%sizeb_global(1),Slk_mat%sizeb_global(2),Slk_mat%buffer_cplx,&
3681 & 1,1,Slk_mat%descript%tab,ipiv,info) ! P * L * U  Factorization.
3682 
3683  if (info/=0) then
3684    write(msg,'(a,i7)')" PZGETRF returned info= ",info
3685    MSG_ERROR(msg)
3686  end if
3687 
3688 !Get optimal size of workspace for PZGETRI.
3689  lwork=-1; liwork=-1
3690  ABI_MALLOC(work,(1))
3691  ABI_MALLOC(iwork,(1))
3692 
3693  call PZGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx,1,1,Slk_mat%descript%tab,ipiv,&
3694 & work,lwork,iwork,liwork,info)
3695 
3696  ABI_CHECK(info==0,"PZGETRI: Error during compuation of workspace size")
3697 
3698  lwork = NINT(real(work(1))); liwork=iwork(1)
3699  ABI_FREE(work)
3700  ABI_FREE(iwork)
3701 
3702 !Solve the problem.
3703  ABI_MALLOC(work,(lwork))
3704  ABI_MALLOC(iwork,(liwork))
3705 
3706  call PZGETRI(Slk_mat%sizeb_global(1),Slk_mat%buffer_cplx,1,1,Slk_mat%descript%tab,ipiv,&
3707 & work,lwork,iwork,liwork,info)
3708 
3709  if (info/=0) then
3710    write(msg,'(a,i7)')" PZGETRI returned info= ",info
3711    MSG_ERROR(msg)
3712  end if
3713 
3714  ABI_FREE(work)
3715  ABI_FREE(iwork)
3716  ABI_FREE(ipiv)
3717 
3718 end subroutine slk_zinvert