diff --git a/CHANGELOG.md b/CHANGELOG.md
index 98cd0f40..cabfe501 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,8 @@
 
 - Reworked templates to include example tests.
 
+- Removed unnecessary full facotorization in the examples and made the input 1 based.
+
 ## Changes to Re::Solve in release 0.99.2
 
 ### Major Features
diff --git a/examples/gluRefactor.cpp b/examples/gluRefactor.cpp
index 6ba0f752..42c47f2a 100644
--- a/examples/gluRefactor.cpp
+++ b/examples/gluRefactor.cpp
@@ -163,7 +163,7 @@ int gluRefactor(int argc, char* argv[])
   vector_type* vec_x   = nullptr;
 
   RESOLVE_RANGE_PUSH(__FUNCTION__);
-  for (int i = 0; i < num_systems; ++i)
+  for (int i = 1; i <= num_systems; ++i)
   {
     std::cout << "System " << i << ":\n";
 
@@ -189,7 +189,7 @@ int gluRefactor(int argc, char* argv[])
       return -1;
     }
     bool is_expand_symmetric = true;
-    if (i == 0)
+    if (i == 1)
     {
       A       = io::createCsrFromFile(mat_file, is_expand_symmetric);
       vec_rhs = io::createVectorFromFile(rhs_file);
@@ -216,7 +216,7 @@ int gluRefactor(int argc, char* argv[])
 
     int status = 0;
 
-    if (i < 1)
+    if (i < 2)
     {
       RESOLVE_RANGE_PUSH("KLU");
       // Setup factorization solver
diff --git a/examples/gpuRefactor.cpp b/examples/gpuRefactor.cpp
index f8a4f2a9..ae6c366a 100644
--- a/examples/gpuRefactor.cpp
+++ b/examples/gpuRefactor.cpp
@@ -185,7 +185,7 @@ int gpuRefactor(int argc, char* argv[])
   vector_type* vec_x   = nullptr;
 
   RESOLVE_RANGE_PUSH(__FUNCTION__);
-  for (int i = 0; i < num_systems; ++i)
+  for (int i = 1; i <= num_systems; ++i)
   {
     std::cout << "System " << i << ":\n";
 
@@ -211,7 +211,7 @@ int gpuRefactor(int argc, char* argv[])
       return -1;
     }
     bool is_expand_symmetric = true;
-    if (i == 0)
+    if (i == 1)
     {
       A       = io::createCsrFromFile(mat_file, is_expand_symmetric);
       vec_rhs = io::createVectorFromFile(rhs_file);
@@ -238,7 +238,7 @@ int gpuRefactor(int argc, char* argv[])
 
     int status = 0;
 
-    if (i == 0)
+    if (i == 1)
     {
       RESOLVE_RANGE_PUSH("KLU");
       // Setup factorization solver
@@ -248,10 +248,6 @@ int gpuRefactor(int argc, char* argv[])
       // Analysis (symbolic factorization)
       status = KLU.analyze();
       std::cout << "KLU analysis status: " << status << std::endl;
-    }
-
-    if (i < 2)
-    {
       // Numeric factorization
       status = KLU.factorize();
       std::cout << "KLU factorization status: " << status << std::endl;
@@ -262,26 +258,22 @@ int gpuRefactor(int argc, char* argv[])
 
       // Print summary of results
       helper.printShortSummary(A, vec_rhs, vec_x);
-
-      if (i == 1)
+      // Extract factors and configure refactorization solver
+      matrix::Csr* L = (matrix::Csr*) KLU.getLFactor();
+      matrix::Csr* U = (matrix::Csr*) KLU.getUFactor();
+      if (L == nullptr || U == nullptr)
       {
-        // Extract factors and configure refactorization solver
-        matrix::Csr* L = (matrix::Csr*) KLU.getLFactor();
-        matrix::Csr* U = (matrix::Csr*) KLU.getUFactor();
-        if (L == nullptr || U == nullptr)
-        {
-          std::cout << "Factor extraction from KLU failed!\n";
-        }
-        index_type* P = KLU.getPOrdering();
-        index_type* Q = KLU.getQOrdering();
+        std::cout << "Factor extraction from KLU failed!\n";
+      }
+      index_type* P = KLU.getPOrdering();
+      index_type* Q = KLU.getQOrdering();
 
-        Rf.setup(A, L, U, P, Q, vec_rhs);
+      Rf.setup(A, L, U, P, Q, vec_rhs);
 
-        // Setup iterative refinement solver
-        if (is_iterative_refinement)
-        {
-          FGMRES.setup(A);
-        }
+      // Setup iterative refinement solver
+      if (is_iterative_refinement)
+      {
+        FGMRES.setup(A);
       }
       RESOLVE_RANGE_POP("KLU");
     }
diff --git a/examples/kluRefactor.cpp b/examples/kluRefactor.cpp
index 6834af11..5f68efde 100644
--- a/examples/kluRefactor.cpp
+++ b/examples/kluRefactor.cpp
@@ -130,7 +130,7 @@ int main(int argc, char* argv[])
   LinSolverDirectKLU*      KLU = new LinSolverDirectKLU;
   GramSchmidt              GS(&vector_handler, GramSchmidt::CGS2);
   LinSolverIterativeFGMRES FGMRES(&matrix_handler, &vector_handler, &GS);
-  for (int i = 0; i < num_systems; ++i)
+  for (int i = 1; i <= num_systems; ++i)
   {
     std::cout << "System " << i << ":\n";
 
@@ -153,7 +153,7 @@ int main(int argc, char* argv[])
       return 1;
     }
     bool is_expand_symmetric = true;
-    if (i == 0)
+    if (i == 1)
     {
       A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric);
 
@@ -172,15 +172,12 @@ int main(int argc, char* argv[])
     std::cout << "COO to CSR completed. Expanded NNZ: " << A->getNnz() << std::endl;
     // Now call direct solver
     int status;
-    if (i == 0)
+    if (i == 1)
     {
       vec_rhs->setDataUpdated(ReSolve::memory::HOST);
       KLU->setup(A);
       status = KLU->analyze();
       std::cout << "KLU analysis status: " << status << std::endl;
-    }
-    if (i < 2)
-    {
       status = KLU->factorize();
       std::cout << "KLU factorization status: " << status << std::endl;
     }
diff --git a/examples/sysRefactor.cpp b/examples/sysRefactor.cpp
index 8ce61918..ea09c583 100644
--- a/examples/sysRefactor.cpp
+++ b/examples/sysRefactor.cpp
@@ -223,7 +223,7 @@ int sysRefactor(int argc, char* argv[])
   }
 
   RESOLVE_RANGE_PUSH(__FUNCTION__);
-  for (int i = 0; i < num_systems; ++i)
+  for (int i = 1; i <= num_systems; ++i)
   {
     std::cout << "System " << i << ":\n";
     RESOLVE_RANGE_PUSH("File input");
@@ -250,7 +250,7 @@ int sysRefactor(int argc, char* argv[])
 
     // Refactorization is LU-based, so need to expand symmetric matrices
     bool is_expand_symmetric = true;
-    if (i == 0)
+    if (i == 1)
     {
       A       = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric);
       vec_rhs = ReSolve::io::createVectorFromFile(rhs_file);
@@ -281,7 +281,7 @@ int sysRefactor(int argc, char* argv[])
     printSystemInfo(matrix_pathname_full, A);
 
     // Now call direct solver
-    if (i == 0)
+    if (i == 1)
     {
       // Set matrix in solver after the initial matrix is loaded
       status = solver.setMatrix(A);
@@ -298,13 +298,6 @@ int sysRefactor(int argc, char* argv[])
       // Numeric factorization on the host
       status = solver.factorize();
       std::cout << "Numeric factorization on the host status: " << status << std::endl;
-    }
-    else if (i == 1)
-    {
-      // Numeric factorization on the host
-      status = solver.factorize();
-      std::cout << "Numeric factorization on the host status: " << status << std::endl;
-
       // Set up refactorization solver
       status = solver.refactorizationSetup();
       std::cout << "Refactorization setup status: " << status << std::endl;