document QR algorithms

numerical_methods/qr_decompose.h

@@ -1,8 +1,9 @@
 /**
  * @file
  *
- * Library functions to compute QR decomposition of a
- * given matrix.
+ * \brief Library functions to compute [QR
+ * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of a given
+ * matrix.
  */
 
 #ifndef QR_DECOMPOSE_H
@@ -11,6 +12,9 @@
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
 
 /**
  * function to display matrix on stdout
@@ -36,11 +40,16 @@ void print_matrix(double **A, /**< matrix to print */
  * \f$\vec{a}\cdot\vec{b}=\displaystyle\sum_{i=0}^L a_i\times b_i\f$
  *
  * \returns \f$\vec{a}\cdot\vec{b}\f$
- **/
+ */
 double vector_dot(double *a, double *b, int L)
 {
     double mag = 0.f;
-    for (int i = 0; i < L; i++)
+    int i;
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp parallel for reduction(+ : mag)
+#endif
+    for (i = 0; i < L; i++)
         mag += a[i] * b[i];
 
     return mag;
@@ -53,7 +62,7 @@ double vector_dot(double *a, double *b, int L)
  * \f$\left|\vec{a}\right|=\sqrt{\displaystyle\sum_{i=0}^L a_i^2}\f$
  *
  * \returns \f$\left|\vec{a}\right|\f$
- **/
+ */
 double vector_mag(double *vector, int L)
 {
     double dot = vector_dot(vector, vector, L);
@@ -65,7 +74,7 @@ double vector_mag(double *vector, int L)
  * \f[\text{proj}_\vec{b}\vec{a}=\frac{\vec{a}\cdot\vec{b}}{\left|\vec{b}\right|^2}\vec{b}\f]
  *
  * \returns NULL if error, otherwise pointer to output
- **/
+ */
 double *vector_proj(double *a, double *b, double *out, int L)
 {
     const double num = vector_dot(a, b, L);
@@ -74,7 +83,12 @@ double *vector_proj(double *a, double *b, double *out, int L)
         return NULL;
 
     const double scalar = num / deno;
-    for (int i = 0; i < L; i++)
+    int i;
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+    for (i = 0; i < L; i++)
         out[i] = scalar * b[i];
 
     return out;
@@ -86,14 +100,19 @@ double *vector_proj(double *a, double *b, double *out, int L)
  * \f$\vec{c}=\vec{a}-\vec{b}\f$
  *
  * \returns pointer to output vector
- **/
+ */
 double *vector_sub(double *a,   /**< minuend */
                    double *b,   /**< subtrahend */
                    double *out, /**< resultant vector */
                    int L        /**< length of vectors */
 )
 {
-    for (int i = 0; i < L; i++)
+    int i;
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+    for (i = 0; i < L; i++)
         out[i] = a[i] - b[i];
 
     return out;
@@ -123,7 +142,7 @@ double *vector_sub(double *a,   /**< minuend */
  *                  \vdots & \vdots & \vdots & \ddots
  *      \end{bmatrix}\\
  * \f}
- **/
+ */
 void qr_decompose(double **A, /**< input matrix to decompose */
                   double **Q, /**< output decomposed matrix */
                   double **R, /**< output decomposed matrix */
@@ -137,16 +156,25 @@ void qr_decompose(double **A, /**< input matrix to decompose */
     for (int i = 0; i < N;
          i++) /* for each column => R is a square matrix of NxN */
     {
-        for (int j = 0; j < i; j++) /* second dimension of column */
-            R[i][j] = 0.;           /* make R upper triangular */
+        int j;
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (j = 0; j < i; j++) /* second dimension of column */
+            R[i][j] = 0.;       /* make R upper triangular */
 
-        /* get corresponding Q vector */
-        for (int j = 0; j < M; j++)
+            /* get corresponding Q vector */
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (j = 0; j < M; j++)
         {
             tmp_vector[j] = A[j][i]; /* accumulator for uk */
             col_vector[j] = A[j][i];
         }
-        for (int j = 0; j < i; j++)
+        for (j = 0; j < i; j++)
         {
             for (int k = 0; k < M; k++)
                 col_vector2[k] = Q[k][j];
@@ -154,7 +182,12 @@ void qr_decompose(double **A, /**< input matrix to decompose */
             vector_sub(tmp_vector, col_vector2, tmp_vector, M);
         }
         double mag = vector_mag(tmp_vector, M);
-        for (int j = 0; j < M; j++)
+
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (j = 0; j < M; j++)
             Q[j][i] = tmp_vector[j] / mag;
 
         /* compute upper triangular values of R */

numerical_methods/qr_decomposition.c

@@ -1,7 +1,8 @@
 /**
  * @file
- * Program to compute the QR decomposition of a
- * given matrix.
+ * \brief Program to compute the [QR
+ * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of a given
+ * matrix.
  */
 
 #include "qr_decompose.h"

numerical_methods/qr_eigen_values.c

@@ -1,7 +1,8 @@
 /**
  * @file
- * Compute real eigen values and eigen vectors of a symmetric matrix using QR
- * decomposition method.
+ * \brief Compute real eigen values and eigen vectors of a symmetric matrix
+ * using [QR decomposition](https://en.wikipedia.org/wiki/QR_decomposition)
+ * method.
  */
 #include "qr_decompose.h"
 #include <math.h>