16 #include <shogun/lib/config.h> 45 num_rows(nrows), num_cols(ncols), gpu_ptr(nullptr)
47 m_on_gpu.store(
false, std::memory_order_release);
53 num_rows(nrows), num_cols(ncols)
55 m_on_gpu.store(
false, std::memory_order_release);
60 :
SGReferencedData(ref_counting), num_rows(nrows), num_cols(ncols), gpu_ptr(nullptr)
62 matrix=SG_CALLOC(T, ((int64_t) nrows)*ncols);
63 m_on_gpu.store(
false, std::memory_order_release);
74 m_on_gpu.store(vec.
on_gpu(), std::memory_order_release);
82 REQUIRE(nrows>0,
"Number of rows (%d) has to be a positive integer!\n", nrows);
83 REQUIRE(ncols>0,
"Number of cols (%d) has to be a positive integer!\n", ncols);
84 REQUIRE(vec.
vlen==nrows*ncols,
"Number of elements in the matrix (%d) must " 85 "be the same as the number of elements in the vector (%d)!\n",
86 nrows*ncols, vec.
vlen);
92 m_on_gpu.store(vec.
on_gpu(), std::memory_order_release);
100 m_on_gpu.store(
true, std::memory_order_release);
112 num_rows(mat.rows()), num_cols(mat.cols()), gpu_ptr(nullptr)
114 m_on_gpu.store(
false, std::memory_order_release);
132 copy_refcount(other);
152 if (matrix==
nullptr || other.
matrix==
nullptr)
158 return std::equal(matrix, matrix+size(), other.
matrix);
162 #define REAL_EQUALS(real_t) \ 164 bool SGMatrix<real_t>::equals(const SGMatrix<real_t>& other) const \ 169 if (matrix==nullptr || other.matrix==nullptr) \ 172 if (num_rows!=other.num_rows || num_cols!=other.num_cols) \ 175 return std::equal(matrix, matrix+size(), other.matrix, \ 176 [](const real_t& a, const real_t& b) \ 178 return CMath::fequals<real_t>(a, b, std::numeric_limits<real_t>::epsilon()); \ 186 #endif // REAL_EQUALS 194 if (matrix==
nullptr || other.
matrix==
nullptr)
200 return std::equal(matrix, matrix+size(), other.
matrix,
203 return CMath::fequals<float64_t>(a.real(), b.real(), LDBL_EPSILON) &&
204 CMath::fequals<float64_t>(a.imag(), b.imag(), LDBL_EPSILON);
213 REQUIRE(matrix!=
nullptr,
"The underlying matrix is not allocated!\n");
214 REQUIRE(num_rows>0,
"Number of rows (%d) has to be positive!\n", num_rows);
215 REQUIRE(num_cols>0,
"Number of cols (%d) has to be positive!\n", num_cols);
217 std::fill(matrix, matrix+size(), const_elem);
231 REQUIRE(matrix!=
nullptr,
"The underlying matrix is not allocated!\n");
232 REQUIRE(num_rows>0,
"Number of rows (%d) has to be positive!\n", num_rows);
233 REQUIRE(num_cols>0,
"Number of cols (%d) has to be positive!\n", num_cols);
235 if (num_rows!=num_cols)
238 for (
index_t i=0; i<num_rows; ++i)
240 for (
index_t j=i+1; j<num_cols; ++j)
242 if (matrix[j*num_rows+i]!=matrix[i*num_rows+j])
250 #ifndef REAL_IS_SYMMETRIC 251 #define REAL_IS_SYMMETRIC(real_t) \ 253 bool SGMatrix<real_t>::is_symmetric() const \ 257 REQUIRE(matrix!=nullptr, "The underlying matrix is not allocated!\n"); \ 258 REQUIRE(num_rows>0, "Number of rows (%d) has to be positive!\n", num_rows); \ 259 REQUIRE(num_cols>0, "Number of cols (%d) has to be positive!\n", num_cols); \ 261 if (num_rows!=num_cols) \ 264 for (index_t i=0; i<num_rows; ++i) \ 266 for (index_t j=i+1; j<num_cols; ++j) \ 268 if (!CMath::fequals<real_t>(matrix[j*num_rows+i], \ 269 matrix[i*num_rows+j], std::numeric_limits<real_t>::epsilon())) \ 280 #undef REAL_IS_SYMMETRIC 281 #endif // REAL_IS_SYMMETRIC 288 REQUIRE(matrix!=
nullptr,
"The underlying matrix is not allocated!\n");
289 REQUIRE(num_rows>0,
"Number of rows (%d) has to be positive!\n", num_rows);
290 REQUIRE(num_cols>0,
"Number of cols (%d) has to be positive!\n", num_cols);
292 if (num_rows!=num_cols)
295 for (
index_t i=0; i<num_rows; ++i)
297 for (
index_t j=i+1; j<num_cols; ++j)
299 if (!(CMath::fequals<float64_t>(matrix[j*num_rows+i].real(),
300 matrix[i*num_rows+j].real(), DBL_EPSILON) &&
301 CMath::fequals<float64_t>(matrix[j*num_rows+i].imag(),
302 matrix[i*num_rows+j].imag(), DBL_EPSILON)))
315 REQUIRE(matrix!=
nullptr,
"The underlying matrix is not allocated!\n");
316 REQUIRE(num_rows>0,
"Number of rows (%d) has to be positive!\n", num_rows);
317 REQUIRE(num_cols>0,
"Number of cols (%d) has to be positive!\n", num_cols);
319 return *std::max_element(matrix, matrix+size());
325 SG_SERROR(
"SGMatrix::max_single():: Not supported for complex128_t\n");
334 return SGMatrix<T>(gpu_ptr->clone_vector(gpu_ptr.get(),
335 num_rows*num_cols), num_rows, num_cols);
339 return SGMatrix<T>(clone_matrix(matrix, num_rows, num_cols),
347 REQUIRE(matrix!=
nullptr,
"The underlying matrix is not allocated!\n");
348 REQUIRE(nrows>0,
"Number of rows (%d) has to be positive!\n", nrows);
349 REQUIRE(ncols>0,
"Number of cols (%d) has to be positive!\n", ncols);
351 auto size=int64_t(nrows)*ncols;
352 T* result=SG_MALLOC(T, size);
353 sg_memcpy(result, matrix, size*
sizeof(T));
361 T* transposed=SG_MALLOC(T, int64_t(num_vec)*num_feat);
362 for (int64_t i=0; i<num_vec; i++)
364 for (int64_t j=0; j<num_feat; j++)
365 transposed[i+j*num_vec]=matrix[i*num_feat+j];
371 CMath::swap(num_feat, num_vec);
378 for(int64_t i=0;i<size;i++)
380 for(int64_t j=0;j<size;j++)
383 matrix[j*size+i]=v[i];
392 float64_t* mat, int32_t cols, int32_t rows)
395 for (int64_t i=0; i<rows; i++)
396 trace+=mat[i*cols+i];
404 T* rowsums=SG_CALLOC(T, n);
406 for (int64_t i=0; i<n; i++)
408 for (int64_t j=0; j<m; j++)
409 rowsums[i]+=matrix[j+i*m];
418 T* colsums=SG_CALLOC(T, m);
420 for (int64_t i=0; i<n; i++)
422 for (int64_t j=0; j<m; j++)
423 colsums[j]+=matrix[j+i*m];
432 center_matrix(matrix, num_rows, num_cols);
440 T* colsums=get_column_sum(matrix, m,n);
441 T* rowsums=get_row_sum(matrix, m,n);
443 for (int32_t i=0; i<m; i++)
444 colsums[i]/=num_data;
445 for (int32_t j=0; j<n; j++)
446 rowsums[j]/=num_data;
450 for (int64_t i=0; i<n; i++)
452 for (int64_t j=0; j<m; j++)
453 matrix[i*m+j]+=s-colsums[j]-rowsums[i];
467 T* means=get_row_sum(matrix, num_rows, num_cols);
470 for (int64_t i=0; i<num_cols; ++i)
473 for (int64_t j=0; j<num_rows; ++j)
474 matrix[i*num_rows+j]-=means[i];
483 display_matrix(matrix, num_rows, num_cols, name);
496 const bool* matrix, int32_t rows, int32_t cols,
const char* name,
499 ASSERT(rows>=0 && cols>=0)
501 for (int64_t i=0; i<rows; i++)
504 for (int64_t j=0; j<cols; j++)
505 SG_SPRINT(
"%s\t%d%s", prefix, matrix[j*rows+i] ? 1 : 0,
506 j==cols-1?
"" :
",");
507 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
514 const char* matrix, int32_t rows, int32_t cols,
const char* name,
517 ASSERT(rows>=0 && cols>=0)
519 for (int64_t i=0; i<rows; i++)
522 for (int64_t j=0; j<cols; j++)
523 SG_SPRINT(
"%s\t%c%s", prefix, matrix[j*rows+i],
524 j==cols-1?
"" :
",");
525 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
532 const int8_t* matrix, int32_t rows, int32_t cols,
const char* name,
535 ASSERT(rows>=0 && cols>=0)
537 for (int64_t i=0; i<rows; i++)
540 for (int64_t j=0; j<cols; j++)
541 SG_SPRINT(
"%s\t%d%s", prefix, matrix[j*rows+i],
542 j==cols-1?
"" :
",");
543 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
550 const uint8_t* matrix, int32_t rows, int32_t cols,
const char* name,
553 ASSERT(rows>=0 && cols>=0)
555 for (int64_t i=0; i<rows; i++)
558 for (int64_t j=0; j<cols; j++)
559 SG_SPRINT(
"%s\t%d%s", prefix, matrix[j*rows+i],
560 j==cols-1?
"" :
",");
561 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
568 const int16_t* matrix, int32_t rows, int32_t cols,
const char* name,
571 ASSERT(rows>=0 && cols>=0)
573 for (int64_t i=0; i<rows; i++)
576 for (int64_t j=0; j<cols; j++)
577 SG_SPRINT(
"%s\t%d%s", prefix, matrix[j*rows+i],
578 j==cols-1?
"" :
",");
579 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
586 const uint16_t* matrix, int32_t rows, int32_t cols,
const char* name,
589 ASSERT(rows>=0 && cols>=0)
591 for (int64_t i=0; i<rows; i++)
594 for (int64_t j=0; j<cols; j++)
595 SG_SPRINT(
"%s\t%d%s", prefix, matrix[j*rows+i],
596 j==cols-1?
"" :
",");
597 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
605 const int32_t* matrix, int32_t rows, int32_t cols,
const char* name,
608 ASSERT(rows>=0 && cols>=0)
610 for (int64_t i=0; i<rows; i++)
613 for (int64_t j=0; j<cols; j++)
614 SG_SPRINT(
"%s\t%d%s", prefix, matrix[j*rows+i],
615 j==cols-1?
"" :
",");
616 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
623 const uint32_t* matrix, int32_t rows, int32_t cols,
const char* name,
626 ASSERT(rows>=0 && cols>=0)
628 for (int64_t i=0; i<rows; i++)
631 for (int64_t j=0; j<cols; j++)
632 SG_SPRINT(
"%s\t%d%s", prefix, matrix[j*rows+i],
633 j==cols-1?
"" :
",");
634 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
640 const int64_t* matrix, int32_t rows, int32_t cols,
const char* name,
643 ASSERT(rows>=0 && cols>=0)
645 for (int64_t i=0; i<rows; i++)
648 for (int64_t j=0; j<cols; j++)
649 SG_SPRINT(
"%s\t%d%s", prefix, matrix[j*rows+i],
650 j==cols-1?
"" :
",");
651 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
658 const uint64_t* matrix, int32_t rows, int32_t cols,
const char* name,
661 ASSERT(rows>=0 && cols>=0)
663 for (int64_t i=0; i<rows; i++)
666 for (int64_t j=0; j<cols; j++)
667 SG_SPRINT(
"%s\t%d%s", prefix, matrix[j*rows+i],
668 j==cols-1?
"" :
",");
669 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
676 const float32_t* matrix, int32_t rows, int32_t cols,
const char* name,
679 ASSERT(rows>=0 && cols>=0)
681 for (int64_t i=0; i<rows; i++)
684 for (int64_t j=0; j<cols; j++)
685 SG_SPRINT(
"%s\t%.18g%s", prefix, (
float) matrix[j*rows+i],
686 j==cols-1?
"" :
",");
687 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
694 const float64_t* matrix, int32_t rows, int32_t cols,
const char* name,
697 ASSERT(rows>=0 && cols>=0)
699 for (int64_t i=0; i<rows; i++)
702 for (int64_t j=0; j<cols; j++)
703 SG_SPRINT(
"%s\t%.18g%s", prefix, (
double) matrix[j*rows+i],
704 j==cols-1?
"" :
",");
705 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
712 const floatmax_t* matrix, int32_t rows, int32_t cols,
const char* name,
715 ASSERT(rows>=0 && cols>=0)
717 for (int64_t i=0; i<rows; i++)
720 for (int64_t j=0; j<cols; j++)
721 SG_SPRINT(
"%s\t%.18g%s", prefix, (
double) matrix[j*rows+i],
722 j==cols-1?
"" :
",");
723 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
730 const complex128_t* matrix, int32_t rows, int32_t cols,
const char* name,
733 ASSERT(rows>=0 && cols>=0)
735 for (int64_t i=0; i<rows; i++)
738 for (int64_t j=0; j<cols; j++)
739 SG_SPRINT(
"%s\t(%.18g+i%.18g)%s", prefix, matrix[j*rows+i].real(),
740 matrix[j*rows+i].imag(), j==cols-1?
"" :
",");
741 SG_SPRINT(
"%s]%s\n", prefix, i==rows-1?
"" :
",")
760 I(i,j)=i==j ?
scale : 0.0;
773 I(i,j)=i==j ?
scale : 0.0;
799 I(i,j)=i==j ?
scale : 0.0;
812 I(i,j)=i==j ?
scale : 0.0;
825 I(i,j)=i==j ?
scale : 0.0;
838 I(i,j)=i==j ?
scale : 0.0;
851 I(i,j)=i==j ?
scale : 0.0;
864 I(i,j)=i==j ?
scale : 0.0;
877 I(i,j)=i==j ?
scale : 0.0;
890 I(i,j)=i==j ?
scale : 0.0;
903 I(i,j)=i==j ?
scale : 0.0;
947 int32_t lsize=CMath::min((int32_t) m, (int32_t) n);
948 double* s=SG_MALLOC(
double, lsize);
949 double* u=SG_MALLOC(
double, m*m);
950 double* vt=SG_MALLOC(
double, n*n);
952 wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
955 for (int64_t i=0; i<n; i++)
957 for (int64_t j=0; j<lsize; j++)
958 vt[i*n+j]=vt[i*n+j]/s[j];
961 cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
974 REQUIRE(!matrix.
on_gpu(),
"Operation is not possible when data is in GPU memory.\n");
976 int32_t* ipiv = SG_MALLOC(int32_t, matrix.
num_cols);
985 REQUIRE(!matrix.
on_gpu(),
"Operation is not possible when data is in GPU memory.\n");
988 SG_SERROR(
"SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix" 989 " rows and columns are not equal!\n");
1010 double* eigenvalues=SG_CALLOC(
float64_t, n+1);
1014 eigenvalues, &info);
1017 SG_SERROR(
"DSYEV failed with code %d\n", info)
1024 int n,
int il,
int iu)
1026 eigenvalues = SG_MALLOC(
double, n);
1027 eigenvectors = SG_MALLOC(
double, (iu-il+1)*n);
1029 wrap_dsyevr(
'V',
'U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status);
1033 #endif //HAVE_LAPACK 1042 "Operation is not possible when data is in GPU memory.\n");
1053 SG_SERROR(
"SGMatrix::matrix_multiply(): Dimension mismatch: " 1054 "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
1062 cblas_dgemm(CblasColMajor,
1063 transpose_A ? CblasTrans : CblasNoTrans,
1064 transpose_B ? CblasTrans : CblasNoTrans,
1065 rows_A, cols_B, cols_A, scale,
1070 for (int32_t i=0; i<rows_A; i++)
1072 for (int32_t j=0; j<cols_B; j++)
1074 for (int32_t k=0; k<cols_A; k++)
1076 float64_t x1=transpose_A ? A(k,i):A(i,k);
1077 float64_t x2=transpose_B ? B(j,k):B(k,j);
1084 #endif //HAVE_LAPACK 1098 result=pre_allocated;
1101 if (pre_allocated.
num_rows!=num_rows ||
1104 SG_SERROR(
"SGMatrix<T>::get_allocated_matrix(). Provided target" 1105 "matrix dimensions (%dx%d) do not match passed data " 1106 "dimensions (%dx%d)!\n", pre_allocated.
num_rows,
1107 pre_allocated.
num_cols, num_rows, num_cols);
1122 gpu_ptr=((
SGMatrix*)(&orig))->gpu_ptr;
1123 matrix=((
SGMatrix*)(&orig))->matrix;
1124 num_rows=((
SGMatrix*)(&orig))->num_rows;
1125 num_cols=((
SGMatrix*)(&orig))->num_cols;
1126 m_on_gpu.store(((
SGMatrix*)(&orig))->m_on_gpu.load(
1127 std::memory_order_acquire), std::memory_order_release);
1137 m_on_gpu.store(
false, std::memory_order_release);
1168 SG_SERROR(
"SGMatrix::load():: Not supported for complex128_t\n");
1177 writer->
set_matrix(matrix, num_rows, num_cols);
1184 SG_SERROR(
"SGMatrix::save():: Not supported for complex128_t\n");
1192 for (int64_t i = 0; i < num_cols; i++)
1194 rowv[i] = matrix[i*num_rows+row];
1203 index_t diag_vlen=CMath::min(num_cols, num_rows);
1206 for (int64_t i=0; i<diag_vlen; i++)
1208 diag[i]=matrix[i*num_rows+i];
std::shared_ptr< GPUMemoryBase< T > > gpu_ptr
#define REAL_EQUALS(real_t)
void wrap_dsyevr(char jobz, char uplo, int n, double *a, int lda, int il, int iu, double *eigenvalues, double *eigenvectors, int *info)
std::complex< float64_t > complex128_t
virtual void set_matrix(const bool *matrix, int32_t num_feat, int32_t num_vec)
void scale(SGVector< T > &a, SGVector< T > &result, T alpha=1)
#define SG_SNOTIMPLEMENTED
void set_const(Container< T > &a, T value)
void display_matrix(const char *name="matrix") const
void wrap_dgesvd(char jobu, char jobvt, int m, int n, double *a, int lda, double *sing, double *u, int ldu, double *vt, int ldvt, int *info)
std::shared_ptr< GPUMemoryBase< T > > gpu_ptr
virtual void get_matrix(bool *&matrix, int32_t &num_feat, int32_t &num_vec)
shogun reference count managed data
A File access base class.
void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info)
all of classes and functions are contained in the shogun namespace
T sum(const Container< T > &a, bool no_diag=false)
Interface for GPU memory libraries.
#define REAL_IS_SYMMETRIC(real_t)