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) 2011-2013 Heiko Strathmann 00008 * Written (W) 2012 Fernando Jose Iglesias Garcia 00009 * Written (W) 2010,2012 Soeren Sonnenburg 00010 * Copyright (C) 2010 Berlin Institute of Technology 00011 * Copyright (C) 2012 Soeren Sonnenburg 00012 */ 00013 00014 #include <shogun/lib/config.h> 00015 #include <shogun/io/SGIO.h> 00016 #include <shogun/io/File.h> 00017 #include <shogun/lib/SGMatrix.h> 00018 #include <shogun/lib/SGVector.h> 00019 #include <shogun/mathematics/Math.h> 00020 #include <shogun/mathematics/lapack.h> 00021 #include <shogun/lib/SGMatrixList.h> 00022 00023 namespace shogun { 00024 00025 template <class T> 00026 SGMatrix<T>::SGMatrix() : SGReferencedData() 00027 { 00028 init_data(); 00029 } 00030 00031 template <class T> 00032 SGMatrix<T>::SGMatrix(bool ref_counting) : SGReferencedData(ref_counting) 00033 { 00034 init_data(); 00035 } 00036 00037 template <class T> 00038 SGMatrix<T>::SGMatrix(T* m, index_t nrows, index_t ncols, bool ref_counting) 00039 : SGReferencedData(ref_counting), matrix(m), 00040 num_rows(nrows), num_cols(ncols) { } 00041 00042 template <class T> 00043 SGMatrix<T>::SGMatrix(index_t nrows, index_t ncols, bool ref_counting) 00044 : SGReferencedData(ref_counting), num_rows(nrows), num_cols(ncols) 00045 { 00046 matrix=SG_MALLOC(T, ((int64_t) nrows)*ncols); 00047 } 00048 00049 template <class T> 00050 SGMatrix<T>::SGMatrix(const SGMatrix &orig) : SGReferencedData(orig) 00051 { 00052 copy_data(orig); 00053 } 00054 00055 template <class T> 00056 SGMatrix<T>::~SGMatrix() 00057 { 00058 unref(); 00059 } 00060 00061 template <class T> 00062 bool SGMatrix<T>::operator==(SGMatrix<T>& other) 00063 { 00064 if (num_rows!=other.num_rows || num_cols!=other.num_cols) 00065 return false; 00066 00067 if (matrix!=other.matrix) 00068 return false; 00069 00070 return true; 00071 } 00072 00073 template <class T> 00074 bool SGMatrix<T>::equals(SGMatrix<T>& other) 00075 { 00076 if (num_rows!=other.num_rows || num_cols!=other.num_cols) 00077 return false; 00078 00079 for (int64_t i=0; i<int64_t(num_rows)*num_cols; ++i) 00080 { 00081 if (matrix[i]!=other.matrix[i]) 00082 return false; 00083 } 00084 00085 return true; 00086 } 00087 00088 template <class T> 00089 void SGMatrix<T>::set_const(T const_elem) 00090 { 00091 for (int64_t i=0; i<int64_t(num_rows)*num_cols; i++) 00092 matrix[i]=const_elem ; 00093 } 00094 00095 template <class T> 00096 void SGMatrix<T>::zero() 00097 { 00098 if (matrix && (int64_t(num_rows)*num_cols)) 00099 set_const(0); 00100 } 00101 00102 template <> 00103 void SGMatrix<complex128_t>::zero() 00104 { 00105 if (matrix && (int64_t(num_rows)*num_cols)) 00106 set_const(complex128_t(0.0)); 00107 } 00108 00109 template <class T> 00110 T SGMatrix<T>::max_single() 00111 { 00112 T max=matrix[0]; 00113 for (int64_t i=1; i<int64_t(num_rows)*num_cols; ++i) 00114 { 00115 if (matrix[i]>max) 00116 max=matrix[i]; 00117 } 00118 00119 return max; 00120 } 00121 00122 template <> 00123 complex128_t SGMatrix<complex128_t>::max_single() 00124 { 00125 SG_SERROR("SGMatrix::max_single():: Not supported for complex128_t\n"); 00126 return complex128_t(0.0); 00127 } 00128 00129 template <class T> 00130 SGMatrix<T> SGMatrix<T>::clone() 00131 { 00132 return SGMatrix<T>(clone_matrix(matrix, num_rows, num_cols), 00133 num_rows, num_cols); 00134 } 00135 00136 template <class T> 00137 T* SGMatrix<T>::clone_matrix(const T* matrix, int32_t nrows, int32_t ncols) 00138 { 00139 T* result = SG_MALLOC(T, int64_t(nrows)*ncols); 00140 for (int64_t i=0; i<int64_t(nrows)*ncols; i++) 00141 result[i]=matrix[i]; 00142 00143 return result; 00144 } 00145 00146 template <class T> 00147 void SGMatrix<T>::transpose_matrix( 00148 T*& matrix, int32_t& num_feat, int32_t& num_vec) 00149 { 00150 /* this should be done in-place! Heiko */ 00151 T* transposed=SG_MALLOC(T, int64_t(num_vec)*num_feat); 00152 for (int64_t i=0; i<num_vec; i++) 00153 { 00154 for (int64_t j=0; j<num_feat; j++) 00155 transposed[i+j*num_vec]=matrix[i*num_feat+j]; 00156 } 00157 00158 SG_FREE(matrix); 00159 matrix=transposed; 00160 00161 CMath::swap(num_feat, num_vec); 00162 } 00163 00164 template <class T> 00165 void SGMatrix<T>::create_diagonal_matrix(T* matrix, T* v,int32_t size) 00166 { 00167 for(int64_t i=0;i<size;i++) 00168 { 00169 for(int64_t j=0;j<size;j++) 00170 { 00171 if(i==j) 00172 matrix[j*size+i]=v[i]; 00173 else 00174 matrix[j*size+i]=0; 00175 } 00176 } 00177 } 00178 00179 template <class T> 00180 float64_t SGMatrix<T>::trace( 00181 float64_t* mat, int32_t cols, int32_t rows) 00182 { 00183 float64_t trace=0; 00184 for (int64_t i=0; i<rows; i++) 00185 trace+=mat[i*cols+i]; 00186 return trace; 00187 } 00188 00189 template <class T> 00190 T* SGMatrix<T>::get_row_sum(T* matrix, int32_t m, int32_t n) 00191 { 00192 T* rowsums=SG_CALLOC(T, n); 00193 00194 for (int64_t i=0; i<n; i++) 00195 { 00196 for (int64_t j=0; j<m; j++) 00197 rowsums[i]+=matrix[j+i*m]; 00198 } 00199 return rowsums; 00200 } 00201 00202 template <class T> 00203 T* SGMatrix<T>::get_column_sum(T* matrix, int32_t m, int32_t n) 00204 { 00205 T* colsums=SG_CALLOC(T, m); 00206 00207 for (int64_t i=0; i<n; i++) 00208 { 00209 for (int64_t j=0; j<m; j++) 00210 colsums[j]+=matrix[j+i*m]; 00211 } 00212 return colsums; 00213 } 00214 00215 template <class T> 00216 void SGMatrix<T>::center() 00217 { 00218 center_matrix(matrix, num_rows, num_cols); 00219 } 00220 00221 template <class T> 00222 void SGMatrix<T>::center_matrix(T* matrix, int32_t m, int32_t n) 00223 { 00224 float64_t num_data=n; 00225 00226 T* colsums=get_column_sum(matrix, m,n); 00227 T* rowsums=get_row_sum(matrix, m,n); 00228 00229 for (int32_t i=0; i<m; i++) 00230 colsums[i]/=num_data; 00231 for (int32_t j=0; j<n; j++) 00232 rowsums[j]/=num_data; 00233 00234 T s=SGVector<T>::sum(rowsums, n)/num_data; 00235 00236 for (int64_t i=0; i<n; i++) 00237 { 00238 for (int64_t j=0; j<m; j++) 00239 matrix[i*m+j]+=s-colsums[j]-rowsums[i]; 00240 } 00241 00242 SG_FREE(rowsums); 00243 SG_FREE(colsums); 00244 } 00245 00246 template <class T> 00247 void SGMatrix<T>::remove_column_mean() 00248 { 00249 /* compute "row" sums (which I would call column sums), i.e. sum of all 00250 * elements in a fixed column */ 00251 T* means=get_row_sum(matrix, num_rows, num_cols); 00252 00253 /* substract column mean from every corresponding entry */ 00254 for (int64_t i=0; i<num_cols; ++i) 00255 { 00256 means[i]/=num_rows; 00257 for (int64_t j=0; j<num_rows; ++j) 00258 matrix[i*num_rows+j]-=means[i]; 00259 } 00260 00261 SG_FREE(means); 00262 } 00263 00264 template<class T> void SGMatrix<T>::display_matrix(const char* name) const 00265 { 00266 display_matrix(matrix, num_rows, num_cols, name); 00267 } 00268 00269 template <class T> 00270 void SGMatrix<T>::display_matrix( 00271 const SGMatrix<T> matrix, const char* name, 00272 const char* prefix) 00273 { 00274 matrix.display_matrix(); 00275 } 00276 00277 template <> 00278 void SGMatrix<bool>::display_matrix( 00279 const bool* matrix, int32_t rows, int32_t cols, const char* name, 00280 const char* prefix) 00281 { 00282 ASSERT(rows>=0 && cols>=0) 00283 SG_SPRINT("%s%s=[\n", prefix, name) 00284 for (int64_t i=0; i<rows; i++) 00285 { 00286 SG_SPRINT("%s[", prefix) 00287 for (int64_t j=0; j<cols; j++) 00288 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i] ? 1 : 0, 00289 j==cols-1? "" : ","); 00290 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00291 } 00292 SG_SPRINT("%s]\n", prefix) 00293 } 00294 00295 template <> 00296 void SGMatrix<char>::display_matrix( 00297 const char* matrix, int32_t rows, int32_t cols, const char* name, 00298 const char* prefix) 00299 { 00300 ASSERT(rows>=0 && cols>=0) 00301 SG_SPRINT("%s%s=[\n", prefix, name) 00302 for (int64_t i=0; i<rows; i++) 00303 { 00304 SG_SPRINT("%s[", prefix) 00305 for (int64_t j=0; j<cols; j++) 00306 SG_SPRINT("%s\t%c%s", prefix, matrix[j*rows+i], 00307 j==cols-1? "" : ","); 00308 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00309 } 00310 SG_SPRINT("%s]\n", prefix) 00311 } 00312 00313 template <> 00314 void SGMatrix<int8_t>::display_matrix( 00315 const int8_t* matrix, int32_t rows, int32_t cols, const char* name, 00316 const char* prefix) 00317 { 00318 ASSERT(rows>=0 && cols>=0) 00319 SG_SPRINT("%s%s=[\n", prefix, name) 00320 for (int64_t i=0; i<rows; i++) 00321 { 00322 SG_SPRINT("%s[", prefix) 00323 for (int64_t j=0; j<cols; j++) 00324 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00325 j==cols-1? "" : ","); 00326 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00327 } 00328 SG_SPRINT("%s]\n", prefix) 00329 } 00330 00331 template <> 00332 void SGMatrix<uint8_t>::display_matrix( 00333 const uint8_t* matrix, int32_t rows, int32_t cols, const char* name, 00334 const char* prefix) 00335 { 00336 ASSERT(rows>=0 && cols>=0) 00337 SG_SPRINT("%s%s=[\n", prefix, name) 00338 for (int64_t i=0; i<rows; i++) 00339 { 00340 SG_SPRINT("%s[", prefix) 00341 for (int64_t j=0; j<cols; j++) 00342 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00343 j==cols-1? "" : ","); 00344 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00345 } 00346 SG_SPRINT("%s]\n", prefix) 00347 } 00348 00349 template <> 00350 void SGMatrix<int16_t>::display_matrix( 00351 const int16_t* matrix, int32_t rows, int32_t cols, const char* name, 00352 const char* prefix) 00353 { 00354 ASSERT(rows>=0 && cols>=0) 00355 SG_SPRINT("%s%s=[\n", prefix, name) 00356 for (int64_t i=0; i<rows; i++) 00357 { 00358 SG_SPRINT("%s[", prefix) 00359 for (int64_t j=0; j<cols; j++) 00360 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00361 j==cols-1? "" : ","); 00362 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00363 } 00364 SG_SPRINT("%s]\n", prefix) 00365 } 00366 00367 template <> 00368 void SGMatrix<uint16_t>::display_matrix( 00369 const uint16_t* matrix, int32_t rows, int32_t cols, const char* name, 00370 const char* prefix) 00371 { 00372 ASSERT(rows>=0 && cols>=0) 00373 SG_SPRINT("%s%s=[\n", prefix, name) 00374 for (int64_t i=0; i<rows; i++) 00375 { 00376 SG_SPRINT("%s[", prefix) 00377 for (int64_t j=0; j<cols; j++) 00378 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00379 j==cols-1? "" : ","); 00380 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00381 } 00382 SG_SPRINT("%s]\n", prefix) 00383 } 00384 00385 00386 template <> 00387 void SGMatrix<int32_t>::display_matrix( 00388 const int32_t* matrix, int32_t rows, int32_t cols, const char* name, 00389 const char* prefix) 00390 { 00391 ASSERT(rows>=0 && cols>=0) 00392 SG_SPRINT("%s%s=[\n", prefix, name) 00393 for (int64_t i=0; i<rows; i++) 00394 { 00395 SG_SPRINT("%s[", prefix) 00396 for (int64_t j=0; j<cols; j++) 00397 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00398 j==cols-1? "" : ","); 00399 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00400 } 00401 SG_SPRINT("%s]\n", prefix) 00402 } 00403 00404 template <> 00405 void SGMatrix<uint32_t>::display_matrix( 00406 const uint32_t* matrix, int32_t rows, int32_t cols, const char* name, 00407 const char* prefix) 00408 { 00409 ASSERT(rows>=0 && cols>=0) 00410 SG_SPRINT("%s%s=[\n", prefix, name) 00411 for (int64_t i=0; i<rows; i++) 00412 { 00413 SG_SPRINT("%s[", prefix) 00414 for (int64_t j=0; j<cols; j++) 00415 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00416 j==cols-1? "" : ","); 00417 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00418 } 00419 SG_SPRINT("%s]\n", prefix) 00420 } 00421 template <> 00422 void SGMatrix<int64_t>::display_matrix( 00423 const int64_t* matrix, int32_t rows, int32_t cols, const char* name, 00424 const char* prefix) 00425 { 00426 ASSERT(rows>=0 && cols>=0) 00427 SG_SPRINT("%s%s=[\n", prefix, name) 00428 for (int64_t i=0; i<rows; i++) 00429 { 00430 SG_SPRINT("%s[", prefix) 00431 for (int64_t j=0; j<cols; j++) 00432 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00433 j==cols-1? "" : ","); 00434 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00435 } 00436 SG_SPRINT("%s]\n", prefix) 00437 } 00438 00439 template <> 00440 void SGMatrix<uint64_t>::display_matrix( 00441 const uint64_t* matrix, int32_t rows, int32_t cols, const char* name, 00442 const char* prefix) 00443 { 00444 ASSERT(rows>=0 && cols>=0) 00445 SG_SPRINT("%s%s=[\n", prefix, name) 00446 for (int64_t i=0; i<rows; i++) 00447 { 00448 SG_SPRINT("%s[", prefix) 00449 for (int64_t j=0; j<cols; j++) 00450 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00451 j==cols-1? "" : ","); 00452 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00453 } 00454 SG_SPRINT("%s]\n", prefix) 00455 } 00456 00457 template <> 00458 void SGMatrix<float32_t>::display_matrix( 00459 const float32_t* matrix, int32_t rows, int32_t cols, const char* name, 00460 const char* prefix) 00461 { 00462 ASSERT(rows>=0 && cols>=0) 00463 SG_SPRINT("%s%s=[\n", prefix, name) 00464 for (int64_t i=0; i<rows; i++) 00465 { 00466 SG_SPRINT("%s[", prefix) 00467 for (int64_t j=0; j<cols; j++) 00468 SG_SPRINT("%s\t%.18g%s", prefix, (float) matrix[j*rows+i], 00469 j==cols-1? "" : ","); 00470 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00471 } 00472 SG_SPRINT("%s]\n", prefix) 00473 } 00474 00475 template <> 00476 void SGMatrix<float64_t>::display_matrix( 00477 const float64_t* matrix, int32_t rows, int32_t cols, const char* name, 00478 const char* prefix) 00479 { 00480 ASSERT(rows>=0 && cols>=0) 00481 SG_SPRINT("%s%s=[\n", prefix, name) 00482 for (int64_t i=0; i<rows; i++) 00483 { 00484 SG_SPRINT("%s[", prefix) 00485 for (int64_t j=0; j<cols; j++) 00486 SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i], 00487 j==cols-1? "" : ","); 00488 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00489 } 00490 SG_SPRINT("%s]\n", prefix) 00491 } 00492 00493 template <> 00494 void SGMatrix<floatmax_t>::display_matrix( 00495 const floatmax_t* matrix, int32_t rows, int32_t cols, const char* name, 00496 const char* prefix) 00497 { 00498 ASSERT(rows>=0 && cols>=0) 00499 SG_SPRINT("%s%s=[\n", prefix, name) 00500 for (int64_t i=0; i<rows; i++) 00501 { 00502 SG_SPRINT("%s[", prefix) 00503 for (int64_t j=0; j<cols; j++) 00504 SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i], 00505 j==cols-1? "" : ","); 00506 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00507 } 00508 SG_SPRINT("%s]\n", prefix) 00509 } 00510 00511 template <> 00512 void SGMatrix<complex128_t>::display_matrix( 00513 const complex128_t* matrix, int32_t rows, int32_t cols, const char* name, 00514 const char* prefix) 00515 { 00516 ASSERT(rows>=0 && cols>=0) 00517 SG_SPRINT("%s%s=[\n", prefix, name) 00518 for (int64_t i=0; i<rows; i++) 00519 { 00520 SG_SPRINT("%s[", prefix) 00521 for (int64_t j=0; j<cols; j++) 00522 SG_SPRINT("%s\t(%.18g+i%.18g)%s", prefix, matrix[j*rows+i].real(), 00523 matrix[j*rows+i].imag(), j==cols-1? "" : ","); 00524 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",") 00525 } 00526 SG_SPRINT("%s]\n", prefix) 00527 } 00528 00529 template <> 00530 SGMatrix<char> SGMatrix<char>::create_identity_matrix(index_t size, char scale) 00531 { 00532 SG_SNOTIMPLEMENTED 00533 return SGMatrix<char>(); 00534 } 00535 00536 template <> 00537 SGMatrix<int8_t> SGMatrix<int8_t>::create_identity_matrix(index_t size, int8_t scale) 00538 { 00539 SGMatrix<int8_t> I(size, size); 00540 for (index_t i=0; i<size; ++i) 00541 { 00542 for (index_t j=0; j<size; ++j) 00543 I(i,j)=i==j ? scale : 0.0; 00544 } 00545 00546 return I; 00547 } 00548 00549 template <> 00550 SGMatrix<uint8_t> SGMatrix<uint8_t>::create_identity_matrix(index_t size, uint8_t scale) 00551 { 00552 SGMatrix<uint8_t> I(size, size); 00553 for (index_t i=0; i<size; ++i) 00554 { 00555 for (index_t j=0; j<size; ++j) 00556 I(i,j)=i==j ? scale : 0.0; 00557 } 00558 00559 return I; 00560 } 00561 00562 template <> 00563 SGMatrix<bool> SGMatrix<bool>::create_identity_matrix(index_t size, bool scale) 00564 { 00565 SGMatrix<bool> I(size, size); 00566 for (index_t i=0; i<size; ++i) 00567 { 00568 for (index_t j=0; j<size; ++j) 00569 I(i,j)=i==j ? scale : (!scale); 00570 } 00571 00572 return I; 00573 } 00574 00575 template <> 00576 SGMatrix<int16_t> SGMatrix<int16_t>::create_identity_matrix(index_t size, int16_t scale) 00577 { 00578 SGMatrix<int16_t> I(size, size); 00579 for (index_t i=0; i<size; ++i) 00580 { 00581 for (index_t j=0; j<size; ++j) 00582 I(i,j)=i==j ? scale : 0.0; 00583 } 00584 00585 return I; 00586 } 00587 00588 template <> 00589 SGMatrix<uint16_t> SGMatrix<uint16_t>::create_identity_matrix(index_t size, uint16_t scale) 00590 { 00591 SGMatrix<uint16_t> I(size, size); 00592 for (index_t i=0; i<size; ++i) 00593 { 00594 for (index_t j=0; j<size; ++j) 00595 I(i,j)=i==j ? scale : 0.0; 00596 } 00597 00598 return I; 00599 } 00600 00601 template <> 00602 SGMatrix<int32_t> SGMatrix<int32_t>::create_identity_matrix(index_t size, int32_t scale) 00603 { 00604 SGMatrix<int32_t> I(size, size); 00605 for (index_t i=0; i<size; ++i) 00606 { 00607 for (index_t j=0; j<size; ++j) 00608 I(i,j)=i==j ? scale : 0.0; 00609 } 00610 00611 return I; 00612 } 00613 00614 template <> 00615 SGMatrix<uint32_t> SGMatrix<uint32_t>::create_identity_matrix(index_t size, uint32_t scale) 00616 { 00617 SGMatrix<uint32_t> I(size, size); 00618 for (index_t i=0; i<size; ++i) 00619 { 00620 for (index_t j=0; j<size; ++j) 00621 I(i,j)=i==j ? scale : 0.0; 00622 } 00623 00624 return I; 00625 } 00626 00627 template <> 00628 SGMatrix<int64_t> SGMatrix<int64_t>::create_identity_matrix(index_t size, int64_t scale) 00629 { 00630 SGMatrix<int64_t> I(size, size); 00631 for (index_t i=0; i<size; ++i) 00632 { 00633 for (index_t j=0; j<size; ++j) 00634 I(i,j)=i==j ? scale : 0.0; 00635 } 00636 00637 return I; 00638 } 00639 00640 template <> 00641 SGMatrix<uint64_t> SGMatrix<uint64_t>::create_identity_matrix(index_t size, uint64_t scale) 00642 { 00643 SGMatrix<uint64_t> I(size, size); 00644 for (index_t i=0; i<size; ++i) 00645 { 00646 for (index_t j=0; j<size; ++j) 00647 I(i,j)=i==j ? scale : 0.0; 00648 } 00649 00650 return I; 00651 } 00652 00653 template <> 00654 SGMatrix<float32_t> SGMatrix<float32_t>::create_identity_matrix(index_t size, float32_t scale) 00655 { 00656 SGMatrix<float32_t> I(size, size); 00657 for (index_t i=0; i<size; ++i) 00658 { 00659 for (index_t j=0; j<size; ++j) 00660 I(i,j)=i==j ? scale : 0.0; 00661 } 00662 00663 return I; 00664 } 00665 00666 template <> 00667 SGMatrix<float64_t> SGMatrix<float64_t>::create_identity_matrix(index_t size, float64_t scale) 00668 { 00669 SGMatrix<float64_t> I(size, size); 00670 for (index_t i=0; i<size; ++i) 00671 { 00672 for (index_t j=0; j<size; ++j) 00673 I(i,j)=i==j ? scale : 0.0; 00674 } 00675 00676 return I; 00677 } 00678 00679 template <> 00680 SGMatrix<floatmax_t> SGMatrix<floatmax_t>::create_identity_matrix(index_t size, floatmax_t scale) 00681 { 00682 SGMatrix<floatmax_t> I(size, size); 00683 for (index_t i=0; i<size; ++i) 00684 { 00685 for (index_t j=0; j<size; ++j) 00686 I(i,j)=i==j ? scale : 0.0; 00687 } 00688 00689 return I; 00690 } 00691 00692 template <> 00693 SGMatrix<complex128_t> SGMatrix<complex128_t>::create_identity_matrix(index_t size, complex128_t scale) 00694 { 00695 SGMatrix<complex128_t> I(size, size); 00696 for (index_t i=0; i<size; ++i) 00697 { 00698 for (index_t j=0; j<size; ++j) 00699 I(i,j)=i==j ? scale : complex128_t(0.0); 00700 } 00701 00702 return I; 00703 } 00704 00705 00706 template <class T> 00707 SGMatrix<float64_t> SGMatrix<T>::create_centering_matrix(index_t size) 00708 { 00709 SGMatrix<float64_t> H=SGMatrix<float64_t>::create_identity_matrix(size, 1.0); 00710 00711 float64_t subtract=1.0/size; 00712 for (index_t i=0; i<size; ++i) 00713 { 00714 for (index_t j=0; j<0; ++j) 00715 H(i,j)-=subtract; 00716 } 00717 00718 return H; 00719 } 00720 00721 //Howto construct the pseudo inverse (from "The Matrix Cookbook") 00722 // 00723 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m). 00724 // 00725 //The matrix A+ known as the pseudo inverse is unique and does always exist. 00726 // 00727 //The pseudo inverse A+ can be constructed from the singular value 00728 //decomposition A = UDV^T , by A^+ = V(D+)U^T. 00729 00730 #ifdef HAVE_LAPACK 00731 template <class T> 00732 float64_t* SGMatrix<T>::pinv( 00733 float64_t* matrix, int32_t rows, int32_t cols, float64_t* target) 00734 { 00735 if (!target) 00736 target=SG_MALLOC(float64_t, rows*cols); 00737 00738 char jobu='A'; 00739 char jobvt='A'; 00740 int m=rows; /* for calling external lib */ 00741 int n=cols; /* for calling external lib */ 00742 int lda=m; /* for calling external lib */ 00743 int ldu=m; /* for calling external lib */ 00744 int ldvt=n; /* for calling external lib */ 00745 int info=-1; /* for calling external lib */ 00746 int32_t lsize=CMath::min((int32_t) m, (int32_t) n); 00747 double* s=SG_MALLOC(double, lsize); 00748 double* u=SG_MALLOC(double, m*m); 00749 double* vt=SG_MALLOC(double, n*n); 00750 00751 wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info); 00752 ASSERT(info==0) 00753 00754 for (int64_t i=0; i<n; i++) 00755 { 00756 for (int64_t j=0; j<lsize; j++) 00757 vt[i*n+j]=vt[i*n+j]/s[j]; 00758 } 00759 00760 cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m); 00761 00762 SG_FREE(u); 00763 SG_FREE(vt); 00764 SG_FREE(s); 00765 00766 return target; 00767 } 00768 00770 template <class T> 00771 void SGMatrix<T>::inverse(SGMatrix<float64_t> matrix) 00772 { 00773 ASSERT(matrix.num_cols==matrix.num_rows) 00774 int32_t* ipiv = SG_MALLOC(int32_t, matrix.num_cols); 00775 clapack_dgetrf(CblasColMajor,matrix.num_cols,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv); 00776 clapack_dgetri(CblasColMajor,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv); 00777 SG_FREE(ipiv); 00778 } 00779 00780 template <class T> 00781 SGVector<float64_t> SGMatrix<T>::compute_eigenvectors(SGMatrix<float64_t> matrix) 00782 { 00783 if (matrix.num_rows!=matrix.num_rows) 00784 { 00785 SG_SERROR("SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix" 00786 " rows and columns are not equal!\n"); 00787 } 00788 00789 /* use reference counting for SGVector */ 00790 SGVector<float64_t> result(NULL, 0, true); 00791 result.vlen=matrix.num_rows; 00792 result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows, 00793 matrix.num_rows); 00794 return result; 00795 } 00796 00797 template <class T> 00798 double* SGMatrix<T>::compute_eigenvectors(double* matrix, int n, int m) 00799 { 00800 ASSERT(n == m) 00801 00802 char V='V'; 00803 char U='U'; 00804 int info; 00805 int ord=n; 00806 int lda=n; 00807 double* eigenvalues=SG_CALLOC(float64_t, n+1); 00808 00809 // lapack sym matrix eigenvalues+vectors 00810 wrap_dsyev(V, U, ord, matrix, lda, 00811 eigenvalues, &info); 00812 00813 if (info!=0) 00814 SG_SERROR("DSYEV failed with code %d\n", info) 00815 00816 return eigenvalues; 00817 } 00818 00819 template <class T> 00820 void SGMatrix<T>::compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors, 00821 int n, int il, int iu) 00822 { 00823 eigenvalues = SG_MALLOC(double, n); 00824 eigenvectors = SG_MALLOC(double, (iu-il+1)*n); 00825 int status = 0; 00826 wrap_dsyevr('V','U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status); 00827 ASSERT(status==0) 00828 } 00829 00830 #endif //HAVE_LAPACK 00831 00832 template <class T> 00833 SGMatrix<float64_t> SGMatrix<T>::matrix_multiply( 00834 SGMatrix<float64_t> A, SGMatrix<float64_t> B, 00835 bool transpose_A, bool transpose_B, float64_t scale) 00836 { 00837 /* these variables store size of transposed matrices*/ 00838 index_t cols_A=transpose_A ? A.num_rows : A.num_cols; 00839 index_t rows_A=transpose_A ? A.num_cols : A.num_rows; 00840 index_t rows_B=transpose_B ? B.num_cols : B.num_rows; 00841 index_t cols_B=transpose_B ? B.num_rows : B.num_cols; 00842 00843 /* do a dimension check */ 00844 if (cols_A!=rows_B) 00845 { 00846 SG_SERROR("SGMatrix::matrix_multiply(): Dimension mismatch: " 00847 "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B); 00848 } 00849 00850 /* allocate result matrix */ 00851 SGMatrix<float64_t> C(rows_A, cols_B); 00852 C.zero(); 00853 #ifdef HAVE_LAPACK 00854 /* multiply */ 00855 cblas_dgemm(CblasColMajor, 00856 transpose_A ? CblasTrans : CblasNoTrans, 00857 transpose_B ? CblasTrans : CblasNoTrans, 00858 rows_A, cols_B, cols_A, scale, 00859 A.matrix, A.num_rows, B.matrix, B.num_rows, 00860 0.0, C.matrix, C.num_rows); 00861 #else 00862 for (int32_t i=0; i<rows_A; i++) 00863 { 00864 for (int32_t j=0; j<cols_B; j++) 00865 { 00866 for (int32_t k=0; k<cols_A; k++) 00867 C(i,j) += A(i,k)*B(k,j); 00868 } 00869 } 00870 #endif //HAVE_LAPACK 00871 00872 return C; 00873 } 00874 00875 template<class T> 00876 SGMatrix<T> SGMatrix<T>::get_allocated_matrix(index_t num_rows, 00877 index_t num_cols, SGMatrix<T> pre_allocated) 00878 { 00879 SGMatrix<T> result; 00880 00881 /* evtl use pre-allocated space */ 00882 if (pre_allocated.matrix) 00883 { 00884 result=pre_allocated; 00885 00886 /* check dimension consistency */ 00887 if (pre_allocated.num_rows!=num_rows || 00888 pre_allocated.num_cols!=num_cols) 00889 { 00890 SG_SERROR("SGMatrix<T>::get_allocated_matrix(). Provided target" 00891 "matrix dimensions (%dx%d) do not match passed data " 00892 "dimensions (%dx%d)!\n", pre_allocated.num_rows, 00893 pre_allocated.num_cols, num_rows, num_cols); 00894 } 00895 } 00896 else 00897 { 00898 /* otherwise, allocate space */ 00899 result=SGMatrix<T>(num_rows, num_cols); 00900 } 00901 00902 return result; 00903 } 00904 00905 template<class T> 00906 void SGMatrix<T>::copy_data(const SGReferencedData &orig) 00907 { 00908 matrix=((SGMatrix*)(&orig))->matrix; 00909 num_rows=((SGMatrix*)(&orig))->num_rows; 00910 num_cols=((SGMatrix*)(&orig))->num_cols; 00911 } 00912 00913 template<class T> 00914 void SGMatrix<T>::init_data() 00915 { 00916 matrix=NULL; 00917 num_rows=0; 00918 num_cols=0; 00919 } 00920 00921 template<class T> 00922 void SGMatrix<T>::free_data() 00923 { 00924 SG_FREE(matrix); 00925 matrix=NULL; 00926 num_rows=0; 00927 num_cols=0; 00928 } 00929 00930 template<class T> 00931 void SGMatrix<T>::load(CFile* loader) 00932 { 00933 ASSERT(loader) 00934 unref(); 00935 00936 SG_SET_LOCALE_C; 00937 SGMatrix<T> mat; 00938 loader->get_matrix(mat.matrix, mat.num_rows, mat.num_cols); 00939 copy_data(mat); 00940 copy_refcount(mat); 00941 ref(); 00942 SG_RESET_LOCALE; 00943 } 00944 00945 template<> 00946 void SGMatrix<complex128_t>::load(CFile* loader) 00947 { 00948 SG_SERROR("SGMatrix::load():: Not supported for complex128_t\n"); 00949 } 00950 00951 template<class T> 00952 void SGMatrix<T>::save(CFile* writer) 00953 { 00954 ASSERT(writer) 00955 SG_SET_LOCALE_C; 00956 writer->set_matrix(matrix, num_rows, num_cols); 00957 SG_RESET_LOCALE; 00958 } 00959 00960 template<> 00961 void SGMatrix<complex128_t>::save(CFile* saver) 00962 { 00963 SG_SERROR("SGMatrix::save():: Not supported for complex128_t\n"); 00964 } 00965 00966 template<class T> 00967 SGVector<T> SGMatrix<T>::get_row_vector(index_t row) const 00968 { 00969 SGVector<T> rowv(num_cols); 00970 for (int64_t i = 0; i < num_cols; i++) 00971 { 00972 rowv[i] = matrix[i*num_rows+row]; 00973 } 00974 return rowv; 00975 } 00976 00977 template<class T> 00978 SGVector<T> SGMatrix<T>::get_diagonal_vector() const 00979 { 00980 index_t diag_vlen=CMath::min(num_cols, num_rows); 00981 SGVector<T> diag(diag_vlen); 00982 00983 for (int64_t i=0; i<diag_vlen; i++) 00984 { 00985 diag[i]=matrix[i*num_rows+i]; 00986 } 00987 00988 return diag; 00989 } 00990 00991 template class SGMatrix<bool>; 00992 template class SGMatrix<char>; 00993 template class SGMatrix<int8_t>; 00994 template class SGMatrix<uint8_t>; 00995 template class SGMatrix<int16_t>; 00996 template class SGMatrix<uint16_t>; 00997 template class SGMatrix<int32_t>; 00998 template class SGMatrix<uint32_t>; 00999 template class SGMatrix<int64_t>; 01000 template class SGMatrix<uint64_t>; 01001 template class SGMatrix<float32_t>; 01002 template class SGMatrix<float64_t>; 01003 template class SGMatrix<floatmax_t>; 01004 template class SGMatrix<complex128_t>; 01005 }