qr_decompose.h File Reference

Library functions to compute QR decomposition of a given matrix. More...

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

Include dependency graph for qr_decompose.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
void	print_matrix (double **A, int M, int N)
	function to display matrix on stdout
double	vector_dot (double *a, double *b, int L)
	Compute dot product of two vectors of equal lengths.
double	vector_mag (double *vector, int L)
	Compute magnitude of vector.
double *	vector_proj (double *a, double *b, double *out, int L)
	Compute projection of vector \(\vec{a}\) on \(\vec{b}\) defined as.
double *	vector_sub (double *a, double *b, double *out, int L)
	Compute vector subtraction.
void	qr_decompose (double **A, double **Q, double **R, int M, int N)
	Decompose matrix \(A\) using Gram-Schmidt process.

Detailed Description

Library functions to compute QR decomposition of a given matrix.

Author

Krishna Vedala

Function Documentation

print_matrix()

void print_matrix	(	double **	A,
		int	M,
		int	N
	)

function to display matrix on stdout

Parameters

A	matrix to print
M	number of rows of matrix
N	number of columns of matrix

{
    for (int row = 0; row < M; row++)
    {
        for (int col = 0; col < N; col++) printf("% 9.3g\t", A[row][col]);
        putchar('\n');
    }
    putchar('\n');
}

qr_decompose()

void qr_decompose	(	double **	A,
		double **	Q,
		double **	R,
		int	M,
		int	N
	)

Decompose matrix \(A\) using Gram-Schmidt process.

\begin{eqnarray*} \text{given that}\quad A &=& *\left[\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_{N-1},\right]\\ \text{where}\quad\mathbf{a}_i &=& *\left[a_{0i},a_{1i},a_{2i},\ldots,a_{(M-1)i}\right]^T\quad\ldots\mbox{(column *vectors)}\\ \text{then}\quad\mathbf{u}_i &=& \mathbf{a}_i *-\sum_{j=0}^{i-1}\text{proj}_{\mathbf{u}_j}\mathbf{a}_i\\ \mathbf{e}_i &=&\frac{\mathbf{u}_i}{\left|\mathbf{u}_i\right|}\\ Q &=& \begin{bmatrix}\mathbf{e}_0 & \mathbf{e}_1 & \mathbf{e}_2 & \dots & *\mathbf{e}_{N-1}\end{bmatrix}\\ R &=& \begin{bmatrix}\langle\mathbf{e}_0\,,\mathbf{a}_0\rangle & *\langle\mathbf{e}_1\,,\mathbf{a}_1\rangle & *\langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots \\ 0 & \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle & *\langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots\\ 0 & 0 & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots\\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}\\ \end{eqnarray*}

Parameters

A	input matrix to decompose
Q	output decomposed matrix
R	output decomposed matrix
M	number of rows of matrix A
N	number of columns of matrix A

{
    double *col_vector = (double *)malloc(M * sizeof(double));
    double *col_vector2 = (double *)malloc(M * sizeof(double));
    double *tmp_vector = (double *)malloc(M * sizeof(double));
    for (int i = 0; i < N;
         i++) /* for each column => R is a square matrix of NxN */
    {
        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 */
#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 (j = 0; j < i; j++)
        {
            for (int k = 0; k < M; k++) col_vector2[k] = Q[k][j];
            vector_proj(col_vector, col_vector2, col_vector2, M);
            vector_sub(tmp_vector, col_vector2, tmp_vector, M);
        }
        double mag = vector_mag(tmp_vector, M);
 
#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 */
        for (int kk = 0; kk < M; kk++) col_vector[kk] = Q[kk][i];
        for (int k = i; k < N; k++)
        {
            for (int kk = 0; kk < M; kk++) col_vector2[kk] = A[kk][k];
            R[i][k] = vector_dot(col_vector, col_vector2, M);
        }
    }
 
    free(col_vector);
    free(col_vector2);
    free(tmp_vector);
}

Here is the call graph for this function:

vector_dot()

double vector_dot	(	double *	a,
		double *	b,
		int	L
	)

Compute dot product of two vectors of equal lengths.

If \(\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\) and \(\vec{b}=\left[b_0,b_1,b_1,...,b_L\right]\) then \(\vec{a}\cdot\vec{b}=\displaystyle\sum_{i=0}^L a_i\times b_i\)

Returns

\(\vec{a}\cdot\vec{b}\)

{
    double mag = 0.f;
    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;
}

vector_mag()

double vector_mag	(	double *	vector,
		int	L
	)

Compute magnitude of vector.

If \(\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\) then \(\left|\vec{a}\right|=\sqrt{\displaystyle\sum_{i=0}^L a_i^2}\)

Returns

\(\left|\vec{a}\right|\)

{
    double dot = vector_dot(vector, vector, L);
    return sqrt(dot);
}

Here is the call graph for this function:

vector_proj()

double * vector_proj	(	double *	a,
		double *	b,
		double *	out,
		int	L
	)

Compute projection of vector \(\vec{a}\) on \(\vec{b}\) defined as.

\[\text{proj}_\vec{b}\vec{a}=\frac{\vec{a}\cdot\vec{b}}{\left|\vec{b}\right|^2}\vec{b}\]

Returns

NULL if error, otherwise pointer to output

check for division by zero

{
    const double num = vector_dot(a, b, L);
    const double deno = vector_dot(b, b, L);
    if (deno == 0) /*! check for division by zero */
        return NULL;
 
    const double scalar = num / deno;
    int i;
#ifdef _OPENMP
// parallelize on threads
#pragma omp for
#endif
    for (i = 0; i < L; i++) out[i] = scalar * b[i];
 
    return out;
}

Here is the call graph for this function:

vector_sub()

double * vector_sub	(	double *	a,
		double *	b,
		double *	out,
		int	L
	)

Compute vector subtraction.

\(\vec{c}=\vec{a}-\vec{b}\)

Returns

pointer to output vector

Parameters

a	minuend
b	subtrahend
out	resultant vector
L	length of vectors

{
    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;
}