SHOGUN  6.1.3
SGMatrix.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2016-2017 Pan Deng
8  * Written (W) 2013-2017 Soumyajit De
9  * Written (W) 2011-2013 Heiko Strathmann
10  * Written (W) 2012 Fernando Jose Iglesias Garcia
11  * Written (W) 2010,2012 Soeren Sonnenburg
12  * Copyright (C) 2010 Berlin Institute of Technology
13  * Copyright (C) 2012 Soeren Sonnenburg
14  */
15 
16 #include <shogun/lib/config.h>
17 #include <shogun/io/SGIO.h>
18 #include <shogun/io/File.h>
19 #include <shogun/lib/SGMatrix.h>
20 #include <shogun/lib/SGVector.h>
24 #include <limits>
25 #include <algorithm>
26 
27 namespace shogun
28 {
29 
30 template <class T>
32 {
33  init_data();
34 }
35 
36 template <class T>
37 SGMatrix<T>::SGMatrix(bool ref_counting) : SGReferencedData(ref_counting)
38 {
39  init_data();
40 }
41 
42 template <class T>
43 SGMatrix<T>::SGMatrix(T* m, index_t nrows, index_t ncols, bool ref_counting)
44  : SGReferencedData(ref_counting), matrix(m),
45  num_rows(nrows), num_cols(ncols), gpu_ptr(nullptr)
46 {
47  m_on_gpu.store(false, std::memory_order_release);
48 }
49 
50 template <class T>
51 SGMatrix<T>::SGMatrix(T* m, index_t nrows, index_t ncols, index_t offset)
52  : SGReferencedData(false), matrix(m+offset),
53  num_rows(nrows), num_cols(ncols)
54 {
55  m_on_gpu.store(false, std::memory_order_release);
56 }
57 
58 template <class T>
59 SGMatrix<T>::SGMatrix(index_t nrows, index_t ncols, bool ref_counting)
60  : SGReferencedData(ref_counting), num_rows(nrows), num_cols(ncols), gpu_ptr(nullptr)
61 {
62  matrix=SG_CALLOC(T, ((int64_t) nrows)*ncols);
63  m_on_gpu.store(false, std::memory_order_release);
64 }
65 
66 template <class T>
68 {
69  REQUIRE((vec.vector || vec.gpu_ptr), "Vector not initialized!\n");
70  matrix=vec.vector;
71  num_rows=vec.vlen;
72  num_cols=1;
73  gpu_ptr = vec.gpu_ptr;
74  m_on_gpu.store(vec.on_gpu(), std::memory_order_release);
75 }
76 
77 template <class T>
79 : SGReferencedData(vec)
80 {
81  REQUIRE((vec.vector || vec.gpu_ptr), "Vector not initialized!\n");
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);
87 
88  matrix=vec.vector;
89  num_rows=nrows;
90  num_cols=ncols;
91  gpu_ptr = vec.gpu_ptr;
92  m_on_gpu.store(vec.on_gpu(), std::memory_order_release);
93 }
94 
95 template<class T>
97  : SGReferencedData(true), matrix(NULL), num_rows(nrows), num_cols(ncols),
98  gpu_ptr(std::shared_ptr<GPUMemoryBase<T>>(mat))
99 {
100  m_on_gpu.store(true, std::memory_order_release);
101 }
102 
103 template <class T>
105 {
106  copy_data(orig);
107 }
108 
109 template <class T>
111 : SGReferencedData(false), matrix(mat.data()),
112  num_rows(mat.rows()), num_cols(mat.cols()), gpu_ptr(nullptr)
113 {
114  m_on_gpu.store(false, std::memory_order_release);
115 }
116 
117 template <class T>
119 {
120  assert_on_cpu();
122 }
123 
124 template<class T>
126 {
127  if(&other == this)
128  return *this;
129 
130  unref();
131  copy_data(other);
132  copy_refcount(other);
133  ref();
134  return *this;
135 }
136 
137 template <class T>
139 {
140  unref();
141 }
142 
143 template <class T>
144 bool SGMatrix<T>::equals(const SGMatrix<T>& other) const
145 {
146  // avoid comparing elements when both are same.
147  // the case where both matrices are uninitialized is handled here as well.
148  if (*this==other)
149  return true;
150 
151  // avoid uninitialized memory read in case the matrices are not initialized
152  if (matrix==nullptr || other.matrix==nullptr)
153  return false;
154 
155  if (num_rows!=other.num_rows || num_cols!=other.num_cols)
156  return false;
157 
158  return std::equal(matrix, matrix+size(), other.matrix);
159 }
160 
161 #ifndef REAL_EQUALS
162 #define REAL_EQUALS(real_t) \
163 template <> \
164 bool SGMatrix<real_t>::equals(const SGMatrix<real_t>& other) const \
165 { \
166  if (*this==other) \
167  return true; \
168  \
169  if (matrix==nullptr || other.matrix==nullptr) \
170  return false; \
171  \
172  if (num_rows!=other.num_rows || num_cols!=other.num_cols) \
173  return false; \
174  \
175  return std::equal(matrix, matrix+size(), other.matrix, \
176  [](const real_t& a, const real_t& b) \
177  { \
178  return CMath::fequals<real_t>(a, b, std::numeric_limits<real_t>::epsilon()); \
179  }); \
180 }
181 
185 #undef REAL_EQUALS
186 #endif // REAL_EQUALS
187 
188 template <>
190 {
191  if (*this==other)
192  return true;
193 
194  if (matrix==nullptr || other.matrix==nullptr)
195  return false;
196 
197  if (num_rows!=other.num_rows || num_cols!=other.num_cols)
198  return false;
199 
200  return std::equal(matrix, matrix+size(), other.matrix,
201  [](const complex128_t& a, const complex128_t& b)
202  {
203  return CMath::fequals<float64_t>(a.real(), b.real(), LDBL_EPSILON) &&
204  CMath::fequals<float64_t>(a.imag(), b.imag(), LDBL_EPSILON);
205  });
206 }
207 
208 template <class T>
209 void SGMatrix<T>::set_const(T const_elem)
210 {
211  assert_on_cpu();
212 
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);
216 
217  std::fill(matrix, matrix+size(), const_elem);
218 }
219 
220 template <class T>
222 {
223  set_const(static_cast<T>(0));
224 }
225 
226 template <class T>
228 {
229  assert_on_cpu();
230 
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);
234 
235  if (num_rows!=num_cols)
236  return false;
237 
238  for (index_t i=0; i<num_rows; ++i)
239  {
240  for (index_t j=i+1; j<num_cols; ++j)
241  {
242  if (matrix[j*num_rows+i]!=matrix[i*num_rows+j])
243  return false;
244  }
245  }
246 
247  return true;
248 }
249 
250 #ifndef REAL_IS_SYMMETRIC
251 #define REAL_IS_SYMMETRIC(real_t) \
252 template <> \
253 bool SGMatrix<real_t>::is_symmetric() const \
254 { \
255  assert_on_cpu(); \
256  \
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); \
260  \
261  if (num_rows!=num_cols) \
262  return false; \
263  \
264  for (index_t i=0; i<num_rows; ++i) \
265  { \
266  for (index_t j=i+1; j<num_cols; ++j) \
267  { \
268  if (!CMath::fequals<real_t>(matrix[j*num_rows+i], \
269  matrix[i*num_rows+j], std::numeric_limits<real_t>::epsilon())) \
270  return false; \
271  } \
272  } \
273  \
274  return true; \
275 }
276 
280 #undef REAL_IS_SYMMETRIC
281 #endif // REAL_IS_SYMMETRIC
282 
283 template <>
285 {
286  assert_on_cpu();
287 
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);
291 
292  if (num_rows!=num_cols)
293  return false;
294 
295  for (index_t i=0; i<num_rows; ++i)
296  {
297  for (index_t j=i+1; j<num_cols; ++j)
298  {
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)))
303  return false;
304  }
305  }
306 
307  return true;
308 }
309 
310 template <class T>
312 {
313  assert_on_cpu();
314 
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);
318 
319  return *std::max_element(matrix, matrix+size());
320 }
321 
322 template <>
324 {
325  SG_SERROR("SGMatrix::max_single():: Not supported for complex128_t\n");
326  return complex128_t(0.0);
327 }
328 
329 template <class T>
331 {
332  if (on_gpu())
333  {
334  return SGMatrix<T>(gpu_ptr->clone_vector(gpu_ptr.get(),
336  }
337  else
338  {
340  num_rows, num_cols);
341  }
342 }
343 
344 template <class T>
345 T* SGMatrix<T>::clone_matrix(const T* matrix, int32_t nrows, int32_t ncols)
346 {
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);
350 
351  auto size=int64_t(nrows)*ncols;
352  T* result=SG_MALLOC(T, size);
353  sg_memcpy(result, matrix, size*sizeof(T));
354  return result;
355 }
356 
357 template <class T>
358 void SGMatrix<T>::transpose_matrix(T*& matrix, int32_t& num_feat, int32_t& num_vec)
359 {
360  /* this should be done in-place! Heiko */
361  T* transposed=SG_MALLOC(T, int64_t(num_vec)*num_feat);
362  for (int64_t i=0; i<num_vec; i++)
363  {
364  for (int64_t j=0; j<num_feat; j++)
365  transposed[i+j*num_vec]=matrix[i*num_feat+j];
366  }
367 
368  SG_FREE(matrix);
369  matrix=transposed;
370 
371  CMath::swap(num_feat, num_vec);
372 }
373 
374 template <class T>
376 {
377  /* Need assert v.size() */
378  for(int64_t i=0;i<size;i++)
379  {
380  for(int64_t j=0;j<size;j++)
381  {
382  if(i==j)
383  matrix[j*size+i]=v[i];
384  else
385  matrix[j*size+i]=0;
386  }
387  }
388 }
389 
390 template <class T>
392 {
393  assert_on_cpu();
394  return SGMatrix<T>(
395  get_column_vector(col_start), num_rows, col_end - col_start, false);
396 }
397 
398 template <class T>
400 {
401  assert_on_cpu();
402  return SGVector<T>(get_column_vector(col), num_rows, false);
403 }
404 
405 template <class T>
407 {
408  assert_on_cpu();
409  ASSERT(!vec.on_gpu())
410  ASSERT(vec.vlen == num_rows)
411  sg_memcpy(&matrix[num_rows * col], vec.vector, sizeof(T) * num_rows);
412 }
413 
414 template <class T>
415 float64_t SGMatrix<T>::trace(float64_t* mat, int32_t cols, int32_t rows)
416 {
417  float64_t trace=0;
418  for (int64_t i=0; i<rows; i++)
419  trace+=mat[i*cols+i];
420  return trace;
421 }
422 
423 /* Already in linalg */
424 template <class T>
425 T* SGMatrix<T>::get_row_sum(T* matrix, int32_t m, int32_t n)
426 {
427  T* rowsums=SG_CALLOC(T, n);
428 
429  for (int64_t i=0; i<n; i++)
430  {
431  for (int64_t j=0; j<m; j++)
432  rowsums[i]+=matrix[j+i*m];
433  }
434  return rowsums;
435 }
436 
437 /* Already in linalg */
438 template <class T>
439 T* SGMatrix<T>::get_column_sum(T* matrix, int32_t m, int32_t n)
440 {
441  T* colsums=SG_CALLOC(T, m);
442 
443  for (int64_t i=0; i<n; i++)
444  {
445  for (int64_t j=0; j<m; j++)
446  colsums[j]+=matrix[j+i*m];
447  }
448  return colsums;
449 }
450 
451 template <class T>
453 {
454  assert_on_cpu();
456 }
457 
458 template <class T>
459 void SGMatrix<T>::center_matrix(T* matrix, int32_t m, int32_t n)
460 {
461  float64_t num_data=n;
462 
463  T* colsums=get_column_sum(matrix, m,n);
464  T* rowsums=get_row_sum(matrix, m,n);
465 
466  for (int32_t i=0; i<m; i++)
467  colsums[i]/=num_data;
468  for (int32_t j=0; j<n; j++)
469  rowsums[j]/=num_data;
470 
471  T s=SGVector<T>::sum(rowsums, n)/num_data;
472 
473  for (int64_t i=0; i<n; i++)
474  {
475  for (int64_t j=0; j<m; j++)
476  matrix[i*m+j]+=s-colsums[j]-rowsums[i];
477  }
478 
479  SG_FREE(rowsums);
480  SG_FREE(colsums);
481 }
482 
483 template <class T>
485 {
486  assert_on_cpu();
487 
488  /* compute "row" sums (which I would call column sums), i.e. sum of all
489  * elements in a fixed column */
490  T* means=get_row_sum(matrix, num_rows, num_cols);
491 
492  /* substract column mean from every corresponding entry */
493  for (int64_t i=0; i<num_cols; ++i)
494  {
495  means[i]/=num_rows;
496  for (int64_t j=0; j<num_rows; ++j)
497  matrix[i*num_rows+j]-=means[i];
498  }
499 
500  SG_FREE(means);
501 }
502 
503 template<class T> void SGMatrix<T>::display_matrix(const char* name) const
504 {
505  assert_on_cpu();
507 }
508 
509 template <class T>
511  const SGMatrix<T> matrix, const char* name,
512  const char* prefix)
513 {
514  matrix.display_matrix();
515 }
516 
517 template <>
519  const bool* matrix, int32_t rows, int32_t cols, const char* name,
520  const char* prefix)
521 {
522  ASSERT(rows>=0 && cols>=0)
523  SG_SPRINT("%s%s=[\n", prefix, name)
524  for (int64_t i=0; i<rows; i++)
525  {
526  SG_SPRINT("%s[", prefix)
527  for (int64_t j=0; j<cols; j++)
528  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i] ? 1 : 0,
529  j==cols-1? "" : ",");
530  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
531  }
532  SG_SPRINT("%s]\n", prefix)
533 }
534 
535 template <>
537  const char* matrix, int32_t rows, int32_t cols, const char* name,
538  const char* prefix)
539 {
540  ASSERT(rows>=0 && cols>=0)
541  SG_SPRINT("%s%s=[\n", prefix, name)
542  for (int64_t i=0; i<rows; i++)
543  {
544  SG_SPRINT("%s[", prefix)
545  for (int64_t j=0; j<cols; j++)
546  SG_SPRINT("%s\t%c%s", prefix, matrix[j*rows+i],
547  j==cols-1? "" : ",");
548  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
549  }
550  SG_SPRINT("%s]\n", prefix)
551 }
552 
553 template <>
555  const int8_t* matrix, int32_t rows, int32_t cols, const char* name,
556  const char* prefix)
557 {
558  ASSERT(rows>=0 && cols>=0)
559  SG_SPRINT("%s%s=[\n", prefix, name)
560  for (int64_t i=0; i<rows; i++)
561  {
562  SG_SPRINT("%s[", prefix)
563  for (int64_t j=0; j<cols; j++)
564  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
565  j==cols-1? "" : ",");
566  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
567  }
568  SG_SPRINT("%s]\n", prefix)
569 }
570 
571 template <>
573  const uint8_t* matrix, int32_t rows, int32_t cols, const char* name,
574  const char* prefix)
575 {
576  ASSERT(rows>=0 && cols>=0)
577  SG_SPRINT("%s%s=[\n", prefix, name)
578  for (int64_t i=0; i<rows; i++)
579  {
580  SG_SPRINT("%s[", prefix)
581  for (int64_t j=0; j<cols; j++)
582  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
583  j==cols-1? "" : ",");
584  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
585  }
586  SG_SPRINT("%s]\n", prefix)
587 }
588 
589 template <>
591  const int16_t* matrix, int32_t rows, int32_t cols, const char* name,
592  const char* prefix)
593 {
594  ASSERT(rows>=0 && cols>=0)
595  SG_SPRINT("%s%s=[\n", prefix, name)
596  for (int64_t i=0; i<rows; i++)
597  {
598  SG_SPRINT("%s[", prefix)
599  for (int64_t j=0; j<cols; j++)
600  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
601  j==cols-1? "" : ",");
602  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
603  }
604  SG_SPRINT("%s]\n", prefix)
605 }
606 
607 template <>
609  const uint16_t* matrix, int32_t rows, int32_t cols, const char* name,
610  const char* prefix)
611 {
612  ASSERT(rows>=0 && cols>=0)
613  SG_SPRINT("%s%s=[\n", prefix, name)
614  for (int64_t i=0; i<rows; i++)
615  {
616  SG_SPRINT("%s[", prefix)
617  for (int64_t j=0; j<cols; j++)
618  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
619  j==cols-1? "" : ",");
620  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
621  }
622  SG_SPRINT("%s]\n", prefix)
623 }
624 
625 
626 template <>
628  const int32_t* matrix, int32_t rows, int32_t cols, const char* name,
629  const char* prefix)
630 {
631  ASSERT(rows>=0 && cols>=0)
632  SG_SPRINT("%s%s=[\n", prefix, name)
633  for (int64_t i=0; i<rows; i++)
634  {
635  SG_SPRINT("%s[", prefix)
636  for (int64_t j=0; j<cols; j++)
637  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
638  j==cols-1? "" : ",");
639  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
640  }
641  SG_SPRINT("%s]\n", prefix)
642 }
643 
644 template <>
646  const uint32_t* matrix, int32_t rows, int32_t cols, const char* name,
647  const char* prefix)
648 {
649  ASSERT(rows>=0 && cols>=0)
650  SG_SPRINT("%s%s=[\n", prefix, name)
651  for (int64_t i=0; i<rows; i++)
652  {
653  SG_SPRINT("%s[", prefix)
654  for (int64_t j=0; j<cols; j++)
655  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
656  j==cols-1? "" : ",");
657  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
658  }
659  SG_SPRINT("%s]\n", prefix)
660 }
661 template <>
663  const int64_t* matrix, int32_t rows, int32_t cols, const char* name,
664  const char* prefix)
665 {
666  ASSERT(rows>=0 && cols>=0)
667  SG_SPRINT("%s%s=[\n", prefix, name)
668  for (int64_t i=0; i<rows; i++)
669  {
670  SG_SPRINT("%s[", prefix)
671  for (int64_t j=0; j<cols; j++)
672  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
673  j==cols-1? "" : ",");
674  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
675  }
676  SG_SPRINT("%s]\n", prefix)
677 }
678 
679 template <>
681  const uint64_t* matrix, int32_t rows, int32_t cols, const char* name,
682  const char* prefix)
683 {
684  ASSERT(rows>=0 && cols>=0)
685  SG_SPRINT("%s%s=[\n", prefix, name)
686  for (int64_t i=0; i<rows; i++)
687  {
688  SG_SPRINT("%s[", prefix)
689  for (int64_t j=0; j<cols; j++)
690  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
691  j==cols-1? "" : ",");
692  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
693  }
694  SG_SPRINT("%s]\n", prefix)
695 }
696 
697 template <>
699  const float32_t* matrix, int32_t rows, int32_t cols, const char* name,
700  const char* prefix)
701 {
702  ASSERT(rows>=0 && cols>=0)
703  SG_SPRINT("%s%s=[\n", prefix, name)
704  for (int64_t i=0; i<rows; i++)
705  {
706  SG_SPRINT("%s[", prefix)
707  for (int64_t j=0; j<cols; j++)
708  SG_SPRINT("%s\t%.18g%s", prefix, (float) matrix[j*rows+i],
709  j==cols-1? "" : ",");
710  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
711  }
712  SG_SPRINT("%s]\n", prefix)
713 }
714 
715 template <>
717  const float64_t* matrix, int32_t rows, int32_t cols, const char* name,
718  const char* prefix)
719 {
720  ASSERT(rows>=0 && cols>=0)
721  SG_SPRINT("%s%s=[\n", prefix, name)
722  for (int64_t i=0; i<rows; i++)
723  {
724  SG_SPRINT("%s[", prefix)
725  for (int64_t j=0; j<cols; j++)
726  SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i],
727  j==cols-1? "" : ",");
728  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
729  }
730  SG_SPRINT("%s]\n", prefix)
731 }
732 
733 template <>
735  const floatmax_t* matrix, int32_t rows, int32_t cols, const char* name,
736  const char* prefix)
737 {
738  ASSERT(rows>=0 && cols>=0)
739  SG_SPRINT("%s%s=[\n", prefix, name)
740  for (int64_t i=0; i<rows; i++)
741  {
742  SG_SPRINT("%s[", prefix)
743  for (int64_t j=0; j<cols; j++)
744  SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i],
745  j==cols-1? "" : ",");
746  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
747  }
748  SG_SPRINT("%s]\n", prefix)
749 }
750 
751 template <>
753  const complex128_t* matrix, int32_t rows, int32_t cols, const char* name,
754  const char* prefix)
755 {
756  ASSERT(rows>=0 && cols>=0)
757  SG_SPRINT("%s%s=[\n", prefix, name)
758  for (int64_t i=0; i<rows; i++)
759  {
760  SG_SPRINT("%s[", prefix)
761  for (int64_t j=0; j<cols; j++)
762  SG_SPRINT("%s\t(%.18g+i%.18g)%s", prefix, matrix[j*rows+i].real(),
763  matrix[j*rows+i].imag(), j==cols-1? "" : ",");
764  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
765  }
766  SG_SPRINT("%s]\n", prefix)
767 }
768 
769 template <>
771 {
773  return SGMatrix<char>();
774 }
775 
776 template <>
778 {
779  SGMatrix<int8_t> identity_matrix(size, size);
780  for (index_t i=0; i<size; ++i)
781  {
782  for (index_t j=0; j<size; ++j)
783  identity_matrix(i,j)=i==j ? scale : 0.0;
784  }
785 
786  return identity_matrix;
787 }
788 
789 template <>
791 {
792  SGMatrix<uint8_t> identity_matrix(size, size);
793  for (index_t i=0; i<size; ++i)
794  {
795  for (index_t j=0; j<size; ++j)
796  identity_matrix(i,j)=i==j ? scale : 0.0;
797  }
798 
799  return identity_matrix;
800 }
801 
802 template <>
804 {
805  SGMatrix<bool> identity_matrix(size, size);
806  for (index_t i=0; i<size; ++i)
807  {
808  for (index_t j=0; j<size; ++j)
809  identity_matrix(i,j)=i==j ? scale : (!scale);
810  }
811 
812  return identity_matrix;
813 }
814 
815 template <>
817 {
818  SGMatrix<int16_t> identity_matrix(size, size);
819  for (index_t i=0; i<size; ++i)
820  {
821  for (index_t j=0; j<size; ++j)
822  identity_matrix(i,j)=i==j ? scale : 0.0;
823  }
824 
825  return identity_matrix;
826 }
827 
828 template <>
830 {
831  SGMatrix<uint16_t> identity_matrix(size, size);
832  for (index_t i=0; i<size; ++i)
833  {
834  for (index_t j=0; j<size; ++j)
835  identity_matrix(i,j)=i==j ? scale : 0.0;
836  }
837 
838  return identity_matrix;
839 }
840 
841 template <>
843 {
844  SGMatrix<int32_t> identity_matrix(size, size);
845  for (index_t i=0; i<size; ++i)
846  {
847  for (index_t j=0; j<size; ++j)
848  identity_matrix(i,j)=i==j ? scale : 0.0;
849  }
850 
851  return identity_matrix;
852 }
853 
854 template <>
856 {
857  SGMatrix<uint32_t> identity_matrix(size, size);
858  for (index_t i=0; i<size; ++i)
859  {
860  for (index_t j=0; j<size; ++j)
861  identity_matrix(i,j)=i==j ? scale : 0.0;
862  }
863 
864  return identity_matrix;
865 }
866 
867 template <>
869 {
870  SGMatrix<int64_t> identity_matrix(size, size);
871  for (index_t i=0; i<size; ++i)
872  {
873  for (index_t j=0; j<size; ++j)
874  identity_matrix(i,j)=i==j ? scale : 0.0;
875  }
876 
877  return identity_matrix;
878 }
879 
880 template <>
882 {
883  SGMatrix<uint64_t> identity_matrix(size, size);
884  for (index_t i=0; i<size; ++i)
885  {
886  for (index_t j=0; j<size; ++j)
887  identity_matrix(i,j)=i==j ? scale : 0.0;
888  }
889 
890  return identity_matrix;
891 }
892 
893 template <>
895 {
896  SGMatrix<float32_t> identity_matrix(size, size);
897  for (index_t i=0; i<size; ++i)
898  {
899  for (index_t j=0; j<size; ++j)
900  identity_matrix(i,j)=i==j ? scale : 0.0;
901  }
902 
903  return identity_matrix;
904 }
905 
906 template <>
908 {
909  SGMatrix<float64_t> identity_matrix(size, size);
910  for (index_t i=0; i<size; ++i)
911  {
912  for (index_t j=0; j<size; ++j)
913  identity_matrix(i,j)=i==j ? scale : 0.0;
914  }
915 
916  return identity_matrix;
917 }
918 
919 template <>
921 {
922  SGMatrix<floatmax_t> identity_matrix(size, size);
923  for (index_t i=0; i<size; ++i)
924  {
925  for (index_t j=0; j<size; ++j)
926  identity_matrix(i,j)=i==j ? scale : 0.0;
927  }
928 
929  return identity_matrix;
930 }
931 
932 template <>
934 {
935  SGMatrix<complex128_t> identity_matrix(size, size);
936  for (index_t i=0; i<size; ++i)
937  {
938  for (index_t j=0; j<size; ++j)
939  identity_matrix(i,j)=i==j ? scale : complex128_t(0.0);
940  }
941 
942  return identity_matrix;
943 }
944 
945 //Howto construct the pseudo inverse (from "The Matrix Cookbook")
946 //
947 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
948 //
949 //The matrix A+ known as the pseudo inverse is unique and does always exist.
950 //
951 //The pseudo inverse A+ can be constructed from the singular value
952 //decomposition A = UDV^T , by A^+ = V(D+)U^T.
953 
954 #ifdef HAVE_LAPACK
955 template <class T>
957  float64_t* matrix, int32_t rows, int32_t cols, float64_t* target)
958 {
959  if (!target)
960  target=SG_MALLOC(float64_t, rows*cols);
961 
962  char jobu='A';
963  char jobvt='A';
964  int m=rows; /* for calling external lib */
965  int n=cols; /* for calling external lib */
966  int lda=m; /* for calling external lib */
967  int ldu=m; /* for calling external lib */
968  int ldvt=n; /* for calling external lib */
969  int info=-1; /* for calling external lib */
970  int32_t lsize=CMath::min((int32_t) m, (int32_t) n);
971  double* s=SG_MALLOC(double, lsize);
972  double* u=SG_MALLOC(double, m*m);
973  double* vt=SG_MALLOC(double, n*n);
974 
975  wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
976  ASSERT(info==0)
977 
978  for (int64_t i=0; i<n; i++)
979  {
980  for (int64_t j=0; j<lsize; j++)
981  vt[i*n+j]=vt[i*n+j]/s[j];
982  }
983 
984  cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
985 
986  SG_FREE(u);
987  SG_FREE(vt);
988  SG_FREE(s);
989 
990  return target;
991 }
992 
994 template <class T>
996 {
997  REQUIRE(!matrix.on_gpu(), "Operation is not possible when data is in GPU memory.\n");
998  ASSERT(matrix.num_cols==matrix.num_rows)
999  int32_t* ipiv = SG_MALLOC(int32_t, matrix.num_cols);
1000  clapack_dgetrf(CblasColMajor,matrix.num_cols,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
1001  clapack_dgetri(CblasColMajor,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
1002  SG_FREE(ipiv);
1003 }
1004 
1005 template <class T>
1007 {
1008  REQUIRE(!matrix.on_gpu(), "Operation is not possible when data is in GPU memory.\n");
1009  if (matrix.num_rows!=matrix.num_cols)
1010  {
1011  SG_SERROR("SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix"
1012  " rows and columns are not equal!\n");
1013  }
1014 
1015  /* use reference counting for SGVector */
1016  SGVector<float64_t> result(NULL, 0, true);
1017  result.vlen=matrix.num_rows;
1018  result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows,
1019  matrix.num_rows);
1020  return result;
1021 }
1022 
1023 template <class T>
1024 double* SGMatrix<T>::compute_eigenvectors(double* matrix, int n, int m)
1025 {
1026  ASSERT(n == m)
1027 
1028  char V='V';
1029  char U='U';
1030  int info;
1031  int ord=n;
1032  int lda=n;
1033  double* eigenvalues=SG_CALLOC(float64_t, n+1);
1034 
1035  // lapack sym matrix eigenvalues+vectors
1036  wrap_dsyev(V, U, ord, matrix, lda,
1037  eigenvalues, &info);
1038 
1039  if (info!=0)
1040  SG_SERROR("DSYEV failed with code %d\n", info)
1041 
1042  return eigenvalues;
1043 }
1044 
1045 template <class T>
1046 void SGMatrix<T>::compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors,
1047  int n, int il, int iu)
1048 {
1049  eigenvalues = SG_MALLOC(double, n);
1050  eigenvectors = SG_MALLOC(double, (iu-il+1)*n);
1051  int status = 0;
1052  wrap_dsyevr('V','U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status);
1053  ASSERT(status==0)
1054 }
1055 
1056 #endif //HAVE_LAPACK
1057 
1058 /* Already in linalg */
1059 template <class T>
1062  bool transpose_A, bool transpose_B, float64_t scale)
1063 {
1064  REQUIRE((!A.on_gpu()) && (!B.on_gpu()),
1065  "Operation is not possible when data is in GPU memory.\n");
1066 
1067  /* these variables store size of transposed matrices*/
1068  index_t cols_A=transpose_A ? A.num_rows : A.num_cols;
1069  index_t rows_A=transpose_A ? A.num_cols : A.num_rows;
1070  index_t rows_B=transpose_B ? B.num_cols : B.num_rows;
1071  index_t cols_B=transpose_B ? B.num_rows : B.num_cols;
1072 
1073  /* do a dimension check */
1074  if (cols_A!=rows_B)
1075  {
1076  SG_SERROR("SGMatrix::matrix_multiply(): Dimension mismatch: "
1077  "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
1078  }
1079 
1080  /* allocate result matrix */
1081  SGMatrix<float64_t> C(rows_A, cols_B);
1082  C.zero();
1083 #ifdef HAVE_LAPACK
1084  /* multiply */
1085  cblas_dgemm(CblasColMajor,
1086  transpose_A ? CblasTrans : CblasNoTrans,
1087  transpose_B ? CblasTrans : CblasNoTrans,
1088  rows_A, cols_B, cols_A, scale,
1089  A.matrix, A.num_rows, B.matrix, B.num_rows,
1090  0.0, C.matrix, C.num_rows);
1091 #else
1092  /* C(i,j) = scale * \Sigma A(i,k)*B(k,j) */
1093  for (int32_t i=0; i<rows_A; i++)
1094  {
1095  for (int32_t j=0; j<cols_B; j++)
1096  {
1097  for (int32_t k=0; k<cols_A; k++)
1098  {
1099  float64_t x1=transpose_A ? A(k,i):A(i,k);
1100  float64_t x2=transpose_B ? B(j,k):B(k,j);
1101  C(i,j)+=x1*x2;
1102  }
1103 
1104  C(i,j)*=scale;
1105  }
1106  }
1107 #endif //HAVE_LAPACK
1108 
1109  return C;
1110 }
1111 
1112 template<class T>
1114  index_t num_cols, SGMatrix<T> pre_allocated)
1115 {
1116  SGMatrix<T> result;
1117 
1118  /* evtl use pre-allocated space */
1119  if (pre_allocated.matrix || pre_allocated.gpu_ptr)
1120  {
1121  result=pre_allocated;
1122 
1123  /* check dimension consistency */
1124  if (pre_allocated.num_rows!=num_rows ||
1125  pre_allocated.num_cols!=num_cols)
1126  {
1127  SG_SERROR("SGMatrix<T>::get_allocated_matrix(). Provided target"
1128  "matrix dimensions (%dx%d) do not match passed data "
1129  "dimensions (%dx%d)!\n", pre_allocated.num_rows,
1130  pre_allocated.num_cols, num_rows, num_cols);
1131  }
1132  }
1133  else
1134  {
1135  /* otherwise, allocate space */
1136  result=SGMatrix<T>(num_rows, num_cols);
1137  }
1138 
1139  return result;
1140 }
1141 
1142 template<class T>
1144 {
1145  gpu_ptr=((SGMatrix*)(&orig))->gpu_ptr;
1146  matrix=((SGMatrix*)(&orig))->matrix;
1147  num_rows=((SGMatrix*)(&orig))->num_rows;
1148  num_cols=((SGMatrix*)(&orig))->num_cols;
1149  m_on_gpu.store(((SGMatrix*)(&orig))->m_on_gpu.load(
1150  std::memory_order_acquire), std::memory_order_release);
1151 }
1152 
1153 template<class T>
1155 {
1156  matrix=NULL;
1157  num_rows=0;
1158  num_cols=0;
1159  gpu_ptr=nullptr;
1160  m_on_gpu.store(false, std::memory_order_release);
1161 }
1162 
1163 template<class T>
1165 {
1166  SG_FREE(matrix);
1167  matrix=NULL;
1168  num_rows=0;
1169  num_cols=0;
1170 }
1171 
1172 template<class T>
1174 {
1175  ASSERT(loader)
1176  unref();
1177 
1179  SGMatrix<T> mat;
1180  loader->get_matrix(mat.matrix, mat.num_rows, mat.num_cols);
1181  mat.gpu_ptr = nullptr;
1182  copy_data(mat);
1183  copy_refcount(mat);
1184  ref();
1186 }
1187 
1188 template<>
1190 {
1191  SG_SERROR("SGMatrix::load():: Not supported for complex128_t\n");
1192 }
1193 
1194 template<class T>
1196 {
1197  assert_on_cpu();
1198  ASSERT(writer)
1200  writer->set_matrix(matrix, num_rows, num_cols);
1202 }
1203 
1204 template<>
1206 {
1207  SG_SERROR("SGMatrix::save():: Not supported for complex128_t\n");
1208 }
1209 
1210 template<class T>
1212 {
1213  assert_on_cpu();
1214  SGVector<T> rowv(num_cols);
1215  for (int64_t i = 0; i < num_cols; i++)
1216  {
1217  rowv[i] = matrix[i*num_rows+row];
1218  }
1219  return rowv;
1220 }
1221 
1222 template<class T>
1224 {
1225  assert_on_cpu();
1226  index_t diag_vlen=CMath::min(num_cols, num_rows);
1227  SGVector<T> diag(diag_vlen);
1228 
1229  for (int64_t i=0; i<diag_vlen; i++)
1230  {
1231  diag[i]=matrix[i*num_rows+i];
1232  }
1233 
1234  return diag;
1235 }
1236 
1237 template class SGMatrix<bool>;
1238 template class SGMatrix<char>;
1239 template class SGMatrix<int8_t>;
1240 template class SGMatrix<uint8_t>;
1241 template class SGMatrix<int16_t>;
1242 template class SGMatrix<uint16_t>;
1243 template class SGMatrix<int32_t>;
1244 template class SGMatrix<uint32_t>;
1245 template class SGMatrix<int64_t>;
1246 template class SGMatrix<uint64_t>;
1247 template class SGMatrix<float32_t>;
1248 template class SGMatrix<float64_t>;
1249 template class SGMatrix<floatmax_t>;
1250 template class SGMatrix<complex128_t>;
1251 }
std::shared_ptr< GPUMemoryBase< T > > gpu_ptr
Definition: SGVector.h:573
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)
Definition: SGMatrix.cpp:1060
#define REAL_EQUALS(real_t)
Definition: SGMatrix.cpp:162
void wrap_dsyevr(char jobz, char uplo, int n, double *a, int lda, int il, int iu, double *eigenvalues, double *eigenvectors, int *info)
Definition: lapack.cpp:310
#define SG_RESET_LOCALE
Definition: SGIO.h:85
SGVector< T > get_row_vector(index_t row) const
Definition: SGMatrix.cpp:1211
std::complex< float64_t > complex128_t
Definition: common.h:77
virtual void set_matrix(const bool *matrix, int32_t num_feat, int32_t num_vec)
Definition: File.cpp:126
SGMatrix< T > & operator=(const SGMatrix< T > &)
Definition: SGMatrix.cpp:125
SGVector< T > get_column(index_t col) const
Definition: SGMatrix.cpp:399
bool is_symmetric() const
Definition: SGMatrix.cpp:227
Eigen::Map< EigenMatrixXt, 0, Eigen::Stride< 0, 0 > > EigenMatrixXtMap
Definition: SGMatrix.h:48
static T * get_row_sum(T *matrix, int32_t m, int32_t n)
Definition: SGMatrix.cpp:425
void display_matrix(const char *name="matrix") const
Definition: SGMatrix.cpp:503
int32_t index_t
Definition: common.h:72
bool on_gpu() const
Definition: SGMatrix.h:83
static T sum(T *vec, int32_t len)
Return sum(vec)
Definition: SGVector.h:418
static void transpose_matrix(T *&matrix, int32_t &num_feat, int32_t &num_vec)
Definition: SGMatrix.cpp:358
void scale(SGVector< T > &a, SGVector< T > &result, T alpha=1)
void load(CFile *loader)
Definition: SGMatrix.cpp:1173
static float64_t * pinv(float64_t *matrix, int32_t rows, int32_t cols, float64_t *target=NULL)
Definition: SGMatrix.cpp:956
Definition: basetag.h:132
#define REQUIRE(x,...)
Definition: SGIO.h:181
static T * get_column_sum(T *matrix, int32_t m, int32_t n)
Definition: SGMatrix.cpp:439
#define SG_SNOTIMPLEMENTED
Definition: SGIO.h:172
SGMatrix< T > submatrix(index_t col_start, index_t col_end) const
Definition: SGMatrix.cpp:391
#define SG_SET_LOCALE_C
Definition: SGIO.h:84
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)
Definition: lapack.cpp:256
SGMatrix< T > clone() const
Definition: SGMatrix.cpp:330
void save(CFile *saver)
Definition: SGMatrix.cpp:1195
#define SG_SPRINT(...)
Definition: SGIO.h:165
#define ASSERT(x)
Definition: SGIO.h:176
virtual void init_data()
Definition: SGMatrix.cpp:1154
static SGVector< float64_t > compute_eigenvectors(SGMatrix< float64_t > matrix)
Definition: SGMatrix.cpp:1006
void copy_refcount(const SGReferencedData &orig)
std::shared_ptr< GPUMemoryBase< T > > gpu_ptr
Definition: SGMatrix.h:499
virtual void get_matrix(bool *&matrix, int32_t &num_feat, int32_t &num_vec)
Definition: File.cpp:109
static T * clone_matrix(const T *matrix, int32_t nrows, int32_t ncols)
Definition: SGMatrix.cpp:345
shogun reference count managed data
double float64_t
Definition: common.h:60
long double floatmax_t
Definition: common.h:61
static SGMatrix< T > create_identity_matrix(index_t size, T scale)
A File access base class.
Definition: File.h:34
T max_single() const
Definition: SGMatrix.cpp:311
static void create_diagonal_matrix(T *matrix, T *v, int32_t size)
Definition: SGMatrix.cpp:375
index_t num_rows
Definition: SGMatrix.h:495
shogun vector
virtual ~SGMatrix()
Definition: SGMatrix.cpp:138
void remove_column_mean()
Definition: SGMatrix.cpp:484
static float64_t trace(float64_t *mat, int32_t cols, int32_t rows)
Definition: SGMatrix.cpp:415
index_t num_cols
Definition: SGMatrix.h:497
void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info)
Definition: lapack.cpp:238
static void inverse(SGMatrix< float64_t > matrix)
inverses square matrix in-place
Definition: SGMatrix.cpp:995
float float32_t
Definition: common.h:59
shogun matrix
void compute_few_eigenvectors(double *matrix_, double *&eigenvalues, double *&eigenvectors, int n, int il, int iu)
Definition: SGMatrix.cpp:1046
bool equals(const SGMatrix< T > &other) const
Definition: SGMatrix.cpp:144
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
Interface for GPU memory libraries.
Definition: SGMatrix.h:35
#define SG_SERROR(...)
Definition: SGIO.h:164
T * data() const
Definition: SGMatrix.h:269
#define REAL_IS_SYMMETRIC(real_t)
Definition: SGMatrix.cpp:251
static void swap(T &a, T &b)
Definition: Math.h:406
void set_const(T const_elem)
Definition: SGMatrix.cpp:209
static SGMatrix< T > get_allocated_matrix(index_t num_rows, index_t num_cols, SGMatrix< T > pre_allocated=SGMatrix< T >())
Definition: SGMatrix.cpp:1113
SGVector< T > get_diagonal_vector() const
Definition: SGMatrix.cpp:1223
bool on_gpu() const
Definition: SGVector.h:91
virtual void copy_data(const SGReferencedData &orig)
Definition: SGMatrix.cpp:1143
static T min(T a, T b)
Definition: Math.h:138
void set_column(index_t col, const SGVector< T > vec)
Definition: SGMatrix.cpp:406
T * get_column_vector(index_t col) const
Definition: SGMatrix.h:144
static void center_matrix(T *matrix, int32_t m, int32_t n)
Definition: SGMatrix.cpp:459
virtual void free_data()
Definition: SGMatrix.cpp:1164
int64_t size() const
Definition: SGMatrix.h:275
index_t vlen
Definition: SGVector.h:571

SHOGUN Machine Learning Toolbox - Documentation