Exercise 11 problem 1

ex11/problem1/main.cpp

+/* 
+ * Programming techniques for scientific simulation: Problem 11.1
+ * 
+ * */
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include "main.hpp"
+
+int main() {
+  // Open file
+  std::ofstream myfile;
+  // Initialize inputs to LAPACK function
+  int const N = 10;
+  int const K = 2;
+  int const m = 4;
+  matrix_t A = buildMatrix(N, K, m);
+  double w[N];
+  double wkopt;
+  int lwork = -1;
+  int info;
+  // Call function to get optimal lwork (written inside wkopt)
+  dsyev_('V', 'U', N, A.data(), N, w, &wkopt, lwork, info);
+  lwork = (int)wkopt;  // cast value to integer
+  matrix_t work(lwork);  // create array with length lwork
+  dsyev_('V', 'U', N, A.data(), N, w, work.data(), lwork, info);
+  // Compute theoretical eigenvalues
+  matrix_t ev = computeEV(N, K, m);
+  // Open file
+  myfile.open("results.txt");
+  // Print results
+  if (info  > 0) {
+    myfile << "The algorithm failed to compute eigenvalues\n";
+  }
+  else if (info < 0) {
+    myfile << "The input " << info << " has an illegal value\n";
+  }
+  myfile << "The eigenvalues are: \n";
+  for (int i = 0; i < N; i++) {
+    myfile << std::sqrt(w[i]) << "  ";
+  }
+  myfile << "\n-----------------------------------\n";
+  myfile << "The theoretical eigenvalues are: \n";
+  for (int i = 0; i < N; i++) {
+    myfile << ev[i] << "  ";
+  }
+  myfile << "\n-----------------------------------\n";
+  myfile << "The eigenvectors are: \n";
+  for (int j = 0; j < N; j++) {
+    myfile << "(  ";
+    for (int i = 0; i < N; i++) {
+      myfile << A[i + j*N] << "  ";
+    }
+    myfile << ")\n";
+  }
+  myfile.close();
+  return 0;
+} 

ex11/problem1/main.hpp

+/* HEADER FILE */
+
+#ifndef __MAIN_H__
+#define __MAIN_H__
+
+#include <vector>
+#include <cmath>
+#include <numbers>
+
+using value_t = double;
+using matrix_t = std::vector<value_t>;
+
+matrix_t buildMatrix(
+  const unsigned int N,
+  const double K,
+  const double m
+  ) {
+  // Builds a matrix of size N*N where the
+  // diagonal values are 2 and the values next 
+  // to the diagonal are -1
+  matrix_t M(N * N);
+  double tmp = K/m;
+  for (unsigned int j = 0; j < N; j++) {
+    for (unsigned int i = 0; i < N; i++) {
+      if (i == j) { M[i + j * N] = 2*tmp; }
+      else if ((j == i + 1) | (j == i - 1)) { M[i + j * N] = -tmp; }
+    }
+  }
+  return M;
+}
+
+matrix_t computeEV(int const unsigned N, double const K, double const m) {
+  matrix_t ev(N);
+  for (int unsigned i = 1; i <= N; i++) {
+    ev[i-1] = std::sqrt((K/m)*(2 - 2*std::cos((M_PI * i)/(N + 1))));
+  }
+  return ev;
+}
+
+extern "C" double dsyev_(
+    char const   & JOBZ,
+    char const   & UPLO,
+    int const    & N,
+    double       * A,
+    int const    & LDA,
+    double       * W,
+    double const * WORK,
+    int const    & LWORK,
+    int const    & INFO
+);
+
+
+#endif

ex11/problem1/results.txt

+The eigenvalues are: 
+0.201264  0.39843  0.587486  0.764582  0.926113  1.06879  1.18971  1.28641  1.35693  1.39982  
+-----------------------------------
+The theoretical eigenvalues are: 
+0.201264  0.39843  0.587486  0.764582  0.926113  1.06879  1.18971  1.28641  1.35693  1.39982  
+-----------------------------------
+The eigenvectors are: 
+(  -0.120131  -0.23053  -0.322253  -0.387868  -0.422061  -0.422061  -0.387868  -0.322253  -0.23053  -0.120131  )
+(  0.23053  0.387868  0.422061  0.322253  0.120131  -0.120131  -0.322253  -0.422061  -0.387868  -0.23053  )
+(  -0.322253  -0.422061  -0.23053  0.120131  0.387868  0.387868  0.120131  -0.23053  -0.422061  -0.322253  )
+(  0.387868  0.322253  -0.120131  -0.422061  -0.23053  0.23053  0.422061  0.120131  -0.322253  -0.387868  )
+(  -0.422061  -0.120131  0.387868  0.23053  -0.322253  -0.322253  0.23053  0.387868  -0.120131  -0.422061  )
+(  0.422061  -0.120131  -0.387868  0.23053  0.322253  -0.322253  -0.23053  0.387868  0.120131  -0.422061  )
+(  0.387868  -0.322253  -0.120131  0.422061  -0.23053  -0.23053  0.422061  -0.120131  -0.322253  0.387868  )
+(  0.322253  -0.422061  0.23053  0.120131  -0.387868  0.387868  -0.120131  -0.23053  0.422061  -0.322253  )
+(  0.23053  -0.387868  0.422061  -0.322253  0.120131  0.120131  -0.322253  0.422061  -0.387868  0.23053  )
+(  -0.120131  0.23053  -0.322253  0.387868  -0.422061  0.422061  -0.387868  0.322253  -0.23053  0.120131  )