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