1 /*
2  * Copyright (c) 2007 - 2015 Joseph Gaeddert
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining a copy
5  * of this software and associated documentation files (the "Software"), to deal
6  * in the Software without restriction, including without limitation the rights
7  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8  * copies of the Software, and to permit persons to whom the Software is
9  * furnished to do so, subject to the following conditions:
10  *
11  * The above copyright notice and this permission notice shall be included in
12  * all copies or substantial portions of the Software.
13  *
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20  * THE SOFTWARE.
21  */
22 
23 //
24 // Matrix Q/R decomposition method definitions
25 //
26 
27 #include <math.h>
28 #include "liquid.internal.h"
29 
30 #define DEBUG_MATRIX_QRDECOMP 1
31 
32 // Q/R decomposition using the Gram-Schmidt algorithm
MATRIX(_qrdecomp_gramschmidt)33 void MATRIX(_qrdecomp_gramschmidt)(T *          _x,
34                                    unsigned int _rx,
35                                    unsigned int _cx,
36                                    T *          _Q,
37                                    T *          _R)
38 {
39     // validate input
40     if (_rx != _cx) {
41         fprintf(stderr,"error: matrix_qrdecomp_gramschmidt(), input matrix not square\n");
42         exit(-1);
43     }
44     unsigned int n = _rx;
45 
46     unsigned int i;
47     unsigned int j;
48     unsigned int k;
49 
50     // generate and initialize matrices
51     T e[n*n];   // normalized...
52     for (i=0; i<n*n; i++)
53         e[i] = 0.0f;
54 
55     for (k=0; k<n; k++) {
56         // e(i,k) <- _x(i,k)
57         for (i=0; i<n; i++)
58             matrix_access(e,n,n,i,k) = matrix_access(_x,n,n,i,k);
59 
60         // subtract...
61         for (i=0; i<k; i++) {
62             // compute dot product _x(:,k) * e(:,i)
63             T g = 0;
64             for (j=0; j<n; j++) {
65                 T prod = matrix_access(_x,n,n,j,k) * conj( matrix_access(e,n,n,j,i) );
66                 g += prod;
67             }
68             //printf("  i=%2u, g = %12.4e\n", i, crealf(g));
69             for (j=0; j<n; j++)
70                 matrix_access(e,n,n,j,k) -= matrix_access(e,n,n,j,i) * g;
71         }
72 
73         // compute e_k = e_k / |e_k|
74         float ek = 0.0f;
75         T ak;
76         TP ak2;
77         for (i=0; i<n; i++) {
78             ak  = matrix_access(e,n,n,i,k);
79             ak2 = T_ABS(ak);
80             ak2 = ak2 * ak2;
81             ek += ak2;
82         }
83         ek = sqrtf(ek);
84 
85         // normalize e
86         for (i=0; i<n; i++)
87             matrix_access(e,n,n,i,k) /= ek;
88     }
89 
90     // move Q
91     memmove(_Q, e, n*n*sizeof(T));
92 
93     // compute R
94     // j : row
95     // k : column
96     for (j=0; j<n; j++) {
97         for (k=0; k<n; k++) {
98             if (k < j) {
99                 matrix_access(_R,n,n,j,k) = 0.0f;
100             } else {
101                 // compute dot product between and Q(:,j) and _x(:,k)
102                 T g = 0;
103                 for (i=0; i<n; i++) {
104                     T prod = conj( matrix_access(_Q,n,n,i,j) ) * matrix_access(_x,n,n,i,k);
105                     g += prod;
106                 }
107                 matrix_access(_R,n,n,j,k) = g;
108             }
109         }
110     }
111 }
112 
113