replace einsum by gemm in orth_projection

Author: Cstandardlib

source/module_hsolver/diago_bpcg.cpp

@@ -165,45 +165,45 @@ void DiagoBPCG<T, Device>::orth_projection(
 {
     ct::EinsumOption option(
         /*conj_x=*/false, /*conj_y=*/true, /*alpha=*/1.0, /*beta=*/0.0, /*Tensor out=*/&hsub_in);
-    hsub_in = ct::op::einsum("ij,kj->ik", grad_out, psi_in, option);
+    // hsub_in = ct::op::einsum("ij,kj->ik", grad_out, psi_in, option);
 
     // this->orth_projection(this->psi, this->hsub, this->grad);
-    // gemm: hsub_in(n_band x n_band) = grad_out^T(n_band x n_basis) * psi_in(n_basis x n_band)
-    // gemm_op()(this->ctx,
-    //             'C',
-    //             'N',
-    //             this->n_band,      //m
-    //             this->n_band,      //n
-    //             this->n_dim,       //k
-    //             this->one,
-    //             grad_out.data<T>(),
-    //             this->n_basis,     //lda
-    //             psi_in.data<T>(),
-    //             this->n_basis,     //ldb
-    //             this->zero,
-    //             hsub_in.data<T>(),
-    //             this->n_band);     //ldc
+    // gemm: hsub_in(n_band x n_band) = psi_in^T(n_band x n_basis) * grad_out(n_basis x n_band)
+    gemm_op()(this->ctx,
+                'C',
+                'N',
+                this->n_band,      //m
+                this->n_band,      //n
+                this->n_dim,       //k
+                this->one,
+                psi_in.data<T>(),
+                this->n_basis,     //lda
+                grad_out.data<T>(),
+                this->n_basis,     //ldb
+                this->zero,
+                hsub_in.data<T>(),
+                this->n_band);     //ldc
 
     // set_matrix_op()('L', hsub_in->data<T>(), this->n_band);
     option = ct::EinsumOption(
         /*conj_x=*/false, /*conj_y=*/false, /*alpha=*/-1.0, /*beta=*/1.0, /*Tensor out=*/&grad_out);
-    grad_out = ct::op::einsum("ij,jk->ik", hsub_in, psi_in, option);
+    // grad_out = ct::op::einsum("ij,jk->ik", hsub_in, psi_in, option);
 
     // grad_out(n_basis x n_band) = psi_in(n_basis x n_band) * hsub_in(n_band x n_band)
-    // gemm_op()(this->ctx,
-    //             'N',
-    //             'N',
-    //             this->n_basis,    //m
-    //             this->n_band,     //n
-    //             this->n_band,     //k
-    //             this->one,
-    //             psi_in.data<T>(),
-    //             this->n_basis,    //lda
-    //             hsub_in.data<T>(),
-    //             this->n_band,     //ldb
-    //             this->zero,
-    //             grad_out.data<T>(),
-    //             this->n_basis);  //ldc
+    gemm_op()(this->ctx,
+                'N',
+                'N',
+                this->n_dim,    //m
+                this->n_band,     //n
+                this->n_band,     //k
+                this->neg_one,
+                psi_in.data<T>(),
+                this->n_basis,    //lda
+                hsub_in.data<T>(),
+                this->n_band,     //ldb
+                this->one,
+                grad_out.data<T>(),
+                this->n_basis);  //ldc
 
     return;
 }

source/module_hsolver/test/diago_bpcg_test.cpp

@@ -249,7 +249,15 @@ TEST(DiagoBPCGTest, readH)
     // read Hamilt matrix from file data-H
     std::vector<std::complex<double>> hm;
     std::ifstream ifs;
-    ifs.open("H-KPoints-Si64.dat");
+    std::string filename = "H-KPoints-Si64.dat"; //"H-small-6x6.dat";
+    std::cout << "Reading file " << filename << std::endl;
+    ifs.open(filename);
+    // open file and check status
+    if (!ifs.is_open())
+    {
+        std::cout << "Error opening file " << filename << std::endl;
+        exit(1);
+    }
     DIAGOTEST::readh(ifs, hm);
     ifs.close();
     int dim = DIAGOTEST::npw;