naromero

Index: hs_operators.py===================================================================--- hs_operators.py	(revision 5859)+++ 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 5859)+++ test/parallel/ut_hsops.py	(working copy)@@ -42,7 +42,7 @@     Setup a simple band parallel calculation."""      # Number of bands-    nbands = 36+    nbands = 12      # Spin-paired, single kpoint     nspins = 1@@ -434,7 +434,7 @@      # ================================= -    def test_wavefunction_content(self):+    def Xtest_wavefunction_content(self):         # Integrate diagonal brakets of pseudo wavefunctions         gpts_c = self.gd.get_size_of_global_array() @@ -453,7 +453,7 @@         intpsit_n = self.bd.collect(intpsit_myn, broadcast=True)         self.assertAlmostEqual(np.abs(intpsit_n-np.arange(self.nbands)).max(), 0, 9) -    def test_projection_content(self):+    def Xtest_projection_content(self):         # Distribute inverse effective charges to everybody in domain         all_Qeff_a = np.empty(len(self.atoms), dtype=float)         for a,rank in enumerate(self.rank_a):@@ -477,7 +477,7 @@         if all_fingerprints.ptp(0).any():             raise RuntimeError('Distributed eff. charges are not identical!') -    def test_overlaps_hermitian(self):+    def Xtest_overlaps_hermitian(self):         # Set up Hermitian overlap operator:         S = lambda x: x         dS = lambda a, P_ni: np.dot(P_ni, self.setups[a].dO_ii)@@ -499,7 +499,7 @@          self.check_and_plot(S_nn, self.S0_nn, 9, 'overlaps,hermitian') -    def test_overlaps_nonhermitian(self):+    def Xtest_overlaps_nonhermitian(self):         alpha = np.random.normal(size=1).astype(self.dtype)         if self.dtype == complex:             alpha += 1j*np.random.normal(size=1)@@ -522,10 +522,71 @@          self.check_and_plot(S_nn, alpha*self.S0_nn, 9, 'overlaps,nonhermitian') -    def test_multiply_orthonormal(self):+    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+        from gpaw.blacs import parallelprint+        # parallelprint(world, S_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, True) # 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)+        print 'after calc matrix element'+        parallelprint(world, S_Nn)+        # if self.bd.comm.rank == 0:+        # self.gd.comm.broadcast(S_Nn, 0)++        print 'after broadcast'+        parallelprint(world, S_Nn)+        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, 6+        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) +        diagonalizer.diagonalize(S_Nn, C_nN, eps_n)++        if self.gd.comm.rank == 0:+            print world.rank, 'eps_n:', eps_n, 'ref_n:', W_n++    def Xtest_multiply_orthonormal(self):+        # 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)@@ -584,7 +645,7 @@         self.assertAlmostEqual(np.abs(D0_nn-np.diag(W_n)).max(), 0, 9)         self.check_and_plot(D_nn, D0_nn, 9, 'multiply,orthonormal') -    def test_multiply_randomized(self):+    def Xtest_multiply_randomized(self):         # Known starting point of S_nn = <psit_m|S|psit_n>         S_nn = self.S0_nn @@ -620,7 +681,7 @@         D0_nn = np.dot(C_nn.T.conj(), np.dot(S_nn, C_nn))         self.check_and_plot(D_nn, D0_nn, 9, 'multiply,randomized') -    def test_multiply_nonhermitian(self):+    def Xtest_multiply_nonhermitian(self):         alpha = np.random.normal(size=1).astype(self.dtype)         if self.dtype == complex:             alpha += 1j*np.random.normal(size=1)@@ -715,10 +776,10 @@     parinfo = np.unique(np.sort(parinfo)).tolist()      testcases = []-    for dtype in [float, complex]:-        for parstride_bands in [False, True]:-            for blocking in ['fast', 'best']: # 'light'-                for async in [False, True]:+    for dtype in [float]:+        for parstride_bands in [False]: #[False, True]:+            for blocking in ['fast']: # 'light'+                for async in [False]:                     testcases.append(UTConstantWavefunctionFactory(dtype, \                         parstride_bands, blocking, async))