Merge pull request #14 from kvedala/gitpod

Added Gitpod support - online IDE to create, edit and pull programs, build and test run

.github/pull_request_template.md

@@ -3,17 +3,17 @@
 Thank you for your Pull Request. Please provide a description above and review
 the requirements below.
 
-Contributors guide: https://github.com/TheAlgorithms/C-Plus-Plus/CONTRIBUTION.md
+Contributors guide: https://github.com/TheAlgorithms/C-Plus-Plus/CONTRIBUTING.md
 -->
 
 #### Checklist
 <!-- Remove items that do not apply. For completed items, change [ ] to [x]. -->
 
 - [ ] Added description of change
-- [ ] Added file name matches [File name guidelines](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTION.md#New-File-Name-guidelines)
+- [ ] Added file name matches [File name guidelines](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTING.md#New-File-Name-guidelines)
 - [ ] Added tests and example, test must pass
 - [ ] Relevant documentation/comments is changed or added
-- [ ] PR title follows semantic [commit guidelines](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTION.md#Commit-Guidelines)
+- [ ] PR title follows semantic [commit guidelines](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTING.md#Commit-Guidelines)
 - [ ] Search previous suggestions before making a new one, as yours may be a duplicate.
 - [ ] Sort by alphabetical order
 - [ ] I acknowledge that all my contributions will be made under the project's license.

.gitpod.yml

@@ -0,0 +1,6 @@
+tasks:
+- init: >
+  pip install cpplint
+  cmake -S . -B build -DCMAKE_BUILD_TYPE=Debug -DCMAKE_EXPORT_COMPILE_COMMANDS=ON &&
+  cmake --build build -t all 
+  

README.md

@@ -1,5 +1,7 @@
 # The Algorithms - C # {#mainpage}
-[![contributions welcome](https://img.shields.io/static/v1.svg?label=Contributions&message=Welcome&color=0059b3&style=flat-square)](https://github.com/kvedala/C-Plus-Plus/blob/master/CONTRIBUTION.md) 
+[![Gitpod Ready-to-Code](https://img.shields.io/badge/Gitpod-Ready--to--Code-blue?logo=gitpod)](https://gitpod.io/#https://github.com/TheAlgorithms/C-Plus-Plus) 
+[![Gitter chat](https://img.shields.io/badge/Chat-Gitter-ff69b4.svg?label=Chat&logo=gitter&style=flat-square)](https://gitter.im/TheAlgorithms)
+[![contributions welcome](https://img.shields.io/static/v1.svg?label=Contributions&message=Welcome&color=0059b3&style=flat-square)](https://github.com/kvedala/C-Plus-Plus/blob/master/CONTRIBUTING.md) 
 ![GitHub repo size](https://img.shields.io/github/repo-size/kvedala/C-Plus-Plus?color=red&style=flat-square)
 ![GitHub closed pull requests](https://img.shields.io/github/issues-pr-closed/kvedala/C?color=green&style=flat-square)
 ![Doxygen CI](https://github.com/kvedala/C/workflows/Doxygen%20CI/badge.svg)
@@ -15,4 +17,4 @@ All the code can be executed and tested online: [![using Google Colab Notebook](
 The implementations are for learning purpose. They may be less efficient than the implementation in the standard library.
 
 ### Contribute Guidelines
-Read our [Contribution Guidelines](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTION.md) before you contribute.
+Read our [Contribution Guidelines](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTING.md) before you contribute.

numerical_methods/qr_eigen_values.c

@@ -6,12 +6,17 @@
  * \author [Krishna Vedala](https://github.com/kvedala)
  */
 #include "qr_decompose.h"
+#include <assert.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <time.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
 
-#define LIMS 9 /**< limit of range of matrix values */
+#define LIMS 9        /**< limit of range of matrix values */
+#define EPSILON 1e-10 /**< accuracy tolerance limit */
 
 /**
  * create a square matrix of given size with random elements
@@ -21,8 +26,10 @@
 void create_matrix(double **A, int N)
 {
     int i, j, tmp, lim2 = LIMS >> 1;
-    srand(time(NULL));
 
+#ifdef _OPENMP
+#pragma omp for
+#endif
     for (i = 0; i < N; i++)
     {
         A[i][i] = (rand() % LIMS) - lim2;
@@ -56,7 +63,12 @@ double **mat_mul(double **A, double **B, double **OUT, int R1, int C1, int R2,
         perror("Matrix dimensions mismatch!");
         return OUT;
     }
-    for (int i = 0; i < R1; i++)
+
+    int i;
+#ifdef _OPENMP
+#pragma omp for
+#endif
+    for (i = 0; i < R1; i++)
         for (int j = 0; j < C2; j++)
         {
             OUT[i][j] = 0.f;
@@ -66,14 +78,240 @@ double **mat_mul(double **A, double **B, double **OUT, int R1, int C1, int R2,
     return OUT;
 }
 
+/** Compute eigen values using iterative shifted QR decomposition algorithm as
+ * follows:
+ * 1. Use last diagonal element of A as eigen value approximation \f$c\f$
+ * 2. Shift diagonals of matrix \f$A' = A - cI\f$
+ * 3. Decompose matrix \f$A'=QR\f$
+ * 4. Compute next approximation \f$A'_1 = RQ \f$
+ * 5. Shift diagonals back \f$A_1 = A'_1 + cI\f$
+ * 6. Termination condition check: last element below diagonal is almost 0
+ *   1. If not 0, go back to step 1 with the new approximation \f$A_1\f$
+ *   2. If 0, continue to step 7
+ * 7. Save last known \f$c\f$ as the eigen value.
+ * 8. Are all eigen values found?
+ *   1. If not, remove last row and column of \f$A_1\f$ and go back to step 1.
+ *   2. If yes, stop.
+ *
+ * \note The matrix \f$A\f$ gets modified
+ *
+ * \param[in,out] A matrix to compute eigen values for
+ * \param[out] eig_vals resultant vector containing computed eigen values
+ * \param[in] mat_size matrix size
+ * \param[in] debug_print 1 to print intermediate Q & R matrices, 0 for not to
+ * \returns time for computation in seconds
+ */
+double eigen_values(double **A, double *eigen_vals, int mat_size,
+                    char debug_print)
+{
+    double **R = (double **)malloc(sizeof(double *) * mat_size);
+    double **Q = (double **)malloc(sizeof(double *) * mat_size);
+
+    if (!eigen_vals)
+    {
+        perror("Output eigen value vector cannot be NULL!");
+        return -1;
+    }
+    else if (!Q || !R)
+    {
+        perror("Unable to allocate memory for Q & R!");
+        return -1;
+    }
+
+    /* allocate dynamic memory for matrices */
+    for (int i = 0; i < mat_size; i++)
+    {
+        R[i] = (double *)malloc(sizeof(double) * mat_size);
+        Q[i] = (double *)malloc(sizeof(double) * mat_size);
+        if (!Q[i] || !R[i])
+        {
+            perror("Unable to allocate memory for Q & R.");
+            return -1;
+        }
+    }
+
+    if (debug_print)
+        print_matrix(A, mat_size, mat_size);
+
+    int rows = mat_size, columns = mat_size;
+    int counter = 0, num_eigs = rows - 1;
+    double last_eig = 0;
+
+    clock_t t1 = clock();
+    while (num_eigs > 0) /* continue till all eigen values are found */
+    {
+        /* iterate with QR decomposition */
+        while (fabs(A[num_eigs][num_eigs - 1]) > EPSILON)
+        {
+            last_eig = A[num_eigs][num_eigs];
+            for (int i = 0; i < rows; i++)
+                A[i][i] -= last_eig; /* A - cI */
+            qr_decompose(A, Q, R, rows, columns);
+
+            if (debug_print)
+            {
+                print_matrix(A, rows, columns);
+                print_matrix(Q, rows, columns);
+                print_matrix(R, columns, columns);
+                printf("-------------------- %d ---------------------\n",
+                       ++counter);
+            }
+
+            mat_mul(R, Q, A, columns, columns, rows, columns);
+            for (int i = 0; i < rows; i++)
+                A[i][i] += last_eig; /* A + cI */
+        }
+
+        /* store the converged eigen value */
+        eigen_vals[num_eigs] = last_eig;
+
+        if (debug_print)
+        {
+            printf("========================\n");
+            printf("Eigen value: % g,\n", last_eig);
+            printf("========================\n");
+        }
+
+        num_eigs--;
+        rows--;
+        columns--;
+    }
+    eigen_vals[0] = A[0][0];
+    double dtime = (double)(clock() - t1) / CLOCKS_PER_SEC;
+
+    if (debug_print)
+    {
+        print_matrix(R, mat_size, mat_size);
+        print_matrix(Q, mat_size, mat_size);
+    }
+
+    /* cleanup dynamic memory */
+    for (int i = 0; i < mat_size; i++)
+    {
+        free(R[i]);
+        free(Q[i]);
+    }
+    free(R);
+    free(Q);
+
+    return dtime;
+}
+
+/**
+ * test function to compute eigen values of a 2x2 matrix
+ * \f[\begin{bmatrix}
+ * 5 & 7\\
+ * 7 & 11
+ * \end{bmatrix}\f]
+ * which are approximately, {15.56158, 0.384227}
+ */
+void test1()
+{
+    int mat_size = 2;
+    double X[][2] = {{5, 7}, {7, 11}};
+    double y[] = {15.56158, 0.384227}; // corresponding y-values
+    double eig_vals[2];
+
+    // The following steps are to convert a "double[][]" to "double **"
+    double **A = (double **)malloc(mat_size * sizeof(double *));
+    for (int i = 0; i < mat_size; i++)
+        A[i] = X[i];
+
+    printf("------- Test 1 -------\n");
+
+    double dtime = eigen_values(A, eig_vals, mat_size, 0);
+
+    for (int i = 0; i < mat_size; i++)
+    {
+        printf("%d/5 Checking for %.3g --> ", i + 1, y[i]);
+        char result = 0;
+        for (int j = 0; j < mat_size && !result; j++)
+        {
+            if (fabs(y[i] - eig_vals[j]) < 0.1)
+            {
+                result = 1;
+                printf("(%.3g) ", eig_vals[j]);
+            }
+        }
+
+        // ensure that i^th expected eigen value was computed
+        assert(result != 0);
+        printf("found\n");
+    }
+    printf("Test 1 Passed in %.3g sec\n\n", dtime);
+    free(A);
+}
+
+/**
+ * test function to compute eigen values of a 2x2 matrix
+ * \f[\begin{bmatrix}
+ * -4& 4& 2& 0& -3\\
+ * 4& -4& 4& -3& -1\\
+ * 2& 4& 4& 3& -3\\
+ * 0& -3& 3& -1&-1\\
+ * -3& -1& -3& -3& 0
+ * \end{bmatrix}\f]
+ * which are approximately, {9.27648, -9.26948, 2.0181, -1.03516, -5.98994}
+ */
+void test2()
+{
+    int mat_size = 5;
+    double X[][5] = {{-4, 4, 2, 0, -3},
+                     {4, -4, 4, -3, -1},
+                     {2, 4, 4, 3, -3},
+                     {0, -3, 3, -1, -3},
+                     {-3, -1, -3, -3, 0}};
+    double y[] = {9.27648, -9.26948, 2.0181, -1.03516,
+                  -5.98994}; // corresponding y-values
+    double eig_vals[5];
+
+    // The following steps are to convert a "double[][]" to "double **"
+    double **A = (double **)malloc(mat_size * sizeof(double *));
+    for (int i = 0; i < mat_size; i++)
+        A[i] = X[i];
+
+    printf("------- Test 2 -------\n");
+
+    double dtime = eigen_values(A, eig_vals, mat_size, 0);
+
+    for (int i = 0; i < mat_size; i++)
+    {
+        printf("%d/5 Checking for %.3g --> ", i + 1, y[i]);
+        char result = 0;
+        for (int j = 0; j < mat_size && !result; j++)
+        {
+            if (fabs(y[i] - eig_vals[j]) < 0.1)
+            {
+                result = 1;
+                printf("(%.3g) ", eig_vals[j]);
+            }
+        }
+
+        // ensure that i^th expected eigen value was computed
+        assert(result != 0);
+        printf("found\n");
+    }
+    printf("Test 2 Passed in %.3g sec\n\n", dtime);
+    free(A);
+}
+
 /**
  * main function
  */
 int main(int argc, char **argv)
 {
+    srand(time(NULL));
+
     int mat_size = 5;
     if (argc == 2)
         mat_size = atoi(argv[1]);
+    else
+    { // if invalid input argument is given run tests
+        test1();
+        test2();
+        printf("Usage: ./qr_eigen_values [mat_size]\n");
+        return 0;
+    }
 
     if (mat_size < 2)
     {
@@ -84,26 +322,16 @@ int main(int argc, char **argv)
     int i, rows = mat_size, columns = mat_size;
 
     double **A = (double **)malloc(sizeof(double *) * mat_size);
-    double **R = (double **)malloc(sizeof(double *) * mat_size);
-    double **Q = (double **)malloc(sizeof(double *) * mat_size);
-
     /* number of eigen values = matrix size */
     double *eigen_vals = (double *)malloc(sizeof(double) * mat_size);
-    if (!Q || !R || !eigen_vals)
+    if (!eigen_vals)
     {
-        perror("Unable to allocate memory for Q & R!");
+        perror("Unable to allocate memory for eigen values!");
         return -1;
     }
     for (i = 0; i < mat_size; i++)
     {
         A[i] = (double *)malloc(sizeof(double) * mat_size);
-        R[i] = (double *)malloc(sizeof(double) * mat_size);
-        Q[i] = (double *)malloc(sizeof(double) * mat_size);
-        if (!Q[i] || !R[i])
-        {
-            perror("Unable to allocate memory for Q & R.");
-            return -1;
-        }
     }
 
     /* create a random matrix */
@@ -111,64 +339,15 @@ int main(int argc, char **argv)
 
     print_matrix(A, mat_size, mat_size);
 
-    int counter = 0, num_eigs = rows - 1;
-    double last_eig = 0;
-
-    clock_t t1 = clock();
-    while (num_eigs > 0) /* continue till all eigen values are found */
-    {
-        /* iterate with QR decomposition */
-        while (fabs(A[num_eigs][num_eigs - 1]) > 1e-10)
-        {
-            last_eig = A[num_eigs][num_eigs];
-            for (int i = 0; i < rows; i++)
-                A[i][i] -= last_eig; /* A - cI */
-            qr_decompose(A, Q, R, rows, columns);
-
-#if defined(DEBUG) || !defined(NDEBUG)
-            print_matrix(A, rows, columns);
-            print_matrix(Q, rows, columns);
-            print_matrix(R, columns, columns);
-            printf("-------------------- %d ---------------------\n",
-                   ++counter);
-#endif
-            mat_mul(R, Q, A, columns, columns, rows, columns);
-            for (int i = 0; i < rows; i++)
-                A[i][i] += last_eig; /* A + cI */
-        }
-
-        /* store the converged eigen value */
-        eigen_vals[num_eigs] = A[num_eigs][num_eigs];
-#if defined(DEBUG) || !defined(NDEBUG)
-        printf("========================\n");
-        printf("Eigen value: % g,\n", last_eig);
-        printf("========================\n");
-#endif
-        num_eigs--;
-        rows--;
-        columns--;
-    }
-    eigen_vals[0] = A[0][0];
-    double dtime = (double)(clock() - t1) / CLOCKS_PER_SEC;
-
-#if defined(DEBUG) || !defined(NDEBUG)
-    print_matrix(R, mat_size, mat_size);
-    print_matrix(Q, mat_size, mat_size);
-#endif
+    double dtime = eigen_values(A, eigen_vals, mat_size, 0);
     printf("Eigen vals: ");
     for (i = 0; i < mat_size; i++)
         printf("% 9.4g\t", eigen_vals[i]);
     printf("\nTime taken to compute: % .4g sec\n", dtime);
 
     for (int i = 0; i < mat_size; i++)
-    {
         free(A[i]);
-        free(R[i]);
-        free(Q[i]);
-    }
     free(A);
-    free(R);
-    free(Q);
     free(eigen_vals);
     return 0;
 }