00001
00006 #include "matrix_ex.h"
00007 #include <gsl/gsl_linalg.h>
00008 #include <gsl/gsl_blas.h>
00009
00016 void MatrixEx::matrix_pseudo_inverse(gsl_matrix *dest, const gsl_matrix *src)
00017 {
00018 if (src->size1 < src->size2)
00019 pseudo_inverse_right(dest, src);
00020 else if (src->size1 > src->size2)
00021 pseudo_inverse_left(dest, src);
00022 else
00023 matrix_inverse(dest, src);
00024 }
00025
00032 void MatrixEx::pseudo_inverse_right(gsl_matrix *dest, const gsl_matrix *src)
00033 {
00034 gsl_matrix *trans, *prod, *inv;
00035
00036 trans = gsl_matrix_alloc(src->size2, src->size1);
00037 prod = gsl_matrix_alloc(src->size1, src->size1);
00038 inv = gsl_matrix_alloc(src->size1, src->size1);
00039
00040 gsl_matrix_transpose_memcpy(trans, src);
00041 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, src, trans, 0.0, prod);
00042 matrix_inverse(inv, prod);
00043 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, trans, inv, 0.0, dest);
00044
00045 gsl_matrix_free(trans);
00046 gsl_matrix_free(prod);
00047 gsl_matrix_free(inv);
00048 }
00049
00056 void MatrixEx::pseudo_inverse_left(gsl_matrix *dest, const gsl_matrix *src)
00057 {
00058 gsl_matrix *trans, *prod, *inv;
00059
00060 trans = gsl_matrix_alloc(src->size2, src->size1);
00061 prod = gsl_matrix_alloc(src->size2, src->size2);
00062 inv = gsl_matrix_alloc(src->size2, src->size2);
00063
00064 gsl_matrix_transpose_memcpy(trans, src);
00065 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, trans, src, 0.0, prod);
00066 matrix_inverse(inv, prod);
00067 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv, trans, 0.0, dest);
00068
00069 gsl_matrix_free(trans);
00070 gsl_matrix_free(prod);
00071 gsl_matrix_free(inv);
00072 }
00073
00080 void MatrixEx::matrix_inverse(gsl_matrix *dest, const gsl_matrix *src)
00081 {
00082 gsl_matrix *tmp;
00083 int signum;
00084 gsl_permutation *p = gsl_permutation_alloc(src->size1);
00085 tmp = gsl_matrix_alloc(src->size1, src->size2);
00086 gsl_matrix_memcpy(tmp, src);
00087
00088 gsl_linalg_LU_decomp(tmp, p, &signum);
00089 gsl_linalg_LU_invert(tmp, p, dest);
00090
00091 gsl_permutation_free(p);
00092 gsl_matrix_free(tmp);
00093 }
