The Tor Browser / annotate

media/libopus/silk/float/solve_LS_FLP.c@b8a032363ba2 (annotated)

media/libopus/silk/float/solve_LS_FLP.c

Thu, 22 Jan 2015 13:21:57 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Thu, 22 Jan 2015 13:21:57 +0100
branch
TOR_BUG_9701
changeset 15
b8a032363ba2
permissions
-rw-r--r--

Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6

michael@0 1 /***********************************************************************
michael@0 2 Copyright (c) 2006-2011, Skype Limited. All rights reserved.
michael@0 3 Redistribution and use in source and binary forms, with or without
michael@0 4 modification, are permitted provided that the following conditions
michael@0 5 are met:
michael@0 6 - Redistributions of source code must retain the above copyright notice,
michael@0 7 this list of conditions and the following disclaimer.
michael@0 8 - Redistributions in binary form must reproduce the above copyright
michael@0 9 notice, this list of conditions and the following disclaimer in the
michael@0 10 documentation and/or other materials provided with the distribution.
michael@0 11 - Neither the name of Internet Society, IETF or IETF Trust, nor the
michael@0 12 names of specific contributors, may be used to endorse or promote
michael@0 13 products derived from this software without specific prior written
michael@0 14 permission.
michael@0 15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
michael@0 16 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
michael@0 17 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
michael@0 18 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
michael@0 19 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
michael@0 20 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
michael@0 21 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
michael@0 22 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
michael@0 23 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
michael@0 24 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
michael@0 25 POSSIBILITY OF SUCH DAMAGE.
michael@0 26 ***********************************************************************/
michael@0 27
michael@0 28 #ifdef HAVE_CONFIG_H
michael@0 29 #include "config.h"
michael@0 30 #endif
michael@0 31
michael@0 32 #include "main_FLP.h"
michael@0 33 #include "tuning_parameters.h"
michael@0 34
michael@0 35 /**********************************************************************
michael@0 36 * LDL Factorisation. Finds the upper triangular matrix L and the diagonal
michael@0 37 * Matrix D (only the diagonal elements returned in a vector)such that
michael@0 38 * the symmetric matric A is given by A = L*D*L'.
michael@0 39 **********************************************************************/
michael@0 40 static OPUS_INLINE void silk_LDL_FLP(
michael@0 41 silk_float *A, /* I/O Pointer to Symetric Square Matrix */
michael@0 42 opus_int M, /* I Size of Matrix */
michael@0 43 silk_float *L, /* I/O Pointer to Square Upper triangular Matrix */
michael@0 44 silk_float *Dinv /* I/O Pointer to vector holding the inverse diagonal elements of D */
michael@0 45 );
michael@0 46
michael@0 47 /**********************************************************************
michael@0 48 * Function to solve linear equation Ax = b, when A is a MxM lower
michael@0 49 * triangular matrix, with ones on the diagonal.
michael@0 50 **********************************************************************/
michael@0 51 static OPUS_INLINE void silk_SolveWithLowerTriangularWdiagOnes_FLP(
michael@0 52 const silk_float *L, /* I Pointer to Lower Triangular Matrix */
michael@0 53 opus_int M, /* I Dim of Matrix equation */
michael@0 54 const silk_float *b, /* I b Vector */
michael@0 55 silk_float *x /* O x Vector */
michael@0 56 );
michael@0 57
michael@0 58 /**********************************************************************
michael@0 59 * Function to solve linear equation (A^T)x = b, when A is a MxM lower
michael@0 60 * triangular, with ones on the diagonal. (ie then A^T is upper triangular)
michael@0 61 **********************************************************************/
michael@0 62 static OPUS_INLINE void silk_SolveWithUpperTriangularFromLowerWdiagOnes_FLP(
michael@0 63 const silk_float *L, /* I Pointer to Lower Triangular Matrix */
michael@0 64 opus_int M, /* I Dim of Matrix equation */
michael@0 65 const silk_float *b, /* I b Vector */
michael@0 66 silk_float *x /* O x Vector */
michael@0 67 );
michael@0 68
michael@0 69 /**********************************************************************
michael@0 70 * Function to solve linear equation Ax = b, when A is a MxM
michael@0 71 * symmetric square matrix - using LDL factorisation
michael@0 72 **********************************************************************/
michael@0 73 void silk_solve_LDL_FLP(
michael@0 74 silk_float *A, /* I/O Symmetric square matrix, out: reg. */
michael@0 75 const opus_int M, /* I Size of matrix */
michael@0 76 const silk_float *b, /* I Pointer to b vector */
michael@0 77 silk_float *x /* O Pointer to x solution vector */
michael@0 78 )
michael@0 79 {
michael@0 80 opus_int i;
michael@0 81 silk_float L[ MAX_MATRIX_SIZE ][ MAX_MATRIX_SIZE ];
michael@0 82 silk_float T[ MAX_MATRIX_SIZE ];
michael@0 83 silk_float Dinv[ MAX_MATRIX_SIZE ]; /* inverse diagonal elements of D*/
michael@0 84
michael@0 85 silk_assert( M <= MAX_MATRIX_SIZE );
michael@0 86
michael@0 87 /***************************************************
michael@0 88 Factorize A by LDL such that A = L*D*(L^T),
michael@0 89 where L is lower triangular with ones on diagonal
michael@0 90 ****************************************************/
michael@0 91 silk_LDL_FLP( A, M, &L[ 0 ][ 0 ], Dinv );
michael@0 92
michael@0 93 /****************************************************
michael@0 94 * substitute D*(L^T) = T. ie:
michael@0 95 L*D*(L^T)*x = b => L*T = b <=> T = inv(L)*b
michael@0 96 ******************************************************/
michael@0 97 silk_SolveWithLowerTriangularWdiagOnes_FLP( &L[ 0 ][ 0 ], M, b, T );
michael@0 98
michael@0 99 /****************************************************
michael@0 100 D*(L^T)*x = T <=> (L^T)*x = inv(D)*T, because D is
michael@0 101 diagonal just multiply with 1/d_i
michael@0 102 ****************************************************/
michael@0 103 for( i = 0; i < M; i++ ) {
michael@0 104 T[ i ] = T[ i ] * Dinv[ i ];
michael@0 105 }
michael@0 106 /****************************************************
michael@0 107 x = inv(L') * inv(D) * T
michael@0 108 *****************************************************/
michael@0 109 silk_SolveWithUpperTriangularFromLowerWdiagOnes_FLP( &L[ 0 ][ 0 ], M, T, x );
michael@0 110 }
michael@0 111
michael@0 112 static OPUS_INLINE void silk_SolveWithUpperTriangularFromLowerWdiagOnes_FLP(
michael@0 113 const silk_float *L, /* I Pointer to Lower Triangular Matrix */
michael@0 114 opus_int M, /* I Dim of Matrix equation */
michael@0 115 const silk_float *b, /* I b Vector */
michael@0 116 silk_float *x /* O x Vector */
michael@0 117 )
michael@0 118 {
michael@0 119 opus_int i, j;
michael@0 120 silk_float temp;
michael@0 121 const silk_float *ptr1;
michael@0 122
michael@0 123 for( i = M - 1; i >= 0; i-- ) {
michael@0 124 ptr1 = matrix_adr( L, 0, i, M );
michael@0 125 temp = 0;
michael@0 126 for( j = M - 1; j > i ; j-- ) {
michael@0 127 temp += ptr1[ j * M ] * x[ j ];
michael@0 128 }
michael@0 129 temp = b[ i ] - temp;
michael@0 130 x[ i ] = temp;
michael@0 131 }
michael@0 132 }
michael@0 133
michael@0 134 static OPUS_INLINE void silk_SolveWithLowerTriangularWdiagOnes_FLP(
michael@0 135 const silk_float *L, /* I Pointer to Lower Triangular Matrix */
michael@0 136 opus_int M, /* I Dim of Matrix equation */
michael@0 137 const silk_float *b, /* I b Vector */
michael@0 138 silk_float *x /* O x Vector */
michael@0 139 )
michael@0 140 {
michael@0 141 opus_int i, j;
michael@0 142 silk_float temp;
michael@0 143 const silk_float *ptr1;
michael@0 144
michael@0 145 for( i = 0; i < M; i++ ) {
michael@0 146 ptr1 = matrix_adr( L, i, 0, M );
michael@0 147 temp = 0;
michael@0 148 for( j = 0; j < i; j++ ) {
michael@0 149 temp += ptr1[ j ] * x[ j ];
michael@0 150 }
michael@0 151 temp = b[ i ] - temp;
michael@0 152 x[ i ] = temp;
michael@0 153 }
michael@0 154 }
michael@0 155
michael@0 156 static OPUS_INLINE void silk_LDL_FLP(
michael@0 157 silk_float *A, /* I/O Pointer to Symetric Square Matrix */
michael@0 158 opus_int M, /* I Size of Matrix */
michael@0 159 silk_float *L, /* I/O Pointer to Square Upper triangular Matrix */
michael@0 160 silk_float *Dinv /* I/O Pointer to vector holding the inverse diagonal elements of D */
michael@0 161 )
michael@0 162 {
michael@0 163 opus_int i, j, k, loop_count, err = 1;
michael@0 164 silk_float *ptr1, *ptr2;
michael@0 165 double temp, diag_min_value;
michael@0 166 silk_float v[ MAX_MATRIX_SIZE ], D[ MAX_MATRIX_SIZE ]; /* temp arrays*/
michael@0 167
michael@0 168 silk_assert( M <= MAX_MATRIX_SIZE );
michael@0 169
michael@0 170 diag_min_value = FIND_LTP_COND_FAC * 0.5f * ( A[ 0 ] + A[ M * M - 1 ] );
michael@0 171 for( loop_count = 0; loop_count < M && err == 1; loop_count++ ) {
michael@0 172 err = 0;
michael@0 173 for( j = 0; j < M; j++ ) {
michael@0 174 ptr1 = matrix_adr( L, j, 0, M );
michael@0 175 temp = matrix_ptr( A, j, j, M ); /* element in row j column j*/
michael@0 176 for( i = 0; i < j; i++ ) {
michael@0 177 v[ i ] = ptr1[ i ] * D[ i ];
michael@0 178 temp -= ptr1[ i ] * v[ i ];
michael@0 179 }
michael@0 180 if( temp < diag_min_value ) {
michael@0 181 /* Badly conditioned matrix: add white noise and run again */
michael@0 182 temp = ( loop_count + 1 ) * diag_min_value - temp;
michael@0 183 for( i = 0; i < M; i++ ) {
michael@0 184 matrix_ptr( A, i, i, M ) += ( silk_float )temp;
michael@0 185 }
michael@0 186 err = 1;
michael@0 187 break;
michael@0 188 }
michael@0 189 D[ j ] = ( silk_float )temp;
michael@0 190 Dinv[ j ] = ( silk_float )( 1.0f / temp );
michael@0 191 matrix_ptr( L, j, j, M ) = 1.0f;
michael@0 192
michael@0 193 ptr1 = matrix_adr( A, j, 0, M );
michael@0 194 ptr2 = matrix_adr( L, j + 1, 0, M);
michael@0 195 for( i = j + 1; i < M; i++ ) {
michael@0 196 temp = 0.0;
michael@0 197 for( k = 0; k < j; k++ ) {
michael@0 198 temp += ptr2[ k ] * v[ k ];
michael@0 199 }
michael@0 200 matrix_ptr( L, i, j, M ) = ( silk_float )( ( ptr1[ i ] - temp ) * Dinv[ j ] );
michael@0 201 ptr2 += M; /* go to next column*/
michael@0 202 }
michael@0 203 }
michael@0 204 }
michael@0 205 silk_assert( err == 0 );
michael@0 206 }
michael@0 207

Mercurial Repository: The Tor Browser

Europalab SCM Repository Server

mercurial