From b575cda8d69b39a691c4204653b2f6fa714358b3 Mon Sep 17 00:00:00 2001
From: Pablo Oliveira <pablo@sifflez.org>
Date: Tue, 9 Feb 2021 14:44:54 +0100
Subject: [PATCH] Optimize out Identity matrix

---
 SM_MaponiA3.cpp | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/SM_MaponiA3.cpp b/SM_MaponiA3.cpp
index 7b6a11d..6e77d95 100644
--- a/SM_MaponiA3.cpp
+++ b/SM_MaponiA3.cpp
@@ -8,7 +8,6 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int
 
     unsigned int k, l, lbar, i, j, tmp = M;
     unsigned int *p = new unsigned int[M+1];
-    double *Id = new double[M*M];
     double alpha, beta;
     double *breakdown = new double[M+1];
     double *Al = new double[M*M];
@@ -26,14 +25,6 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int
         }
     }
 
-    // Initialize identity matrix
-    for (i = 0; i < M; i++) {
-        for (j = 0; j < M; j++) {
-            if (i != j) Id[i*M+j] = 0;
-            else Id[i*M+j] = 1;
-        }
-    }
-
     // Initialize ylk with zeros
     for (l = 0; l < M; l++) {
         for (k = 0; k < M+1; k++) {
@@ -80,14 +71,21 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int
     // pointer 'copy' that points to whereever 'Slater_inv' points to now.
     double *copy  = Slater_inv;
 
+    double *U = new double[M*M];
     // Construct A-inverse from A0-inverse and the ylk
     for (l = 0; l < M; l++) {
         k = l+1;
-        double * U = outProd(ylk[l][p[k]], (Id + (p[k]-1)*M), M);
+        for (unsigned int i = 0; i < M; i++) {
+          for (unsigned int j = 0; j < M; j++) {
+            U[i*M+j] = ylk[l][p[k]][i+1] * ((p[k]-1) == j);
+          }
+        }
+
+
         beta = 1 + ylk[l][p[k]][p[k]];
         for (i = 0; i < M; i++) {
             for (j = 0; j < M; j++) {
-                Al[i*M+j] = Id[i*M+j] - U[i*M+j] / beta;
+                Al[i*M+j] = (i == j) - U[i*M+j] / beta;
             }
         }
         Slater_inv = matMul(Al, Slater_inv, M);
@@ -106,7 +104,7 @@ void MaponiA3(double *Slater0, double *Slater_inv, unsigned int M, unsigned int
         }
         delete [] ylk[l];
     }
-    delete [] Id, Al;
+    delete [] Al, U;
     delete [] p, breakdown;
 }