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.2

1.2     ! takayama    1: /* $OpenXM: OpenXM/src/ox_gsl/call_gsl_eigen.c,v 1.1 2019/10/23 07:00:43 takayama Exp $ */
1.1       takayama    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:
1.2     ! takayama  118:   push(gsl_vector_complex_to_cmo_array(eval,n));
        !           119:   push(gsl_matrix_complex_to_cmo_array_transposed(evec,n,n));
        !           120:
1.1       takayama  121:   gsl_vector_complex_free(eval);
                    122:   gsl_matrix_complex_free(evec);
                    123:
                    124:   return;
                    125: }