gpc_utils.c
Go to the documentation of this file.
00001 #include "gpc.h"
00002 #include "gpc_utils.h"
00003 
00004 void m_trans(const gsl_matrix*A, gsl_matrix*A_t){
00005         gsl_matrix_transpose_memcpy(A_t,A);
00006 }
00007 
00008 void m_mult(const gsl_matrix*A, const gsl_matrix*B, gsl_matrix*AB){
00009         gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,A,B,0.0,AB);
00010 }
00011 
00012 void m_add_to(const gsl_matrix*A, gsl_matrix*B){
00013         gsl_matrix_add(B, A);
00014 }
00015 
00016 void m_scale(double m, gsl_matrix*A){
00017         gsl_matrix_scale(A,m);
00018 }
00019 
00020 void m_add (const gsl_matrix*A, const gsl_matrix*B, gsl_matrix*ApB){
00021         gsl_matrix_memcpy(ApB,A);
00022         gsl_matrix_add(ApB,B);  
00023 }
00024 
00025 void m_inv(const gsl_matrix*A, gsl_matrix*invA) {
00026         unsigned int n = A->size1;
00027         gsl_matrix * m = gsl_matrix_alloc(n,n);
00028         gsl_matrix_memcpy(m,A);
00029         gsl_permutation * perm = gsl_permutation_alloc (n);
00030         /* Make LU decomposition of matrix m */
00031         int s;
00032         gsl_linalg_LU_decomp (m, perm, &s);
00033         /* Invert the matrix m */
00034         gsl_linalg_LU_invert (m, perm, invA);
00035         gsl_permutation_free(perm);
00036         gsl_matrix_free(m);
00037 }
00038 
00039 double m_det(const gsl_matrix*A) {
00040         unsigned int n = A->size1;
00041         gsl_matrix * m = gsl_matrix_alloc(n,n);
00042         gsl_matrix_memcpy(m,A);
00043         gsl_permutation * perm = gsl_permutation_alloc (n);
00044         int sign;
00045         gsl_linalg_LU_decomp (m, perm, &sign);
00046         double det = gsl_linalg_LU_det(m, sign);
00047         
00048         gsl_permutation_free(perm);
00049         gsl_matrix_free(m);
00050         return det;
00051 }
00052 
00053 double m_dot(const gsl_matrix*A,const gsl_matrix*B) {
00054         double sum = 0;
00055         unsigned int j;
00056         for(j=0;j<A->size2;j++)
00057                 sum += gmg(A,0,j)*gmg(B,j,0);
00058         return sum;
00059 }
00060 
00061 int poly_real_roots(unsigned int n, const double*a, double *roots) {
00062         double z[(n-1)*2];
00063         gsl_poly_complex_workspace * w  = gsl_poly_complex_workspace_alloc(n);
00064         if(GSL_SUCCESS != gsl_poly_complex_solve (a, n, w, z)) {
00065                 return 0;
00066         }
00067         gsl_poly_complex_workspace_free (w);
00068         
00069         unsigned int i=0;
00070         for(i=0;i<n-1;i++) {
00071                 roots[i] = z[2*i];
00072         }       
00073         return 1;
00074 }
00075 
00076 
00077 int poly_greatest_real_root(unsigned int n, const double*a, double *root) {
00078         double z[(n-1)*2];
00079         gsl_poly_complex_workspace * w  = gsl_poly_complex_workspace_alloc(n);
00080         if(GSL_SUCCESS != gsl_poly_complex_solve (a, n, w, z)) {
00081                 return 0;
00082         }
00083         gsl_poly_complex_workspace_free (w);
00084         if(TRACE_ALGO) {
00085                 printf("Solving the equation\n a = [");
00086                 for(unsigned int i=0;i<n;i++) {
00087                         printf("%lf ", a[i]);
00088                 }
00089                 printf("]\n");
00090         }
00091 
00092         unsigned int i;
00093         double lambda = 0; int assigned = 0;
00094         for (i = 0; i < n-1; i++) {
00095                 if(TRACE_ALGO) {
00096                         fprintf (stderr, "root z%d = %+.18f + %+.18f i \n", i, z[2*i], z[2*i+1]);
00097                 }
00098 /*               XXX ==0 is bad */
00099                 if( z[2*i+1]==0 ) /* real root */
00100                         if(!assigned || (z[2*i]>lambda)) {
00101                                 assigned = 1;
00102                                 lambda = z[2*i];
00103                         }
00104         }
00105         if(TRACE_ALGO)
00106                 fprintf (stderr, "lambda = %+.18f \n", lambda); 
00107         if(!assigned) {
00108                 fprintf(stderr, "poly_greatest_real_root: Could not find real root for polynomial.\n");
00109                 fprintf(stderr, "polynomial coefficients : ");
00110                 for (i = 0; i < n; i++)
00111                         fprintf(stderr, " %lf ", a[i]);
00112                 fprintf(stderr, "\nRoots:\n");
00113                 
00114                 for (i = 0; i < n-1; i++)
00115                         fprintf (stderr, "root z%d = %+.18f + %+.18f i \n", i, z[2*i], z[2*i+1]);
00116 
00117                 return 0;
00118         }
00119         
00120         *root = lambda;
00121         return 1;
00122 }
00123 
00124 void m_display(const char*str, gsl_matrix*m) {
00125         printf("%s= \n", str);
00126         unsigned int i,j;
00127         for(i=0;i<m->size1;i++) {
00128                 printf("   ");
00129                 for(j=0;j<m->size2;j++)
00130                         printf("%e ", gmg(m,i,j));
00131                 printf("\n");
00132         }
00133 }
00134 


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