16 #include <shogun/lib/config.h> 47 m_on_gpu.store(
false, std::memory_order_release);
55 m_on_gpu.store(
false, std::memory_order_release);
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);
114 m_on_gpu.store(
false, std::memory_order_release);
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 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");
231 REQUIRE(
matrix!=
nullptr,
"The underlying matrix is not allocated!\n");
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");
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");
325 SG_SERROR(
"SGMatrix::max_single():: Not supported for complex128_t\n");
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];
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];
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];
474 matrix[i*num_rows+j]-=means[i];
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?
"" :
",")
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?
"" :
",")
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?
"" :
",")
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?
"" :
",")
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;
786 I(i,j)=i==j ? scale : (!
scale);
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);
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");
1184 SG_SERROR(
"SGMatrix::save():: Not supported for complex128_t\n");
1192 for (int64_t i = 0; i <
num_cols; i++)
1206 for (int64_t i=0; i<diag_vlen; i++)
std::shared_ptr< GPUMemoryBase< T > > gpu_ptr
static SGMatrix< float64_t > matrix_multiply(SGMatrix< float64_t > A, SGMatrix< float64_t > B, bool transpose_A=false, bool transpose_B=false, float64_t scale=1.0)
#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)
SGVector< T > get_row_vector(index_t row) const
std::complex< float64_t > complex128_t
virtual void set_matrix(const bool *matrix, int32_t num_feat, int32_t num_vec)
SGMatrix< T > & operator=(const SGMatrix< T > &)
bool is_symmetric() const
Eigen::Map< EigenMatrixXt, 0, Eigen::Stride< 0, 0 > > EigenMatrixXtMap
static T * get_row_sum(T *matrix, int32_t m, int32_t n)
void display_matrix(const char *name="matrix") const
static T sum(T *vec, int32_t len)
Return sum(vec)
static void transpose_matrix(T *&matrix, int32_t &num_feat, int32_t &num_vec)
void scale(SGVector< T > &a, SGVector< T > &result, T alpha=1)
static float64_t * pinv(float64_t *matrix, int32_t rows, int32_t cols, float64_t *target=NULL)
static T * get_column_sum(T *matrix, int32_t m, int32_t n)
#define SG_SNOTIMPLEMENTED
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)
SGMatrix< T > clone() const
static SGVector< float64_t > compute_eigenvectors(SGMatrix< float64_t > matrix)
void copy_refcount(const SGReferencedData &orig)
std::shared_ptr< GPUMemoryBase< T > > gpu_ptr
virtual void get_matrix(bool *&matrix, int32_t &num_feat, int32_t &num_vec)
static T * clone_matrix(const T *matrix, int32_t nrows, int32_t ncols)
shogun reference count managed data
static SGMatrix< T > create_identity_matrix(index_t size, T scale)
A File access base class.
static void create_diagonal_matrix(T *matrix, T *v, int32_t size)
void remove_column_mean()
static float64_t trace(float64_t *mat, int32_t cols, int32_t rows)
void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info)
static void inverse(SGMatrix< float64_t > matrix)
inverses square matrix in-place
void compute_few_eigenvectors(double *matrix_, double *&eigenvalues, double *&eigenvectors, int n, int il, int iu)
bool equals(const SGMatrix< T > &other) const
all of classes and functions are contained in the shogun namespace
Interface for GPU memory libraries.
#define REAL_IS_SYMMETRIC(real_t)
static void swap(T &a, T &b)
void set_const(T const_elem)
static SGMatrix< T > get_allocated_matrix(index_t num_rows, index_t num_cols, SGMatrix< T > pre_allocated=SGMatrix< T >())
SGVector< T > get_diagonal_vector() const
virtual void copy_data(const SGReferencedData &orig)
static void center_matrix(T *matrix, int32_t m, int32_t n)