egsl_ops.c
Go to the documentation of this file.
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 /*      egsl_print("m1",v1);
00056         egsl_print("m2",v2); */
00057         gsl_matrix * m1 = egsl_gslm(v1); 
00058         gsl_matrix * m2 = egsl_gslm(v2);
00059         
00060 /*      printf("m1 size = %d,%d j = %d\n",m1->size1,m1->size2,j); */
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         /* Make LU decomposition of matrix m */
00124         int s;
00125         gsl_linalg_LU_decomp (m, perm, &s);
00126         /* Invert the matrix m */
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         /* Check for v to be square */
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 


csm
Author(s): Andrea Censi
autogenerated on Mon Jan 16 2017 03:48:29