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

SHOGUN Machine Learning Toolbox - Documentation