Return to gsl.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

Fixed a bug.

/*$OpenXM: OpenXM/src/asir-contrib/packages/src/gsl.rr,v 1.2 2019/10/23 07:28:11 takayama Exp $ */

module gsl;

/* ------------ list of local functions ---------- */
localf init_gsl_proc$
localf which_ox_gsl$
localf get_gsl_proc$
localf set_gsl_proc$
localf find_proc$
localf start$
localf start_unix$
localf start_windows$

localf integration_qags$
localf test1$
localf to_double$
localf eigen_nonsymmv$
localf test2$
localf list_to_complex$
/*  ------------ static variables ---------------- */
static GSL_proc$
def init_gsl_proc() {
  GSL_proc = -1$
}

/* Functions to get and set static variables out of the module */
def get_gsl_proc() {
  return GSL_proc;
}

def set_gsl_proc(A) {
  GSL_proc = A;
}

#define    GSL_FIND_PROC(P)  P = getopt(proc);\
                          if (type(P) == -1) {\
                             P = find_proc();\
                          }

def which_ox_gsl() {
  L = ox_get_serverinfo();
  for (I=0; I<length(L); I++) {  
    if (base_position("ox_gsl",L[I][1][0])>0) return L[I][0];
  }
  return -1;
}

def find_proc() {
  /*! extern GSL_proc; */
  if (GSL_proc == -1) {
     GSL_proc = which_ox_gsl();
     if (GSL_proc < 0) GSL_proc = start();
  }
  return(GSL_proc);
}
/* Search : oxPrintMessage, see cmoDebugCMO, too. */
/************** end of configure *******************************/

def start() {
  extern Xm_unix;
  if (ox_ostype()[0] == "windows" && Xm_unix == 0)
    return start_windows(0);
  else
    return start_unix();
}

def start_unix() {
 extern Xm_noX;
 if (Xm_noX) {
   P = ox_launch_nox(0,getenv("OpenXM_HOME")+"/src/ox_gsl/ox_gsl");
 }else{
   P = ox_launch(0,getenv("OpenXM_HOME")+"/src/ox_gsl/ox_gsl");
 }
// ox_check_errors(P);  //workaround to avoid mathcap error.
 GSL_proc = P;
 return(P);
}

def start_windows(U) {
  error("Not implemeted.");
}

def integration_qags(F,Low,Top) {
  GSL_FIND_PROC(Pid);
  ox_cmo_rpc(Pid,"gsl_integration_qags",F,deval(Low*1.0),deval(Top*1.0));
  Ans=ox_pop_cmo(Pid);
  return Ans;
}
def test1() {
  integration_qags(quote(log(x)/x^(1/2)),0.01,1);
}

def list_to_complex(L,Level) {
  if (type(L)<4) return L;
  if (Level == 0) return L[0]+@i*L[1];
  if (Level >= 1) return map(list_to_complex,L,Level-1);
}
def to_double(L) {
  if (L==0) return(0.0);
  if (type(L)<4) return(deval(1.0*L));
  return map(to_double,L);
}
def eigen_nonsymmv(Mat) {
  GSL_FIND_PROC(Pid);
  Mat_double = matrix_list_to_matrix(Mat);
  M=size(Mat_double)[0]; N=size(Mat_double)[1];
  if (M != N) error("eigen_nonsymmv accepts only square matrix.");
  Mat_double=to_double(Mat_double);
  print(Mat_double);
  ox_cmo_rpc(Pid,"gsl_eigen_nonsymmv",base_flatten(Mat_double));
  Ans2=ox_pop_cmo(Pid);
  Ans1=ox_pop_cmo(Pid);
  return [list_to_complex(Ans1,1),list_to_complex(Ans2,2)];
}
def test2() {
  M=[[6.79502,-43.2617,20.7214,0,0,0],[12.7671,-67.4408,21.0846,0,0,0],[45.5035,-358.137,226.992,0,0,0],[0,0,0,6.79502,-43.2617,20.7214],[0,0,0,12.7671,-67.4408,21.0846],[0,0,0,45.5035,-358.137,226.992]];
  // M=[[1,0],[2,3]];
  M = matrix_list_to_matrix(M);
  Ans=eigen_nonsymmv(M);
  N = size(M)[0];
  printf("eigenvalues=%a\neigenvectors=%a\n",Ans[0],Ans[1]);
  for (I=0; I<N; I++) {
    printf("Is it zero? %a\n",(M-Ans[0][I]*matrix_identity_matrix(N))*matrix_list_to_matrix(Ans[1][I]));
  }
  return Ans;
}

endmodule;

gsl.init_gsl_proc()$

end$