FortranZT
FortranZT

Reputation: 1

Using ScaLAPACK subroutine PDGEMR2D gives errors

I am working on a problem where I need the full spectrum of a Hamiltonian matrix. I order to have a better performance, I am using ScaLAPACK routine PDSYEV to compute the eigenstates and eigenvalues of very large matrices (N > 15000), and the results are very impresive. But, after computing the eigenstates, I need to calculate some other dynamical properties, thus I need to have all the eigenstates into a single matrix in a single process. Using pdgemr2d works fine for N up to 20000, but for larger values I get the error : xxmr2d core aborted. My cluster has 128GB RAM , 32 CPU/64 threads.

This is how I organized the code.

   .... 
   call MPI_INIT(info)
   call MPI_COMM_RANK(MPI_COMM_WORLD, myrank, info)
   call MPI_COMM_RANK(MPI_COMM_WORLD, nprocs, info)
   call MPI_COMM_SIZE(MPI_COMM_WORLD, nsize, info)
   nprow = 4
   npcol = 4
   block = 6
   call  blacs_get(-1,0,ictxt)
   call  blacs_gridinit(ictxt,"R",nprow, npcol)
   call  blacs_gridinfo(ictxt, nprow, npcol, myrow, mycol)
  ! print*, "myrank", myrank, myrow,mycol
   np = numroc(nst,block,myrow,0,nprow)
   nq = numroc(nst,block, mycol, 0,npcol)
   lld = max(1,np)
   call descinit(descAloc,nst,nst,block,block, 0,0, ictxt,np,info)
  ! print*, "INFO descinit A", info
  ! print*, "descAloc", descAloc
   call DESCINIT(descZ, nst, nst, block, block, 0, 0, ictxt,np, info)
   call DESCINIT(descZglobal, nst, nst, nst, nst, 0,0,ictxt,nst,info)
   call DESCINIT(descdumm, nst,nst,block,block,0,0,ictxt,np,info)
   !print*, "info descinit Z", info
   pi=acos(-1.d0) ! value of pi
   fi = 0.0/N_site
   phase =  dcmplx(cos(fi),sin(fi))
  
   !if  we  eigenvalues only 
   ! lworkk = 5*nst + max(3*block,block*(np+1))+1 
   !  if also eigenvectors 
    lworkk =  max( 1 + 6*nst + 2*np*nq, 3*nst + max( block*( np+1 ), 3*block)) + 2*nst
   allocate(Z(lld,lld), Aloc(lld,lld),workk(lworkk), wrk(np)) ! allocate memory to construct the full Hamiltonian matrix, j operator matrix for the calculated number of states in the Sz=0 sector 
      do local_i = 1 , np
        do local_j = 1 , nq
           ! Determine the global indicies (i,j) from the local indices   (local_i,local_j)
           i = indxl2g(local_i,block,myrow,0,nprow)
       j = indxl2g(local_j,block,mycol,0,npcol)
           ! Compute Hmatc(i,j) directly
           if (i == j) then 
           Aloc(local_i,local_j) = real(AD(i))
       else 
       Aloc(local_i, local_j) = 0.d0  ! Initialize to zero 
           do ih = row_ptr(i),row_ptr(i+1) - 1
          ij = col_ind(ih)
          isig = sign(1,col_ind_oper(ih))
              if (ij == j) then 
          if (isig == 1) then 
          Aloc(local_i,local_j) = Aloc(local_i,local_j) + 0.5*J_int
          else 
          Aloc(local_i,local_j) = Aloc(local_i,local_j) - 0.5*isig*J_int
              end if 
          end if 
           end do 
           end if 
        end do 
      end do 
   call pdsyev("v","u",nst,Aloc,1,1,descAloc,W,Z,1,1,descZ,workk,lworkk,info)
   if(myrank ==0) then 
   write(0,*) "info pdsyev", info 
   end if 
   if(myrank == 0) then
        write(0,*) "First and last eigenvalues :"
        write(0,"(6f5.2)") W(1),W(2),W(3),W(nst-2),W(nst-1),W(nst)
   write(0,*) "info pdsyev", info 
   end if     
   allocate(Zglobal(nst,nst)) ! global matrix to store all eigenvectors 
   call pdgemr2d(nst,nst, Z,1,1,descz, Zglobal, 1,1, descZglobal,ictxt)


   if(myrank == 0) then 
   deallocate(Zglobal)
   end if
    deallocate (Z)
    call blacs_gridexit(ictxt)
    call MPI_FINALIZE(info)

I compile the code as in following:

 mpiifx problem.f90 -qmkl=cluster
 mpirun -np 16 ./a.out

the error is : cannot allocate memory, aborted xxmr2d: out of memory

Upvotes: 0

Views: 115

Answers (0)

Related Questions