Skip to content

Commit

Permalink
-added convienience methods for display_matrix/vector based on SGVector
Browse files Browse the repository at this point in the history 
-added matrix_multiply based on LAPACK
  • Loading branch information
@karlnapf
karlnapf committed Jun 3, 2012
1 parent 4343cf8 commit e04a0e6
Show file tree
Hide file tree
Showing 3 changed files with 180 additions and 24 deletions.
102 changes: 79 additions & 23 deletions examples/undocumented/libshogun/mathematics_lapack.cpp
View file
Expand Up @@ -23,8 +23,6 @@ bool is_equal(float64_t a, float64_t b, float64_t eps)

void test_ev()
{
#ifdef HAVE_LAPACK

SGMatrix<float64_t> A(3,3);
A(0,0)=0;
A(0,1)=1;
Expand Down Expand Up @@ -60,13 +58,72 @@ void test_ev()
ASSERT(is_equal(ev[0], -sqrt22*2, eps));
ASSERT(is_equal(ev[1], 0, eps));
ASSERT(is_equal(ev[2], sqrt22*2, eps));
}

#endif /* HAVE_LAPACK */
void test_matrix_multiply()
{
index_t n=10;
SGMatrix<float64_t> I(n, n);
CMath::fill_vector(I.matrix, n*n, 0.0);
for (index_t i=0; i<n; ++i)
I(i,i)=1.0;

index_t m=4;
SGMatrix<float64_t> A(n, m);
CMath::range_fill_vector(A.matrix, m*n);
CMath::display_matrix(I, "I");
CMath::transpose_matrix(A.matrix, A.num_rows, A.num_cols);
CMath::display_matrix(A, "A transposed");
CMath::transpose_matrix(A.matrix, A.num_rows, A.num_cols);
CMath::display_matrix(A, "A");

SG_SPRINT("multiply A by I and check result\n");
SGMatrix<float64_t> A2=CMath::matrix_multiply(I, A);
ASSERT(A2.num_rows==A.num_rows);
ASSERT(A2.num_cols==A.num_cols);
CMath::display_matrix(A2);
for (index_t i=0; i<A2.num_rows; ++i)
{
for (index_t j=0; j<A2.num_cols; ++j)
ASSERT(A(i,j)==A2(i,j));
}

SG_SPRINT("multiply A by transposed I and check result\n");
SGMatrix<float64_t> A3=CMath::matrix_multiply(I, A, true);
ASSERT(A3.num_rows==I.num_rows);
ASSERT(A3.num_cols==A.num_cols);
CMath::display_matrix(A3);
for (index_t i=0; i<A2.num_rows; ++i)
{
for (index_t j=0; j<A2.num_cols; ++j)
ASSERT(A(i,j)==A3(i,j));
}

SG_SPRINT("multiply transposed A by I and check result\n");
SGMatrix<float64_t> A4=CMath::matrix_multiply(A, I, true, false);
ASSERT(A4.num_rows==A.num_cols);
ASSERT(A4.num_cols==I.num_cols);
CMath::display_matrix(A4);
for (index_t i=0; i<A.num_rows; ++i)
{
for (index_t j=0; j<A.num_cols; ++j)
ASSERT(A(i,j)==A4(j,i));
}

SG_SPRINT("multiply A by scaled I and check result\n");
SGMatrix<float64_t> A5=CMath::matrix_multiply(I, A, false, false, n);
ASSERT(A5.num_rows==I.num_rows);
ASSERT(A5.num_cols==A.num_cols);
CMath::display_matrix(A5);
for (index_t i=0; i<A2.num_rows; ++i)
{
for (index_t j=0; j<A2.num_cols; ++j)
ASSERT(n*A(i,j)==A5(i,j));
}
}

void test_lapack()
{
#ifdef HAVE_LAPACK
// size of square matrix
int N = 100;

Expand All @@ -93,10 +150,8 @@ void test_lapack()
status = 0;
wrap_dsygvx(1,'V','U',N,double_matrix,N,double_matrix,N,1,3,double_eigenvalues,double_eigenvectors,&status);
if (status!=0)
{
printf("DSYGVX/SSYGVX failed with code %d\n",status);
ASSERT(false);
}
SG_SERROR("DSYGVX/SSYGVX failed with code %d\n",status);

delete[] double_eigenvectors;

// DGEQRF+DORGQR
Expand All @@ -105,10 +160,8 @@ void test_lapack()
wrap_dgeqrf(N,N,double_matrix,N,double_tau,&status);
wrap_dorgqr(N,N,N,double_matrix,N,double_tau,&status);
if (status!=0)
{
printf("DGEQRF/DORGQR failed with code %d\n",status);
ASSERT(false);
}
SG_SERROR("DGEQRF/DORGQR failed with code %d\n",status);

delete[] double_tau;

// DGESVD
Expand All @@ -120,10 +173,8 @@ void test_lapack()
status = 0;
wrap_dgesvd('A','A',N,N,double_matrix,N,double_s,double_U,N,double_Vt,N,&status);
if (status!=0)
{
printf("DGESVD failed with code %d\n",status);
ASSERT(false);
}
SG_SERROR("DGESVD failed with code %d\n",status);

delete[] double_s;
delete[] double_U;
delete[] double_Vt;
Expand All @@ -132,24 +183,29 @@ void test_lapack()
status = 0;
wrap_dsyev('V','U',N,double_matrix,N,double_eigenvalues,&status);
if (status!=0)
{
printf("DSYEV failed with code %d\n",status);
ASSERT(false);
}
SG_SERROR("DSYEV failed with code %d\n",status);

delete[] double_eigenvalues;
delete[] double_matrix;

#endif // HAVE_LAPACK
}

