00001 #include <gsl/gsl_matrix.h>
00002 #include <gsl/gsl_blas.h>
00003 #include <gsl/gsl_linalg.h>
00004 #include <gsl/gsl_eigen.h>
00005
00006 #include "egsl.h"
00007
00008
00009
00010 val egsl_sub(val v1,val v2){
00011 return egsl_sum(v1, egsl_scale(-1.0,v2));
00012 }
00013
00014 val egsl_compose_col(val v1, val v2){
00015 gsl_matrix *m1 = egsl_gslm(v1);
00016 gsl_matrix *m2 = egsl_gslm(v2);
00017 egsl_expect_size(v2, 0, m1->size2);
00018 val v3 = egsl_alloc(m1->size1+m2->size1,m1->size2);
00019 gsl_matrix *m3 = egsl_gslm(v3);
00020 size_t i,j;
00021 for(j=0;j<m1->size2;j++) {
00022 for(i=0;i<m1->size1;i++)
00023 gsl_matrix_set(m3, i, j, gsl_matrix_get(m1,i,j));
00024
00025 for(i=0;i<m2->size1;i++)
00026 gsl_matrix_set(m3, m1->size1+i, j, gsl_matrix_get(m2,i,j));
00027 }
00028 return v3;
00029 }
00030
00031 val egsl_compose_row(val v1, val v2){
00032 gsl_matrix *m1 = egsl_gslm(v1);
00033 gsl_matrix *m2 = egsl_gslm(v2);
00034 egsl_expect_size(v2, m1->size1, 0);
00035 val v3 = egsl_alloc(m1->size1, m1->size2 + m2->size2);
00036 gsl_matrix *m3 = egsl_gslm(v3);
00037 size_t i,j;
00038 for(i=0;i<m1->size1;i++) {
00039 for(j=0;j<m1->size2;j++)
00040 gsl_matrix_set(m3, i, j, gsl_matrix_get(m1,i,j));
00041
00042 for(j=0;j<m2->size2;j++)
00043 gsl_matrix_set(m3, i, m1->size2+j, gsl_matrix_get(m2,i,j));
00044 }
00045 return v3;
00046 }
00047
00048 void egsl_add_to(val v1, val v2) {
00049 gsl_matrix * m1 = egsl_gslm(v1);
00050 gsl_matrix * m2 = egsl_gslm(v2);
00051 gsl_matrix_add(m1,m2);
00052 }
00053
00054 void egsl_add_to_col(val v1, size_t j, val v2) {
00055
00056
00057 gsl_matrix * m1 = egsl_gslm(v1);
00058 gsl_matrix * m2 = egsl_gslm(v2);
00059
00060
00061 egsl_expect_size(v2, m1->size1, 1);
00062 size_t i;
00063 for(i=0;i<m1->size1;i++) {
00064 *gsl_matrix_ptr(m1, i, j) += gsl_matrix_get(m2,i,0);
00065 }
00066 }
00067
00068
00069 val egsl_copy_val(val v1) {
00070 gsl_matrix * m1 = egsl_gslm(v1);
00071 val v2 = egsl_alloc(m1->size1,m1->size2);
00072 gsl_matrix * m2 = egsl_gslm(v2);
00073 gsl_matrix_memcpy(m2,m1);
00074 return v2;
00075 }
00076
00077 val egsl_scale(double s, val v1){
00078 val v2 = egsl_copy_val(v1);
00079 gsl_matrix * m2 = egsl_gslm(v2);
00080 gsl_matrix_scale(m2, s);
00081 return v2;
00082 }
00083
00084 val egsl_sum(val v1, val v2){
00085 gsl_matrix * m1 = egsl_gslm(v1);
00086 gsl_matrix * m2 = egsl_gslm(v2);
00087 val v3 = egsl_alloc(m1->size1,m1->size2);
00088 gsl_matrix * m3 = egsl_gslm(v3);
00089 gsl_matrix_memcpy(m3,m1);
00090 gsl_matrix_add(m3,m2);
00091 return v3;
00092 }
00093
00094 val egsl_sum3(val v1, val v2, val v3){
00095 return egsl_sum(v1, egsl_sum(v2,v3));
00096 }
00097
00098 val egsl_mult(val v1, val v2){
00099 gsl_matrix * a = egsl_gslm(v1);
00100 gsl_matrix * b = egsl_gslm(v2);
00101 val v = egsl_alloc(a->size1,b->size2);
00102 gsl_matrix * ab = egsl_gslm(v);
00103 gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,a,b,0.0,ab);
00104 return v;
00105 }
00106
00107 val egsl_transpose(val v1){
00108 gsl_matrix * m1 = egsl_gslm(v1);
00109 val v2 = egsl_alloc(m1->size2,m1->size1);
00110 gsl_matrix * m2 = egsl_gslm(v2);
00111 gsl_matrix_transpose_memcpy(m2,m1);
00112 return v2;
00113 }
00114
00115 val egsl_inverse(val v1){
00116 gsl_matrix*A = egsl_gslm(v1);
00117 val v2 = egsl_alloc(A->size1,A->size1);
00118 gsl_matrix*invA = egsl_gslm(v2);
00119 size_t n = A->size1;
00120 gsl_matrix * m = gsl_matrix_alloc(n,n);
00121 gsl_matrix_memcpy(m,A);
00122 gsl_permutation * perm = gsl_permutation_alloc (n);
00123
00124 int s;
00125 gsl_linalg_LU_decomp (m, perm, &s);
00126
00127 gsl_linalg_LU_invert (m, perm, invA);
00128 gsl_permutation_free(perm);
00129 gsl_matrix_free(m);
00130 return v2;
00131 }
00132
00133 void egsl_symm_eig(val v, double* eigenvalues, val* eigenvectors) {
00134 gsl_matrix *m = egsl_gslm(v);
00135 size_t N = m->size1;
00136
00137
00138 gsl_matrix *A = gsl_matrix_alloc(N,N);
00139 gsl_matrix_memcpy(A, m);
00140
00141 gsl_vector *eval = gsl_vector_alloc(N);
00142 gsl_matrix *evec = gsl_matrix_alloc(N,N);
00143
00144 gsl_eigen_symmv_workspace * ws = gsl_eigen_symmv_alloc(N);
00145 gsl_eigen_symmv(A, eval, evec, ws);
00146 gsl_eigen_symmv_free(ws);
00147
00148
00149 gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);
00150
00151 size_t j;
00152 for(j=0;j<N;j++) {
00153 eigenvalues[j] = gsl_vector_get(eval, j);
00154 eigenvectors[j] = egsl_alloc(N,1);
00155 size_t i;
00156 for(i=0;i<N;i++)
00157 *egsl_atmp(eigenvectors[j],i,0) = gsl_matrix_get(evec,i,j);
00158 }
00159
00160
00161 gsl_vector_free(eval);
00162 gsl_matrix_free(evec);
00163 gsl_matrix_free(A);
00164 }
00165
00166
00167
00168
00169
00170
00171
00172