SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2013 Thoralf Klein 00008 * Written (W) 2011-2013 Heiko Strathmann 00009 * Written (W) 2013 Soumyajit De 00010 * Written (W) 2012 Fernando José Iglesias García 00011 * Written (W) 2010,2012 Soeren Sonnenburg 00012 * Copyright (C) 2010 Berlin Institute of Technology 00013 * Copyright (C) 2012 Soeren Sonnenburg 00014 */ 00015 00016 #include <shogun/lib/config.h> 00017 #include <shogun/lib/SGVector.h> 00018 #include <shogun/lib/SGMatrix.h> 00019 #include <shogun/lib/SGSparseVector.h> 00020 #include <shogun/lib/SGReferencedData.h> 00021 #include <shogun/io/File.h> 00022 00023 #include <shogun/mathematics/Math.h> 00024 #include <shogun/mathematics/lapack.h> 00025 #include <algorithm> 00026 00027 #include <shogun/mathematics/eigen3.h> 00028 00029 #define COMPLEX128_ERROR_NOARG(function) \ 00030 template <> \ 00031 void SGVector<complex128_t>::function() \ 00032 { \ 00033 SG_SERROR("SGVector::%s():: Not supported for complex128_t\n",\ 00034 #function);\ 00035 } 00036 00037 #define BOOL_ERROR_ONEARG(function) \ 00038 template <> \ 00039 void SGVector<bool>::function(bool a) \ 00040 { \ 00041 SG_SERROR("SGVector::%s():: Not supported for bool\n",\ 00042 #function);\ 00043 } 00044 00045 #define COMPLEX128_ERROR_ONEARG(function) \ 00046 template <> \ 00047 void SGVector<complex128_t>::function(complex128_t a) \ 00048 { \ 00049 SG_SERROR("SGVector::%s():: Not supported for complex128_t\n",\ 00050 #function);\ 00051 } 00052 00053 #define COMPLEX128_ERROR_TWOARGS(function) \ 00054 template <> \ 00055 void SGVector<complex128_t>::function(complex128_t a, complex128_t b) \ 00056 { \ 00057 SG_SERROR("SGVector::%s():: Not supported for complex128_t\n",\ 00058 #function);\ 00059 } 00060 00061 #define COMPLEX128_ERROR_THREEARGS(function) \ 00062 template <> \ 00063 void SGVector<complex128_t>::function(complex128_t a, complex128_t b,\ 00064 complex128_t c) \ 00065 { \ 00066 SG_SERROR("SGVector::%s():: Not supported for complex128_t\n",\ 00067 #function);\ 00068 } 00069 00070 namespace shogun { 00071 00072 template<class T> 00073 SGVector<T>::SGVector() : SGReferencedData() 00074 { 00075 init_data(); 00076 } 00077 00078 template<class T> 00079 SGVector<T>::SGVector(T* v, index_t len, bool ref_counting) 00080 : SGReferencedData(ref_counting), vector(v), vlen(len) 00081 { 00082 } 00083 00084 template<class T> 00085 SGVector<T>::SGVector(index_t len, bool ref_counting) 00086 : SGReferencedData(ref_counting), vlen(len) 00087 { 00088 vector=SG_MALLOC(T, len); 00089 } 00090 00091 template<class T> 00092 SGVector<T>::SGVector(const SGVector &orig) : SGReferencedData(orig) 00093 { 00094 copy_data(orig); 00095 } 00096 00097 template<class T> 00098 void SGVector<T>::set(SGVector<T> orig) 00099 { 00100 *this = SGVector<T>(orig); 00101 } 00102 00103 template<class T> 00104 SGVector<T>::~SGVector() 00105 { 00106 unref(); 00107 } 00108 00109 template<class T> 00110 void SGVector<T>::zero() 00111 { 00112 if (vector && vlen) 00113 set_const(0); 00114 } 00115 00116 template <> 00117 void SGVector<complex128_t>::zero() 00118 { 00119 if (vector && vlen) 00120 set_const(complex128_t(0.0)); 00121 } 00122 00123 template<class T> 00124 void SGVector<T>::set_const(T const_elem) 00125 { 00126 for (index_t i=0; i<vlen; i++) 00127 vector[i]=const_elem ; 00128 } 00129 00130 #if HAVE_CATLAS 00131 template<> 00132 void SGVector<float64_t>::set_const(float64_t const_elem) 00133 { 00134 catlas_dset(vlen, const_elem, vector, 1); 00135 } 00136 00137 template<> 00138 void SGVector<float32_t>::set_const(float32_t const_elem) 00139 { 00140 catlas_sset(vlen, const_elem, vector, 1); 00141 } 00142 #endif // HAVE_CATLAS 00143 00144 template<class T> 00145 void SGVector<T>::range_fill(T start) 00146 { 00147 range_fill_vector(vector, vlen, start); 00148 } 00149 00150 COMPLEX128_ERROR_ONEARG(range_fill) 00151 00152 template<class T> 00153 void SGVector<T>::random(T min_value, T max_value) 00154 { 00155 random_vector(vector, vlen, min_value, max_value); 00156 } 00157 00158 COMPLEX128_ERROR_TWOARGS(random) 00159 00160 template<class T> 00161 void SGVector<T>::randperm() 00162 { 00163 randperm(vector, vlen); 00164 } 00165 00166 COMPLEX128_ERROR_NOARG(randperm) 00167 00168 template <class T> 00169 void SGVector<T>::qsort() 00170 { 00171 CMath::qsort<T>(vector, vlen); 00172 } 00173 00174 COMPLEX128_ERROR_NOARG(qsort) 00175 00176 00177 template<class T> 00178 struct IndexSorter 00179 { 00181 IndexSorter(const SGVector<T> *vec) { data = vec->vector; } 00182 00184 bool operator() (index_t i, index_t j) const 00185 { 00186 return data[i] < data[j]; 00187 } 00188 00190 const T* data; 00191 }; 00192 00193 template<class T> 00194 SGVector<index_t> SGVector<T>::argsort() 00195 { 00196 IndexSorter<T> cmp(this); 00197 SGVector<index_t> idx(vlen); 00198 for (index_t i=0; i < vlen; ++i) 00199 idx[i] = i; 00200 00201 std::sort(idx.vector, idx.vector+vlen, cmp); 00202 00203 return idx; 00204 } 00205 00206 template <> 00207 SGVector<index_t> SGVector<complex128_t>::argsort() 00208 { 00209 SG_SERROR("SGVector::argsort():: Not supported for complex128_t\n"); 00210 SGVector<index_t> idx(vlen); 00211 return idx; 00212 } 00213 00214 template <class T> 00215 bool SGVector<T>::is_sorted() const 00216 { 00217 if (vlen < 2) 00218 return true; 00219 00220 for(int32_t i=1; i<vlen; i++) 00221 { 00222 if (vector[i-1] > vector[i]) 00223 return false; 00224 } 00225 00226 return true; 00227 } 00228 00229 template <> 00230 bool SGVector<complex128_t>::is_sorted() const 00231 { 00232 SG_SERROR("SGVector::is_sorted():: Not supported for complex128_t\n"); 00233 return false; 00234 } 00235 00236 template <class T> 00237 index_t SGVector<T>::find_position_to_insert(T element) 00238 { 00239 index_t i; 00240 for (i=0; i<vlen; ++i) 00241 { 00242 if (vector[i]>element) 00243 break; 00244 } 00245 return i; 00246 } 00247 00248 template <> 00249 index_t SGVector<complex128_t>::find_position_to_insert(complex128_t element) 00250 { 00251 SG_SERROR("SGVector::find_position_to_insert():: \ 00252 Not supported for complex128_t\n"); 00253 return index_t(-1); 00254 } 00255 00256 template<class T> 00257 SGVector<T> SGVector<T>::clone() const 00258 { 00259 return SGVector<T>(clone_vector(vector, vlen), vlen); 00260 } 00261 00262 template<class T> 00263 T* SGVector<T>::clone_vector(const T* vec, int32_t len) 00264 { 00265 T* result = SG_MALLOC(T, len); 00266 memcpy(result, vec, sizeof(T)*len); 00267 return result; 00268 } 00269 00270 template<class T> 00271 void SGVector<T>::fill_vector(T* vec, int32_t len, T value) 00272 { 00273 for (int32_t i=0; i<len; i++) 00274 vec[i]=value; 00275 } 00276 00277 template<class T> 00278 void SGVector<T>::range_fill_vector(T* vec, int32_t len, T start) 00279 { 00280 for (int32_t i=0; i<len; i++) 00281 vec[i]=i+start; 00282 } 00283 00284 template <> 00285 void SGVector<complex128_t>::range_fill_vector(complex128_t* vec, 00286 int32_t len, complex128_t start) 00287 { 00288 SG_SERROR("SGVector::range_fill_vector():: \ 00289 Not supported for complex128_t\n"); 00290 } 00291 00292 template<class T> 00293 const T& SGVector<T>::get_element(index_t index) 00294 { 00295 ASSERT(vector && (index>=0) && (index<vlen)) 00296 return vector[index]; 00297 } 00298 00299 template<class T> 00300 void SGVector<T>::set_element(const T& p_element, index_t index) 00301 { 00302 ASSERT(vector && (index>=0) && (index<vlen)) 00303 vector[index]=p_element; 00304 } 00305 00306 template<class T> 00307 void SGVector<T>::resize_vector(int32_t n) 00308 { 00309 vector=SG_REALLOC(T, vector, vlen, n); 00310 00311 if (n > vlen) 00312 memset(&vector[vlen], 0, (n-vlen)*sizeof(T)); 00313 vlen=n; 00314 } 00315 00317 template<class T> 00318 SGVector<T> SGVector<T>::operator+ (SGVector<T> x) 00319 { 00320 ASSERT(x.vector && vector) 00321 ASSERT(x.vlen == vlen) 00322 00323 SGVector<T> result=clone(); 00324 result.add(x); 00325 return result; 00326 } 00327 00328 template<class T> 00329 void SGVector<T>::add(const SGVector<T> x) 00330 { 00331 ASSERT(x.vector && vector) 00332 ASSERT(x.vlen == vlen) 00333 00334 for (int32_t i=0; i<vlen; i++) 00335 vector[i]+=x.vector[i]; 00336 } 00337 00338 template<class T> 00339 void SGVector<T>::add(const T x) 00340 { 00341 ASSERT(vector) 00342 00343 for (int32_t i=0; i<vlen; i++) 00344 vector[i]+=x; 00345 } 00346 00347 template<class T> 00348 void SGVector<T>::add(const SGSparseVector<T>& x) 00349 { 00350 if (x.features) 00351 { 00352 for (int32_t i=0; i < x.num_feat_entries; i++) 00353 { 00354 index_t idx = x.features[i].feat_index; 00355 ASSERT(idx < vlen) 00356 vector[idx] += x.features[i].entry; 00357 } 00358 } 00359 } 00360 00361 template<class T> 00362 void SGVector<T>::display_size() const 00363 { 00364 SG_SPRINT("SGVector '%p' of size: %d\n", vector, vlen) 00365 } 00366 00367 template<class T> 00368 void SGVector<T>::copy_data(const SGReferencedData &orig) 00369 { 00370 vector=((SGVector*)(&orig))->vector; 00371 vlen=((SGVector*)(&orig))->vlen; 00372 } 00373 00374 template<class T> 00375 void SGVector<T>::init_data() 00376 { 00377 vector=NULL; 00378 vlen=0; 00379 } 00380 00381 template<class T> 00382 void SGVector<T>::free_data() 00383 { 00384 SG_FREE(vector); 00385 vector=NULL; 00386 vlen=0; 00387 } 00388 00389 template<class T> 00390 bool SGVector<T>::equals(SGVector<T>& other) 00391 { 00392 if (other.vlen!=vlen) 00393 return false; 00394 00395 for (index_t i=0; i<vlen; ++i) 00396 { 00397 if (other.vector[i]!=vector[i]) 00398 return false; 00399 } 00400 00401 return true; 00402 } 00403 00404 template<class T> 00405 void SGVector<T>::display_vector(const char* name, 00406 const char* prefix) const 00407 { 00408 display_vector(vector, vlen, name, prefix); 00409 } 00410 00411 template <class T> 00412 void SGVector<T>::display_vector(const SGVector<T> vector, const char* name, 00413 const char* prefix) 00414 { 00415 vector.display_vector(prefix); 00416 } 00417 00418 template <> 00419 void SGVector<bool>::display_vector(const bool* vector, int32_t n, const char* name, 00420 const char* prefix) 00421 { 00422 ASSERT(n>=0) 00423 SG_SPRINT("%s%s=[", prefix, name) 00424 for (int32_t i=0; i<n; i++) 00425 SG_SPRINT("%s%d%s", prefix, vector[i] ? 1 : 0, i==n-1? "" : ",") 00426 SG_SPRINT("%s]\n", prefix) 00427 } 00428 00429 template <> 00430 void SGVector<char>::display_vector(const char* vector, int32_t n, const char* name, 00431 const char* prefix) 00432 { 00433 ASSERT(n>=0) 00434 SG_SPRINT("%s%s=[", prefix, name) 00435 for (int32_t i=0; i<n; i++) 00436 SG_SPRINT("%s%c%s", prefix, vector[i], i==n-1? "" : ",") 00437 SG_SPRINT("%s]\n", prefix) 00438 } 00439 00440 template <> 00441 void SGVector<uint8_t>::display_vector(const uint8_t* vector, int32_t n, const char* name, 00442 const char* prefix) 00443 { 00444 ASSERT(n>=0) 00445 SG_SPRINT("%s%s=[", prefix, name) 00446 for (int32_t i=0; i<n; i++) 00447 SG_SPRINT("%s%u%s", prefix, vector[i], i==n-1? "" : ",") 00448 SG_SPRINT("%s]\n", prefix) 00449 } 00450 00451 template <> 00452 void SGVector<int8_t>::display_vector(const int8_t* vector, int32_t n, const char* name, 00453 const char* prefix) 00454 { 00455 ASSERT(n>=0) 00456 SG_SPRINT("%s%s=[", prefix, name) 00457 for (int32_t i=0; i<n; i++) 00458 SG_SPRINT("%s%d%s", prefix, vector[i], i==n-1? "" : ",") 00459 SG_SPRINT("%s]\n", prefix) 00460 } 00461 00462 template <> 00463 void SGVector<uint16_t>::display_vector(const uint16_t* vector, int32_t n, const char* name, 00464 const char* prefix) 00465 { 00466 ASSERT(n>=0) 00467 SG_SPRINT("%s%s=[", prefix, name) 00468 for (int32_t i=0; i<n; i++) 00469 SG_SPRINT("%s%u%s", prefix, vector[i], i==n-1? "" : ",") 00470 SG_SPRINT("%s]\n", prefix) 00471 } 00472 00473 template <> 00474 void SGVector<int16_t>::display_vector(const int16_t* vector, int32_t n, const char* name, 00475 const char* prefix) 00476 { 00477 ASSERT(n>=0) 00478 SG_SPRINT("%s%s=[", prefix, name) 00479 for (int32_t i=0; i<n; i++) 00480 SG_SPRINT("%s%d%s", prefix, vector[i], i==n-1? "" : ",") 00481 SG_SPRINT("%s]\n", prefix) 00482 } 00483 00484 template <> 00485 void SGVector<int32_t>::display_vector(const int32_t* vector, int32_t n, const char* name, 00486 const char* prefix) 00487 { 00488 ASSERT(n>=0) 00489 SG_SPRINT("%s%s=[", prefix, name) 00490 for (int32_t i=0; i<n; i++) 00491 SG_SPRINT("%s%d%s", prefix, vector[i], i==n-1? "" : ",") 00492 SG_SPRINT("%s]\n", prefix) 00493 } 00494 00495 template <> 00496 void SGVector<uint32_t>::display_vector(const uint32_t* vector, int32_t n, const char* name, 00497 const char* prefix) 00498 { 00499 ASSERT(n>=0) 00500 SG_SPRINT("%s%s=[", prefix, name) 00501 for (int32_t i=0; i<n; i++) 00502 SG_SPRINT("%s%u%s", prefix, vector[i], i==n-1? "" : ",") 00503 SG_SPRINT("%s]\n", prefix) 00504 } 00505 00506 00507 template <> 00508 void SGVector<int64_t>::display_vector(const int64_t* vector, int32_t n, const char* name, 00509 const char* prefix) 00510 { 00511 ASSERT(n>=0) 00512 SG_SPRINT("%s%s=[", prefix, name) 00513 for (int32_t i=0; i<n; i++) 00514 SG_SPRINT("%s%lld%s", prefix, vector[i], i==n-1? "" : ",") 00515 SG_SPRINT("%s]\n", prefix) 00516 } 00517 00518 template <> 00519 void SGVector<uint64_t>::display_vector(const uint64_t* vector, int32_t n, const char* name, 00520 const char* prefix) 00521 { 00522 ASSERT(n>=0) 00523 SG_SPRINT("%s%s=[", prefix, name) 00524 for (int32_t i=0; i<n; i++) 00525 SG_SPRINT("%s%llu%s", prefix, vector[i], i==n-1? "" : ",") 00526 SG_SPRINT("%s]\n", prefix) 00527 } 00528 00529 template <> 00530 void SGVector<float32_t>::display_vector(const float32_t* vector, int32_t n, const char* name, 00531 const char* prefix) 00532 { 00533 ASSERT(n>=0) 00534 SG_SPRINT("%s%s=[", prefix, name) 00535 for (int32_t i=0; i<n; i++) 00536 SG_SPRINT("%s%g%s", prefix, vector[i], i==n-1? "" : ",") 00537 SG_SPRINT("%s]\n", prefix) 00538 } 00539 00540 template <> 00541 void SGVector<float64_t>::display_vector(const float64_t* vector, int32_t n, const char* name, 00542 const char* prefix) 00543 { 00544 ASSERT(n>=0) 00545 SG_SPRINT("%s%s=[", prefix, name) 00546 for (int32_t i=0; i<n; i++) 00547 SG_SPRINT("%s%.18g%s", prefix, vector[i], i==n-1? "" : ",") 00548 SG_SPRINT("%s]\n", prefix) 00549 } 00550 00551 template <> 00552 void SGVector<floatmax_t>::display_vector(const floatmax_t* vector, int32_t n, 00553 const char* name, const char* prefix) 00554 { 00555 ASSERT(n>=0) 00556 SG_SPRINT("%s%s=[", prefix, name) 00557 for (int32_t i=0; i<n; i++) 00558 { 00559 SG_SPRINT("%s%.36Lg%s", prefix, (long double) vector[i], 00560 i==n-1? "" : ","); 00561 } 00562 SG_SPRINT("%s]\n", prefix) 00563 } 00564 00565 template <> 00566 void SGVector<complex128_t>::display_vector(const complex128_t* vector, int32_t n, 00567 const char* name, const char* prefix) 00568 { 00569 ASSERT(n>=0) 00570 SG_SPRINT("%s%s=[", prefix, name) 00571 for (int32_t i=0; i<n; i++) 00572 { 00573 SG_SPRINT("%s(%.36lg+i%.36lg)%s", prefix, vector[i].real(), 00574 vector[i].imag(), i==n-1? "" : ","); 00575 } 00576 SG_SPRINT("%s]\n", prefix) 00577 } 00578 00579 template <class T> 00580 void SGVector<T>::vec1_plus_scalar_times_vec2(T* vec1, 00581 const T scalar, const T* vec2, int32_t n) 00582 { 00583 for (int32_t i=0; i<n; i++) 00584 vec1[i]+=scalar*vec2[i]; 00585 } 00586 00587 template <> 00588 void SGVector<float64_t>::vec1_plus_scalar_times_vec2(float64_t* vec1, 00589 float64_t scalar, const float64_t* vec2, int32_t n) 00590 { 00591 #ifdef HAVE_LAPACK 00592 int32_t skip=1; 00593 cblas_daxpy(n, scalar, vec2, skip, vec1, skip); 00594 #else 00595 for (int32_t i=0; i<n; i++) 00596 vec1[i]+=scalar*vec2[i]; 00597 #endif 00598 } 00599 00600 template <> 00601 void SGVector<float32_t>::vec1_plus_scalar_times_vec2(float32_t* vec1, 00602 float32_t scalar, const float32_t* vec2, int32_t n) 00603 { 00604 #ifdef HAVE_LAPACK 00605 int32_t skip=1; 00606 cblas_saxpy(n, scalar, vec2, skip, vec1, skip); 00607 #else 00608 for (int32_t i=0; i<n; i++) 00609 vec1[i]+=scalar*vec2[i]; 00610 #endif 00611 } 00612 00613 template <class T> 00614 float64_t SGVector<T>::dot(const float64_t* v1, const float64_t* v2, int32_t n) 00615 { 00616 float64_t r=0; 00617 #ifdef HAVE_EIGEN3 00618 Eigen::Map<const Eigen::VectorXd> ev1(v1,n); 00619 Eigen::Map<const Eigen::VectorXd> ev2(v2,n); 00620 r = ev1.dot(ev2); 00621 #elif HAVE_LAPACK 00622 int32_t skip=1; 00623 r = cblas_ddot(n, v1, skip, v2, skip); 00624 #else 00625 for (int32_t i=0; i<n; i++) 00626 r+=v1[i]*v2[i]; 00627 #endif 00628 return r; 00629 } 00630 00631 template <class T> 00632 float32_t SGVector<T>::dot(const float32_t* v1, const float32_t* v2, int32_t n) 00633 { 00634 float32_t r=0; 00635 #ifdef HAVE_EIGEN3 00636 Eigen::Map<const Eigen::VectorXf> ev1(v1,n); 00637 Eigen::Map<const Eigen::VectorXf> ev2(v2,n); 00638 r = ev1.dot(ev2); 00639 #elif HAVE_LAPACK 00640 int32_t skip=1; 00641 r = cblas_sdot(n, v1, skip, v2, skip); 00642 #else 00643 for (int32_t i=0; i<n; i++) 00644 r+=v1[i]*v2[i]; 00645 #endif 00646 return r; 00647 } 00648 00649 template <class T> 00650 void SGVector<T>::random_vector(T* vec, int32_t len, T min_value, T max_value) 00651 { 00652 for (int32_t i=0; i<len; i++) 00653 vec[i]=CMath::random(min_value, max_value); 00654 } 00655 00656 template <> 00657 void SGVector<complex128_t>::random_vector(complex128_t* vec, int32_t len, 00658 complex128_t min_value, complex128_t max_value) 00659 { 00660 SG_SNOTIMPLEMENTED 00661 } 00662 00663 template <class T> 00664 SGVector<T> SGVector<T>::randperm_vec(int32_t n) 00665 { 00666 return SGVector<T>(randperm(n), n); 00667 } 00668 00669 template <> 00670 SGVector<complex128_t> SGVector<complex128_t>::randperm_vec(int32_t n) 00671 { 00672 SG_SNOTIMPLEMENTED 00673 SGVector<complex128_t> perm(n); 00674 return perm; 00675 } 00676 00677 template <class T> 00678 T* SGVector<T>::randperm(int32_t n) 00679 { 00680 T* perm = SG_MALLOC(T, n); 00681 randperm(perm, n); 00682 00683 return perm; 00684 } 00685 00686 template <> 00687 complex128_t* SGVector<complex128_t>::randperm(int32_t n) 00688 { 00689 SG_SNOTIMPLEMENTED 00690 SGVector<complex128_t> perm(n); 00691 return perm.vector; 00692 } 00693 00694 template <class T> 00695 void SGVector<T>::randperm(T* perm, int32_t n) 00696 { 00697 for (int32_t i = 0; i < n; i++) 00698 perm[i] = i; 00699 permute(perm,n); 00700 } 00701 00702 template <> 00703 void SGVector<complex128_t>::randperm(complex128_t* perm, int32_t n) 00704 { 00705 SG_SNOTIMPLEMENTED 00706 } 00707 00709 template <class T> 00710 void SGVector<T>::permute(T* vec, int32_t n) 00711 { 00712 for (int32_t i = 0; i < n; i++) 00713 CMath::swap(vec[i], vec[CMath::random(i, n-1)]); 00714 } 00715 00716 template <class T> 00717 void SGVector<T>::permute(T* vec, int32_t n, CRandom * rand) 00718 { 00719 for (int32_t i = 0; i < n; i++) 00720 CMath::swap(vec[i], vec[rand->random(i, n-1)]); 00721 } 00722 00723 template<class T> 00724 void SGVector<T>::permute() 00725 { 00726 SGVector<T>::permute(vector, vlen); 00727 } 00728 00729 template<class T> 00730 void SGVector<T>::permute(CRandom * rand) 00731 { 00732 SGVector<T>::permute(vector, vlen, rand); 00733 } 00734 00735 template <class T> 00736 void SGVector<T>::permute_vector(SGVector<T> vec) 00737 { 00738 for (index_t i=0; i<vec.vlen; ++i) 00739 { 00740 CMath::swap(vec.vector[i], 00741 vec.vector[CMath::random(i, vec.vlen-1)]); 00742 } 00743 } 00744 00745 template <> 00746 bool SGVector<bool>::twonorm(const bool* x, int32_t len) 00747 { 00748 SG_SNOTIMPLEMENTED 00749 return false; 00750 } 00751 00752 template <> 00753 char SGVector<char>::twonorm(const char* x, int32_t len) 00754 { 00755 SG_SNOTIMPLEMENTED 00756 return '\0'; 00757 } 00758 00759 template <> 00760 int8_t SGVector<int8_t>::twonorm(const int8_t* x, int32_t len) 00761 { 00762 float64_t result=0; 00763 for (int32_t i=0; i<len; i++) 00764 result+=x[i]*x[i]; 00765 00766 return CMath::sqrt(result); 00767 } 00768 00769 template <> 00770 uint8_t SGVector<uint8_t>::twonorm(const uint8_t* x, int32_t len) 00771 { 00772 float64_t result=0; 00773 for (int32_t i=0; i<len; i++) 00774 result+=x[i]*x[i]; 00775 00776 return CMath::sqrt(result); 00777 } 00778 00779 template <> 00780 int16_t SGVector<int16_t>::twonorm(const int16_t* x, int32_t len) 00781 { 00782 float64_t result=0; 00783 for (int32_t i=0; i<len; i++) 00784 result+=x[i]*x[i]; 00785 00786 return CMath::sqrt(result); 00787 } 00788 00789 template <> 00790 uint16_t SGVector<uint16_t>::twonorm(const uint16_t* x, int32_t len) 00791 { 00792 float64_t result=0; 00793 for (int32_t i=0; i<len; i++) 00794 result+=x[i]*x[i]; 00795 00796 return CMath::sqrt(result); 00797 } 00798 00799 template <> 00800 int32_t SGVector<int32_t>::twonorm(const int32_t* x, int32_t len) 00801 { 00802 float64_t result=0; 00803 for (int32_t i=0; i<len; i++) 00804 result+=x[i]*x[i]; 00805 00806 return CMath::sqrt(result); 00807 } 00808 00809 template <> 00810 uint32_t SGVector<uint32_t>::twonorm(const uint32_t* x, int32_t len) 00811 { 00812 float64_t result=0; 00813 for (int32_t i=0; i<len; i++) 00814 result+=x[i]*x[i]; 00815 00816 return CMath::sqrt(result); 00817 } 00818 00819 template <> 00820 int64_t SGVector<int64_t>::twonorm(const int64_t* x, int32_t len) 00821 { 00822 float64_t result=0; 00823 for (int32_t i=0; i<len; i++) 00824 result+=x[i]*x[i]; 00825 00826 return CMath::sqrt(result); 00827 } 00828 00829 template <> 00830 uint64_t SGVector<uint64_t>::twonorm(const uint64_t* x, int32_t len) 00831 { 00832 float64_t result=0; 00833 for (int32_t i=0; i<len; i++) 00834 result+=x[i]*x[i]; 00835 00836 return CMath::sqrt(result); 00837 } 00838 00839 template <> 00840 float32_t SGVector<float32_t>::twonorm(const float32_t* x, int32_t len) 00841 { 00842 float64_t result=0; 00843 for (int32_t i=0; i<len; i++) 00844 result+=x[i]*x[i]; 00845 00846 return CMath::sqrt(result); 00847 } 00848 00849 template <> 00850 float64_t SGVector<float64_t>::twonorm(const float64_t* v, int32_t n) 00851 { 00852 float64_t norm = 0.0; 00853 #ifdef HAVE_LAPACK 00854 norm = cblas_dnrm2(n, v, 1); 00855 #else 00856 norm = CMath::sqrt(SGVector::dot(v, v, n)); 00857 #endif 00858 return norm; 00859 } 00860 00861 template <> 00862 floatmax_t SGVector<floatmax_t>::twonorm(const floatmax_t* x, int32_t len) 00863 { 00864 float64_t result=0; 00865 for (int32_t i=0; i<len; i++) 00866 result+=x[i]*x[i]; 00867 00868 return CMath::sqrt(result); 00869 } 00870 00871 template <> 00872 complex128_t SGVector<complex128_t>::twonorm(const complex128_t* x, int32_t len) 00873 { 00874 complex128_t result(0.0); 00875 for (int32_t i=0; i<len; i++) 00876 result+=x[i]*x[i]; 00877 00878 return CMath::sqrt(result); 00879 } 00880 00881 template <class T> 00882 float64_t SGVector<T>::onenorm(T* x, int32_t len) 00883 { 00884 float64_t result=0; 00885 for (int32_t i=0;i<len; ++i) 00886 result+=CMath::abs(x[i]); 00887 00888 return result; 00889 } 00890 00892 template <class T> 00893 T SGVector<T>::qsq(T* x, int32_t len, float64_t q) 00894 { 00895 float64_t result=0; 00896 for (int32_t i=0; i<len; i++) 00897 result+=CMath::pow(fabs(x[i]), q); 00898 00899 return result; 00900 } 00901 00902 template <> 00903 complex128_t SGVector<complex128_t>::qsq(complex128_t* x, int32_t len, float64_t q) 00904 { 00905 SG_SNOTIMPLEMENTED 00906 return complex128_t(0.0); 00907 } 00908 00910 template <class T> 00911 T SGVector<T>::qnorm(T* x, int32_t len, float64_t q) 00912 { 00913 ASSERT(q!=0) 00914 return CMath::pow((float64_t) qsq(x, len, q), 1.0/q); 00915 } 00916 00917 template <> 00918 complex128_t SGVector<complex128_t>::qnorm(complex128_t* x, int32_t len, float64_t q) 00919 { 00920 SG_SNOTIMPLEMENTED 00921 return complex128_t(0.0); 00922 } 00923 00925 template <class T> 00926 T SGVector<T>::min(T* vec, int32_t len) 00927 { 00928 ASSERT(len>0) 00929 T minv=vec[0]; 00930 00931 for (int32_t i=1; i<len; i++) 00932 minv=CMath::min(vec[i], minv); 00933 00934 return minv; 00935 } 00936 00937 #ifdef HAVE_LAPACK 00938 template <> 00939 float64_t SGVector<float64_t>::max_abs(float64_t* vec, int32_t len) 00940 { 00941 ASSERT(len>0) 00942 int32_t skip = 1; 00943 int32_t idx = cblas_idamax(len, vec, skip); 00944 00945 return CMath::abs(vec[idx]); 00946 } 00947 00948 template <> 00949 float32_t SGVector<float32_t>::max_abs(float32_t* vec, int32_t len) 00950 { 00951 ASSERT(len>0) 00952 int32_t skip = 1; 00953 int32_t idx = cblas_isamax(len, vec, skip); 00954 00955 return CMath::abs(vec[idx]); 00956 } 00957 #endif 00958 00960 template <class T> 00961 T SGVector<T>::max_abs(T* vec, int32_t len) 00962 { 00963 ASSERT(len>0) 00964 T maxv=CMath::abs(vec[0]); 00965 00966 for (int32_t i=1; i<len; i++) 00967 maxv=CMath::max(CMath::abs(vec[i]), maxv); 00968 00969 return maxv; 00970 } 00971 00972 template <> 00973 complex128_t SGVector<complex128_t>::max_abs(complex128_t* vec, int32_t len) 00974 { 00975 SG_SNOTIMPLEMENTED 00976 return complex128_t(0.0); 00977 } 00978 00980 template <class T> 00981 T SGVector<T>::max(T* vec, int32_t len) 00982 { 00983 ASSERT(len>0) 00984 T maxv=vec[0]; 00985 00986 for (int32_t i=1; i<len; i++) 00987 maxv=CMath::max(vec[i], maxv); 00988 00989 return maxv; 00990 } 00991 00992 #ifdef HAVE_LAPACK 00993 template <> 00994 int32_t SGVector<float64_t>::arg_max_abs(float64_t* vec, int32_t inc, int32_t len, float64_t* maxv_ptr) 00995 { 00996 ASSERT(len>0 || inc > 0) 00997 int32_t idx = cblas_idamax(len, vec, inc); 00998 00999 if (maxv_ptr != NULL) 01000 *maxv_ptr = CMath::abs(vec[idx*inc]); 01001 01002 return idx; 01003 } 01004 01005 template <> 01006 int32_t SGVector<float32_t>::arg_max_abs(float32_t* vec, int32_t inc, int32_t len, float32_t* maxv_ptr) 01007 { 01008 ASSERT(len>0 || inc > 0) 01009 int32_t idx = cblas_isamax(len, vec, inc); 01010 01011 if (maxv_ptr != NULL) 01012 *maxv_ptr = CMath::abs(vec[idx*inc]); 01013 01014 return idx; 01015 } 01016 #endif 01017 01018 template <class T> 01019 int32_t SGVector<T>::arg_max_abs(T * vec, int32_t inc, int32_t len, T * maxv_ptr) 01020 { 01021 ASSERT(len > 0 || inc > 0) 01022 01023 T maxv = CMath::abs(vec[0]); 01024 int32_t maxIdx = 0; 01025 01026 for (int32_t i = 1, j = inc ; i < len ; i++, j += inc) 01027 { 01028 if (CMath::abs(vec[j]) > maxv) 01029 maxv = CMath::abs(vec[j]), maxIdx = i; 01030 } 01031 01032 if (maxv_ptr != NULL) 01033 *maxv_ptr = maxv; 01034 01035 return maxIdx; 01036 } 01037 01038 template <> 01039 int32_t SGVector<complex128_t>::arg_max_abs(complex128_t * vec, int32_t inc, 01040 int32_t len, complex128_t * maxv_ptr) 01041 { 01042 int32_t maxIdx = 0; 01043 SG_SERROR("SGVector::arg_max_abs():: Not supported for complex128_t\n"); 01044 return maxIdx; 01045 } 01046 01047 template <class T> 01048 int32_t SGVector<T>::arg_max(T * vec, int32_t inc, int32_t len, T * maxv_ptr) 01049 { 01050 ASSERT(len > 0 || inc > 0) 01051 01052 T maxv = vec[0]; 01053 int32_t maxIdx = 0; 01054 01055 for (int32_t i = 1, j = inc ; i < len ; i++, j += inc) 01056 { 01057 if (vec[j] > maxv) 01058 maxv = vec[j], maxIdx = i; 01059 } 01060 01061 if (maxv_ptr != NULL) 01062 *maxv_ptr = maxv; 01063 01064 return maxIdx; 01065 } 01066 01067 template <> 01068 int32_t SGVector<complex128_t>::arg_max(complex128_t * vec, int32_t inc, 01069 int32_t len, complex128_t * maxv_ptr) 01070 { 01071 int32_t maxIdx=0; 01072 SG_SERROR("SGVector::arg_max():: Not supported for complex128_t\n"); 01073 return maxIdx; 01074 } 01075 01076 01078 template <class T> 01079 int32_t SGVector<T>::arg_min(T * vec, int32_t inc, int32_t len, T * minv_ptr) 01080 { 01081 ASSERT(len > 0 || inc > 0) 01082 01083 T minv = vec[0]; 01084 int32_t minIdx = 0; 01085 01086 for (int32_t i = 1, j = inc ; i < len ; i++, j += inc) 01087 { 01088 if (vec[j] < minv) 01089 minv = vec[j], minIdx = i; 01090 } 01091 01092 if (minv_ptr != NULL) 01093 *minv_ptr = minv; 01094 01095 return minIdx; 01096 } 01097 01098 template <> 01099 int32_t SGVector<complex128_t>::arg_min(complex128_t * vec, int32_t inc, 01100 int32_t len, complex128_t * minv_ptr) 01101 { 01102 int32_t minIdx=0; 01103 SG_SERROR("SGVector::arg_min():: Not supported for complex128_t\n"); 01104 return minIdx; 01105 } 01106 01108 template <class T> 01109 T SGVector<T>::sum_abs(T* vec, int32_t len) 01110 { 01111 T result=0; 01112 for (int32_t i=0; i<len; i++) 01113 result+=CMath::abs(vec[i]); 01114 01115 return result; 01116 } 01117 01118 #if HAVE_LAPACK 01119 template <> 01120 float64_t SGVector<float64_t>::sum_abs(float64_t* vec, int32_t len) 01121 { 01122 float64_t result=0; 01123 result = cblas_dasum(len, vec, 1); 01124 return result; 01125 } 01126 01127 template <> 01128 float32_t SGVector<float32_t>::sum_abs(float32_t* vec, int32_t len) 01129 { 01130 float32_t result=0; 01131 result = cblas_sasum(len, vec, 1); 01132 return result; 01133 } 01134 #endif 01135 01137 template <class T> 01138 bool SGVector<T>::fequal(T x, T y, float64_t precision) 01139 { 01140 return CMath::abs(x-y)<precision; 01141 } 01142 01143 template <class T> 01144 int32_t SGVector<T>::unique(T* output, int32_t size) 01145 { 01146 CMath::qsort<T>(output, size); 01147 int32_t j=0; 01148 01149 for (int32_t i=0; i<size; i++) 01150 { 01151 if (i==0 || output[i]!=output[i-1]) 01152 output[j++]=output[i]; 01153 } 01154 return j; 01155 } 01156 01157 template <> 01158 int32_t SGVector<complex128_t>::unique(complex128_t* output, int32_t size) 01159 { 01160 int32_t j=0; 01161 SG_SERROR("SGVector::unique():: Not supported for complex128_t\n"); 01162 return j; 01163 } 01164 01165 template <class T> 01166 SGVector<index_t> SGVector<T>::find(T elem) 01167 { 01168 SGVector<index_t> idx(vlen); 01169 index_t k=0; 01170 01171 for (index_t i=0; i < vlen; ++i) 01172 if (vector[i] == elem) 01173 idx[k++] = i; 01174 idx.vlen = k; 01175 return idx; 01176 } 01177 01178 template<class T> 01179 void SGVector<T>::scale_vector(T alpha, T* vec, int32_t len) 01180 { 01181 for (int32_t i=0; i<len; i++) 01182 vec[i]*=alpha; 01183 } 01184 01185 #ifdef HAVE_LAPACK 01186 template<> 01187 void SGVector<float64_t>::scale_vector(float64_t alpha, float64_t* vec, int32_t len) 01188 { 01189 cblas_dscal(len, alpha, vec, 1); 01190 } 01191 01192 template<> 01193 void SGVector<float32_t>::scale_vector(float32_t alpha, float32_t* vec, int32_t len) 01194 { 01195 cblas_sscal(len, alpha, vec, 1); 01196 } 01197 #endif 01198 01199 template<class T> 01200 void SGVector<T>::scale(T alpha) 01201 { 01202 scale_vector(alpha, vector, vlen); 01203 } 01204 01205 template<class T> float64_t SGVector<T>::mean() const 01206 { 01207 float64_t cum = 0; 01208 01209 for ( index_t i = 0 ; i < vlen ; ++i ) 01210 cum += vector[i]; 01211 01212 return cum/vlen; 01213 } 01214 01215 template <> 01216 float64_t SGVector<complex128_t>::mean() const 01217 { 01218 SG_SNOTIMPLEMENTED 01219 return float64_t(0.0); 01220 } 01221 01222 template<class T> void SGVector<T>::load(CFile* loader) 01223 { 01224 ASSERT(loader) 01225 unref(); 01226 01227 SG_SET_LOCALE_C; 01228 SGVector<T> vec; 01229 loader->get_vector(vec.vector, vec.vlen); 01230 copy_data(vec); 01231 copy_refcount(vec); 01232 ref(); 01233 SG_RESET_LOCALE; 01234 } 01235 01236 template<> 01237 void SGVector<complex128_t>::load(CFile* loader) 01238 { 01239 SG_SERROR("SGVector::load():: Not supported for complex128_t\n"); 01240 } 01241 01242 template<class T> void SGVector<T>::save(CFile* saver) 01243 { 01244 ASSERT(saver) 01245 01246 SG_SET_LOCALE_C; 01247 saver->set_vector(vector, vlen); 01248 SG_RESET_LOCALE; 01249 } 01250 01251 template<> 01252 void SGVector<complex128_t>::save(CFile* saver) 01253 { 01254 SG_SERROR("SGVector::save():: Not supported for complex128_t\n"); 01255 } 01256 01257 01258 #define MATHOP(op) \ 01259 template <class T> void SGVector<T>::op() \ 01260 { \ 01261 for (int32_t i=0; i<vlen; i++) \ 01262 vector[i]=(T) CMath::op((double) vector[i]); \ 01263 } 01264 01265 MATHOP(abs) 01266 MATHOP(acos) 01267 MATHOP(asin) 01268 MATHOP(atan) 01269 MATHOP(cos) 01270 MATHOP(cosh) 01271 MATHOP(exp) 01272 MATHOP(log) 01273 MATHOP(log10) 01274 MATHOP(sin) 01275 MATHOP(sinh) 01276 MATHOP(sqrt) 01277 MATHOP(tan) 01278 MATHOP(tanh) 01279 #undef MATHOP 01280 01281 #define COMPLEX128_MATHOP(op) \ 01282 template <>\ 01283 void SGVector<complex128_t>::op() \ 01284 { \ 01285 for (int32_t i=0; i<vlen; i++) \ 01286 vector[i]=complex128_t(CMath::op(vector[i])); \ 01287 } 01288 01289 COMPLEX128_MATHOP(abs) 01290 COMPLEX128_MATHOP(sin) 01291 COMPLEX128_MATHOP(cos) 01292 COMPLEX128_MATHOP(tan) 01293 COMPLEX128_MATHOP(sinh) 01294 COMPLEX128_MATHOP(cosh) 01295 COMPLEX128_MATHOP(tanh) 01296 COMPLEX128_MATHOP(exp) 01297 COMPLEX128_MATHOP(log) 01298 COMPLEX128_MATHOP(log10) 01299 COMPLEX128_MATHOP(sqrt) 01300 #undef COMPLEX128_MATHOP 01301 01302 #define COMPLEX128_MATHOP_NOTIMPLEMENTED(op) \ 01303 template <>\ 01304 void SGVector<complex128_t>::op() \ 01305 { \ 01306 SG_SERROR("SGVector::%s():: Not supported for complex128_t\n",#op);\ 01307 } 01308 01309 COMPLEX128_MATHOP_NOTIMPLEMENTED(asin) 01310 COMPLEX128_MATHOP_NOTIMPLEMENTED(acos) 01311 COMPLEX128_MATHOP_NOTIMPLEMENTED(atan) 01312 #undef COMPLEX128_MATHOP_NOTIMPLEMENTED 01313 01314 template <class T> void SGVector<T>::atan2(T x) 01315 { 01316 for (int32_t i=0; i<vlen; i++) 01317 vector[i]=CMath::atan2(vector[i], x); 01318 } 01319 01320 BOOL_ERROR_ONEARG(atan2) 01321 COMPLEX128_ERROR_ONEARG(atan2) 01322 01323 template <class T> void SGVector<T>::pow(T q) 01324 { 01325 for (int32_t i=0; i<vlen; i++) 01326 vector[i]=(T) CMath::pow((double) vector[i], (double) q); 01327 } 01328 01329 template <class T> 01330 SGVector<float64_t> SGVector<T>::linspace_vec(T start, T end, int32_t n) 01331 { 01332 return SGVector<float64_t>(linspace(start, end, n), n); 01333 } 01334 01335 template <class T> 01336 float64_t* SGVector<T>::linspace(T start, T end, int32_t n) 01337 { 01338 float64_t* output = SG_MALLOC(float64_t, n); 01339 CMath::linspace(output, start, end, n); 01340 01341 return output; 01342 } 01343 01344 template <> 01345 float64_t* SGVector<complex128_t>::linspace(complex128_t start, complex128_t end, int32_t n) 01346 { 01347 float64_t* output = SG_MALLOC(float64_t, n); 01348 SG_SERROR("SGVector::linspace():: Not supported for complex128_t\n"); 01349 return output; 01350 } 01351 01352 template <> 01353 void SGVector<complex128_t>::pow(complex128_t q) 01354 { 01355 for (int32_t i=0; i<vlen; i++) 01356 vector[i]=CMath::pow(vector[i], q); 01357 } 01358 01359 template <class T> SGVector<float64_t> SGVector<T>::get_real() 01360 { 01361 SGVector<float64_t> real(vlen); 01362 for (int32_t i=0; i<vlen; i++) 01363 real[i]=CMath::real(vector[i]); 01364 return real; 01365 } 01366 01367 template <class T> SGVector<float64_t> SGVector<T>::get_imag() 01368 { 01369 SGVector<float64_t> imag(vlen); 01370 for (int32_t i=0; i<vlen; i++) 01371 imag[i]=CMath::imag(vector[i]); 01372 return imag; 01373 } 01374 01375 template <class T> 01376 SGMatrix<T> SGVector<T>::convert_to_matrix(SGVector<T> vector, 01377 index_t nrows, index_t ncols, bool fortran_order) 01378 { 01379 if (nrows*ncols>vector.size()) 01380 SG_SERROR("SGVector::convert_to_matrix():: Dimensions mismatch\n"); 01381 01382 T* data=NULL; 01383 SGVector<T>::convert_to_matrix(data, nrows, ncols, vector.vector, vector.vlen, fortran_order); 01384 01385 SGMatrix<T> matrix=SGMatrix<T>(data, nrows, ncols); 01386 return matrix; 01387 } 01388 01389 template <class T> 01390 void SGVector<T>::convert_to_matrix(T*& matrix, index_t nrows, index_t ncols, const T* vector, int32_t vlen, bool fortran_order) 01391 { 01392 if (nrows*ncols>vlen) 01393 SG_SERROR("SGVector::convert_to_matrix():: Dimensions mismatch\n"); 01394 01395 if (matrix!=NULL) 01396 SG_FREE(matrix); 01397 matrix=SG_MALLOC(T, nrows*ncols); 01398 01399 if (fortran_order) 01400 { 01401 for (index_t i=0; i<ncols*nrows; i++) 01402 matrix[i]=vector[i]; 01403 } 01404 else 01405 { 01406 for (index_t i=0; i<nrows; i++) 01407 { 01408 for (index_t j=0; j<ncols; j++) 01409 matrix[i+j*nrows]=vector[j+i*ncols]; 01410 } 01411 } 01412 } 01413 01414 #define UNDEFINED(function, type) \ 01415 template <> \ 01416 SGVector<float64_t> SGVector<type>::function() \ 01417 { \ 01418 SG_SERROR("SGVector::%s():: Not supported for %s\n", \ 01419 #function, #type); \ 01420 SGVector<float64_t> ret(vlen); \ 01421 return ret; \ 01422 } 01423 01424 UNDEFINED(get_real, bool) 01425 UNDEFINED(get_real, char) 01426 UNDEFINED(get_real, int8_t) 01427 UNDEFINED(get_real, uint8_t) 01428 UNDEFINED(get_real, int16_t) 01429 UNDEFINED(get_real, uint16_t) 01430 UNDEFINED(get_real, int32_t) 01431 UNDEFINED(get_real, uint32_t) 01432 UNDEFINED(get_real, int64_t) 01433 UNDEFINED(get_real, uint64_t) 01434 UNDEFINED(get_real, float32_t) 01435 UNDEFINED(get_real, float64_t) 01436 UNDEFINED(get_real, floatmax_t) 01437 UNDEFINED(get_imag, bool) 01438 UNDEFINED(get_imag, char) 01439 UNDEFINED(get_imag, int8_t) 01440 UNDEFINED(get_imag, uint8_t) 01441 UNDEFINED(get_imag, int16_t) 01442 UNDEFINED(get_imag, uint16_t) 01443 UNDEFINED(get_imag, int32_t) 01444 UNDEFINED(get_imag, uint32_t) 01445 UNDEFINED(get_imag, int64_t) 01446 UNDEFINED(get_imag, uint64_t) 01447 UNDEFINED(get_imag, float32_t) 01448 UNDEFINED(get_imag, float64_t) 01449 UNDEFINED(get_imag, floatmax_t) 01450 #undef UNDEFINED 01451 01452 template class SGVector<bool>; 01453 template class SGVector<char>; 01454 template class SGVector<int8_t>; 01455 template class SGVector<uint8_t>; 01456 template class SGVector<int16_t>; 01457 template class SGVector<uint16_t>; 01458 template class SGVector<int32_t>; 01459 template class SGVector<uint32_t>; 01460 template class SGVector<int64_t>; 01461 template class SGVector<uint64_t>; 01462 template class SGVector<float32_t>; 01463 template class SGVector<float64_t>; 01464 template class SGVector<floatmax_t>; 01465 template class SGVector<complex128_t>; 01466 } 01467 01468 #undef COMPLEX128_ERROR_NOARG 01469 #undef COMPLEX128_ERROR_ONEARG 01470 #undef COMPLEX128_ERROR_TWOARGS 01471 #undef COMPLEX128_ERROR_THREEARGS