1c----------------------------------------------------------------------- 2c 3 subroutine wfn1_lindep(g_s,g_s12,g_v,nindep,tol) 4 implicit none 5C> 6C> \brief Diagonalize the overlap matrix and find from the eigenvalues 7C> how many linearly independent vectors there are 8C> 9#include "errquit.fh" 10#include "stdio.fh" 11#include "mafdecls.fh" 12#include "global.fh" 13c 14 integer g_s ! the GA handle for the overlap matrix 15 integer g_s12 ! the GA handle for the S^{-1/2} matrix 16 integer g_v ! the GA handle for the eigenvectors 17c 18 integer, intent(out) :: nindep ! the number of linearly 19 ! independent vectors 20 double precision, intent(in) :: tol ! the eigenvalue cutoff 21c 22 integer itype 23 integer ndim 24 integer idims(ga_max_dim) 25 integer nbf ! the number of basis functions 26 integer ii ! counter 27 integer jj ! counter 28 integer kk ! counter 29 integer nproc ! the number of processor ranks 30 integer iproc ! the current processor rank 31 integer ilo,ihi ! row limits on the local block 32 integer jlo,jhi ! column limits on the local block 33 integer islo,ishi ! row limits on the local block 34 integer jslo,jshi ! column limits on the local block 35 integer k_v ! memory location for the local block of v 36 integer ldv ! leading dimension of the local block of v 37c 38 double precision :: dnrm 39 double precision, allocatable :: eigval(:) 40 double precision, allocatable :: vi(:,:) 41 double precision, allocatable :: vj(:,:) 42 double precision, allocatable :: s12(:,:) 43c 44 character*(12) pname 45 pname = "wfn1_lindep:" 46c 47 nproc = ga_nnodes() 48 iproc = ga_nodeid() 49c 50 call nga_inquire(g_s,itype,ndim,idims) 51 nbf = idims(1) 52 if (ndim.ne.2) then 53 call errquit(pname//" wrong number of dimensions",ndim,UERR) 54 endif 55 if (idims(1).ne.idims(2)) then 56 call errquit(pname//" non-square overlap matrix", 57 & idims(2)-idims(1),UERR) 58 endif 59c 60 allocate(eigval(1:nbf)) 61c 62#ifdef PARALLEL_DIAG 63#ifdef SCALAPACK 64 call ga_pdsyev(g_s, g_v, eigval, 0) 65#else 66 call ga_diag_std(g_s, g_v, eigval) 67#endif 68#else 69 call ga_diag_std_seq(g_s, g_v, eigval) 70#endif 71c 72 nindep = 0 73 do ii = 1, nbf 74 if (eigval(ii).ge.tol) nindep = nindep + 1 75 enddo 76 if (nbf-nindep.gt.0.and.iproc.eq.0) then 77 write(luout,'(" !! The overlap matrix has ",i5, 78 & " vectors deemed linearly dependent with"/ 79 & " eigenvalues:")')nbf-nindep 80 write(luout,'(1p,8d9.2)')(eigval(ii),ii=0,nbf-nindep-1) 81 write(luout,*) 82 endif 83c 84c compute S^{1/2} matrix 85c 86 call ga_distribution(g_s12,iproc,islo,ishi,jslo,jshi) 87 allocate(vi(islo:ishi,nbf)) 88 allocate(vj(jslo:jshi,nbf)) 89 allocate(s12(islo:ishi,jslo:jshi)) 90 call ga_get(g_v,islo,ishi,1,nbf,vi,ishi-islo+1) 91 call ga_get(g_v,jslo,jshi,1,nbf,vj,jshi-jslo+1) 92 s12 = 0.0d0 93 do kk = 1, nbf 94 do jj = jslo, jshi 95 do ii = islo, ishi 96 s12(ii,jj) = s12(ii,jj) 97 & + vi(ii,kk)*vj(jj,kk)*sqrt(eigval(kk)) 98 enddo ! ii 99 enddo ! jj 100 enddo ! kk 101 call ga_put(g_s12,islo,ishi,jslo,jshi,s12,ishi-islo+1) 102 deallocate(vi,vj,s12) 103c 104c scale eigenvectors to obtain S^{-1/2} 105c 106 call ga_distribution(g_v,iproc,ilo,ihi,jlo,jhi) 107 call ga_access(g_v,ilo,ihi,jlo,jhi,k_v,ldv) 108 do jj = 0, jhi-jlo 109 dnrm = 1.0d0/sqrt(eigval(jlo+jj)) 110 do ii = 0, ihi-ilo 111 dbl_mb(k_v+jj*ldv+ii) = dbl_mb(k_v+jj*ldv+ii)*dnrm 112 enddo 113 enddo 114c 115 deallocate(eigval) 116c 117 end 118c 119c----------------------------------------------------------------------- 120