Mercurial > octave-nkf
annotate liboctave/Array-d.cc @ 7961:a5d1e27ee1f4 ss-3-1-51
3.1.51 snapshot
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 22 Jul 2008 11:40:48 -0400 |
| parents | 4976f66d469b |
| children | 7cbe01c21986 |
| rev | line source |
|---|---|
| 757 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1994, 1995, 1996, 1997, 1998, 2000, 2001, 2003, 2004, |
| 4 2005, 2006, 2007 John W. Eaton | |
| 757 | 5 |
| 6 This file is part of Octave. | |
| 7 | |
| 8 Octave is free software; you can redistribute it and/or modify it | |
| 9 under the terms of the GNU General Public License as published by the | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 757 | 12 |
| 13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 16 for more details. | |
| 17 | |
| 18 You should have received a copy of the GNU General Public License | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 757 | 21 |
| 22 */ | |
| 23 | |
| 2006 | 24 #ifdef HAVE_CONFIG_H |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 757 | 28 // Instantiate Arrays of double values. |
| 29 | |
| 30 #include "Array.h" | |
| 31 #include "Array.cc" | |
| 7433 | 32 #include "oct-sort.cc" |
| 33 | |
| 34 #if defined (HAVE_IEEE754_DATA_FORMAT) | |
| 35 | |
| 36 static inline uint64_t | |
| 37 FloatFlip (uint64_t f) | |
| 38 { | |
| 39 uint64_t mask | |
| 40 = -static_cast<int64_t>(f >> 63) | 0x8000000000000000ULL; | |
| 41 | |
| 42 return f ^ mask; | |
| 43 } | |
| 44 | |
| 45 static inline uint64_t | |
| 46 IFloatFlip (uint64_t f) | |
| 47 { | |
| 48 uint64_t mask = ((f >> 63) - 1) | 0x8000000000000000ULL; | |
| 49 | |
| 50 return f ^ mask; | |
| 51 } | |
| 52 | |
| 53 template <> | |
| 54 bool | |
| 55 ascending_compare (double a, double b) | |
| 56 { | |
| 57 return (xisnan (b) || (a < b)); | |
| 58 } | |
| 59 | |
| 60 template <> | |
| 61 bool | |
| 62 ascending_compare (vec_index<double> *a, vec_index<double> *b) | |
| 63 { | |
| 64 return (xisnan (b->vec) || (a->vec < b->vec)); | |
| 65 } | |
| 66 | |
| 67 template <> | |
| 68 bool | |
| 69 descending_compare (double a, double b) | |
| 70 { | |
| 71 return (xisnan (a) || (a > b)); | |
| 72 } | |
| 73 | |
| 74 template <> | |
| 75 bool | |
| 76 descending_compare (vec_index<double> *a, vec_index<double> *b) | |
| 77 { | |
| 78 return (xisnan (b->vec) || (a->vec > b->vec)); | |
| 79 } | |
| 80 | |
| 81 INSTANTIATE_ARRAY_SORT (uint64_t); | |
| 82 | |
| 83 template <> | |
| 84 Array<double> | |
| 85 Array<double>::sort (octave_idx_type dim, sortmode mode) const | |
| 86 { | |
| 87 Array<double> m = *this; | |
| 88 | |
| 89 dim_vector dv = m.dims (); | |
| 90 | |
| 91 if (m.length () < 1) | |
| 92 return m; | |
| 93 | |
| 94 octave_idx_type ns = dv(dim); | |
| 95 octave_idx_type iter = dv.numel () / ns; | |
| 96 octave_idx_type stride = 1; | |
| 97 for (int i = 0; i < dim; i++) | |
| 98 stride *= dv(i); | |
| 99 | |
| 100 double *v = m.fortran_vec (); | |
| 101 | |
| 102 uint64_t *p = reinterpret_cast<uint64_t *> (v); | |
| 103 | |
| 104 octave_sort<uint64_t> lsort; | |
| 105 | |
| 106 if (mode == ASCENDING) | |
| 107 lsort.set_compare (ascending_compare); | |
| 108 else if (mode == DESCENDING) | |
| 109 lsort.set_compare (descending_compare); | |
|
7463
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7443
diff
changeset
|
110 else |
|
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7443
diff
changeset
|
111 abort (); |
| 7433 | 112 |
| 113 if (stride == 1) | |
| 114 { | |
| 115 for (octave_idx_type j = 0; j < iter; j++) | |
| 116 { | |
| 117 // Flip the data in the vector so that int compares on | |
| 118 // IEEE754 give the correct ordering. | |
| 119 | |
| 120 for (octave_idx_type i = 0; i < ns; i++) | |
| 121 p[i] = FloatFlip (p[i]); | |
| 122 | |
| 123 lsort.sort (p, ns); | |
| 124 | |
| 125 // Flip the data out of the vector so that int compares | |
| 126 // on IEEE754 give the correct ordering. | |
| 127 | |
| 128 for (octave_idx_type i = 0; i < ns; i++) | |
| 129 p[i] = IFloatFlip (p[i]); | |
| 130 | |
| 131 // There are two representations of NaN. One will be | |
| 132 // sorted to the beginning of the vector and the other | |
| 133 // to the end. If it will be sorted incorrectly, fix | |
| 134 // things up. | |
| 135 | |
| 136 if (lo_ieee_signbit (octave_NaN)) | |
| 137 { | |
|
7463
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7443
diff
changeset
|
138 if (mode == ASCENDING) |
| 7433 | 139 { |
| 140 octave_idx_type i = 0; | |
| 141 double *vtmp = reinterpret_cast<double *> (p); | |
| 7924 | 142 while (xisnan (vtmp[i++]) && i < ns) |
| 143 /* do nothing */; | |
| 7433 | 144 for (octave_idx_type l = 0; l < ns - i + 1; l++) |
| 145 vtmp[l] = vtmp[l+i-1]; | |
| 146 for (octave_idx_type l = ns - i + 1; l < ns; l++) | |
| 147 vtmp[l] = octave_NaN; | |
| 148 } | |
| 149 else | |
| 150 { | |
| 151 octave_idx_type i = ns; | |
| 152 double *vtmp = reinterpret_cast<double *> (p); | |
| 7924 | 153 while (xisnan (vtmp[--i]) && i > 0) |
| 154 /* do nothing */; | |
| 7433 | 155 for (octave_idx_type l = i; l >= 0; l--) |
| 156 vtmp[l-i+ns-1] = vtmp[l]; | |
| 157 for (octave_idx_type l = 0; l < ns - i - 1; l++) | |
| 158 vtmp[l] = octave_NaN; | |
| 159 } | |
| 160 } | |
| 161 | |
| 162 p += ns; | |
| 163 } | |
| 164 } | |
| 165 else | |
| 166 { | |
| 167 OCTAVE_LOCAL_BUFFER (uint64_t, vi, ns); | |
| 168 | |
| 169 for (octave_idx_type j = 0; j < iter; j++) | |
| 170 { | |
| 171 octave_idx_type offset = j; | |
| 172 octave_idx_type offset2 = 0; | |
| 173 while (offset >= stride) | |
| 174 { | |
| 175 offset -= stride; | |
| 176 offset2++; | |
| 177 } | |
| 178 offset += offset2 * stride * ns; | |
| 179 | |
| 180 // Flip the data in the vector so that int compares on | |
| 181 // IEEE754 give the correct ordering. | |
| 182 | |
| 183 for (octave_idx_type i = 0; i < ns; i++) | |
| 184 vi[i] = FloatFlip (p[i*stride + offset]); | |
| 185 | |
| 186 lsort.sort (vi, ns); | |
| 187 | |
| 188 // Flip the data out of the vector so that int compares | |
| 189 // on IEEE754 give the correct ordering. | |
| 190 | |
| 191 for (octave_idx_type i = 0; i < ns; i++) | |
| 192 p[i*stride + offset] = IFloatFlip (vi[i]); | |
| 193 | |
| 194 // There are two representations of NaN. One will be | |
| 195 // sorted to the beginning of the vector and the other | |
| 196 // to the end. If it will be sorted to the beginning, | |
| 197 // fix things up. | |
| 198 | |
| 199 if (lo_ieee_signbit (octave_NaN)) | |
| 200 { | |
|
7463
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7443
diff
changeset
|
201 if (mode == ASCENDING) |
| 7433 | 202 { |
| 203 octave_idx_type i = 0; | |
| 7924 | 204 while (xisnan (v[i++*stride + offset]) && i < ns) |
| 205 /* do nothing */; | |
| 7433 | 206 for (octave_idx_type l = 0; l < ns - i + 1; l++) |
| 207 v[l*stride + offset] = v[(l+i-1)*stride + offset]; | |
| 208 for (octave_idx_type l = ns - i + 1; l < ns; l++) | |
| 209 v[l*stride + offset] = octave_NaN; | |
| 210 } | |
| 211 else | |
| 212 { | |
| 213 octave_idx_type i = ns; | |
| 7924 | 214 while (xisnan (v[--i*stride + offset]) && i > 0) |
| 215 /* do nothing */; | |
| 7433 | 216 for (octave_idx_type l = i; l >= 0; l--) |
| 217 v[(l-i+ns-1)*stride + offset] = v[l*stride + offset]; | |
| 218 for (octave_idx_type l = 0; l < ns - i - 1; l++) | |
| 219 v[l*stride + offset] = octave_NaN; | |
| 220 } | |
| 221 } | |
| 222 } | |
| 223 } | |
| 224 | |
| 225 return m; | |
| 226 } | |
| 227 | |
| 228 template <> | |
| 229 Array<double> | |
| 230 Array<double>::sort (Array<octave_idx_type> &sidx, octave_idx_type dim, | |
| 231 sortmode mode) const | |
| 232 { | |
| 233 Array<double> m = *this; | |
| 234 | |
| 235 dim_vector dv = m.dims (); | |
| 236 | |
| 237 if (m.length () < 1) | |
| 238 { | |
| 239 sidx = Array<octave_idx_type> (dv); | |
| 240 return m; | |
| 241 } | |
| 242 | |
| 243 octave_idx_type ns = dv(dim); | |
| 244 octave_idx_type iter = dv.numel () / ns; | |
| 245 octave_idx_type stride = 1; | |
| 246 for (int i = 0; i < dim; i++) | |
| 247 stride *= dv(i); | |
| 248 | |
| 249 double *v = m.fortran_vec (); | |
| 250 | |
| 251 uint64_t *p = reinterpret_cast<uint64_t *> (v); | |
| 252 | |
| 253 octave_sort<vec_index<uint64_t> *> indexed_sort; | |
| 254 | |
| 255 if (mode == ASCENDING) | |
| 256 indexed_sort.set_compare (ascending_compare); | |
| 257 else if (mode == DESCENDING) | |
| 258 indexed_sort.set_compare (descending_compare); | |
|
7463
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7443
diff
changeset
|
259 else |
|
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7443
diff
changeset
|
260 abort (); |
| 7433 | 261 |
| 262 OCTAVE_LOCAL_BUFFER (vec_index<uint64_t> *, vi, ns); | |
| 263 OCTAVE_LOCAL_BUFFER (vec_index<uint64_t>, vix, ns); | |
| 264 | |
| 265 for (octave_idx_type i = 0; i < ns; i++) | |
| 266 vi[i] = &vix[i]; | |
| 267 | |
| 268 sidx = Array<octave_idx_type> (dv); | |
| 269 | |
| 270 for (octave_idx_type j = 0; j < iter; j++) | |
| 271 { | |
| 272 octave_idx_type offset = j; | |
| 273 octave_idx_type offset2 = 0; | |
| 274 while (offset >= stride) | |
| 275 { | |
| 276 offset -= stride; | |
| 277 offset2++; | |
| 278 } | |
| 279 offset += offset2 * stride * ns; | |
| 280 | |
| 281 // Flip the data in the vector so that int compares on | |
| 282 // IEEE754 give the correct ordering. | |
| 283 | |
| 284 for (octave_idx_type i = 0; i < ns; i++) | |
| 285 { | |
| 286 vi[i]->vec = FloatFlip (p[i*stride + offset]); | |
| 287 vi[i]->indx = i; | |
| 288 } | |
| 289 | |
| 290 indexed_sort.sort (vi, ns); | |
| 291 | |
| 292 // Flip the data out of the vector so that int compares on | |
| 293 // IEEE754 give the correct ordering | |
| 294 | |
| 295 for (octave_idx_type i = 0; i < ns; i++) | |
| 296 { | |
| 297 p[i*stride + offset] = IFloatFlip (vi[i]->vec); | |
| 298 sidx(i*stride + offset) = vi[i]->indx; | |
| 299 } | |
| 300 | |
| 301 // There are two representations of NaN. One will be sorted | |
| 302 // to the beginning of the vector and the other to the end. | |
| 303 // If it will be sorted to the beginning, fix things up. | |
| 304 | |
| 305 if (lo_ieee_signbit (octave_NaN)) | |
| 306 { | |
|
7463
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7443
diff
changeset
|
307 if (mode == ASCENDING) |
| 7433 | 308 { |
| 309 octave_idx_type i = 0; | |
| 7924 | 310 while (xisnan (v[i++*stride+offset]) && i < ns) |
| 311 /* do nothing */; | |
| 7433 | 312 OCTAVE_LOCAL_BUFFER (double, itmp, i - 1); |
| 313 for (octave_idx_type l = 0; l < i -1; l++) | |
| 314 itmp[l] = sidx(l*stride + offset); | |
| 315 for (octave_idx_type l = 0; l < ns - i + 1; l++) | |
| 316 { | |
| 317 v[l*stride + offset] = v[(l+i-1)*stride + offset]; | |
| 318 sidx(l*stride + offset) = sidx((l+i-1)*stride + offset); | |
| 319 } | |
| 320 for (octave_idx_type k = 0, l = ns - i + 1; l < ns; l++, k++) | |
| 321 { | |
| 322 v[l*stride + offset] = octave_NaN; | |
| 323 sidx(l*stride + offset) = | |
| 324 static_cast<octave_idx_type>(itmp[k]); | |
| 325 } | |
| 326 } | |
| 327 else | |
| 328 { | |
| 329 octave_idx_type i = ns; | |
| 7924 | 330 while (xisnan (v[--i*stride+offset]) && i > 0) |
| 331 /* do nothing */; | |
| 7433 | 332 OCTAVE_LOCAL_BUFFER (double, itmp, ns - i - 1); |
| 333 for (octave_idx_type l = 0; l < ns - i -1; l++) | |
| 334 itmp[l] = sidx((l+i+1)*stride + offset); | |
| 335 for (octave_idx_type l = i; l >= 0; l--) | |
| 336 { | |
| 337 v[(l-i+ns-1)*stride + offset] = v[l*stride + offset]; | |
| 338 sidx((l-i+ns-1)*stride + offset) = sidx(l*stride + offset); | |
| 339 } | |
| 340 for (octave_idx_type k = 0, l = 0; l < ns - i - 1; l++, k++) | |
| 341 { | |
| 342 v[l*stride + offset] = octave_NaN; | |
| 343 sidx(l*stride + offset) = | |
| 344 static_cast<octave_idx_type>(itmp[k]); | |
| 345 } | |
| 346 } | |
| 347 } | |
| 348 } | |
| 349 | |
| 350 return m; | |
| 351 } | |
| 352 | |
| 353 #else | |
| 354 | |
| 355 template <> | |
| 356 bool | |
| 7443 | 357 ascending_compare (double a, double b) |
| 7433 | 358 { |
| 359 return (xisnan (b) || (a < b)); | |
| 360 } | |
| 361 | |
| 362 template <> | |
| 363 bool | |
| 7443 | 364 ascending_compare (vec_index<double> *a, |
| 365 vec_index<double> *b) | |
| 7433 | 366 { |
| 367 return (xisnan (b->vec) || (a->vec < b->vec)); | |
| 368 } | |
| 369 | |
| 370 template <> | |
| 371 bool | |
| 7443 | 372 descending_compare (double a, double b) |
| 7433 | 373 { |
| 374 return (xisnan (a) || (a > b)); | |
| 375 } | |
| 376 | |
| 377 template <> | |
| 378 bool | |
| 7443 | 379 descending_compare (vec_index<double> *a, |
| 380 vec_index<double> *b) | |
| 7433 | 381 { |
| 382 return (xisnan (b->vec) || (a->vec > b->vec)); | |
| 383 } | |
| 384 | |
| 385 INSTANTIATE_ARRAY_SORT (double); | |
| 386 | |
| 387 #endif | |
| 757 | 388 |
| 6708 | 389 INSTANTIATE_ARRAY_AND_ASSIGN (double, OCTAVE_API); |
| 757 | 390 |
| 6708 | 391 INSTANTIATE_ARRAY_ASSIGN (double, int, OCTAVE_API); |
| 392 INSTANTIATE_ARRAY_ASSIGN (double, short, OCTAVE_API); | |
| 393 INSTANTIATE_ARRAY_ASSIGN (double, char, OCTAVE_API); | |
| 3836 | 394 |
| 1989 | 395 #include "Array2.h" |
| 396 | |
| 6153 | 397 template class OCTAVE_API Array2<double>; |
| 1989 | 398 |
| 3665 | 399 #include "ArrayN.h" |
| 400 #include "ArrayN.cc" | |
| 401 | |
| 6153 | 402 template class OCTAVE_API ArrayN<double>; |
| 4505 | 403 |
| 6153 | 404 template OCTAVE_API std::ostream& operator << (std::ostream&, const ArrayN<double>&); |
| 3665 | 405 |
| 1989 | 406 #include "DiagArray2.h" |
| 407 #include "DiagArray2.cc" | |
| 408 | |
| 6153 | 409 template class OCTAVE_API DiagArray2<double>; |
| 1989 | 410 |
| 757 | 411 /* |
| 412 ;;; Local Variables: *** | |
| 413 ;;; mode: C++ *** | |
| 414 ;;; End: *** | |
| 415 */ |