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 Gram-Schmidt Orthonormalization
25 //
26
27 #include <math.h>
28 #include "liquid.internal.h"
29
30 #define DEBUG_MATRIX_GRAMSCHMIDT 0
31
32 // compute projection of _u onto _v, store in _e
MATRIX(_proj)33 void MATRIX(_proj)(T * _u,
34 T * _v,
35 unsigned int _n,
36 T * _e)
37 {
38 // compute dot-product between _u and _v
39 unsigned int i;
40 T uv = 0.;
41 T uu = 0.;
42 for (i=0; i<_n; i++) {
43 uv += _u[i] * _v[i];
44 uu += _u[i] * _u[i];
45 }
46
47 // TODO : check magnitude of _uu
48 T g = uv / uu;
49 for (i=0; i<_n; i++)
50 _e[i] = _u[i] * g;
51 }
52
53 // Orthnormalization using the Gram-Schmidt algorithm
MATRIX(_gramschmidt)54 void MATRIX(_gramschmidt)(T * _x,
55 unsigned int _rx,
56 unsigned int _cx,
57 T * _v)
58 {
59 // validate input
60 if (_rx == 0 || _cx == 0) {
61 fprintf(stderr,"error: matrix_gramschmidt(), input matrix cannot have zero-length dimensions\n");
62 exit(1);
63 }
64
65 unsigned int i;
66 unsigned int j;
67 unsigned int k;
68
69 // copy _x to _u
70 memmove(_v, _x, _rx * _cx * sizeof(T));
71
72 unsigned int n = _rx; // dimensionality of each vector
73 T proj_ij[n];
74 for (j=0; j<_cx; j++) {
75 for (i=0; i<j; i++) {
76 // v_j <- v_j - proj(v_i, v_j)
77
78 #if DEBUG_MATRIX_GRAMSCHMIDT
79 printf("computing proj(v_%u, v_%u)\n", i, j);
80 #endif
81
82 // compute proj(v_i, v_j)
83 T vij = 0.; // dotprod(v_i, v_j)
84 T vii = 0.; // dotprod(v_i, v_i)
85 T ti;
86 T tj;
87 for (k=0; k<n; k++) {
88 ti = matrix_access(_v, _rx, _cx, k, i);
89 tj = matrix_access(_v, _rx, _cx, k, j);
90
91 T prodij = ti * conj(tj);
92 vij += prodij;
93
94 T prodii = ti * conj(ti);
95 vii += prodii;
96 }
97 // TODO : vii should be 1.0 from normalization step below
98 T g = vij / vii;
99
100 // complete projection
101 for (k=0; k<n; k++)
102 proj_ij[k] = matrix_access(_v, _rx, _cx, k, i) * g;
103
104 // subtract projection from v_j
105 for (k=0; k<n; k++)
106 matrix_access(_v, _rx, _cx, k, j) -= proj_ij[k];
107 }
108
109 // normalize v_j
110 T vjj = 0.; // dotprod(v_j, v_j)
111 T tj = 0.;
112 for (k=0; k<n; k++) {
113 tj = matrix_access(_v, _rx, _cx, k, j);
114 T prodjj = tj * conj(tj);
115 vjj += prodjj;
116 }
117 // TODO : check magnitude of vjj
118 T g = 1. / sqrt( creal(vjj) );
119 for (k=0; k<n; k++)
120 matrix_access(_v, _rx, _cx, k, j) *= g;
121
122 #if DEBUG_MATRIX_GRAMSCHMIDT
123 MATRIX(_print)(_v, _rx, _cx);
124 #endif
125 }
126 }
127
128