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

SHOGUN Machine Learning Toolbox - Documentation