SHOGUN  6.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
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>
23 #include <limits>
24 #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();
121  return EigenMatrixXtMap(matrix, num_rows, num_cols);
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(),
335  num_rows*num_cols), num_rows, num_cols);
336  }
337  else
338  {
339  return SGMatrix<T>(clone_matrix(matrix, num_rows, num_cols),
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>
375 void SGMatrix<T>::create_diagonal_matrix(T* matrix, T* v,int32_t size)
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  float64_t* mat, int32_t cols, int32_t rows)
393 {
394  float64_t trace=0;
395  for (int64_t i=0; i<rows; i++)
396  trace+=mat[i*cols+i];
397  return trace;
398 }
399 
400 /* Already in linalg */
401 template <class T>
402 T* SGMatrix<T>::get_row_sum(T* matrix, int32_t m, int32_t n)
403 {
404  T* rowsums=SG_CALLOC(T, n);
405 
406  for (int64_t i=0; i<n; i++)
407  {
408  for (int64_t j=0; j<m; j++)
409  rowsums[i]+=matrix[j+i*m];
410  }
411  return rowsums;
412 }
413 
414 /* Already in linalg */
415 template <class T>
416 T* SGMatrix<T>::get_column_sum(T* matrix, int32_t m, int32_t n)
417 {
418  T* colsums=SG_CALLOC(T, m);
419 
420  for (int64_t i=0; i<n; i++)
421  {
422  for (int64_t j=0; j<m; j++)
423  colsums[j]+=matrix[j+i*m];
424  }
425  return colsums;
426 }
427 
428 template <class T>
430 {
431  assert_on_cpu();
432  center_matrix(matrix, num_rows, num_cols);
433 }
434 
435 template <class T>
436 void SGMatrix<T>::center_matrix(T* matrix, int32_t m, int32_t n)
437 {
438  float64_t num_data=n;
439 
440  T* colsums=get_column_sum(matrix, m,n);
441  T* rowsums=get_row_sum(matrix, m,n);
442 
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;
447 
448  T s=SGVector<T>::sum(rowsums, n)/num_data;
449 
450  for (int64_t i=0; i<n; i++)
451  {
452  for (int64_t j=0; j<m; j++)
453  matrix[i*m+j]+=s-colsums[j]-rowsums[i];
454  }
455 
456  SG_FREE(rowsums);
457  SG_FREE(colsums);
458 }
459 
460 template <class T>
462 {
463  assert_on_cpu();
464 
465  /* compute "row" sums (which I would call column sums), i.e. sum of all
466  * elements in a fixed column */
467  T* means=get_row_sum(matrix, num_rows, num_cols);
468 
469  /* substract column mean from every corresponding entry */
470  for (int64_t i=0; i<num_cols; ++i)
471  {
472  means[i]/=num_rows;
473  for (int64_t j=0; j<num_rows; ++j)
474  matrix[i*num_rows+j]-=means[i];
475  }
476 
477  SG_FREE(means);
478 }
479 
480 template<class T> void SGMatrix<T>::display_matrix(const char* name) const
481 {
482  assert_on_cpu();
483  display_matrix(matrix, num_rows, num_cols, name);
484 }
485 
486 template <class T>
488  const SGMatrix<T> matrix, const char* name,
489  const char* prefix)
490 {
491  matrix.display_matrix();
492 }
493 
494 template <>
496  const bool* matrix, int32_t rows, int32_t cols, const char* name,
497  const char* prefix)
498 {
499  ASSERT(rows>=0 && cols>=0)
500  SG_SPRINT("%s%s=[\n", prefix, name)
501  for (int64_t i=0; i<rows; i++)
502  {
503  SG_SPRINT("%s[", prefix)
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? "" : ",")
508  }
509  SG_SPRINT("%s]\n", prefix)
510 }
511 
512 template <>
514  const char* matrix, int32_t rows, int32_t cols, const char* name,
515  const char* prefix)
516 {
517  ASSERT(rows>=0 && cols>=0)
518  SG_SPRINT("%s%s=[\n", prefix, name)
519  for (int64_t i=0; i<rows; i++)
520  {
521  SG_SPRINT("%s[", prefix)
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? "" : ",")
526  }
527  SG_SPRINT("%s]\n", prefix)
528 }
529 
530 template <>
532  const int8_t* matrix, int32_t rows, int32_t cols, const char* name,
533  const char* prefix)
534 {
535  ASSERT(rows>=0 && cols>=0)
536  SG_SPRINT("%s%s=[\n", prefix, name)
537  for (int64_t i=0; i<rows; i++)
538  {
539  SG_SPRINT("%s[", prefix)
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? "" : ",")
544  }
545  SG_SPRINT("%s]\n", prefix)
546 }
547 
548 template <>
550  const uint8_t* matrix, int32_t rows, int32_t cols, const char* name,
551  const char* prefix)
552 {
553  ASSERT(rows>=0 && cols>=0)
554  SG_SPRINT("%s%s=[\n", prefix, name)
555  for (int64_t i=0; i<rows; i++)
556  {
557  SG_SPRINT("%s[", prefix)
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? "" : ",")
562  }
563  SG_SPRINT("%s]\n", prefix)
564 }
565 
566 template <>
568  const int16_t* matrix, int32_t rows, int32_t cols, const char* name,
569  const char* prefix)
570 {
571  ASSERT(rows>=0 && cols>=0)
572  SG_SPRINT("%s%s=[\n", prefix, name)
573  for (int64_t i=0; i<rows; i++)
574  {
575  SG_SPRINT("%s[", prefix)
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? "" : ",")
580  }
581  SG_SPRINT("%s]\n", prefix)
582 }
583 
584 template <>
586  const uint16_t* matrix, int32_t rows, int32_t cols, const char* name,
587  const char* prefix)
588 {
589  ASSERT(rows>=0 && cols>=0)
590  SG_SPRINT("%s%s=[\n", prefix, name)
591  for (int64_t i=0; i<rows; i++)
592  {
593  SG_SPRINT("%s[", prefix)
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? "" : ",")
598  }
599  SG_SPRINT("%s]\n", prefix)
600 }
601 
602 
603 template <>
605  const int32_t* matrix, int32_t rows, int32_t cols, const char* name,
606  const char* prefix)
607 {
608  ASSERT(rows>=0 && cols>=0)
609  SG_SPRINT("%s%s=[\n", prefix, name)
610  for (int64_t i=0; i<rows; i++)
611  {
612  SG_SPRINT("%s[", prefix)
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? "" : ",")
617  }
618  SG_SPRINT("%s]\n", prefix)
619 }
620 
621 template <>
623  const uint32_t* matrix, int32_t rows, int32_t cols, const char* name,
624  const char* prefix)
625 {
626  ASSERT(rows>=0 && cols>=0)
627  SG_SPRINT("%s%s=[\n", prefix, name)
628  for (int64_t i=0; i<rows; i++)
629  {
630  SG_SPRINT("%s[", prefix)
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? "" : ",")
635  }
636  SG_SPRINT("%s]\n", prefix)
637 }
638 template <>
640  const int64_t* matrix, int32_t rows, int32_t cols, const char* name,
641  const char* prefix)
642 {
643  ASSERT(rows>=0 && cols>=0)
644  SG_SPRINT("%s%s=[\n", prefix, name)
645  for (int64_t i=0; i<rows; i++)
646  {
647  SG_SPRINT("%s[", prefix)
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? "" : ",")
652  }
653  SG_SPRINT("%s]\n", prefix)
654 }
655 
656 template <>
658  const uint64_t* matrix, int32_t rows, int32_t cols, const char* name,
659  const char* prefix)
660 {
661  ASSERT(rows>=0 && cols>=0)
662  SG_SPRINT("%s%s=[\n", prefix, name)
663  for (int64_t i=0; i<rows; i++)
664  {
665  SG_SPRINT("%s[", prefix)
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? "" : ",")
670  }
671  SG_SPRINT("%s]\n", prefix)
672 }
673 
674 template <>
676  const float32_t* matrix, int32_t rows, int32_t cols, const char* name,
677  const char* prefix)
678 {
679  ASSERT(rows>=0 && cols>=0)
680  SG_SPRINT("%s%s=[\n", prefix, name)
681  for (int64_t i=0; i<rows; i++)
682  {
683  SG_SPRINT("%s[", prefix)
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? "" : ",")
688  }
689  SG_SPRINT("%s]\n", prefix)
690 }
691 
692 template <>
694  const float64_t* matrix, int32_t rows, int32_t cols, const char* name,
695  const char* prefix)
696 {
697  ASSERT(rows>=0 && cols>=0)
698  SG_SPRINT("%s%s=[\n", prefix, name)
699  for (int64_t i=0; i<rows; i++)
700  {
701  SG_SPRINT("%s[", prefix)
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? "" : ",")
706  }
707  SG_SPRINT("%s]\n", prefix)
708 }
709 
710 template <>
712  const floatmax_t* matrix, int32_t rows, int32_t cols, const char* name,
713  const char* prefix)
714 {
715  ASSERT(rows>=0 && cols>=0)
716  SG_SPRINT("%s%s=[\n", prefix, name)
717  for (int64_t i=0; i<rows; i++)
718  {
719  SG_SPRINT("%s[", prefix)
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? "" : ",")
724  }
725  SG_SPRINT("%s]\n", prefix)
726 }
727 
728 template <>
730  const complex128_t* matrix, int32_t rows, int32_t cols, const char* name,
731  const char* prefix)
732 {
733  ASSERT(rows>=0 && cols>=0)
734  SG_SPRINT("%s%s=[\n", prefix, name)
735  for (int64_t i=0; i<rows; i++)
736  {
737  SG_SPRINT("%s[", prefix)
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? "" : ",")
742  }
743  SG_SPRINT("%s]\n", prefix)
744 }
745 
746 template <>
748 {
750  return SGMatrix<char>();
751 }
752 
753 template <>
755 {
756  SGMatrix<int8_t> I(size, size);
757  for (index_t i=0; i<size; ++i)
758  {
759  for (index_t j=0; j<size; ++j)
760  I(i,j)=i==j ? scale : 0.0;
761  }
762 
763  return I;
764 }
765 
766 template <>
768 {
769  SGMatrix<uint8_t> I(size, size);
770  for (index_t i=0; i<size; ++i)
771  {
772  for (index_t j=0; j<size; ++j)
773  I(i,j)=i==j ? scale : 0.0;
774  }
775 
776  return I;
777 }
778 
779 template <>
781 {
782  SGMatrix<bool> I(size, size);
783  for (index_t i=0; i<size; ++i)
784  {
785  for (index_t j=0; j<size; ++j)
786  I(i,j)=i==j ? scale : (!scale);
787  }
788 
789  return I;
790 }
791 
792 template <>
794 {
795  SGMatrix<int16_t> I(size, size);
796  for (index_t i=0; i<size; ++i)
797  {
798  for (index_t j=0; j<size; ++j)
799  I(i,j)=i==j ? scale : 0.0;
800  }
801 
802  return I;
803 }
804 
805 template <>
807 {
808  SGMatrix<uint16_t> I(size, size);
809  for (index_t i=0; i<size; ++i)
810  {
811  for (index_t j=0; j<size; ++j)
812  I(i,j)=i==j ? scale : 0.0;
813  }
814 
815  return I;
816 }
817 
818 template <>
820 {
821  SGMatrix<int32_t> I(size, size);
822  for (index_t i=0; i<size; ++i)
823  {
824  for (index_t j=0; j<size; ++j)
825  I(i,j)=i==j ? scale : 0.0;
826  }
827 
828  return I;
829 }
830 
831 template <>
833 {
834  SGMatrix<uint32_t> I(size, size);
835  for (index_t i=0; i<size; ++i)
836  {
837  for (index_t j=0; j<size; ++j)
838  I(i,j)=i==j ? scale : 0.0;
839  }
840 
841  return I;
842 }
843 
844 template <>
846 {
847  SGMatrix<int64_t> I(size, size);
848  for (index_t i=0; i<size; ++i)
849  {
850  for (index_t j=0; j<size; ++j)
851  I(i,j)=i==j ? scale : 0.0;
852  }
853 
854  return I;
855 }
856 
857 template <>
859 {
860  SGMatrix<uint64_t> I(size, size);
861  for (index_t i=0; i<size; ++i)
862  {
863  for (index_t j=0; j<size; ++j)
864  I(i,j)=i==j ? scale : 0.0;
865  }
866 
867  return I;
868 }
869 
870 template <>
872 {
873  SGMatrix<float32_t> I(size, size);
874  for (index_t i=0; i<size; ++i)
875  {
876  for (index_t j=0; j<size; ++j)
877  I(i,j)=i==j ? scale : 0.0;
878  }
879 
880  return I;
881 }
882 
883 template <>
885 {
886  SGMatrix<float64_t> I(size, size);
887  for (index_t i=0; i<size; ++i)
888  {
889  for (index_t j=0; j<size; ++j)
890  I(i,j)=i==j ? scale : 0.0;
891  }
892 
893  return I;
894 }
895 
896 template <>
898 {
899  SGMatrix<floatmax_t> I(size, size);
900  for (index_t i=0; i<size; ++i)
901  {
902  for (index_t j=0; j<size; ++j)
903  I(i,j)=i==j ? scale : 0.0;
904  }
905 
906  return I;
907 }
908 
909 template <>
911 {
912  SGMatrix<complex128_t> I(size, size);
913  for (index_t i=0; i<size; ++i)
914  {
915  for (index_t j=0; j<size; ++j)
916  I(i,j)=i==j ? scale : complex128_t(0.0);
917  }
918 
919  return I;
920 }
921 
922 //Howto construct the pseudo inverse (from "The Matrix Cookbook")
923 //
924 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
925 //
926 //The matrix A+ known as the pseudo inverse is unique and does always exist.
927 //
928 //The pseudo inverse A+ can be constructed from the singular value
929 //decomposition A = UDV^T , by A^+ = V(D+)U^T.
930 
931 #ifdef HAVE_LAPACK
932 template <class T>
934  float64_t* matrix, int32_t rows, int32_t cols, float64_t* target)
935 {
936  if (!target)
937  target=SG_MALLOC(float64_t, rows*cols);
938 
939  char jobu='A';
940  char jobvt='A';
941  int m=rows; /* for calling external lib */
942  int n=cols; /* for calling external lib */
943  int lda=m; /* for calling external lib */
944  int ldu=m; /* for calling external lib */
945  int ldvt=n; /* for calling external lib */
946  int info=-1; /* for calling external lib */
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);
951 
952  wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
953  ASSERT(info==0)
954 
955  for (int64_t i=0; i<n; i++)
956  {
957  for (int64_t j=0; j<lsize; j++)
958  vt[i*n+j]=vt[i*n+j]/s[j];
959  }
960 
961  cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
962 
963  SG_FREE(u);
964  SG_FREE(vt);
965  SG_FREE(s);
966 
967  return target;
968 }
969 
971 template <class T>
973 {
974  REQUIRE(!matrix.on_gpu(), "Operation is not possible when data is in GPU memory.\n");
975  ASSERT(matrix.num_cols==matrix.num_rows)
976  int32_t* ipiv = SG_MALLOC(int32_t, matrix.num_cols);
977  clapack_dgetrf(CblasColMajor,matrix.num_cols,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
978  clapack_dgetri(CblasColMajor,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
979  SG_FREE(ipiv);
980 }
981 
982 template <class T>
984 {
985  REQUIRE(!matrix.on_gpu(), "Operation is not possible when data is in GPU memory.\n");
986  if (matrix.num_rows!=matrix.num_cols)
987  {
988  SG_SERROR("SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix"
989  " rows and columns are not equal!\n");
990  }
991 
992  /* use reference counting for SGVector */
993  SGVector<float64_t> result(NULL, 0, true);
994  result.vlen=matrix.num_rows;
995  result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows,
996  matrix.num_rows);
997  return result;
998 }
999 
1000 template <class T>
1001 double* SGMatrix<T>::compute_eigenvectors(double* matrix, int n, int m)
1002 {
1003  ASSERT(n == m)
1004 
1005  char V='V';
1006  char U='U';
1007  int info;
1008  int ord=n;
1009  int lda=n;
1010  double* eigenvalues=SG_CALLOC(float64_t, n+1);
1011 
1012  // lapack sym matrix eigenvalues+vectors
1013  wrap_dsyev(V, U, ord, matrix, lda,
1014  eigenvalues, &info);
1015 
1016  if (info!=0)
1017  SG_SERROR("DSYEV failed with code %d\n", info)
1018 
1019  return eigenvalues;
1020 }
1021 
1022 template <class T>
1023 void SGMatrix<T>::compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors,
1024  int n, int il, int iu)
1025 {
1026  eigenvalues = SG_MALLOC(double, n);
1027  eigenvectors = SG_MALLOC(double, (iu-il+1)*n);
1028  int status = 0;
1029  wrap_dsyevr('V','U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status);
1030  ASSERT(status==0)
1031 }
1032 
1033 #endif //HAVE_LAPACK
1034 
1035 /* Already in linalg */
1036 template <class T>
1039  bool transpose_A, bool transpose_B, float64_t scale)
1040 {
1041  REQUIRE((!A.on_gpu()) && (!B.on_gpu()),
1042  "Operation is not possible when data is in GPU memory.\n");
1043 
1044  /* these variables store size of transposed matrices*/
1045  index_t cols_A=transpose_A ? A.num_rows : A.num_cols;
1046  index_t rows_A=transpose_A ? A.num_cols : A.num_rows;
1047  index_t rows_B=transpose_B ? B.num_cols : B.num_rows;
1048  index_t cols_B=transpose_B ? B.num_rows : B.num_cols;
1049 
1050  /* do a dimension check */
1051  if (cols_A!=rows_B)
1052  {
1053  SG_SERROR("SGMatrix::matrix_multiply(): Dimension mismatch: "
1054  "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
1055  }
1056 
1057  /* allocate result matrix */
1058  SGMatrix<float64_t> C(rows_A, cols_B);
1059  C.zero();
1060 #ifdef HAVE_LAPACK
1061  /* multiply */
1062  cblas_dgemm(CblasColMajor,
1063  transpose_A ? CblasTrans : CblasNoTrans,
1064  transpose_B ? CblasTrans : CblasNoTrans,
1065  rows_A, cols_B, cols_A, scale,
1066  A.matrix, A.num_rows, B.matrix, B.num_rows,
1067  0.0, C.matrix, C.num_rows);
1068 #else
1069  /* C(i,j) = scale * \Sigma A(i,k)*B(k,j) */
1070  for (int32_t i=0; i<rows_A; i++)
1071  {
1072  for (int32_t j=0; j<cols_B; j++)
1073  {
1074  for (int32_t k=0; k<cols_A; k++)
1075  {
1076  float64_t x1=transpose_A ? A(k,i):A(i,k);
1077  float64_t x2=transpose_B ? B(j,k):B(k,j);
1078  C(i,j)+=x1*x2;
1079  }
1080 
1081  C(i,j)*=scale;
1082  }
1083  }
1084 #endif //HAVE_LAPACK
1085 
1086  return C;
1087 }
1088 
1089 template<class T>
1091  index_t num_cols, SGMatrix<T> pre_allocated)
1092 {
1093  SGMatrix<T> result;
1094 
1095  /* evtl use pre-allocated space */
1096  if (pre_allocated.matrix || pre_allocated.gpu_ptr)
1097  {
1098  result=pre_allocated;
1099 
1100  /* check dimension consistency */
1101  if (pre_allocated.num_rows!=num_rows ||
1102  pre_allocated.num_cols!=num_cols)
1103  {
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);
1108  }
1109  }
1110  else
1111  {
1112  /* otherwise, allocate space */
1113  result=SGMatrix<T>(num_rows, num_cols);
1114  }
1115 
1116  return result;
1117 }
1118 
1119 template<class T>
1121 {
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);
1128 }
1129 
1130 template<class T>
1132 {
1133  matrix=NULL;
1134  num_rows=0;
1135  num_cols=0;
1136  gpu_ptr=nullptr;
1137  m_on_gpu.store(false, std::memory_order_release);
1138 }
1139 
1140 template<class T>
1142 {
1143  SG_FREE(matrix);
1144  matrix=NULL;
1145  num_rows=0;
1146  num_cols=0;
1147 }
1148 
1149 template<class T>
1151 {
1152  ASSERT(loader)
1153  unref();
1154 
1156  SGMatrix<T> mat;
1157  loader->get_matrix(mat.matrix, mat.num_rows, mat.num_cols);
1158  mat.gpu_ptr = nullptr;
1159  copy_data(mat);
1160  copy_refcount(mat);
1161  ref();
1163 }
1164 
1165 template<>
1167 {
1168  SG_SERROR("SGMatrix::load():: Not supported for complex128_t\n");
1169 }
1170 
1171 template<class T>
1173 {
1174  assert_on_cpu();
1175  ASSERT(writer)
1177  writer->set_matrix(matrix, num_rows, num_cols);
1179 }
1180 
1181 template<>
1183 {
1184  SG_SERROR("SGMatrix::save():: Not supported for complex128_t\n");
1185 }
1186 
1187 template<class T>
1189 {
1190  assert_on_cpu();
1191  SGVector<T> rowv(num_cols);
1192  for (int64_t i = 0; i < num_cols; i++)
1193  {
1194  rowv[i] = matrix[i*num_rows+row];
1195  }
1196  return rowv;
1197 }
1198 
1199 template<class T>
1201 {
1202  assert_on_cpu();
1203  index_t diag_vlen=CMath::min(num_cols, num_rows);
1204  SGVector<T> diag(diag_vlen);
1205 
1206  for (int64_t i=0; i<diag_vlen; i++)
1207  {
1208  diag[i]=matrix[i*num_rows+i];
1209  }
1210 
1211  return diag;
1212 }
1213 
1214 template class SGMatrix<bool>;
1215 template class SGMatrix<char>;
1216 template class SGMatrix<int8_t>;
1217 template class SGMatrix<uint8_t>;
1218 template class SGMatrix<int16_t>;
1219 template class SGMatrix<uint16_t>;
1220 template class SGMatrix<int32_t>;
1221 template class SGMatrix<uint32_t>;
1222 template class SGMatrix<int64_t>;
1223 template class SGMatrix<uint64_t>;
1224 template class SGMatrix<float32_t>;
1225 template class SGMatrix<float64_t>;
1226 template class SGMatrix<floatmax_t>;
1227 template class SGMatrix<complex128_t>;
1228 }
void display_matrix(const char *name="matrix") const
Definition: SGMatrix.cpp:480
static void center_matrix(T *matrix, int32_t m, int32_t n)
Definition: SGMatrix.cpp:436
#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:307
#define SG_RESET_LOCALE
Definition: SGIO.h:85
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 > clone() const
Definition: SGMatrix.cpp:330
int32_t index_t
Definition: common.h:72
static T * get_column_sum(T *matrix, int32_t m, int32_t n)
Definition: SGMatrix.cpp:416
void scale(SGVector< T > &a, SGVector< T > &result, T alpha=1)
static void create_diagonal_matrix(T *matrix, T *v, int32_t size)
Definition: SGMatrix.cpp:375
Definition: basetag.h:132
#define REQUIRE(x,...)
Definition: SGIO.h:205
index_t num_cols
Definition: SGMatrix.h:465
#define SG_SNOTIMPLEMENTED
Definition: SGIO.h:197
void set_const(Container< T > &a, T value)
virtual ~SGMatrix()
Definition: SGMatrix.cpp:138
#define SG_SET_LOCALE_C
Definition: SGIO.h:84
index_t num_rows
Definition: SGMatrix.h:463
shogun matrix
static void transpose_matrix(T *&matrix, int32_t &num_feat, int32_t &num_vec)
Definition: SGMatrix.cpp:358
static T * clone_matrix(const T *matrix, int32_t nrows, int32_t ncols)
Definition: SGMatrix.cpp:345
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:253
index_t vlen
Definition: SGVector.h:545
#define SG_SPRINT(...)
Definition: SGIO.h:179
bool on_gpu() const
Definition: SGMatrix.h:79
#define ASSERT(x)
Definition: SGIO.h:200
shogun vector
void remove_column_mean()
Definition: SGMatrix.cpp:461
virtual void get_matrix(bool *&matrix, int32_t &num_feat, int32_t &num_vec)
Definition: File.cpp:109
shogun reference count managed data
double float64_t
Definition: common.h:60
long double floatmax_t
Definition: common.h:61
A File access base class.
Definition: File.h:34
T max_single() const
Definition: SGMatrix.cpp:311
bool equals(const SGMatrix< T > &other) const
Definition: SGMatrix.cpp:144
static T sum(T *vec, int32_t len)
Return sum(vec)
Definition: SGVector.h:392
bool on_gpu() const
Definition: SGVector.h:83
void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info)
Definition: lapack.cpp:235
static float64_t trace(float64_t *mat, int32_t cols, int32_t rows)
Definition: SGMatrix.cpp:391
virtual void copy_data(const SGReferencedData &orig)
Definition: SGMatrix.cpp:1120
SGMatrix< T > & operator=(const SGMatrix< T > &)
Definition: SGMatrix.cpp:125
float float32_t
Definition: common.h:59
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
Interface for GPU memory libraries.
Definition: SGMatrix.h:34
std::shared_ptr< GPUMemoryBase< T > > gpu_ptr
Definition: SGVector.h:547
std::shared_ptr< GPUMemoryBase< T > > gpu_ptr
Definition: SGMatrix.h:467
#define SG_SERROR(...)
Definition: SGIO.h:178
static T * get_row_sum(T *matrix, int32_t m, int32_t n)
Definition: SGMatrix.cpp:402
bool is_symmetric() const
Definition: SGMatrix.cpp:227
#define REAL_IS_SYMMETRIC(real_t)
Definition: SGMatrix.cpp:251
static void swap(T &a, T &b)
Definition: Math.h:433
void set_const(T const_elem)
Definition: SGMatrix.cpp:209
virtual void init_data()
Definition: SGMatrix.cpp:1131

SHOGUN Machine Learning Toolbox - Documentation