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 Cholesky decomposition method definitions
25 //
26
27 #include <math.h>
28 #include "liquid.internal.h"
29
30 #define DEBUG_MATRIX_CHOL 0
31
32 // Compute Cholesky decomposition of a symmetric/Hermitian positive-
33 // definite matrix as A = L * L^T
34 // _A : input square matrix [size: _n x _n]
35 // _n : input matrix dimension
36 // _L : output lower-triangular matrix
MATRIX(_chol)37 void MATRIX(_chol)(T * _A,
38 unsigned int _n,
39 T * _L)
40 {
41 // reset L
42 unsigned int i;
43 for (i=0; i<_n*_n; i++)
44 _L[i] = 0.0;
45
46 unsigned int j;
47 unsigned int k;
48 T A_jj;
49 T L_jj;
50 T L_ik;
51 T L_jk;
52 TP t0;
53 T t1;
54 for (j=0; j<_n; j++) {
55 // assert that A_jj is real, positive
56 A_jj = matrix_access(_A,_n,_n,j,j);
57 if ( creal(A_jj) < 0.0 ) {
58 fprintf(stderr,"warning: matrix_chol(), matrix is not positive definite (real{A[%u,%u]} = %12.4e < 0)\n",j,j,creal(A_jj));
59 return;
60 }
61 #if T_COMPLEX
62 if ( fabs(cimag(A_jj)) > 0.0 ) {
63 fprintf(stderr,"warning: matrix_chol(), matrix is not positive definite (|imag{A[%u,%u]}| = %12.4e > 0)\n",j,j,fabs(cimag(A_jj)));
64 return;
65 }
66 #endif
67
68 // compute L_jj and store it in output matrix
69 t0 = 0.0;
70 for (k=0; k<j; k++) {
71 L_jk = matrix_access(_L,_n,_n,j,k);
72 #if T_COMPLEX
73 t0 += creal( L_jk * conj(L_jk) );
74 #else
75 t0 += L_jk * L_jk;
76 #endif
77 }
78 // test to ensure A_jj > t0
79 if ( creal(A_jj) < t0 ) {
80 fprintf(stderr,"warning: matrix_chol(), matrix is not positive definite (real{A[%u,%u]} = %12.4e < %12.4e)\n",j,j,creal(A_jj),t0);
81 return;
82 }
83 L_jj = sqrt( A_jj - t0 );
84 matrix_access(_L,_n,_n,j,j) = L_jj;
85
86 for (i=j+1; i<_n; i++) {
87 t1 = matrix_access(_A,_n,_n,i,j);
88 for (k=0; k<j; k++) {
89 L_ik = matrix_access(_L,_n,_n,i,k);
90 L_jk = matrix_access(_L,_n,_n,j,k);
91 #if T_COMPLEX
92 t1 -= L_ik * conj(L_jk);
93 #else
94 t1 -= L_ik * L_jk;
95 #endif
96 }
97 // TODO : store inverse of L_jj to reduce number of divisions
98 matrix_access(_L,_n,_n,i,j) = t1 / L_jj;
99 }
100 }
101 }
102
103