Diff

Index: hs_operators.py===================================================================--- hs_operators.py	(revision 5853)+++ hs_operators.py	(working copy)@@ -39,7 +39,7 @@             self.async = async         if hermitian is not None:             self.hermitian = hermitian-        self.bmd = BandMatrixDescriptor(bd, gd) #XXX prefix Blacs for 1D layout+        self.bmd = BlacsBandMatrixDescriptor(bd, gd) #XXX      def allocate_work_arrays(self, dtype):         """This is a little complicated, but let's look at the facts.@@ -74,7 +74,7 @@         ngroups = self.bd.comm.size         mynbands = self.bd.mynbands         nbands = self.bd.nbands-        if ngroups == 1 and self.nblocks == 1:+        if False: #XXX ngroups == 1 and self.nblocks == 1:             self.work1_xG = self.gd.zeros(mynbands, dtype)         else:             assert mynbands % self.nblocks == 0@@ -83,7 +83,7 @@                 X += int(np.ceil(mynbands/self.gd.n_c.prod()))             self.work1_xG = self.gd.zeros(X, dtype)             self.work2_xG = self.gd.zeros(X, dtype)-            if ngroups > 1:+            if True: #XXX ngroups > 1:                 if self.hermitian:                     Q = ngroups // 2 + 1                 else:@@ -260,7 +260,7 @@         N = self.bd.mynbands         B = self.bd.comm.size -        if B == 1 and J == 1:+        if False: #XXX B == 1 and J == 1:             return self.work1_xG         else:             assert N % J == 0, "Can't divide %d bands in %d blocks." % (N,J)@@ -318,7 +318,7 @@                 # dA denotes dA_aii as usual                 dAP_ani[a] = np.dot(P_ni, dA[a])         -        if B == 1 and J == 1:+        if False: #XXX B == 1 and J == 1:             # Simple case:             Apsit_nG = A(psit_nG)             self._pseudo_braket(psit_nG, Apsit_nG, A_NN)@@ -338,7 +338,7 @@         M = N // J          # Buffer for storage of blocks of calculated matrix elements.-        if B == 1:+        if False: #XXX B == 1:             A_qnn = A_NN.reshape((1, N, N))         else:             A_qnn = self.A_qnn@@ -404,7 +404,7 @@          domain_comm.sum(A_qnn, 0) -        if B == 1:+        if False: #XXX B == 1:             return A_NN          if domain_comm.rank == 0:@@ -455,7 +455,7 @@          C_NN = self.bmd.redistribute_input(C_NN) -        if B == 1 and J == 1:+        if False: #XXX B == 1 and J == 1:             # Simple case:             newpsit_nG = self.work1_xG             gemm(1.0, psit_nG, C_NN, 0.0, newpsit_nG)Index: test/parallel/ut_hsops.py===================================================================--- test/parallel/ut_hsops.py	(revision 5853)+++ test/parallel/ut_hsops.py	(working copy)@@ -522,6 +522,62 @@          self.check_and_plot(S_nn, alpha*self.S0_nn, 9, 'overlaps,nonhermitian') +    def test_scalapack_column_layout(self): #XXX XXX XXX+        # Known starting point of S_nn = <psit_m|S|psit_n>+        S_nn = self.S0_nn++        # Eigenvector decomposition S_nn = V_nn * W_nn * V_nn^dag+        # Utilize the fact that they are analytically known (cf. Maple)+        band_indices = np.arange(self.nbands)+        V_nn = np.eye(self.nbands).astype(self.dtype)+        if self.dtype == complex:+            V_nn[1:,1] = np.conj(self.gamma)**band_indices[1:] * band_indices[1:]**0.5+            V_nn[1,2:] = -self.gamma**band_indices[1:-1] * band_indices[2:]**0.5+        else:+            V_nn[2:,1] = band_indices[2:]**0.5+            V_nn[1,2:] = -band_indices[2:]**0.5++        W_n = np.zeros(self.nbands).astype(self.dtype)+        W_n[1] = (1. + self.Qtotal) * self.nbands * (self.nbands - 1) / 2.++        # Find the inverse basis+        Vinv_nn = np.linalg.inv(V_nn)++        # Test analytical eigenvectors for consistency against analytical S_nn+        D_nn = np.dot(Vinv_nn, np.dot(S_nn, V_nn))+        self.assertAlmostEqual(np.abs(D_nn.diagonal()-W_n).max(), 0, 8)+        self.assertAlmostEqual(np.abs(np.tril(D_nn, -1)).max(), 0, 4)+        self.assertAlmostEqual(np.abs(np.triu(D_nn, 1)).max(), 0, 4)+        del Vinv_nn, D_nn++        # Set up Hermitian overlap operator:+        S = lambda x: x+        dS = lambda a, P_ni: np.dot(P_ni, self.setups[a].dO_ii)+        nblocks = self.get_optimal_number_of_blocks(self.blocking)+        overlap = MatrixOperator(self.bd, self.gd, nblocks, self.async, False) # True)+        overlap.bmd.redistribute_output = lambda A_Nn: A_Nn # XXX override full matrix assembly+        S_Nn = overlap.calculate_matrix_elements(self.psit_nG, self.P_ani, S, dS)++        if self.bd.comm.rank == 0:+            self.gd.comm.broadcast(S_Nn, 0)++        self.assertEqual(S_Nn.shape, (self.bd.nbands,self.bd.mynbands))++        from gpaw.utilities.blacs import scalapack_set +        from gpaw.blacs import BlacsGrid, Redistributor, parallelprint, \+            BlacsBandDescriptor++        mcpus, ncpus, blocksize = 2, 2, 1+        bbd = BlacsBandDescriptor(world, self.gd, self.bd, self.kpt_comm, mcpus, ncpus, blocksize)++        # We would create C_nN in the real-space code this way.+        C_nN = np.empty((self.bd.mynbands, self.bd.nbands), dtype=S_Nn.dtype)+        diagonalizer = bbd.get_diagonalizer()+        eps_n = np.zeros(self.bd.mynbands) #, dtype=S_Nn.dtype) # XXX dtype?+        diagonalizer.diagonalize(S_Nn, C_nN, eps_n)+        if self.bd.comm.rank == 0 and self.gd.comm.rank == 0:+            print world.rank, 'eps_n:', eps_n, 'ref_n:', W_n+     def test_multiply_orthonormal(self):         # Known starting point of S_nn = <psit_m|S|psit_n>         S_nn = self.S0_nn@@ -716,7 +772,7 @@      testcases = []     for dtype in [float, complex]:-        for parstride_bands in [False, True]:+        for parstride_bands in [False]: #[False, True]:             for blocking in ['fast', 'best']: # 'light'                 for async in [False, True]:                     testcases.append(UTConstantWavefunctionFactory(dtype, \