SHOGUN
v2.0.0
|
00001 #include <shogun/lib/config.h> 00002 #include <shogun/io/SGIO.h> 00003 #include <shogun/lib/SGMatrix.h> 00004 #include <shogun/lib/SGVector.h> 00005 #include <shogun/mathematics/Math.h> 00006 #include <shogun/mathematics/lapack.h> 00007 #include <shogun/lib/SGMatrixList.cpp> 00008 00009 namespace shogun { 00010 00011 template <class T> 00012 void SGMatrix<T>::transpose_matrix( 00013 T*& matrix, int32_t& num_feat, int32_t& num_vec) 00014 { 00015 /* this should be done in-place! Heiko */ 00016 T* transposed=SG_MALLOC(T, num_vec*num_feat); 00017 for (int32_t i=0; i<num_vec; i++) 00018 { 00019 for (int32_t j=0; j<num_feat; j++) 00020 transposed[i+j*num_vec]=matrix[i*num_feat+j]; 00021 } 00022 00023 SG_FREE(matrix); 00024 matrix=transposed; 00025 00026 CMath::swap(num_feat, num_vec); 00027 } 00028 00029 template <class T> 00030 void SGMatrix<T>::center_matrix(T* matrix, int32_t m, int32_t n) 00031 { 00032 float64_t num_data=n; 00033 00034 T* colsums=get_column_sum(matrix, m,n); 00035 T* rowsums=get_row_sum(matrix, m,n); 00036 00037 for (int32_t i=0; i<m; i++) 00038 colsums[i]/=num_data; 00039 for (int32_t j=0; j<n; j++) 00040 rowsums[j]/=num_data; 00041 00042 T s=SGVector<T>::sum(rowsums, n)/num_data; 00043 00044 for (int32_t i=0; i<n; i++) 00045 { 00046 for (int32_t j=0; j<m; j++) 00047 matrix[int64_t(i)*m+j]+=s-colsums[j]-rowsums[i]; 00048 } 00049 00050 SG_FREE(rowsums); 00051 SG_FREE(colsums); 00052 } 00053 00054 template <class T> 00055 void SGMatrix<T>::remove_column_mean() 00056 { 00057 /* compute "row" sums (which I would call column sums), i.e. sum of all 00058 * elements in a fixed column */ 00059 T* means=get_row_sum(matrix, num_rows, num_cols); 00060 00061 /* substract column mean from every corresponding entry */ 00062 for (index_t i=0; i<num_cols; ++i) 00063 { 00064 means[i]/=num_rows; 00065 for (index_t j=0; j<num_rows; ++j) 00066 matrix[i*num_rows+j]-=means[i]; 00067 } 00068 00069 SG_FREE(means); 00070 } 00071 00072 template<class T> void SGMatrix<T>::display_matrix(const char* name) const 00073 { 00074 display_matrix(matrix, num_rows, num_cols, name); 00075 } 00076 00077 template <class T> 00078 void SGMatrix<T>::display_matrix( 00079 const SGMatrix<T> matrix, const char* name, 00080 const char* prefix) 00081 { 00082 matrix.display_matrix(); 00083 } 00084 00085 template <> 00086 void SGMatrix<bool>::display_matrix( 00087 const bool* matrix, int32_t rows, int32_t cols, const char* name, 00088 const char* prefix) 00089 { 00090 ASSERT(rows>=0 && cols>=0); 00091 SG_SPRINT("%s%s=[\n", prefix, name); 00092 for (int32_t i=0; i<rows; i++) 00093 { 00094 SG_SPRINT("%s[", prefix); 00095 for (int32_t j=0; j<cols; j++) 00096 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i] ? 1 : 0, 00097 j==cols-1? "" : ","); 00098 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00099 } 00100 SG_SPRINT("%s]\n", prefix); 00101 } 00102 00103 template <> 00104 void SGMatrix<char>::display_matrix( 00105 const char* matrix, int32_t rows, int32_t cols, const char* name, 00106 const char* prefix) 00107 { 00108 ASSERT(rows>=0 && cols>=0); 00109 SG_SPRINT("%s%s=[\n", prefix, name); 00110 for (int32_t i=0; i<rows; i++) 00111 { 00112 SG_SPRINT("%s[", prefix); 00113 for (int32_t j=0; j<cols; j++) 00114 SG_SPRINT("%s\t%c%s", prefix, matrix[j*rows+i], 00115 j==cols-1? "" : ","); 00116 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00117 } 00118 SG_SPRINT("%s]\n", prefix); 00119 } 00120 00121 template <> 00122 void SGMatrix<int8_t>::display_matrix( 00123 const int8_t* matrix, int32_t rows, int32_t cols, const char* name, 00124 const char* prefix) 00125 { 00126 ASSERT(rows>=0 && cols>=0); 00127 SG_SPRINT("%s%s=[\n", prefix, name); 00128 for (int32_t i=0; i<rows; i++) 00129 { 00130 SG_SPRINT("%s[", prefix); 00131 for (int32_t j=0; j<cols; j++) 00132 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00133 j==cols-1? "" : ","); 00134 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00135 } 00136 SG_SPRINT("%s]\n", prefix); 00137 } 00138 00139 template <> 00140 void SGMatrix<uint8_t>::display_matrix( 00141 const uint8_t* matrix, int32_t rows, int32_t cols, const char* name, 00142 const char* prefix) 00143 { 00144 ASSERT(rows>=0 && cols>=0); 00145 SG_SPRINT("%s%s=[\n", prefix, name); 00146 for (int32_t i=0; i<rows; i++) 00147 { 00148 SG_SPRINT("%s[", prefix); 00149 for (int32_t j=0; j<cols; j++) 00150 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00151 j==cols-1? "" : ","); 00152 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00153 } 00154 SG_SPRINT("%s]\n", prefix); 00155 } 00156 00157 template <> 00158 void SGMatrix<int16_t>::display_matrix( 00159 const int16_t* matrix, int32_t rows, int32_t cols, const char* name, 00160 const char* prefix) 00161 { 00162 ASSERT(rows>=0 && cols>=0); 00163 SG_SPRINT("%s%s=[\n", prefix, name); 00164 for (int32_t i=0; i<rows; i++) 00165 { 00166 SG_SPRINT("%s[", prefix); 00167 for (int32_t j=0; j<cols; j++) 00168 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00169 j==cols-1? "" : ","); 00170 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00171 } 00172 SG_SPRINT("%s]\n", prefix); 00173 } 00174 00175 template <> 00176 void SGMatrix<uint16_t>::display_matrix( 00177 const uint16_t* matrix, int32_t rows, int32_t cols, const char* name, 00178 const char* prefix) 00179 { 00180 ASSERT(rows>=0 && cols>=0); 00181 SG_SPRINT("%s%s=[\n", prefix, name); 00182 for (int32_t i=0; i<rows; i++) 00183 { 00184 SG_SPRINT("%s[", prefix); 00185 for (int32_t j=0; j<cols; j++) 00186 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00187 j==cols-1? "" : ","); 00188 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00189 } 00190 SG_SPRINT("%s]\n", prefix); 00191 } 00192 00193 00194 template <> 00195 void SGMatrix<int32_t>::display_matrix( 00196 const int32_t* matrix, int32_t rows, int32_t cols, const char* name, 00197 const char* prefix) 00198 { 00199 ASSERT(rows>=0 && cols>=0); 00200 SG_SPRINT("%s%s=[\n", prefix, name); 00201 for (int32_t i=0; i<rows; i++) 00202 { 00203 SG_SPRINT("%s[", prefix); 00204 for (int32_t j=0; j<cols; j++) 00205 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00206 j==cols-1? "" : ","); 00207 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00208 } 00209 SG_SPRINT("%s]\n", prefix); 00210 } 00211 00212 template <> 00213 void SGMatrix<uint32_t>::display_matrix( 00214 const uint32_t* matrix, int32_t rows, int32_t cols, const char* name, 00215 const char* prefix) 00216 { 00217 ASSERT(rows>=0 && cols>=0); 00218 SG_SPRINT("%s%s=[\n", prefix, name); 00219 for (int32_t i=0; i<rows; i++) 00220 { 00221 SG_SPRINT("%s[", prefix); 00222 for (int32_t j=0; j<cols; j++) 00223 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00224 j==cols-1? "" : ","); 00225 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00226 } 00227 SG_SPRINT("%s]\n", prefix); 00228 } 00229 template <> 00230 void SGMatrix<int64_t>::display_matrix( 00231 const int64_t* matrix, int32_t rows, int32_t cols, const char* name, 00232 const char* prefix) 00233 { 00234 ASSERT(rows>=0 && cols>=0); 00235 SG_SPRINT("%s%s=[\n", prefix, name); 00236 for (int32_t i=0; i<rows; i++) 00237 { 00238 SG_SPRINT("%s[", prefix); 00239 for (int32_t j=0; j<cols; j++) 00240 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00241 j==cols-1? "" : ","); 00242 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00243 } 00244 SG_SPRINT("%s]\n", prefix); 00245 } 00246 00247 template <> 00248 void SGMatrix<uint64_t>::display_matrix( 00249 const uint64_t* matrix, int32_t rows, int32_t cols, const char* name, 00250 const char* prefix) 00251 { 00252 ASSERT(rows>=0 && cols>=0); 00253 SG_SPRINT("%s%s=[\n", prefix, name); 00254 for (int32_t i=0; i<rows; i++) 00255 { 00256 SG_SPRINT("%s[", prefix); 00257 for (int32_t j=0; j<cols; j++) 00258 SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i], 00259 j==cols-1? "" : ","); 00260 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00261 } 00262 SG_SPRINT("%s]\n", prefix); 00263 } 00264 00265 template <> 00266 void SGMatrix<float32_t>::display_matrix( 00267 const float32_t* matrix, int32_t rows, int32_t cols, const char* name, 00268 const char* prefix) 00269 { 00270 ASSERT(rows>=0 && cols>=0); 00271 SG_SPRINT("%s%s=[\n", prefix, name); 00272 for (int32_t i=0; i<rows; i++) 00273 { 00274 SG_SPRINT("%s[", prefix); 00275 for (int32_t j=0; j<cols; j++) 00276 SG_SPRINT("%s\t%.18g%s", prefix, (float) matrix[j*rows+i], 00277 j==cols-1? "" : ","); 00278 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00279 } 00280 SG_SPRINT("%s]\n", prefix); 00281 } 00282 00283 template <> 00284 void SGMatrix<float64_t>::display_matrix( 00285 const float64_t* matrix, int32_t rows, int32_t cols, const char* name, 00286 const char* prefix) 00287 { 00288 ASSERT(rows>=0 && cols>=0); 00289 SG_SPRINT("%s%s=[\n", prefix, name); 00290 for (int32_t i=0; i<rows; i++) 00291 { 00292 SG_SPRINT("%s[", prefix); 00293 for (int32_t j=0; j<cols; j++) 00294 SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i], 00295 j==cols-1? "" : ","); 00296 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00297 } 00298 SG_SPRINT("%s]\n", prefix); 00299 } 00300 00301 template <> 00302 void SGMatrix<floatmax_t>::display_matrix( 00303 const floatmax_t* matrix, int32_t rows, int32_t cols, const char* name, 00304 const char* prefix) 00305 { 00306 ASSERT(rows>=0 && cols>=0); 00307 SG_SPRINT("%s%s=[\n", prefix, name); 00308 for (int32_t i=0; i<rows; i++) 00309 { 00310 SG_SPRINT("%s[", prefix); 00311 for (int32_t j=0; j<cols; j++) 00312 SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i], 00313 j==cols-1? "" : ","); 00314 SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ","); 00315 } 00316 SG_SPRINT("%s]\n", prefix); 00317 } 00318 00319 template <> 00320 SGMatrix<char> SGMatrix<char>::create_identity_matrix(index_t size, char scale) 00321 { 00322 SG_SNOTIMPLEMENTED; 00323 return SGMatrix<char>(); 00324 } 00325 00326 template <> 00327 SGMatrix<int8_t> SGMatrix<int8_t>::create_identity_matrix(index_t size, int8_t scale) 00328 { 00329 SGMatrix<int8_t> I(size, size); 00330 for (index_t i=0; i<size; ++i) 00331 { 00332 for (index_t j=0; j<size; ++j) 00333 I(i,j)=i==j ? scale : 0.0; 00334 } 00335 00336 return I; 00337 } 00338 00339 template <> 00340 SGMatrix<uint8_t> SGMatrix<uint8_t>::create_identity_matrix(index_t size, uint8_t scale) 00341 { 00342 SGMatrix<uint8_t> I(size, size); 00343 for (index_t i=0; i<size; ++i) 00344 { 00345 for (index_t j=0; j<size; ++j) 00346 I(i,j)=i==j ? scale : 0.0; 00347 } 00348 00349 return I; 00350 } 00351 00352 template <> 00353 SGMatrix<bool> SGMatrix<bool>::create_identity_matrix(index_t size, bool scale) 00354 { 00355 SGMatrix<bool> I(size, size); 00356 for (index_t i=0; i<size; ++i) 00357 { 00358 for (index_t j=0; j<size; ++j) 00359 I(i,j)=i==j ? scale : (!scale); 00360 } 00361 00362 return I; 00363 } 00364 00365 template <> 00366 SGMatrix<int16_t> SGMatrix<int16_t>::create_identity_matrix(index_t size, int16_t scale) 00367 { 00368 SGMatrix<int16_t> I(size, size); 00369 for (index_t i=0; i<size; ++i) 00370 { 00371 for (index_t j=0; j<size; ++j) 00372 I(i,j)=i==j ? scale : 0.0; 00373 } 00374 00375 return I; 00376 } 00377 00378 template <> 00379 SGMatrix<uint16_t> SGMatrix<uint16_t>::create_identity_matrix(index_t size, uint16_t scale) 00380 { 00381 SGMatrix<uint16_t> I(size, size); 00382 for (index_t i=0; i<size; ++i) 00383 { 00384 for (index_t j=0; j<size; ++j) 00385 I(i,j)=i==j ? scale : 0.0; 00386 } 00387 00388 return I; 00389 } 00390 00391 template <> 00392 SGMatrix<int32_t> SGMatrix<int32_t>::create_identity_matrix(index_t size, int32_t scale) 00393 { 00394 SGMatrix<int32_t> I(size, size); 00395 for (index_t i=0; i<size; ++i) 00396 { 00397 for (index_t j=0; j<size; ++j) 00398 I(i,j)=i==j ? scale : 0.0; 00399 } 00400 00401 return I; 00402 } 00403 00404 template <> 00405 SGMatrix<uint32_t> SGMatrix<uint32_t>::create_identity_matrix(index_t size, uint32_t scale) 00406 { 00407 SGMatrix<uint32_t> I(size, size); 00408 for (index_t i=0; i<size; ++i) 00409 { 00410 for (index_t j=0; j<size; ++j) 00411 I(i,j)=i==j ? scale : 0.0; 00412 } 00413 00414 return I; 00415 } 00416 00417 template <> 00418 SGMatrix<int64_t> SGMatrix<int64_t>::create_identity_matrix(index_t size, int64_t scale) 00419 { 00420 SGMatrix<int64_t> I(size, size); 00421 for (index_t i=0; i<size; ++i) 00422 { 00423 for (index_t j=0; j<size; ++j) 00424 I(i,j)=i==j ? scale : 0.0; 00425 } 00426 00427 return I; 00428 } 00429 00430 template <> 00431 SGMatrix<uint64_t> SGMatrix<uint64_t>::create_identity_matrix(index_t size, uint64_t scale) 00432 { 00433 SGMatrix<uint64_t> I(size, size); 00434 for (index_t i=0; i<size; ++i) 00435 { 00436 for (index_t j=0; j<size; ++j) 00437 I(i,j)=i==j ? scale : 0.0; 00438 } 00439 00440 return I; 00441 } 00442 00443 template <> 00444 SGMatrix<float32_t> SGMatrix<float32_t>::create_identity_matrix(index_t size, float32_t scale) 00445 { 00446 SGMatrix<float32_t> I(size, size); 00447 for (index_t i=0; i<size; ++i) 00448 { 00449 for (index_t j=0; j<size; ++j) 00450 I(i,j)=i==j ? scale : 0.0; 00451 } 00452 00453 return I; 00454 } 00455 00456 template <> 00457 SGMatrix<float64_t> SGMatrix<float64_t>::create_identity_matrix(index_t size, float64_t scale) 00458 { 00459 SGMatrix<float64_t> I(size, size); 00460 for (index_t i=0; i<size; ++i) 00461 { 00462 for (index_t j=0; j<size; ++j) 00463 I(i,j)=i==j ? scale : 0.0; 00464 } 00465 00466 return I; 00467 } 00468 00469 template <> 00470 SGMatrix<floatmax_t> SGMatrix<floatmax_t>::create_identity_matrix(index_t size, floatmax_t scale) 00471 { 00472 SGMatrix<floatmax_t> I(size, size); 00473 for (index_t i=0; i<size; ++i) 00474 { 00475 for (index_t j=0; j<size; ++j) 00476 I(i,j)=i==j ? scale : 0.0; 00477 } 00478 00479 return I; 00480 } 00481 00482 00483 template <class T> 00484 SGMatrix<float64_t> SGMatrix<T>::create_centering_matrix(index_t size) 00485 { 00486 SGMatrix<float64_t> H=SGMatrix<float64_t>::create_identity_matrix(size, 1.0); 00487 00488 float64_t subtract=1.0/size; 00489 for (index_t i=0; i<size; ++i) 00490 { 00491 for (index_t j=0; j<0; ++j) 00492 H(i,j)-=subtract; 00493 } 00494 00495 return H; 00496 } 00497 00498 //Howto construct the pseudo inverse (from "The Matrix Cookbook") 00499 // 00500 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m). 00501 // 00502 //The matrix A+ known as the pseudo inverse is unique and does always exist. 00503 // 00504 //The pseudo inverse A+ can be constructed from the singular value 00505 //decomposition A = UDV^T , by A^+ = V(D+)U^T. 00506 00507 #ifdef HAVE_LAPACK 00508 template <class T> 00509 float64_t* SGMatrix<T>::pinv( 00510 float64_t* matrix, int32_t rows, int32_t cols, float64_t* target) 00511 { 00512 if (!target) 00513 target=SG_MALLOC(float64_t, rows*cols); 00514 00515 char jobu='A'; 00516 char jobvt='A'; 00517 int m=rows; /* for calling external lib */ 00518 int n=cols; /* for calling external lib */ 00519 int lda=m; /* for calling external lib */ 00520 int ldu=m; /* for calling external lib */ 00521 int ldvt=n; /* for calling external lib */ 00522 int info=-1; /* for calling external lib */ 00523 int32_t lsize=CMath::min((int32_t) m, (int32_t) n); 00524 double* s=SG_MALLOC(double, lsize); 00525 double* u=SG_MALLOC(double, m*m); 00526 double* vt=SG_MALLOC(double, n*n); 00527 00528 wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info); 00529 ASSERT(info==0); 00530 00531 for (int32_t i=0; i<n; i++) 00532 { 00533 for (int32_t j=0; j<lsize; j++) 00534 vt[i*n+j]=vt[i*n+j]/s[j]; 00535 } 00536 00537 cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m); 00538 00539 SG_FREE(u); 00540 SG_FREE(vt); 00541 SG_FREE(s); 00542 00543 return target; 00544 } 00545 00547 template <class T> 00548 void SGMatrix<T>::inverse(SGMatrix<float64_t> matrix) 00549 { 00550 ASSERT(matrix.num_cols==matrix.num_rows); 00551 int32_t* ipiv = SG_MALLOC(int32_t, matrix.num_cols); 00552 clapack_dgetrf(CblasColMajor,matrix.num_cols,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv); 00553 clapack_dgetri(CblasColMajor,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv); 00554 SG_FREE(ipiv); 00555 } 00556 00557 template <class T> 00558 SGVector<float64_t> SGMatrix<T>::compute_eigenvectors(SGMatrix<float64_t> matrix) 00559 { 00560 if (matrix.num_rows!=matrix.num_rows) 00561 { 00562 SG_SERROR("SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix" 00563 " rows and columns are not equal!\n"); 00564 } 00565 00566 /* use reference counting for SGVector */ 00567 SGVector<float64_t> result(NULL, 0, true); 00568 result.vlen=matrix.num_rows; 00569 result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows, 00570 matrix.num_rows); 00571 return result; 00572 } 00573 00574 template <class T> 00575 double* SGMatrix<T>::compute_eigenvectors(double* matrix, int n, int m) 00576 { 00577 ASSERT(n == m); 00578 00579 char V='V'; 00580 char U='U'; 00581 int info; 00582 int ord=n; 00583 int lda=n; 00584 double* eigenvalues=SG_CALLOC(float64_t, n+1); 00585 00586 // lapack sym matrix eigenvalues+vectors 00587 wrap_dsyev(V, U, ord, matrix, lda, 00588 eigenvalues, &info); 00589 00590 if (info!=0) 00591 SG_SERROR("DSYEV failed with code %d\n", info); 00592 00593 return eigenvalues; 00594 } 00595 00596 template <class T> 00597 void SGMatrix<T>::compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors, 00598 int n, int il, int iu) 00599 { 00600 eigenvalues = SG_MALLOC(double, n); 00601 eigenvectors = SG_MALLOC(double, (iu-il+1)*n); 00602 int status = 0; 00603 wrap_dsyevr('V','U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status); 00604 ASSERT(status==0); 00605 } 00606 00607 #endif //HAVE_LAPACK 00608 00609 template <class T> 00610 SGMatrix<float64_t> SGMatrix<T>::matrix_multiply( 00611 SGMatrix<float64_t> A, SGMatrix<float64_t> B, 00612 bool transpose_A, bool transpose_B, float64_t scale) 00613 { 00614 /* these variables store size of transposed matrices*/ 00615 index_t cols_A=transpose_A ? A.num_rows : A.num_cols; 00616 index_t rows_A=transpose_A ? A.num_cols : A.num_rows; 00617 index_t rows_B=transpose_B ? B.num_cols : B.num_rows; 00618 index_t cols_B=transpose_B ? B.num_rows : B.num_cols; 00619 00620 /* do a dimension check */ 00621 if (cols_A!=rows_B) 00622 { 00623 SG_SERROR("SGMatrix::matrix_multiply(): Dimension mismatch: " 00624 "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B); 00625 } 00626 00627 /* allocate result matrix */ 00628 SGMatrix<float64_t> C(rows_A, cols_B); 00629 C.zero(); 00630 #ifdef HAVE_LAPACK 00631 /* multiply */ 00632 cblas_dgemm(CblasColMajor, 00633 transpose_A ? CblasTrans : CblasNoTrans, 00634 transpose_B ? CblasTrans : CblasNoTrans, 00635 rows_A, cols_B, cols_A, scale, 00636 A.matrix, A.num_rows, B.matrix, B.num_rows, 00637 0.0, C.matrix, C.num_rows); 00638 #else 00639 for (int32_t i=0; i<rows_A; i++) 00640 { 00641 for (int32_t j=0; j<cols_B; j++) 00642 { 00643 for (int32_t k=0; k<cols_A; k++) 00644 C(i,j) += A(i,k)*B(k,j); 00645 } 00646 } 00647 #endif //HAVE_LAPACK 00648 00649 return C; 00650 } 00651 00652 template<class T> 00653 SGMatrix<T> SGMatrix<T>::get_allocated_matrix(index_t num_rows, 00654 index_t num_cols, SGMatrix<T> pre_allocated) 00655 { 00656 SGMatrix<T> result; 00657 00658 /* evtl use pre-allocated space */ 00659 if (pre_allocated.matrix) 00660 { 00661 result=pre_allocated; 00662 00663 /* check dimension consistency */ 00664 if (pre_allocated.num_rows!=num_rows || 00665 pre_allocated.num_cols!=num_cols) 00666 { 00667 SG_SERROR("SGMatrix<T>::get_allocated_matrix(). Provided target" 00668 "matrix dimensions (%dx%d) do not match passed data " 00669 "dimensions (%dx%d)!\n", pre_allocated.num_rows, 00670 pre_allocated.num_cols, num_rows, num_cols); 00671 } 00672 } 00673 else 00674 { 00675 /* otherwise, allocate space */ 00676 result=SGMatrix<T>(num_rows, num_cols); 00677 } 00678 00679 return result; 00680 } 00681 00682 template class SGMatrix<bool>; 00683 template class SGMatrix<char>; 00684 template class SGMatrix<int8_t>; 00685 template class SGMatrix<uint8_t>; 00686 template class SGMatrix<int16_t>; 00687 template class SGMatrix<uint16_t>; 00688 template class SGMatrix<int32_t>; 00689 template class SGMatrix<uint32_t>; 00690 template class SGMatrix<int64_t>; 00691 template class SGMatrix<uint64_t>; 00692 template class SGMatrix<float32_t>; 00693 template class SGMatrix<float64_t>; 00694 template class SGMatrix<floatmax_t>; 00695 }