1 /*
2  * Copyright (C) 2011 The Android Open Source Project
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *      http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 /* $Id: db_utilities_indexing.h,v 1.3 2011/06/17 14:03:31 mbansal Exp $ */
18 
19 #ifndef DB_UTILITIES_INDEXING
20 #define DB_UTILITIES_INDEXING
21 
22 
23 
24 /*****************************************************************
25 *    Lean and mean begins here                                   *
26 *****************************************************************/
27 
28 #include "db_utilities.h"
29 
30 /*!
31  * \defgroup LMIndexing (LM) Indexing Utilities (Order Statistics, Matrix Operations)
32  */
33 /*\{*/
34 
db_SetupMatrixRefs(double ** ar,long rows,long cols,double * a)35 inline void db_SetupMatrixRefs(double **ar,long rows,long cols,double *a)
36 {
37     long i;
38     for(i=0;i<rows;i++) ar[i]=&a[i*cols];
39 }
40 
db_SymmetricExtendUpperToLower(double ** A,int rows,int cols)41 inline void db_SymmetricExtendUpperToLower(double **A,int rows,int cols)
42 {
43     int i,j;
44     for(i=1;i<rows;i++) for(j=0;j<i;j++) A[i][j]=A[j][i];
45 }
46 
db_MultiplyMatrixVectorAtb(double * c,const double * const * At,const double * b,int arows,int acols)47 void inline db_MultiplyMatrixVectorAtb(double *c,const double * const *At,const double *b,int arows,int acols)
48 {
49     int i,j;
50     double acc;
51 
52     for(i=0;i<arows;i++)
53     {
54         acc=0;
55         for(j=0;j<acols;j++) acc+=At[j][i]*b[j];
56         c[i]=acc;
57     }
58 }
59 
db_MultiplyMatricesAB(double ** C,const double * const * A,const double * const * B,int arows,int acols,int bcols)60 inline void db_MultiplyMatricesAB(double **C,const double * const *A,const double * const *B,int arows,int acols,int bcols)
61 {
62     int i,j,k;
63     double acc;
64 
65     for(i=0;i<arows;i++) for(j=0;j<bcols;j++)
66     {
67         acc=0;
68         for(k=0;k<acols;k++) acc+=A[i][k]*B[k][j];
69         C[i][j]=acc;
70     }
71 }
72 
db_UpperMultiplyMatricesAtB(double ** Cu,const double * const * At,const double * const * B,int arows,int acols,int bcols)73 inline void db_UpperMultiplyMatricesAtB(double **Cu,const double * const *At,const double * const *B,int arows,int acols,int bcols)
74 {
75     int i,j,k;
76     double acc;
77 
78     for(i=0;i<arows;i++) for(j=i;j<bcols;j++)
79     {
80         acc=0;
81         for(k=0;k<acols;k++) acc+=At[k][i]*B[k][j];
82         Cu[i][j]=acc;
83     }
84 }
85 
86 DB_API void db_Zero(double *d,long nr);
87 
db_MaxIndex2(double s[2])88 inline int db_MaxIndex2(double s[2])
89 {
90     if(s[0]>=s[1]) return(0);
91     return(1);
92 }
93 
db_MaxIndex3(const double s[3])94 inline int db_MaxIndex3(const double s[3])
95 {
96     double best;
97     int pos;
98 
99     best=s[0];pos=0;
100     if(s[1]>best){best=s[1];pos=1;}
101     if(s[2]>best){best=s[2];pos=2;}
102     return(pos);
103 }
104 
db_MaxIndex4(const double s[4])105 inline int db_MaxIndex4(const double s[4])
106 {
107     double best;
108     int pos;
109 
110     best=s[0];pos=0;
111     if(s[1]>best){best=s[1];pos=1;}
112     if(s[2]>best){best=s[2];pos=2;}
113     if(s[3]>best){best=s[3];pos=3;}
114     return(pos);
115 }
116 
db_MaxIndex5(const double s[5])117 inline int db_MaxIndex5(const double s[5])
118 {
119     double best;
120     int pos;
121 
122     best=s[0];pos=0;
123     if(s[1]>best){best=s[1];pos=1;}
124     if(s[2]>best){best=s[2];pos=2;}
125     if(s[3]>best){best=s[3];pos=3;}
126     if(s[4]>best){best=s[4];pos=4;}
127     return(pos);
128 }
129 
db_MaxIndex6(const double s[6])130 inline int db_MaxIndex6(const double s[6])
131 {
132     double best;
133     int pos;
134 
135     best=s[0];pos=0;
136     if(s[1]>best){best=s[1];pos=1;}
137     if(s[2]>best){best=s[2];pos=2;}
138     if(s[3]>best){best=s[3];pos=3;}
139     if(s[4]>best){best=s[4];pos=4;}
140     if(s[5]>best){best=s[5];pos=5;}
141     return(pos);
142 }
143 
db_MaxIndex7(const double s[7])144 inline int db_MaxIndex7(const double s[7])
145 {
146     double best;
147     int pos;
148 
149     best=s[0];pos=0;
150     if(s[1]>best){best=s[1];pos=1;}
151     if(s[2]>best){best=s[2];pos=2;}
152     if(s[3]>best){best=s[3];pos=3;}
153     if(s[4]>best){best=s[4];pos=4;}
154     if(s[5]>best){best=s[5];pos=5;}
155     if(s[6]>best){best=s[6];pos=6;}
156     return(pos);
157 }
158 
db_MinIndex7(const double s[7])159 inline int db_MinIndex7(const double s[7])
160 {
161     double best;
162     int pos;
163 
164     best=s[0];pos=0;
165     if(s[1]<best){best=s[1];pos=1;}
166     if(s[2]<best){best=s[2];pos=2;}
167     if(s[3]<best){best=s[3];pos=3;}
168     if(s[4]<best){best=s[4];pos=4;}
169     if(s[5]<best){best=s[5];pos=5;}
170     if(s[6]<best){best=s[6];pos=6;}
171     return(pos);
172 }
173 
db_MinIndex9(const double s[9])174 inline int db_MinIndex9(const double s[9])
175 {
176     double best;
177     int pos;
178 
179     best=s[0];pos=0;
180     if(s[1]<best){best=s[1];pos=1;}
181     if(s[2]<best){best=s[2];pos=2;}
182     if(s[3]<best){best=s[3];pos=3;}
183     if(s[4]<best){best=s[4];pos=4;}
184     if(s[5]<best){best=s[5];pos=5;}
185     if(s[6]<best){best=s[6];pos=6;}
186     if(s[7]<best){best=s[7];pos=7;}
187     if(s[8]<best){best=s[8];pos=8;}
188     return(pos);
189 }
190 
db_MaxAbsIndex3(const double * s)191 inline int db_MaxAbsIndex3(const double *s)
192 {
193     double t,best;
194     int pos;
195 
196     best=fabs(s[0]);pos=0;
197     t=fabs(s[1]);if(t>best){best=t;pos=1;}
198     t=fabs(s[2]);if(t>best){pos=2;}
199     return(pos);
200 }
201 
db_MaxAbsIndex9(const double * s)202 inline int db_MaxAbsIndex9(const double *s)
203 {
204     double t,best;
205     int pos;
206 
207     best=fabs(s[0]);pos=0;
208     t=fabs(s[1]);if(t>best){best=t;pos=1;}
209     t=fabs(s[2]);if(t>best){best=t;pos=2;}
210     t=fabs(s[3]);if(t>best){best=t;pos=3;}
211     t=fabs(s[4]);if(t>best){best=t;pos=4;}
212     t=fabs(s[5]);if(t>best){best=t;pos=5;}
213     t=fabs(s[6]);if(t>best){best=t;pos=6;}
214     t=fabs(s[7]);if(t>best){best=t;pos=7;}
215     t=fabs(s[8]);if(t>best){best=t;pos=8;}
216     return(pos);
217 }
218 
219 
220 /*!
221 Select ordinal pos (zero based) out of nr_elements in s.
222 temp should point to alloced memory of at least nr_elements*2
223 Optimized runtimes on 450MHz:
224 \code
225   30 with   3 microsecs
226  100 with  11 microsecs
227  300 with  30 microsecs
228  500 with  40 microsecs
229 1000 with 100 microsecs
230 5000 with 540 microsecs
231 \endcode
232 so the expected runtime is around
233 (nr_elements/10) microseconds
234 The total quickselect cost of splitting 500 hypotheses recursively
235 is thus around 100 microseconds
236 
237 Does the same operation as std::nth_element().
238 */
239 DB_API double db_LeanQuickSelect(const double *s,long nr_elements,long pos,double *temp);
240 
241 /*!
242  Median of 3 doubles
243  */
db_TripleMedian(double a,double b,double c)244 inline double db_TripleMedian(double a,double b,double c)
245 {
246     if(a>b)
247     {
248         if(c>a) return(a);
249         else if(c>b) return(c);
250         else return(b);
251     }
252     else
253     {
254         if(c>b) return(b);
255         else if(c>a) return(c);
256         else return(a);
257     }
258 }
259 
260 /*!
261 Align float pointer to nr_bytes by moving forward
262 */
263 DB_API float* db_AlignPointer_f(float *p,unsigned long nr_bytes);
264 
265 /*!
266 Align short pointer to nr_bytes by moving forward
267 */
268 DB_API short* db_AlignPointer_s(short *p,unsigned long nr_bytes);
269 
270 #endif /* DB_UTILITIES_INDEXING */
271