Return to call_gsl_eigen.c CVS log

Up to [local] / OpenXM / src / ox_gsl

Annotation of OpenXM/src/ox_gsl/call_gsl_eigen.c, Revision 1.1

1.1     ! takayama    1: /* $OpenXM$ */
        !             2: #include <stdio.h>
        !             3: #include <gsl/gsl_math.h>
        !             4: #include <gsl/gsl_eigen.h>
        !             5: #include "ox_gsl.h"
        !             6: extern int Debug;
        !             7: cmo *gsl_complex_to_cmo_array(gsl_complex z) {
        !             8:   cmo **cr;
        !             9:   cr = (cmo **) GC_malloc(sizeof(cmo *)*2);
        !            10:   cr[0] = (cmo *) new_cmo_double(GSL_REAL(z));
        !            11:   cr[1] = (cmo *) new_cmo_double(GSL_IMAG(z));
        !            12:   return (cmo *)new_cmo_list_array((void *)cr,2);
        !            13: }
        !            14:
        !            15: cmo *gsl_vector_complex_to_cmo_array(gsl_vector_complex * v,int n) {
        !            16:   cmo **cr;
        !            17:   int i;
        !            18:   if (n <= 0) return NULL;
        !            19:   cr = (cmo **) GC_malloc(sizeof(cmo *)*n);
        !            20:   for (i=0; i<n; i++) {
        !            21:     cr[i] = gsl_complex_to_cmo_array(gsl_vector_complex_get(v,i));
        !            22:   }
        !            23:   return (cmo *)new_cmo_list_array((void *)cr,n);
        !            24: }
        !            25:
        !            26: cmo *gsl_matrix_complex_to_cmo_array_transposed(gsl_matrix_complex *m,int row_m, int col_n) {
        !            27: /* it returns the transpose of m */
        !            28:   cmo **cr;
        !            29:   int i;
        !            30:   gsl_vector_complex_view evec_i;
        !            31:   if (col_n <= 0) return NULL;
        !            32:   cr = (cmo **) GC_malloc(sizeof(cmo *)*col_n);
        !            33:   for (i=0; i<col_n; i++) {
        !            34:     evec_i=gsl_matrix_complex_column(m,i);
        !            35:     cr[i] = gsl_vector_complex_to_cmo_array(&evec_i.vector,row_m);
        !            36:   }
        !            37:   return (cmo *)new_cmo_list_array((void *)cr,col_n);
        !            38: }
        !            39:
        !            40: /* copied from GSL manual modified.
        !            41: */
        !            42: void call_gsl_eigen_nonsymmv() {
        !            43:   int argc;
        !            44:   int nn,n0,n,ii,len;
        !            45:   cmo *mat;
        !            46:   double *data;
        !            47:
        !            48:   gsl_matrix_view m;
        !            49:   gsl_vector_complex *eval;
        !            50:   gsl_matrix_complex *evec;
        !            51:
        !            52:   gsl_eigen_nonsymmv_workspace * w;
        !            53:
        !            54:   gsl_set_error_handler((gsl_error_handler_t *)myhandler);
        !            55:
        !            56:   argc = get_i();
        !            57:   if (argc != 1) {
        !            58:     push(make_error2("call_gsl_eigen_nonsymmv accepts one arg.",NULL,0,-1));
        !            59:     return ;
        !            60:   }
        !            61:   mat = pop();
        !            62:   nn = get_length(mat);
        !            63:   n0 = (int) floor(sqrt(nn));
        !            64:   if (n0 <= 0) {
        !            65:     push(make_error2("call_gsl_eigen_nonsymmv: give a non-empty matrix.",NULL,0,-1));
        !            66:     return ;
        !            67:   }
        !            68:   for (n=n0-1; n<=n0+1; n++) {
        !            69:     if (n*n-nn == 0) break;
        !            70:   }
        !            71:   if (n > n0+1) {
        !            72:     push(make_error2("call_gsl_eigen_nonsymmv: matrix size must be n^2",NULL,0,-1));
        !            73:     return ;
        !            74:   }
        !            75:
        !            76:   data = cmo2double_list(&len, mat);
        !            77:   for (ii=0; ii<len; ii++) printf("%lg, ",data[ii]);
        !            78:
        !            79:   m
        !            80:     = gsl_matrix_view_array (data, n, n);
        !            81:
        !            82:   eval = gsl_vector_complex_alloc (n);
        !            83:   evec = gsl_matrix_complex_alloc (n, n);
        !            84:
        !            85:   w =
        !            86:     gsl_eigen_nonsymmv_alloc (n);
        !            87:
        !            88:   gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);
        !            89:
        !            90:   gsl_eigen_nonsymmv_free (w);
        !            91:
        !            92:   gsl_eigen_nonsymmv_sort (eval, evec,
        !            93:                            GSL_EIGEN_SORT_ABS_DESC);
        !            94:
        !            95:   /* show result on the server for debug */
        !            96:   {
        !            97:     int i, j;
        !            98:
        !            99:     for (i = 0; i < n; i++)
        !           100:       {
        !           101:         gsl_complex eval_i
        !           102:            = gsl_vector_complex_get (eval, i);
        !           103:         gsl_vector_complex_view evec_i
        !           104:            = gsl_matrix_complex_column (evec, i);
        !           105:
        !           106:         printf ("eigenvalue = %g + %gi\n",
        !           107:                 GSL_REAL(eval_i), GSL_IMAG(eval_i));
        !           108:         printf ("eigenvector = \n");
        !           109:         for (j = 0; j < n; ++j)
        !           110:           {
        !           111:             gsl_complex z =
        !           112:               gsl_vector_complex_get(&evec_i.vector, j);
        !           113:             printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
        !           114:           }
        !           115:       }
        !           116:   }
        !           117:
        !           118:   gsl_vector_complex_free(eval);
        !           119:   gsl_matrix_complex_free(evec);
        !           120:
        !           121:   push(gsl_vector_complex_to_cmo_array(eval,n));
        !           122:   push(gsl_matrix_complex_to_cmo_array_transposed(evec,n,n));
        !           123:
        !           124:   return;
        !           125: }