SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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) 2011-2013 Heiko Strathmann
8  * Written (W) 2012 Fernando Jose Iglesias Garcia
9  * Written (W) 2010,2012 Soeren Sonnenburg
10  * Copyright (C) 2010 Berlin Institute of Technology
11  * Copyright (C) 2012 Soeren Sonnenburg
12  */
13 
14 #include <shogun/lib/config.h>
15 #include <shogun/io/SGIO.h>
16 #include <shogun/io/File.h>
17 #include <shogun/lib/SGMatrix.h>
18 #include <shogun/lib/SGVector.h>
21 #include <limits>
22 
23 #ifdef HAVE_EIGEN3
25 #endif
26 
27 namespace shogun {
28 
29 template <class T>
31 {
32  init_data();
33 }
34 
35 template <class T>
36 SGMatrix<T>::SGMatrix(bool ref_counting) : SGReferencedData(ref_counting)
37 {
38  init_data();
39 }
40 
41 template <class T>
42 SGMatrix<T>::SGMatrix(T* m, index_t nrows, index_t ncols, bool ref_counting)
43  : SGReferencedData(ref_counting), matrix(m),
44  num_rows(nrows), num_cols(ncols) { }
45 
46 template <class T>
47 SGMatrix<T>::SGMatrix(index_t nrows, index_t ncols, bool ref_counting)
48  : SGReferencedData(ref_counting), num_rows(nrows), num_cols(ncols)
49 {
50  matrix=SG_MALLOC(T, ((int64_t) nrows)*ncols);
51  zero();
52 }
53 
54 template <class T>
56 {
57  copy_data(orig);
58 }
59 
60 #ifdef HAVE_EIGEN3
61 template <class T>
62 SGMatrix<T>::SGMatrix(EigenMatrixXt& mat)
63 : SGReferencedData(false), matrix(mat.data()),
64  num_rows(mat.rows()), num_cols(mat.cols())
65 {
66 
67 }
68 
69 template <class T>
70 SGMatrix<T>::operator EigenMatrixXtMap() const
71 {
72  return EigenMatrixXtMap(matrix, num_rows, num_cols);
73 }
74 #endif
75 
76 template <class T>
78 {
79  unref();
80 }
81 
82 template <class T>
84 {
85  if (num_rows!=other.num_rows || num_cols!=other.num_cols)
86  return false;
87 
88  if (matrix!=other.matrix)
89  return false;
90 
91  return true;
92 }
93 
94 template <class T>
96 {
97  if (num_rows!=other.num_rows || num_cols!=other.num_cols)
98  return false;
99 
100  for (int64_t i=0; i<int64_t(num_rows)*num_cols; ++i)
101  {
102  if (matrix[i]!=other.matrix[i])
103  return false;
104  }
105 
106  return true;
107 }
108 
109 template <class T>
110 void SGMatrix<T>::set_const(T const_elem)
111 {
112  for (int64_t i=0; i<int64_t(num_rows)*num_cols; i++)
113  matrix[i]=const_elem ;
114 }
115 
116 template <class T>
118 {
119  if (matrix && (int64_t(num_rows)*num_cols))
120  set_const(0);
121 }
122 
123 template <>
125 {
126  if (matrix && (int64_t(num_rows)*num_cols))
127  set_const(complex128_t(0.0));
128 }
129 
130 template <class T>
132 {
133  if (num_rows!=num_cols)
134  return false;
135  for (int i=0; i<num_rows; ++i)
136  {
137  for (int j=i+1; j<num_cols; ++j)
138  {
139  if (matrix[j*num_rows+i]!=matrix[i*num_rows+j])
140  return false;
141  }
142  }
143  return true;
144 }
145 
146 template <>
148 {
149  if (num_rows!=num_cols)
150  return false;
151  for (int i=0; i<num_rows; ++i)
152  {
153  for (int j=i+1; j<num_cols; ++j)
154  {
155  if (!CMath::fequals<float32_t>(matrix[j*num_rows+i],
156  matrix[i*num_rows+j], FLT_EPSILON))
157  return false;
158  }
159  }
160  return true;
161 }
162 
163 template <>
165 {
166  if (num_rows!=num_cols)
167  return false;
168  for (int i=0; i<num_rows; ++i)
169  {
170  for (int j=i+1; j<num_cols; ++j)
171  {
172  if (!CMath::fequals<float64_t>(matrix[j*num_rows+i],
173  matrix[i*num_rows+j], DBL_EPSILON))
174  return false;
175  }
176  }
177  return true;
178 }
179 
180 template <>
182 {
183  if (num_rows!=num_cols)
184  return false;
185  for (int i=0; i<num_rows; ++i)
186  {
187  for (int j=i+1; j<num_cols; ++j)
188  {
189  if (!CMath::fequals<floatmax_t>(matrix[j*num_rows+i],
190  matrix[i*num_rows+j], LDBL_EPSILON))
191  return false;
192  }
193  }
194  return true;
195 }
196 
197 template <>
199 {
200  if (num_rows!=num_cols)
201  return false;
202  for (int i=0; i<num_rows; ++i)
203  {
204  for (int j=i+1; j<num_cols; ++j)
205  {
206  if (!(CMath::fequals<float64_t>(matrix[j*num_rows+i].real(),
207  matrix[i*num_rows+j].real(), DBL_EPSILON) &&
208  CMath::fequals<float64_t>(matrix[j*num_rows+i].imag(),
209  matrix[i*num_rows+j].imag(), DBL_EPSILON)))
210  return false;
211  }
212  }
213  return true;
214 }
215 
216 template <class T>
218 {
219  T max=matrix[0];
220  for (int64_t i=1; i<int64_t(num_rows)*num_cols; ++i)
221  {
222  if (matrix[i]>max)
223  max=matrix[i];
224  }
225 
226  return max;
227 }
228 
229 template <>
231 {
232  SG_SERROR("SGMatrix::max_single():: Not supported for complex128_t\n");
233  return complex128_t(0.0);
234 }
235 
236 template <class T>
238 {
239  return SGMatrix<T>(clone_matrix(matrix, num_rows, num_cols),
240  num_rows, num_cols);
241 }
242 
243 template <class T>
244 T* SGMatrix<T>::clone_matrix(const T* matrix, int32_t nrows, int32_t ncols)
245 {
246  T* result = SG_MALLOC(T, int64_t(nrows)*ncols);
247  for (int64_t i=0; i<int64_t(nrows)*ncols; i++)
248  result[i]=matrix[i];
249 
250  return result;
251 }
252 
253 template <class T>
255  T*& matrix, int32_t& num_feat, int32_t& num_vec)
256 {
257  /* this should be done in-place! Heiko */
258  T* transposed=SG_MALLOC(T, int64_t(num_vec)*num_feat);
259  for (int64_t i=0; i<num_vec; i++)
260  {
261  for (int64_t j=0; j<num_feat; j++)
262  transposed[i+j*num_vec]=matrix[i*num_feat+j];
263  }
264 
265  SG_FREE(matrix);
266  matrix=transposed;
267 
268  CMath::swap(num_feat, num_vec);
269 }
270 
271 template <class T>
272 void SGMatrix<T>::create_diagonal_matrix(T* matrix, T* v,int32_t size)
273 {
274  for(int64_t i=0;i<size;i++)
275  {
276  for(int64_t j=0;j<size;j++)
277  {
278  if(i==j)
279  matrix[j*size+i]=v[i];
280  else
281  matrix[j*size+i]=0;
282  }
283  }
284 }
285 
286 template <class T>
288  float64_t* mat, int32_t cols, int32_t rows)
289 {
290  float64_t trace=0;
291  for (int64_t i=0; i<rows; i++)
292  trace+=mat[i*cols+i];
293  return trace;
294 }
295 
296 template <class T>
297 T* SGMatrix<T>::get_row_sum(T* matrix, int32_t m, int32_t n)
298 {
299  T* rowsums=SG_CALLOC(T, n);
300 
301  for (int64_t i=0; i<n; i++)
302  {
303  for (int64_t j=0; j<m; j++)
304  rowsums[i]+=matrix[j+i*m];
305  }
306  return rowsums;
307 }
308 
309 template <class T>
310 T* SGMatrix<T>::get_column_sum(T* matrix, int32_t m, int32_t n)
311 {
312  T* colsums=SG_CALLOC(T, m);
313 
314  for (int64_t i=0; i<n; i++)
315  {
316  for (int64_t j=0; j<m; j++)
317  colsums[j]+=matrix[j+i*m];
318  }
319  return colsums;
320 }
321 
322 template <class T>
324 {
325  center_matrix(matrix, num_rows, num_cols);
326 }
327 
328 template <class T>
329 void SGMatrix<T>::center_matrix(T* matrix, int32_t m, int32_t n)
330 {
331  float64_t num_data=n;
332 
333  T* colsums=get_column_sum(matrix, m,n);
334  T* rowsums=get_row_sum(matrix, m,n);
335 
336  for (int32_t i=0; i<m; i++)
337  colsums[i]/=num_data;
338  for (int32_t j=0; j<n; j++)
339  rowsums[j]/=num_data;
340 
341  T s=SGVector<T>::sum(rowsums, n)/num_data;
342 
343  for (int64_t i=0; i<n; i++)
344  {
345  for (int64_t j=0; j<m; j++)
346  matrix[i*m+j]+=s-colsums[j]-rowsums[i];
347  }
348 
349  SG_FREE(rowsums);
350  SG_FREE(colsums);
351 }
352 
353 template <class T>
355 {
356  /* compute "row" sums (which I would call column sums), i.e. sum of all
357  * elements in a fixed column */
358  T* means=get_row_sum(matrix, num_rows, num_cols);
359 
360  /* substract column mean from every corresponding entry */
361  for (int64_t i=0; i<num_cols; ++i)
362  {
363  means[i]/=num_rows;
364  for (int64_t j=0; j<num_rows; ++j)
365  matrix[i*num_rows+j]-=means[i];
366  }
367 
368  SG_FREE(means);
369 }
370 
371 template<class T> void SGMatrix<T>::display_matrix(const char* name) const
372 {
373  display_matrix(matrix, num_rows, num_cols, name);
374 }
375 
376 template <class T>
378  const SGMatrix<T> matrix, const char* name,
379  const char* prefix)
380 {
381  matrix.display_matrix();
382 }
383 
384 template <>
386  const bool* matrix, int32_t rows, int32_t cols, const char* name,
387  const char* prefix)
388 {
389  ASSERT(rows>=0 && cols>=0)
390  SG_SPRINT("%s%s=[\n", prefix, name)
391  for (int64_t i=0; i<rows; i++)
392  {
393  SG_SPRINT("%s[", prefix)
394  for (int64_t j=0; j<cols; j++)
395  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i] ? 1 : 0,
396  j==cols-1? "" : ",");
397  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
398  }
399  SG_SPRINT("%s]\n", prefix)
400 }
401 
402 template <>
404  const char* matrix, int32_t rows, int32_t cols, const char* name,
405  const char* prefix)
406 {
407  ASSERT(rows>=0 && cols>=0)
408  SG_SPRINT("%s%s=[\n", prefix, name)
409  for (int64_t i=0; i<rows; i++)
410  {
411  SG_SPRINT("%s[", prefix)
412  for (int64_t j=0; j<cols; j++)
413  SG_SPRINT("%s\t%c%s", prefix, matrix[j*rows+i],
414  j==cols-1? "" : ",");
415  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
416  }
417  SG_SPRINT("%s]\n", prefix)
418 }
419 
420 template <>
422  const int8_t* matrix, int32_t rows, int32_t cols, const char* name,
423  const char* prefix)
424 {
425  ASSERT(rows>=0 && cols>=0)
426  SG_SPRINT("%s%s=[\n", prefix, name)
427  for (int64_t i=0; i<rows; i++)
428  {
429  SG_SPRINT("%s[", prefix)
430  for (int64_t j=0; j<cols; j++)
431  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
432  j==cols-1? "" : ",");
433  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
434  }
435  SG_SPRINT("%s]\n", prefix)
436 }
437 
438 template <>
440  const uint8_t* matrix, int32_t rows, int32_t cols, const char* name,
441  const char* prefix)
442 {
443  ASSERT(rows>=0 && cols>=0)
444  SG_SPRINT("%s%s=[\n", prefix, name)
445  for (int64_t i=0; i<rows; i++)
446  {
447  SG_SPRINT("%s[", prefix)
448  for (int64_t j=0; j<cols; j++)
449  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
450  j==cols-1? "" : ",");
451  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
452  }
453  SG_SPRINT("%s]\n", prefix)
454 }
455 
456 template <>
458  const int16_t* matrix, int32_t rows, int32_t cols, const char* name,
459  const char* prefix)
460 {
461  ASSERT(rows>=0 && cols>=0)
462  SG_SPRINT("%s%s=[\n", prefix, name)
463  for (int64_t i=0; i<rows; i++)
464  {
465  SG_SPRINT("%s[", prefix)
466  for (int64_t j=0; j<cols; j++)
467  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
468  j==cols-1? "" : ",");
469  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
470  }
471  SG_SPRINT("%s]\n", prefix)
472 }
473 
474 template <>
476  const uint16_t* matrix, int32_t rows, int32_t cols, const char* name,
477  const char* prefix)
478 {
479  ASSERT(rows>=0 && cols>=0)
480  SG_SPRINT("%s%s=[\n", prefix, name)
481  for (int64_t i=0; i<rows; i++)
482  {
483  SG_SPRINT("%s[", prefix)
484  for (int64_t j=0; j<cols; j++)
485  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
486  j==cols-1? "" : ",");
487  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
488  }
489  SG_SPRINT("%s]\n", prefix)
490 }
491 
492 
493 template <>
495  const int32_t* matrix, int32_t rows, int32_t cols, const char* name,
496  const char* prefix)
497 {
498  ASSERT(rows>=0 && cols>=0)
499  SG_SPRINT("%s%s=[\n", prefix, name)
500  for (int64_t i=0; i<rows; i++)
501  {
502  SG_SPRINT("%s[", prefix)
503  for (int64_t j=0; j<cols; j++)
504  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
505  j==cols-1? "" : ",");
506  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
507  }
508  SG_SPRINT("%s]\n", prefix)
509 }
510 
511 template <>
513  const uint32_t* matrix, int32_t rows, int32_t cols, const char* name,
514  const char* prefix)
515 {
516  ASSERT(rows>=0 && cols>=0)
517  SG_SPRINT("%s%s=[\n", prefix, name)
518  for (int64_t i=0; i<rows; i++)
519  {
520  SG_SPRINT("%s[", prefix)
521  for (int64_t j=0; j<cols; j++)
522  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
523  j==cols-1? "" : ",");
524  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
525  }
526  SG_SPRINT("%s]\n", prefix)
527 }
528 template <>
530  const int64_t* matrix, int32_t rows, int32_t cols, const char* name,
531  const char* prefix)
532 {
533  ASSERT(rows>=0 && cols>=0)
534  SG_SPRINT("%s%s=[\n", prefix, name)
535  for (int64_t i=0; i<rows; i++)
536  {
537  SG_SPRINT("%s[", prefix)
538  for (int64_t j=0; j<cols; j++)
539  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
540  j==cols-1? "" : ",");
541  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
542  }
543  SG_SPRINT("%s]\n", prefix)
544 }
545 
546 template <>
548  const uint64_t* matrix, int32_t rows, int32_t cols, const char* name,
549  const char* prefix)
550 {
551  ASSERT(rows>=0 && cols>=0)
552  SG_SPRINT("%s%s=[\n", prefix, name)
553  for (int64_t i=0; i<rows; i++)
554  {
555  SG_SPRINT("%s[", prefix)
556  for (int64_t j=0; j<cols; j++)
557  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
558  j==cols-1? "" : ",");
559  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
560  }
561  SG_SPRINT("%s]\n", prefix)
562 }
563 
564 template <>
566  const float32_t* matrix, int32_t rows, int32_t cols, const char* name,
567  const char* prefix)
568 {
569  ASSERT(rows>=0 && cols>=0)
570  SG_SPRINT("%s%s=[\n", prefix, name)
571  for (int64_t i=0; i<rows; i++)
572  {
573  SG_SPRINT("%s[", prefix)
574  for (int64_t j=0; j<cols; j++)
575  SG_SPRINT("%s\t%.18g%s", prefix, (float) matrix[j*rows+i],
576  j==cols-1? "" : ",");
577  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
578  }
579  SG_SPRINT("%s]\n", prefix)
580 }
581 
582 template <>
584  const float64_t* matrix, int32_t rows, int32_t cols, const char* name,
585  const char* prefix)
586 {
587  ASSERT(rows>=0 && cols>=0)
588  SG_SPRINT("%s%s=[\n", prefix, name)
589  for (int64_t i=0; i<rows; i++)
590  {
591  SG_SPRINT("%s[", prefix)
592  for (int64_t j=0; j<cols; j++)
593  SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i],
594  j==cols-1? "" : ",");
595  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
596  }
597  SG_SPRINT("%s]\n", prefix)
598 }
599 
600 template <>
602  const floatmax_t* matrix, int32_t rows, int32_t cols, const char* name,
603  const char* prefix)
604 {
605  ASSERT(rows>=0 && cols>=0)
606  SG_SPRINT("%s%s=[\n", prefix, name)
607  for (int64_t i=0; i<rows; i++)
608  {
609  SG_SPRINT("%s[", prefix)
610  for (int64_t j=0; j<cols; j++)
611  SG_SPRINT("%s\t%.18g%s", prefix, (double) matrix[j*rows+i],
612  j==cols-1? "" : ",");
613  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
614  }
615  SG_SPRINT("%s]\n", prefix)
616 }
617 
618 template <>
620  const complex128_t* matrix, int32_t rows, int32_t cols, const char* name,
621  const char* prefix)
622 {
623  ASSERT(rows>=0 && cols>=0)
624  SG_SPRINT("%s%s=[\n", prefix, name)
625  for (int64_t i=0; i<rows; i++)
626  {
627  SG_SPRINT("%s[", prefix)
628  for (int64_t j=0; j<cols; j++)
629  SG_SPRINT("%s\t(%.18g+i%.18g)%s", prefix, matrix[j*rows+i].real(),
630  matrix[j*rows+i].imag(), j==cols-1? "" : ",");
631  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
632  }
633  SG_SPRINT("%s]\n", prefix)
634 }
635 
636 template <>
638 {
640  return SGMatrix<char>();
641 }
642 
643 template <>
645 {
646  SGMatrix<int8_t> I(size, size);
647  for (index_t i=0; i<size; ++i)
648  {
649  for (index_t j=0; j<size; ++j)
650  I(i,j)=i==j ? scale : 0.0;
651  }
652 
653  return I;
654 }
655 
656 template <>
658 {
659  SGMatrix<uint8_t> I(size, size);
660  for (index_t i=0; i<size; ++i)
661  {
662  for (index_t j=0; j<size; ++j)
663  I(i,j)=i==j ? scale : 0.0;
664  }
665 
666  return I;
667 }
668 
669 template <>
671 {
672  SGMatrix<bool> I(size, size);
673  for (index_t i=0; i<size; ++i)
674  {
675  for (index_t j=0; j<size; ++j)
676  I(i,j)=i==j ? scale : (!scale);
677  }
678 
679  return I;
680 }
681 
682 template <>
684 {
685  SGMatrix<int16_t> I(size, size);
686  for (index_t i=0; i<size; ++i)
687  {
688  for (index_t j=0; j<size; ++j)
689  I(i,j)=i==j ? scale : 0.0;
690  }
691 
692  return I;
693 }
694 
695 template <>
697 {
698  SGMatrix<uint16_t> I(size, size);
699  for (index_t i=0; i<size; ++i)
700  {
701  for (index_t j=0; j<size; ++j)
702  I(i,j)=i==j ? scale : 0.0;
703  }
704 
705  return I;
706 }
707 
708 template <>
710 {
711  SGMatrix<int32_t> I(size, size);
712  for (index_t i=0; i<size; ++i)
713  {
714  for (index_t j=0; j<size; ++j)
715  I(i,j)=i==j ? scale : 0.0;
716  }
717 
718  return I;
719 }
720 
721 template <>
723 {
724  SGMatrix<uint32_t> I(size, size);
725  for (index_t i=0; i<size; ++i)
726  {
727  for (index_t j=0; j<size; ++j)
728  I(i,j)=i==j ? scale : 0.0;
729  }
730 
731  return I;
732 }
733 
734 template <>
736 {
737  SGMatrix<int64_t> I(size, size);
738  for (index_t i=0; i<size; ++i)
739  {
740  for (index_t j=0; j<size; ++j)
741  I(i,j)=i==j ? scale : 0.0;
742  }
743 
744  return I;
745 }
746 
747 template <>
749 {
750  SGMatrix<uint64_t> I(size, size);
751  for (index_t i=0; i<size; ++i)
752  {
753  for (index_t j=0; j<size; ++j)
754  I(i,j)=i==j ? scale : 0.0;
755  }
756 
757  return I;
758 }
759 
760 template <>
762 {
763  SGMatrix<float32_t> I(size, size);
764  for (index_t i=0; i<size; ++i)
765  {
766  for (index_t j=0; j<size; ++j)
767  I(i,j)=i==j ? scale : 0.0;
768  }
769 
770  return I;
771 }
772 
773 template <>
775 {
776  SGMatrix<float64_t> I(size, size);
777  for (index_t i=0; i<size; ++i)
778  {
779  for (index_t j=0; j<size; ++j)
780  I(i,j)=i==j ? scale : 0.0;
781  }
782 
783  return I;
784 }
785 
786 template <>
788 {
789  SGMatrix<floatmax_t> I(size, size);
790  for (index_t i=0; i<size; ++i)
791  {
792  for (index_t j=0; j<size; ++j)
793  I(i,j)=i==j ? scale : 0.0;
794  }
795 
796  return I;
797 }
798 
799 template <>
801 {
802  SGMatrix<complex128_t> I(size, size);
803  for (index_t i=0; i<size; ++i)
804  {
805  for (index_t j=0; j<size; ++j)
806  I(i,j)=i==j ? scale : complex128_t(0.0);
807  }
808 
809  return I;
810 }
811 
812 //Howto construct the pseudo inverse (from "The Matrix Cookbook")
813 //
814 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
815 //
816 //The matrix A+ known as the pseudo inverse is unique and does always exist.
817 //
818 //The pseudo inverse A+ can be constructed from the singular value
819 //decomposition A = UDV^T , by A^+ = V(D+)U^T.
820 
821 #ifdef HAVE_LAPACK
822 template <class T>
824  float64_t* matrix, int32_t rows, int32_t cols, float64_t* target)
825 {
826  if (!target)
827  target=SG_MALLOC(float64_t, rows*cols);
828 
829  char jobu='A';
830  char jobvt='A';
831  int m=rows; /* for calling external lib */
832  int n=cols; /* for calling external lib */
833  int lda=m; /* for calling external lib */
834  int ldu=m; /* for calling external lib */
835  int ldvt=n; /* for calling external lib */
836  int info=-1; /* for calling external lib */
837  int32_t lsize=CMath::min((int32_t) m, (int32_t) n);
838  double* s=SG_MALLOC(double, lsize);
839  double* u=SG_MALLOC(double, m*m);
840  double* vt=SG_MALLOC(double, n*n);
841 
842  wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
843  ASSERT(info==0)
844 
845  for (int64_t i=0; i<n; i++)
846  {
847  for (int64_t j=0; j<lsize; j++)
848  vt[i*n+j]=vt[i*n+j]/s[j];
849  }
850 
851  cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
852 
853  SG_FREE(u);
854  SG_FREE(vt);
855  SG_FREE(s);
856 
857  return target;
858 }
859 
861 template <class T>
863 {
864  ASSERT(matrix.num_cols==matrix.num_rows)
865  int32_t* ipiv = SG_MALLOC(int32_t, matrix.num_cols);
866  clapack_dgetrf(CblasColMajor,matrix.num_cols,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
867  clapack_dgetri(CblasColMajor,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
868  SG_FREE(ipiv);
869 }
870 
871 template <class T>
873 {
874  if (matrix.num_rows!=matrix.num_cols)
875  {
876  SG_SERROR("SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix"
877  " rows and columns are not equal!\n");
878  }
879 
880  /* use reference counting for SGVector */
881  SGVector<float64_t> result(NULL, 0, true);
882  result.vlen=matrix.num_rows;
883  result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows,
884  matrix.num_rows);
885  return result;
886 }
887 
888 template <class T>
889 double* SGMatrix<T>::compute_eigenvectors(double* matrix, int n, int m)
890 {
891  ASSERT(n == m)
892 
893  char V='V';
894  char U='U';
895  int info;
896  int ord=n;
897  int lda=n;
898  double* eigenvalues=SG_CALLOC(float64_t, n+1);
899 
900  // lapack sym matrix eigenvalues+vectors
901  wrap_dsyev(V, U, ord, matrix, lda,
902  eigenvalues, &info);
903 
904  if (info!=0)
905  SG_SERROR("DSYEV failed with code %d\n", info)
906 
907  return eigenvalues;
908 }
909 
910 template <class T>
911 void SGMatrix<T>::compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors,
912  int n, int il, int iu)
913 {
914  eigenvalues = SG_MALLOC(double, n);
915  eigenvectors = SG_MALLOC(double, (iu-il+1)*n);
916  int status = 0;
917  wrap_dsyevr('V','U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status);
918  ASSERT(status==0)
919 }
920 
921 #endif //HAVE_LAPACK
922 
923 template <class T>
926  bool transpose_A, bool transpose_B, float64_t scale)
927 {
928  /* these variables store size of transposed matrices*/
929  index_t cols_A=transpose_A ? A.num_rows : A.num_cols;
930  index_t rows_A=transpose_A ? A.num_cols : A.num_rows;
931  index_t rows_B=transpose_B ? B.num_cols : B.num_rows;
932  index_t cols_B=transpose_B ? B.num_rows : B.num_cols;
933 
934  /* do a dimension check */
935  if (cols_A!=rows_B)
936  {
937  SG_SERROR("SGMatrix::matrix_multiply(): Dimension mismatch: "
938  "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
939  }
940 
941  /* allocate result matrix */
942  SGMatrix<float64_t> C(rows_A, cols_B);
943  C.zero();
944 #ifdef HAVE_LAPACK
945  /* multiply */
946  cblas_dgemm(CblasColMajor,
947  transpose_A ? CblasTrans : CblasNoTrans,
948  transpose_B ? CblasTrans : CblasNoTrans,
949  rows_A, cols_B, cols_A, scale,
950  A.matrix, A.num_rows, B.matrix, B.num_rows,
951  0.0, C.matrix, C.num_rows);
952 #else
953  /* C(i,j) = scale * \Sigma A(i,k)*B(k,j) */
954  for (int32_t i=0; i<rows_A; i++)
955  {
956  for (int32_t j=0; j<cols_B; j++)
957  {
958  for (int32_t k=0; k<cols_A; k++)
959  {
960  float64_t x1=transpose_A ? A(k,i):A(i,k);
961  float64_t x2=transpose_B ? B(j,k):B(k,j);
962  C(i,j)+=x1*x2;
963  }
964 
965  C(i,j)*=scale;
966  }
967  }
968 #endif //HAVE_LAPACK
969 
970  return C;
971 }
972 
973 template<class T>
975  index_t num_cols, SGMatrix<T> pre_allocated)
976 {
977  SGMatrix<T> result;
978 
979  /* evtl use pre-allocated space */
980  if (pre_allocated.matrix)
981  {
982  result=pre_allocated;
983 
984  /* check dimension consistency */
985  if (pre_allocated.num_rows!=num_rows ||
986  pre_allocated.num_cols!=num_cols)
987  {
988  SG_SERROR("SGMatrix<T>::get_allocated_matrix(). Provided target"
989  "matrix dimensions (%dx%d) do not match passed data "
990  "dimensions (%dx%d)!\n", pre_allocated.num_rows,
991  pre_allocated.num_cols, num_rows, num_cols);
992  }
993  }
994  else
995  {
996  /* otherwise, allocate space */
997  result=SGMatrix<T>(num_rows, num_cols);
998  }
999 
1000  return result;
1001 }
1002 
1003 template<class T>
1005 {
1006  matrix=((SGMatrix*)(&orig))->matrix;
1007  num_rows=((SGMatrix*)(&orig))->num_rows;
1008  num_cols=((SGMatrix*)(&orig))->num_cols;
1009 }
1010 
1011 template<class T>
1013 {
1014  matrix=NULL;
1015  num_rows=0;
1016  num_cols=0;
1017 }
1018 
1019 template<class T>
1021 {
1022  SG_FREE(matrix);
1023  matrix=NULL;
1024  num_rows=0;
1025  num_cols=0;
1026 }
1027 
1028 template<class T>
1030 {
1031  ASSERT(loader)
1032  unref();
1033 
1035  SGMatrix<T> mat;
1036  loader->get_matrix(mat.matrix, mat.num_rows, mat.num_cols);
1037  copy_data(mat);
1038  copy_refcount(mat);
1039  ref();
1041 }
1042 
1043 template<>
1045 {
1046  SG_SERROR("SGMatrix::load():: Not supported for complex128_t\n");
1047 }
1048 
1049 template<class T>
1051 {
1052  ASSERT(writer)
1054  writer->set_matrix(matrix, num_rows, num_cols);
1056 }
1057 
1058 template<>
1060 {
1061  SG_SERROR("SGMatrix::save():: Not supported for complex128_t\n");
1062 }
1063 
1064 template<class T>
1066 {
1067  SGVector<T> rowv(num_cols);
1068  for (int64_t i = 0; i < num_cols; i++)
1069  {
1070  rowv[i] = matrix[i*num_rows+row];
1071  }
1072  return rowv;
1073 }
1074 
1075 template<class T>
1077 {
1078  index_t diag_vlen=CMath::min(num_cols, num_rows);
1079  SGVector<T> diag(diag_vlen);
1080 
1081  for (int64_t i=0; i<diag_vlen; i++)
1082  {
1083  diag[i]=matrix[i*num_rows+i];
1084  }
1085 
1086  return diag;
1087 }
1088 
1089 template class SGMatrix<bool>;
1090 template class SGMatrix<char>;
1091 template class SGMatrix<int8_t>;
1092 template class SGMatrix<uint8_t>;
1093 template class SGMatrix<int16_t>;
1094 template class SGMatrix<uint16_t>;
1095 template class SGMatrix<int32_t>;
1096 template class SGMatrix<uint32_t>;
1097 template class SGMatrix<int64_t>;
1098 template class SGMatrix<uint64_t>;
1099 template class SGMatrix<float32_t>;
1100 template class SGMatrix<float64_t>;
1101 template class SGMatrix<floatmax_t>;
1102 template class SGMatrix<complex128_t>;
1103 }

SHOGUN Machine Learning Toolbox - Documentation