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
35inline 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
41inline 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
47void 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
60inline 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
73inline 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
86DB_API void db_Zero(double *d,long nr);
87
88inline int db_MaxIndex2(double s[2])
89{
90    if(s[0]>=s[1]) return(0);
91    return(1);
92}
93
94inline 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
105inline 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
117inline 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
130inline 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
144inline 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
159inline 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
174inline 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
191inline 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
202inline 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/*!
221Select ordinal pos (zero based) out of nr_elements in s.
222temp should point to alloced memory of at least nr_elements*2
223Optimized 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
2291000 with 100 microsecs
2305000 with 540 microsecs
231\endcode
232so the expected runtime is around
233(nr_elements/10) microseconds
234The total quickselect cost of splitting 500 hypotheses recursively
235is thus around 100 microseconds
236
237Does the same operation as std::nth_element().
238*/
239DB_API double db_LeanQuickSelect(const double *s,long nr_elements,long pos,double *temp);
240
241/*!
242 Median of 3 doubles
243 */
244inline 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/*!
261Align float pointer to nr_bytes by moving forward
262*/
263DB_API float* db_AlignPointer_f(float *p,unsigned long nr_bytes);
264
265/*!
266Align short pointer to nr_bytes by moving forward
267*/
268DB_API short* db_AlignPointer_s(short *p,unsigned long nr_bytes);
269
270#endif /* DB_UTILITIES_INDEXING */
271