Bolander-Linear-Algebra-Library

Repository containing homework assignments and software for Math 5610: Computational Linear Algebra and Solution of Systems of Equations

Find the Orthogonal QR Factorization of a Matrix (Modified Gram-Schmidt)

This is a part of the student software manual project for Math 5610: Computational Linear Algebra and Solution of Systems of Equations.

Routine Name: qr_factor_modgs

Author: Christian Bolander

Language: Fortran. This code can be compiled using the GNU Fortran compiler by

$ gfortran -c qr_factor_modgs.f90

and can be added to a program using

$ gfortran program.f90 qr_factor_modgs.o

Description/Purpose: This subroutine finds the orthogonal QR factorization of a matrix using the modified Gram-Schmidt orthogonality algorithm to be used in solving least squares problems. That is, it decomposes A such that

$\begin{pmatrix} a_{11} &a_{12} \\ a_{21} &a_{22} \\ a_{31} &a_{32} \end{pmatrix} = \begin{pmatrix} q_{11} &q_{12} \\ q_{21} &q_{22} \\ q_{31} &q_{32} \end{pmatrix}\begin{pmatrix} r_{11} &r_{12} \\ 0 &r_{22} \end{pmatrix}$

and returns the Q and R matrices. Since Q is orthogonal, $Q^TQ$ should produce the identity matrix and Q x R should reproduce A.

Input:

A : REAL - a coefficient matrix of size m x n

m : INTEGER - the number of rows in A and Q

n : INTEGER - the number of columns in A and the size of the square matrix R

Output:

Q : REAL - the orthogonal matrix ‘Q’ of size m x n of the QR factorization

R : REAL - the upper triangular ‘R’ factor of size n x n of the QR factorization

Usage/Example:

This routine can be implemented in a program as follows (which follows example 6.1 in Ascher)

INTEGER :: m, n, i
REAL*8, ALLOCATABLE :: A(:, :), Q(:, :), R(:, :), QTQ(:, :)

m = 5
n = 3
ALLOCATE(A(1:m, 1:n), Q(1:m, 1:n), R(1:n, 1:n), QTQ(1:n, 1:n))
A = RESHAPE((/1.D0, 0.0D0, 1.0D0, &
			& 2.D0, 3.0D0, 5.0D0, &
			& 5.D0, 3.0D0, -2.0D0, &
			& 3.D0, 5.0D0, 4.0D0, &
			& -1.D0, 6.0D0, 3.0D0/), (/m, n/), ORDER=(/2, 1/))
CALL qr_factor_modgs(A, m, n, Q, R)
WRITE(*,*) "A"
DO i = 1, m
	WRITE(*,*) A(i, :)
END DO
WRITE(*,*) "Q"
DO i = 1, m
	WRITE(*,*) Q(i, :)
END DO
WRITE(*,*) "R"
DO i = 1, n
	WRITE(*,*) R(i, :)
END DO
CALL mat_prod(TRANSPOSE(Q), Q, n, m, n, QTQ)
WRITE(*,*) "QTQ"
DO i = 1, n
	WRITE(*,*) QTQ(i, :)
END DO
CALL mat_prod(Q, R, m, n, n, A)
WRITE(*,*) "A = QR"
DO i = 1, m
	WRITE(*,*) A(i, :)
END DO

The outputs from the above code:

 A
   1.0000000000000000        0.0000000000000000        1.0000000000000000     
   2.0000000000000000        3.0000000000000000        5.0000000000000000     
   5.0000000000000000        3.0000000000000000       -2.0000000000000000     
   3.0000000000000000        5.0000000000000000        4.0000000000000000     
  -1.0000000000000000        6.0000000000000000        3.0000000000000000     
 Q
  0.15811388300841897       -9.9778515785660896E-002  0.25545570859468658     
  0.31622776601683794       0.19955703157132179       0.69185921077727641     
  0.79056941504209477       -9.9778515785660840E-002 -0.54639137671641314     
  0.47434164902525688       0.36585455788075660       0.26609969645279863     
 -0.15811388300841897       0.89800664207094805      -0.29448366407443027     
 R
   6.3245553203367590        4.7434164902525691        1.5811388300841895     
   0.0000000000000000        7.5166481891864541        5.2550018313781397     
   0.0000000000000000        0.0000000000000000        4.9884823095017978     
 QTQ
  0.99999999999999989        2.7755575615628914E-017  -6.9388939039072284E-018
   2.7755575615628914E-017  0.99999999999999989        1.1102230246251565E-016
  -6.9388939039072284E-018   1.1102230246251565E-016  0.99999999999999989     
 A = QR
   1.0000000000000000        0.0000000000000000       0.99999999999999989     
   2.0000000000000000        3.0000000000000000        5.0000000000000000     
   5.0000000000000000        3.0000000000000000       -2.0000000000000000     
   3.0000000000000000        5.0000000000000000        4.0000000000000000     
  -1.0000000000000000        6.0000000000000000        3.0000000000000000 

Implementation/Code: The code for qr_factor_modgs can be seen below.

SUBROUTINE qr_factor_modgs(A, m, n, Q, R)
	IMPLICIT NONE
	
	! Takes as a input a coefficient matrix `A` of size `m` x `n` and
	! finds the QR factorization of that matrix using the modified
	! Gram-Schmidt orthogonalization algorithm. Outputs `Q` and `R`,
	! which are of size `m` x `n` and `n` x `n` respectively. Uses the
	! `vec_dot_prod` and `l2_vec_norm` subroutines.
	INTEGER, INTENT(IN) :: m, n
	REAL*8, INTENT(IN) :: A(1:m, 1:n)
	REAL*8, INTENT(OUT) :: Q(1:m, 1:n), R(1:n, 1:n)
	
	INTEGER i, j
	
	! Iterates over each column to perform the QR factorization
	DO j = 1, n
		Q(:, j) = A(:, j)
		
		! Calculates the dot product of the columns to the left of the
		! current column to find R(i, j).
		DO i = 1, j - 1
			! Using Q in this line instead of A is what makes this the
			! modified Gram-Schmidt.
			CALL vec_dot_prod(Q(:, j), Q(:, i), m, R(i, j))
			
			! Updates Q(j) using the previous columns and R(i, j)
			Q(:, j) = Q(:, j) - R(i, j)*Q(:, i)
		END DO
		
		! Finds the new column of Q using R.
		CALL l2_vec_norm(Q(:, j), m, R(j, j))
		Q(:, j) = Q(:, j)/R(j, j)
	END DO
	
END SUBROUTINE