int main(int argc, char** argv)
{
init_shogun();
init_shogun_with_defaults();

#ifdef HAVE_LAPACK
SG_SPRINT("checking lapack\n");
test_lapack();

SG_SPRINT("compute_eigenvectors\n");
test_ev();

SG_SPRINT("matrix_multiply\n");
test_matrix_multiply();
#endif

exit_shogun();
return 0;
}


81 changes: 80 additions & 1 deletion src/shogun/mathematics/Math.cpp
View file
Expand Up @@ -164,6 +164,27 @@ void CMath::display_vector(const floatmax_t* vector, int32_t n,
SG_SPRINT("%s]\n", prefix);
}

template <>
void CMath::display_vector(const SGVector<float64_t> vector, const char* name,
const char* prefix)
{
CMath::display_vector(vector.vector, vector.vlen, name, prefix);
}

template <>
void CMath::display_vector(const SGVector<float32_t> vector, const char* name,
const char* prefix)
{
CMath::display_vector(vector.vector, vector.vlen, name, prefix);
}

template <>
void CMath::display_vector(const SGVector<int32_t> vector, const char* name,
const char* prefix)
{
CMath::display_vector(vector.vector, vector.vlen, name, prefix);
}

template <>
void CMath::display_matrix(
const int32_t* matrix, int32_t rows, int32_t cols, const char* name,
Expand Down Expand Up @@ -218,6 +239,33 @@ void CMath::display_matrix(
SG_SPRINT("%s]\n", prefix);
}

template <>
void CMath::display_matrix(
const SGMatrix<float64_t> matrix, const char* name,
const char* prefix)
{
CMath::display_matrix(matrix.matrix, matrix.num_rows, matrix.num_cols,
name, prefix);
}

template <>
void CMath::display_matrix(
const SGMatrix<float32_t> matrix, const char* name,
const char* prefix)
{
CMath::display_matrix(matrix.matrix, matrix.num_rows, matrix.num_cols,
name, prefix);
}

template <>
void CMath::display_matrix(
const SGMatrix<int32_t> matrix, const char* name,
const char* prefix)
{
CMath::display_matrix(matrix.matrix, matrix.num_rows, matrix.num_cols,
name, prefix);
}

}

SGVector<float64_t> CMath::fishers_exact_test_for_multiple_2x3_tables(SGMatrix<float64_t> tables)
Expand Down Expand Up @@ -488,7 +536,6 @@ float64_t CMath::entropy(float64_t* p, int32_t len)
return (float64_t) e;
}


//Howto construct the pseudo inverse (from "The Matrix Cookbook")
//
//Assume A does not have full rank, i.e. A is n \times m and rank(A) = r < min(n;m).
Expand Down Expand Up @@ -582,6 +629,38 @@ double* CMath::compute_eigenvectors(double* matrix, int n, int m)

return eigenvalues;
}

SGMatrix<float64_t> CMath::matrix_multiply(
SGMatrix<float64_t> A, SGMatrix<float64_t> B,
bool transpose_A, bool transpose_B, float64_t scale)
{
/* these variables store size of transposed matrices*/
index_t cols_A=transpose_A ? A.num_rows : A.num_cols;
index_t rows_A=transpose_A ? A.num_cols : A.num_rows;
index_t rows_B=transpose_B ? B.num_cols : B.num_rows;
index_t cols_B=transpose_B ? B.num_rows : B.num_cols;

/* do a dimension check */
if (cols_A!=rows_B)
{
SG_SERROR("CMath::matrix_multiply(): Dimension mismatch: "
"A(%dx%d)*B(%dx%D)\n", rows_A, cols_A, rows_B, cols_B);
}

/* allocate result matrix */
SGMatrix<float64_t> C(rows_A, cols_B);

/* multiply */
cblas_dgemm(CblasColMajor,
transpose_A ? CblasTrans : CblasNoTrans,
transpose_B ? CblasTrans : CblasNoTrans,
rows_A, cols_B, cols_A, scale,
A.matrix, A.num_rows, B.matrix, B.num_rows,
0.0, C.matrix, C.num_rows);

return C;
}

#endif //HAVE_LAPACK

float64_t CMath::dot(const float64_t* v1, const float64_t* v2, int32_t n)
Expand Down
21 changes: 21 additions & 0 deletions src/shogun/mathematics/Math.h
View file
Expand Up @@ -514,6 +514,19 @@ class CMath : public CSGObject
* */
static double* compute_eigenvectors(double* matrix, int n, int m);

/* Computes scale* A*B, where A and B may be transposed.
* Asserts for matching inner dimensions.
* @param A matrix A
* @param transpose_A optional whether A should be transposed before
* @param B matrix B
* @param transpose_B optional whether B should be transposed before
* @param scale optional scaling factor for result
*/
static SGMatrix<float64_t> matrix_multiply(
SGMatrix<float64_t> A, SGMatrix<float64_t> B,
bool transpose_A=false, bool transpose_B=false,
float64_t scale=1.0);

/// inverses square matrix in-place
static void inverse(SGMatrix<float64_t> matrix);

Expand Down Expand Up @@ -1260,11 +1273,19 @@ class CMath : public CSGObject
const T* vector, int32_t n, const char* name="vector",
const char* prefix="");

template <class T> static void display_vector(
const SGVector<T>, const char* name="vector",
const char* prefix="");

/// display matrix (useful for debugging)
template <class T> static void display_matrix(
const T* matrix, int32_t rows, int32_t cols,
const char* name="matrix", const char* prefix="");

template <class T> static void display_matrix(
const SGMatrix<T> matrix, const char* name="matrix",
const char* prefix="");

/** performs a quicksort on an array output of length size
* it is sorted in ascending order
* (for type T1) and returns the index (type T2)
Expand Down

0 comments on commit e04a0e6

Please sign in to comment.