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 
22 namespace shogun {
23 
24 template <class T>
26 {
27  init_data();
28 }
29 
30 template <class T>
31 SGMatrix<T>::SGMatrix(bool ref_counting) : SGReferencedData(ref_counting)
32 {
33  init_data();
34 }
35 
36 template <class T>
37 SGMatrix<T>::SGMatrix(T* m, index_t nrows, index_t ncols, bool ref_counting)
38  : SGReferencedData(ref_counting), matrix(m),
39  num_rows(nrows), num_cols(ncols) { }
40 
41 template <class T>
42 SGMatrix<T>::SGMatrix(index_t nrows, index_t ncols, bool ref_counting)
43  : SGReferencedData(ref_counting), num_rows(nrows), num_cols(ncols)
44 {
45  matrix=SG_MALLOC(T, ((int64_t) nrows)*ncols);
46 }
47 
48 template <class T>
50 {
51  copy_data(orig);
52 }
53 
54 template <class T>
56 {
57  unref();
58 }
59 
60 template <class T>
62 {
63  if (num_rows!=other.num_rows || num_cols!=other.num_cols)
64  return false;
65 
66  if (matrix!=other.matrix)
67  return false;
68 
69  return true;
70 }
71 
72 template <class T>
74 {
75  if (num_rows!=other.num_rows || num_cols!=other.num_cols)
76  return false;
77 
78  for (int64_t i=0; i<int64_t(num_rows)*num_cols; ++i)
79  {
80  if (matrix[i]!=other.matrix[i])
81  return false;
82  }
83 
84  return true;
85 }
86 
87 template <class T>
88 void SGMatrix<T>::set_const(T const_elem)
89 {
90  for (int64_t i=0; i<int64_t(num_rows)*num_cols; i++)
91  matrix[i]=const_elem ;
92 }
93 
94 template <class T>
96 {
97  if (matrix && (int64_t(num_rows)*num_cols))
98  set_const(0);
99 }
100 
101 template <>
103 {
104  if (matrix && (int64_t(num_rows)*num_cols))
105  set_const(complex128_t(0.0));
106 }
107 
108 template <class T>
110 {
111  T max=matrix[0];
112  for (int64_t i=1; i<int64_t(num_rows)*num_cols; ++i)
113  {
114  if (matrix[i]>max)
115  max=matrix[i];
116  }
117 
118  return max;
119 }
120 
121 template <>
123 {
124  SG_SERROR("SGMatrix::max_single():: Not supported for complex128_t\n");
125  return complex128_t(0.0);
126 }
127 
128 template <class T>
130 {
131  return SGMatrix<T>(clone_matrix(matrix, num_rows, num_cols),
132  num_rows, num_cols);
133 }
134 
135 template <class T>
136 T* SGMatrix<T>::clone_matrix(const T* matrix, int32_t nrows, int32_t ncols)
137 {
138  T* result = SG_MALLOC(T, int64_t(nrows)*ncols);
139  for (int64_t i=0; i<int64_t(nrows)*ncols; i++)
140  result[i]=matrix[i];
141 
142  return result;
143 }
144 
145 template <class T>
147  T*& matrix, int32_t& num_feat, int32_t& num_vec)
148 {
149  /* this should be done in-place! Heiko */
150  T* transposed=SG_MALLOC(T, int64_t(num_vec)*num_feat);
151  for (int64_t i=0; i<num_vec; i++)
152  {
153  for (int64_t j=0; j<num_feat; j++)
154  transposed[i+j*num_vec]=matrix[i*num_feat+j];
155  }
156 
157  SG_FREE(matrix);
158  matrix=transposed;
159 
160  CMath::swap(num_feat, num_vec);
161 }
162 
163 template <class T>
164 void SGMatrix<T>::create_diagonal_matrix(T* matrix, T* v,int32_t size)
165 {
166  for(int64_t i=0;i<size;i++)
167  {
168  for(int64_t j=0;j<size;j++)
169  {
170  if(i==j)
171  matrix[j*size+i]=v[i];
172  else
173  matrix[j*size+i]=0;
174  }
175  }
176 }
177 
178 template <class T>
180  float64_t* mat, int32_t cols, int32_t rows)
181 {
182  float64_t trace=0;
183  for (int64_t i=0; i<rows; i++)
184  trace+=mat[i*cols+i];
185  return trace;
186 }
187 
188 template <class T>
189 T* SGMatrix<T>::get_row_sum(T* matrix, int32_t m, int32_t n)
190 {
191  T* rowsums=SG_CALLOC(T, n);
192 
193  for (int64_t i=0; i<n; i++)
194  {
195  for (int64_t j=0; j<m; j++)
196  rowsums[i]+=matrix[j+i*m];
197  }
198  return rowsums;
199 }
200 
201 template <class T>
202 T* SGMatrix<T>::get_column_sum(T* matrix, int32_t m, int32_t n)
203 {
204  T* colsums=SG_CALLOC(T, m);
205 
206  for (int64_t i=0; i<n; i++)
207  {
208  for (int64_t j=0; j<m; j++)
209  colsums[j]+=matrix[j+i*m];
210  }
211  return colsums;
212 }
213 
214 template <class T>
216 {
217  center_matrix(matrix, num_rows, num_cols);
218 }
219 
220 template <class T>
221 void SGMatrix<T>::center_matrix(T* matrix, int32_t m, int32_t n)
222 {
223  float64_t num_data=n;
224 
225  T* colsums=get_column_sum(matrix, m,n);
226  T* rowsums=get_row_sum(matrix, m,n);
227 
228  for (int32_t i=0; i<m; i++)
229  colsums[i]/=num_data;
230  for (int32_t j=0; j<n; j++)
231  rowsums[j]/=num_data;
232 
233  T s=SGVector<T>::sum(rowsums, n)/num_data;
234 
235  for (int64_t i=0; i<n; i++)
236  {
237  for (int64_t j=0; j<m; j++)
238  matrix[i*m+j]+=s-colsums[j]-rowsums[i];
239  }
240 
241  SG_FREE(rowsums);
242  SG_FREE(colsums);
243 }
244 
245 template <class T>
247 {
248  /* compute "row" sums (which I would call column sums), i.e. sum of all
249  * elements in a fixed column */
250  T* means=get_row_sum(matrix, num_rows, num_cols);
251 
252  /* substract column mean from every corresponding entry */
253  for (int64_t i=0; i<num_cols; ++i)
254  {
255  means[i]/=num_rows;
256  for (int64_t j=0; j<num_rows; ++j)
257  matrix[i*num_rows+j]-=means[i];
258  }
259 
260  SG_FREE(means);
261 }
262 
263 template<class T> void SGMatrix<T>::display_matrix(const char* name) const
264 {
265  display_matrix(matrix, num_rows, num_cols, name);
266 }
267 
268 template <class T>
270  const SGMatrix<T> matrix, const char* name,
271  const char* prefix)
272 {
273  matrix.display_matrix();
274 }
275 
276 template <>
278  const bool* matrix, int32_t rows, int32_t cols, const char* name,
279  const char* prefix)
280 {
281  ASSERT(rows>=0 && cols>=0)
282  SG_SPRINT("%s%s=[\n", prefix, name)
283  for (int64_t i=0; i<rows; i++)
284  {
285  SG_SPRINT("%s[", prefix)
286  for (int64_t j=0; j<cols; j++)
287  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i] ? 1 : 0,
288  j==cols-1? "" : ",");
289  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
290  }
291  SG_SPRINT("%s]\n", prefix)
292 }
293 
294 template <>
296  const char* matrix, int32_t rows, int32_t cols, const char* name,
297  const char* prefix)
298 {
299  ASSERT(rows>=0 && cols>=0)
300  SG_SPRINT("%s%s=[\n", prefix, name)
301  for (int64_t i=0; i<rows; i++)
302  {
303  SG_SPRINT("%s[", prefix)
304  for (int64_t j=0; j<cols; j++)
305  SG_SPRINT("%s\t%c%s", prefix, matrix[j*rows+i],
306  j==cols-1? "" : ",");
307  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
308  }
309  SG_SPRINT("%s]\n", prefix)
310 }
311 
312 template <>
314  const int8_t* matrix, int32_t rows, int32_t cols, const char* name,
315  const char* prefix)
316 {
317  ASSERT(rows>=0 && cols>=0)
318  SG_SPRINT("%s%s=[\n", prefix, name)
319  for (int64_t i=0; i<rows; i++)
320  {
321  SG_SPRINT("%s[", prefix)
322  for (int64_t j=0; j<cols; j++)
323  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
324  j==cols-1? "" : ",");
325  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
326  }
327  SG_SPRINT("%s]\n", prefix)
328 }
329 
330 template <>
332  const uint8_t* matrix, int32_t rows, int32_t cols, const char* name,
333  const char* prefix)
334 {
335  ASSERT(rows>=0 && cols>=0)
336  SG_SPRINT("%s%s=[\n", prefix, name)
337  for (int64_t i=0; i<rows; i++)
338  {
339  SG_SPRINT("%s[", prefix)
340  for (int64_t j=0; j<cols; j++)
341  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
342  j==cols-1? "" : ",");
343  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
344  }
345  SG_SPRINT("%s]\n", prefix)
346 }
347 
348 template <>
350  const int16_t* matrix, int32_t rows, int32_t cols, const char* name,
351  const char* prefix)
352 {
353  ASSERT(rows>=0 && cols>=0)
354  SG_SPRINT("%s%s=[\n", prefix, name)
355  for (int64_t i=0; i<rows; i++)
356  {
357  SG_SPRINT("%s[", prefix)
358  for (int64_t j=0; j<cols; j++)
359  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
360  j==cols-1? "" : ",");
361  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
362  }
363  SG_SPRINT("%s]\n", prefix)
364 }
365 
366 template <>
368  const uint16_t* matrix, int32_t rows, int32_t cols, const char* name,
369  const char* prefix)
370 {
371  ASSERT(rows>=0 && cols>=0)
372  SG_SPRINT("%s%s=[\n", prefix, name)
373  for (int64_t i=0; i<rows; i++)
374  {
375  SG_SPRINT("%s[", prefix)
376  for (int64_t j=0; j<cols; j++)
377  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
378  j==cols-1? "" : ",");
379  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
380  }
381  SG_SPRINT("%s]\n", prefix)
382 }
383 
384 
385 template <>
387  const int32_t* matrix, int32_t rows, int32_t cols, const char* name,
388  const char* prefix)
389 {
390  ASSERT(rows>=0 && cols>=0)
391  SG_SPRINT("%s%s=[\n", prefix, name)
392  for (int64_t i=0; i<rows; i++)
393  {
394  SG_SPRINT("%s[", prefix)
395  for (int64_t j=0; j<cols; j++)
396  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
397  j==cols-1? "" : ",");
398  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
399  }
400  SG_SPRINT("%s]\n", prefix)
401 }
402 
403 template <>
405  const uint32_t* matrix, int32_t rows, int32_t cols, const char* name,
406  const char* prefix)
407 {
408  ASSERT(rows>=0 && cols>=0)
409  SG_SPRINT("%s%s=[\n", prefix, name)
410  for (int64_t i=0; i<rows; i++)
411  {
412  SG_SPRINT("%s[", prefix)
413  for (int64_t j=0; j<cols; j++)
414  SG_SPRINT("%s\t%d%s", prefix, matrix[j*rows+i],
415  j==cols-1? "" : ",");
416  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
417  }
418  SG_SPRINT("%s]\n", prefix)
419 }
420 template <>
422  const int64_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 uint64_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 float32_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%.18g%s", prefix, (float) 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 float64_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%.18g%s", prefix, (double) 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 template <>
494  const floatmax_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%.18g%s", prefix, (double) 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 complex128_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(%.18g+i%.18g)%s", prefix, matrix[j*rows+i].real(),
522  matrix[j*rows+i].imag(), j==cols-1? "" : ",");
523  SG_SPRINT("%s]%s\n", prefix, i==rows-1? "" : ",")
524  }
525  SG_SPRINT("%s]\n", prefix)
526 }
527 
528 template <>
530 {
532  return SGMatrix<char>();
533 }
534 
535 template <>
537 {
538  SGMatrix<int8_t> I(size, size);
539  for (index_t i=0; i<size; ++i)
540  {
541  for (index_t j=0; j<size; ++j)
542  I(i,j)=i==j ? scale : 0.0;
543  }
544 
545  return I;
546 }
547 
548 template <>
550 {
551  SGMatrix<uint8_t> I(size, size);
552  for (index_t i=0; i<size; ++i)
553  {
554  for (index_t j=0; j<size; ++j)
555  I(i,j)=i==j ? scale : 0.0;
556  }
557 
558  return I;
559 }
560 
561 template <>
563 {
564  SGMatrix<bool> I(size, size);
565  for (index_t i=0; i<size; ++i)
566  {
567  for (index_t j=0; j<size; ++j)
568  I(i,j)=i==j ? scale : (!scale);
569  }
570 
571  return I;
572 }
573 
574 template <>
576 {
577  SGMatrix<int16_t> I(size, size);
578  for (index_t i=0; i<size; ++i)
579  {
580  for (index_t j=0; j<size; ++j)
581  I(i,j)=i==j ? scale : 0.0;
582  }
583 
584  return I;
585 }
586 
587 template <>
589 {
590  SGMatrix<uint16_t> I(size, size);
591  for (index_t i=0; i<size; ++i)
592  {
593  for (index_t j=0; j<size; ++j)
594  I(i,j)=i==j ? scale : 0.0;
595  }
596 
597  return I;
598 }
599 
600 template <>
602 {
603  SGMatrix<int32_t> I(size, size);
604  for (index_t i=0; i<size; ++i)
605  {
606  for (index_t j=0; j<size; ++j)
607  I(i,j)=i==j ? scale : 0.0;
608  }
609 
610  return I;
611 }
612 
613 template <>
615 {
616  SGMatrix<uint32_t> I(size, size);
617  for (index_t i=0; i<size; ++i)
618  {
619  for (index_t j=0; j<size; ++j)
620  I(i,j)=i==j ? scale : 0.0;
621  }
622 
623  return I;
624 }
625 
626 template <>
628 {
629  SGMatrix<int64_t> I(size, size);
630  for (index_t i=0; i<size; ++i)
631  {
632  for (index_t j=0; j<size; ++j)
633  I(i,j)=i==j ? scale : 0.0;
634  }
635 
636  return I;
637 }
638 
639 template <>
641 {
642  SGMatrix<uint64_t> I(size, size);
643  for (index_t i=0; i<size; ++i)
644  {
645  for (index_t j=0; j<size; ++j)
646  I(i,j)=i==j ? scale : 0.0;
647  }
648 
649  return I;
650 }
651 
652 template <>
654 {
655  SGMatrix<float32_t> I(size, size);
656  for (index_t i=0; i<size; ++i)
657  {
658  for (index_t j=0; j<size; ++j)
659  I(i,j)=i==j ? scale : 0.0;
660  }
661 
662  return I;
663 }
664 
665 template <>
667 {
668  SGMatrix<float64_t> I(size, size);
669  for (index_t i=0; i<size; ++i)
670  {
671  for (index_t j=0; j<size; ++j)
672  I(i,j)=i==j ? scale : 0.0;
673  }
674 
675  return I;
676 }
677 
678 template <>
680 {
681  SGMatrix<floatmax_t> I(size, size);
682  for (index_t i=0; i<size; ++i)
683  {
684  for (index_t j=0; j<size; ++j)
685  I(i,j)=i==j ? scale : 0.0;
686  }
687 
688  return I;
689 }
690 
691 template <>
693 {
694  SGMatrix<complex128_t> I(size, size);
695  for (index_t i=0; i<size; ++i)
696  {
697  for (index_t j=0; j<size; ++j)
698  I(i,j)=i==j ? scale : complex128_t(0.0);
699  }
700 
701  return I;
702 }
703 
704 
705 template <class T>
707 {
709 
710  float64_t subtract=1.0/size;
711  for (index_t i=0; i<size; ++i)
712  {
713  for (index_t j=0; j<0; ++j)
714  H(i,j)-=subtract;
715  }
716 
717  return H;
718 }
719 
720 //Howto construct the pseudo inverse (from "The Matrix Cookbook")
721 //
722 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
723 //
724 //The matrix A+ known as the pseudo inverse is unique and does always exist.
725 //
726 //The pseudo inverse A+ can be constructed from the singular value
727 //decomposition A = UDV^T , by A^+ = V(D+)U^T.
728 
729 #ifdef HAVE_LAPACK
730 template <class T>
732  float64_t* matrix, int32_t rows, int32_t cols, float64_t* target)
733 {
734  if (!target)
735  target=SG_MALLOC(float64_t, rows*cols);
736 
737  char jobu='A';
738  char jobvt='A';
739  int m=rows; /* for calling external lib */
740  int n=cols; /* for calling external lib */
741  int lda=m; /* for calling external lib */
742  int ldu=m; /* for calling external lib */
743  int ldvt=n; /* for calling external lib */
744  int info=-1; /* for calling external lib */
745  int32_t lsize=CMath::min((int32_t) m, (int32_t) n);
746  double* s=SG_MALLOC(double, lsize);
747  double* u=SG_MALLOC(double, m*m);
748  double* vt=SG_MALLOC(double, n*n);
749 
750  wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
751  ASSERT(info==0)
752 
753  for (int64_t i=0; i<n; i++)
754  {
755  for (int64_t j=0; j<lsize; j++)
756  vt[i*n+j]=vt[i*n+j]/s[j];
757  }
758 
759  cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
760 
761  SG_FREE(u);
762  SG_FREE(vt);
763  SG_FREE(s);
764 
765  return target;
766 }
767 
769 template <class T>
771 {
772  ASSERT(matrix.num_cols==matrix.num_rows)
773  int32_t* ipiv = SG_MALLOC(int32_t, matrix.num_cols);
774  clapack_dgetrf(CblasColMajor,matrix.num_cols,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
775  clapack_dgetri(CblasColMajor,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
776  SG_FREE(ipiv);
777 }
778 
779 template <class T>
781 {
782  if (matrix.num_rows!=matrix.num_rows)
783  {
784  SG_SERROR("SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix"
785  " rows and columns are not equal!\n");
786  }
787 
788  /* use reference counting for SGVector */
789  SGVector<float64_t> result(NULL, 0, true);
790  result.vlen=matrix.num_rows;
791  result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows,
792  matrix.num_rows);
793  return result;
794 }
795 
796 template <class T>
797 double* SGMatrix<T>::compute_eigenvectors(double* matrix, int n, int m)
798 {
799  ASSERT(n == m)
800 
801  char V='V';
802  char U='U';
803  int info;
804  int ord=n;
805  int lda=n;
806  double* eigenvalues=SG_CALLOC(float64_t, n+1);
807 
808  // lapack sym matrix eigenvalues+vectors
809  wrap_dsyev(V, U, ord, matrix, lda,
810  eigenvalues, &info);
811 
812  if (info!=0)
813  SG_SERROR("DSYEV failed with code %d\n", info)
814 
815  return eigenvalues;
816 }
817 
818 template <class T>
819 void SGMatrix<T>::compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors,
820  int n, int il, int iu)
821 {
822  eigenvalues = SG_MALLOC(double, n);
823  eigenvectors = SG_MALLOC(double, (iu-il+1)*n);
824  int status = 0;
825  wrap_dsyevr('V','U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status);
826  ASSERT(status==0)
827 }
828 
829 #endif //HAVE_LAPACK
830 
831 template <class T>
834  bool transpose_A, bool transpose_B, float64_t scale)
835 {
836  /* these variables store size of transposed matrices*/
837  index_t cols_A=transpose_A ? A.num_rows : A.num_cols;
838  index_t rows_A=transpose_A ? A.num_cols : A.num_rows;
839  index_t rows_B=transpose_B ? B.num_cols : B.num_rows;
840  index_t cols_B=transpose_B ? B.num_rows : B.num_cols;
841 
842  /* do a dimension check */
843  if (cols_A!=rows_B)
844  {
845  SG_SERROR("SGMatrix::matrix_multiply(): Dimension mismatch: "
846  "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
847  }
848 
849  /* allocate result matrix */
850  SGMatrix<float64_t> C(rows_A, cols_B);
851  C.zero();
852 #ifdef HAVE_LAPACK
853  /* multiply */
854  cblas_dgemm(CblasColMajor,
855  transpose_A ? CblasTrans : CblasNoTrans,
856  transpose_B ? CblasTrans : CblasNoTrans,
857  rows_A, cols_B, cols_A, scale,
858  A.matrix, A.num_rows, B.matrix, B.num_rows,
859  0.0, C.matrix, C.num_rows);
860 #else
861  for (int32_t i=0; i<rows_A; i++)
862  {
863  for (int32_t j=0; j<cols_B; j++)
864  {
865  for (int32_t k=0; k<cols_A; k++)
866  C(i,j) += A(i,k)*B(k,j);
867  }
868  }
869 #endif //HAVE_LAPACK
870 
871  return C;
872 }
873 
874 template<class T>
876  index_t num_cols, SGMatrix<T> pre_allocated)
877 {
878  SGMatrix<T> result;
879 
880  /* evtl use pre-allocated space */
881  if (pre_allocated.matrix)
882  {
883  result=pre_allocated;
884 
885  /* check dimension consistency */
886  if (pre_allocated.num_rows!=num_rows ||
887  pre_allocated.num_cols!=num_cols)
888  {
889  SG_SERROR("SGMatrix<T>::get_allocated_matrix(). Provided target"
890  "matrix dimensions (%dx%d) do not match passed data "
891  "dimensions (%dx%d)!\n", pre_allocated.num_rows,
892  pre_allocated.num_cols, num_rows, num_cols);
893  }
894  }
895  else
896  {
897  /* otherwise, allocate space */
898  result=SGMatrix<T>(num_rows, num_cols);
899  }
900 
901  return result;
902 }
903 
904 template<class T>
906 {
907  matrix=((SGMatrix*)(&orig))->matrix;
908  num_rows=((SGMatrix*)(&orig))->num_rows;
909  num_cols=((SGMatrix*)(&orig))->num_cols;
910 }
911 
912 template<class T>
914 {
915  matrix=NULL;
916  num_rows=0;
917  num_cols=0;
918 }
919 
920 template<class T>
922 {
923  SG_FREE(matrix);
924  matrix=NULL;
925  num_rows=0;
926  num_cols=0;
927 }
928 
929 template<class T>
931 {
932  ASSERT(loader)
933  unref();
934 
936  SGMatrix<T> mat;
937  loader->get_matrix(mat.matrix, mat.num_rows, mat.num_cols);
938  copy_data(mat);
939  copy_refcount(mat);
940  ref();
942 }
943 
944 template<>
946 {
947  SG_SERROR("SGMatrix::load():: Not supported for complex128_t\n");
948 }
949 
950 template<class T>
952 {
953  ASSERT(writer)
955  writer->set_matrix(matrix, num_rows, num_cols);
957 }
958 
959 template<>
961 {
962  SG_SERROR("SGMatrix::save():: Not supported for complex128_t\n");
963 }
964 
965 template<class T>
967 {
968  SGVector<T> rowv(num_cols);
969  for (int64_t i = 0; i < num_cols; i++)
970  {
971  rowv[i] = matrix[i*num_rows+row];
972  }
973  return rowv;
974 }
975 
976 template<class T>
978 {
979  index_t diag_vlen=CMath::min(num_cols, num_rows);
980  SGVector<T> diag(diag_vlen);
981 
982  for (int64_t i=0; i<diag_vlen; i++)
983  {
984  diag[i]=matrix[i*num_rows+i];
985  }
986 
987  return diag;
988 }
989 
990 template class SGMatrix<bool>;
991 template class SGMatrix<char>;
992 template class SGMatrix<int8_t>;
993 template class SGMatrix<uint8_t>;
994 template class SGMatrix<int16_t>;
995 template class SGMatrix<uint16_t>;
996 template class SGMatrix<int32_t>;
997 template class SGMatrix<uint32_t>;
998 template class SGMatrix<int64_t>;
999 template class SGMatrix<uint64_t>;
1000 template class SGMatrix<float32_t>;
1001 template class SGMatrix<float64_t>;
1002 template class SGMatrix<floatmax_t>;
1003 template class SGMatrix<complex128_t>;
1004 }

SHOGUN Machine Learning Toolbox - Documentation