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 //Howto construct the pseudo inverse (from "The Matrix Cookbook")
705 //
706 //Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
707 //
708 //The matrix A+ known as the pseudo inverse is unique and does always exist.
709 //
710 //The pseudo inverse A+ can be constructed from the singular value
711 //decomposition A = UDV^T , by A^+ = V(D+)U^T.
712 
713 #ifdef HAVE_LAPACK
714 template <class T>
716  float64_t* matrix, int32_t rows, int32_t cols, float64_t* target)
717 {
718  if (!target)
719  target=SG_MALLOC(float64_t, rows*cols);
720 
721  char jobu='A';
722  char jobvt='A';
723  int m=rows; /* for calling external lib */
724  int n=cols; /* for calling external lib */
725  int lda=m; /* for calling external lib */
726  int ldu=m; /* for calling external lib */
727  int ldvt=n; /* for calling external lib */
728  int info=-1; /* for calling external lib */
729  int32_t lsize=CMath::min((int32_t) m, (int32_t) n);
730  double* s=SG_MALLOC(double, lsize);
731  double* u=SG_MALLOC(double, m*m);
732  double* vt=SG_MALLOC(double, n*n);
733 
734  wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
735  ASSERT(info==0)
736 
737  for (int64_t i=0; i<n; i++)
738  {
739  for (int64_t j=0; j<lsize; j++)
740  vt[i*n+j]=vt[i*n+j]/s[j];
741  }
742 
743  cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);
744 
745  SG_FREE(u);
746  SG_FREE(vt);
747  SG_FREE(s);
748 
749  return target;
750 }
751 
753 template <class T>
755 {
756  ASSERT(matrix.num_cols==matrix.num_rows)
757  int32_t* ipiv = SG_MALLOC(int32_t, matrix.num_cols);
758  clapack_dgetrf(CblasColMajor,matrix.num_cols,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
759  clapack_dgetri(CblasColMajor,matrix.num_cols,matrix.matrix,matrix.num_cols,ipiv);
760  SG_FREE(ipiv);
761 }
762 
763 template <class T>
765 {
766  if (matrix.num_rows!=matrix.num_cols)
767  {
768  SG_SERROR("SGMatrix::compute_eigenvectors(SGMatrix<float64_t>): matrix"
769  " rows and columns are not equal!\n");
770  }
771 
772  /* use reference counting for SGVector */
773  SGVector<float64_t> result(NULL, 0, true);
774  result.vlen=matrix.num_rows;
775  result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows,
776  matrix.num_rows);
777  return result;
778 }
779 
780 template <class T>
781 double* SGMatrix<T>::compute_eigenvectors(double* matrix, int n, int m)
782 {
783  ASSERT(n == m)
784 
785  char V='V';
786  char U='U';
787  int info;
788  int ord=n;
789  int lda=n;
790  double* eigenvalues=SG_CALLOC(float64_t, n+1);
791 
792  // lapack sym matrix eigenvalues+vectors
793  wrap_dsyev(V, U, ord, matrix, lda,
794  eigenvalues, &info);
795 
796  if (info!=0)
797  SG_SERROR("DSYEV failed with code %d\n", info)
798 
799  return eigenvalues;
800 }
801 
802 template <class T>
803 void SGMatrix<T>::compute_few_eigenvectors(double* matrix_, double*& eigenvalues, double*& eigenvectors,
804  int n, int il, int iu)
805 {
806  eigenvalues = SG_MALLOC(double, n);
807  eigenvectors = SG_MALLOC(double, (iu-il+1)*n);
808  int status = 0;
809  wrap_dsyevr('V','U',n,matrix_,n,il,iu,eigenvalues,eigenvectors,&status);
810  ASSERT(status==0)
811 }
812 
813 #endif //HAVE_LAPACK
814 
815 template <class T>
818  bool transpose_A, bool transpose_B, float64_t scale)
819 {
820  /* these variables store size of transposed matrices*/
821  index_t cols_A=transpose_A ? A.num_rows : A.num_cols;
822  index_t rows_A=transpose_A ? A.num_cols : A.num_rows;
823  index_t rows_B=transpose_B ? B.num_cols : B.num_rows;
824  index_t cols_B=transpose_B ? B.num_rows : B.num_cols;
825 
826  /* do a dimension check */
827  if (cols_A!=rows_B)
828  {
829  SG_SERROR("SGMatrix::matrix_multiply(): Dimension mismatch: "
830  "A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
831  }
832 
833  /* allocate result matrix */
834  SGMatrix<float64_t> C(rows_A, cols_B);
835  C.zero();
836 #ifdef HAVE_LAPACK
837  /* multiply */
838  cblas_dgemm(CblasColMajor,
839  transpose_A ? CblasTrans : CblasNoTrans,
840  transpose_B ? CblasTrans : CblasNoTrans,
841  rows_A, cols_B, cols_A, scale,
842  A.matrix, A.num_rows, B.matrix, B.num_rows,
843  0.0, C.matrix, C.num_rows);
844 #else
845  for (int32_t i=0; i<rows_A; i++)
846  {
847  for (int32_t j=0; j<cols_B; j++)
848  {
849  for (int32_t k=0; k<cols_A; k++)
850  C(i,j) += A(i,k)*B(k,j);
851  }
852  }
853 #endif //HAVE_LAPACK
854 
855  return C;
856 }
857 
858 template<class T>
860  index_t num_cols, SGMatrix<T> pre_allocated)
861 {
862  SGMatrix<T> result;
863 
864  /* evtl use pre-allocated space */
865  if (pre_allocated.matrix)
866  {
867  result=pre_allocated;
868 
869  /* check dimension consistency */
870  if (pre_allocated.num_rows!=num_rows ||
871  pre_allocated.num_cols!=num_cols)
872  {
873  SG_SERROR("SGMatrix<T>::get_allocated_matrix(). Provided target"
874  "matrix dimensions (%dx%d) do not match passed data "
875  "dimensions (%dx%d)!\n", pre_allocated.num_rows,
876  pre_allocated.num_cols, num_rows, num_cols);
877  }
878  }
879  else
880  {
881  /* otherwise, allocate space */
882  result=SGMatrix<T>(num_rows, num_cols);
883  }
884 
885  return result;
886 }
887 
888 template<class T>
890 {
891  matrix=((SGMatrix*)(&orig))->matrix;
892  num_rows=((SGMatrix*)(&orig))->num_rows;
893  num_cols=((SGMatrix*)(&orig))->num_cols;
894 }
895 
896 template<class T>
898 {
899  matrix=NULL;
900  num_rows=0;
901  num_cols=0;
902 }
903 
904 template<class T>
906 {
907  SG_FREE(matrix);
908  matrix=NULL;
909  num_rows=0;
910  num_cols=0;
911 }
912 
913 template<class T>
915 {
916  ASSERT(loader)
917  unref();
918 
920  SGMatrix<T> mat;
921  loader->get_matrix(mat.matrix, mat.num_rows, mat.num_cols);
922  copy_data(mat);
923  copy_refcount(mat);
924  ref();
926 }
927 
928 template<>
930 {
931  SG_SERROR("SGMatrix::load():: Not supported for complex128_t\n");
932 }
933 
934 template<class T>
936 {
937  ASSERT(writer)
939  writer->set_matrix(matrix, num_rows, num_cols);
941 }
942 
943 template<>
945 {
946  SG_SERROR("SGMatrix::save():: Not supported for complex128_t\n");
947 }
948 
949 template<class T>
951 {
952  SGVector<T> rowv(num_cols);
953  for (int64_t i = 0; i < num_cols; i++)
954  {
955  rowv[i] = matrix[i*num_rows+row];
956  }
957  return rowv;
958 }
959 
960 template<class T>
962 {
963  index_t diag_vlen=CMath::min(num_cols, num_rows);
964  SGVector<T> diag(diag_vlen);
965 
966  for (int64_t i=0; i<diag_vlen; i++)
967  {
968  diag[i]=matrix[i*num_rows+i];
969  }
970 
971  return diag;
972 }
973 
974 template class SGMatrix<bool>;
975 template class SGMatrix<char>;
976 template class SGMatrix<int8_t>;
977 template class SGMatrix<uint8_t>;
978 template class SGMatrix<int16_t>;
979 template class SGMatrix<uint16_t>;
980 template class SGMatrix<int32_t>;
981 template class SGMatrix<uint32_t>;
982 template class SGMatrix<int64_t>;
983 template class SGMatrix<uint64_t>;
984 template class SGMatrix<float32_t>;
985 template class SGMatrix<float64_t>;
986 template class SGMatrix<floatmax_t>;
987 template class SGMatrix<complex128_t>;
988 }

SHOGUN Machine Learning Toolbox - Documentation