Annotation of OpenXM/src/hgm/mh/src/rk.c, Revision 1.4
1.1 takayama 1: /*
2: License: LGPL
1.2 takayama 3: Ref: Copied from this11/misc-2011/A1/wishart/Prog
1.4 ! takayama 4: $OpenXM: OpenXM/src/hgm/mh/src/rk.c,v 1.3 2013/02/20 01:06:38 takayama Exp $
1.1 takayama 5: */
6: #include <stdio.h>
7: #include <stdlib.h>
8: #include <math.h>
9: #include <string.h>
1.4 ! takayama 10: #include "sfile.h"
1.1 takayama 11:
12: extern int MH_RANK;
13: extern int MH_M;
14:
15: char *MH_Gfname;
16: char *MH_Dfname;
17: extern int MH_Mg;
1.2 takayama 18: extern double *MH_Beta; extern double *MH_Ng; extern double MH_X0g;
19: extern int MH_RawName;
20: double MH_Hg = 0.001;
1.4 ! takayama 21: static struct SFILE *Gf =NULL;
! 22: static struct SFILE *Df =NULL;
1.2 takayama 23: int MH_P95=0;
24: int MH_Verbose=0;
1.1 takayama 25: static mypower(int x,int n) {
26: int a,i;
27: a = 1;
28: for (i=0; i<n; i++) a = a*x;
29: return(a);
30: }
1.2 takayama 31: mh_gopen_file() {
1.1 takayama 32: FILE *fp;
33: char fname[1024];
34: int i;
1.4 ! takayama 35: extern int MH_byFile;
1.1 takayama 36: Gf=NULL;
37: Df=NULL;
1.2 takayama 38: if (MH_Verbose) {
1.1 takayama 39: fprintf(stderr,"MH_Gfname=%s\n",MH_Gfname);
40: fprintf(stderr,"MH_Dfname=%s\n",MH_Dfname);
41: }
42: if (MH_RANK != mypower(2,MH_Mg)) {
43: fprintf(stderr,"rk.c MH_RANK=%d is not equal to 2^MH_Mg=%d. input_data_file is broken.\n",MH_RANK,mypower(2,MH_Mg));
1.2 takayama 44: mh_usage();
45: mh_exit(-1);
1.1 takayama 46: }
47: if (MH_Dfname != NULL) {
1.2 takayama 48: if (!MH_RawName) {
49: sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"-m%d-n%lf-x0%lf-b",MH_Mg,*MH_Ng,MH_X0g);
1.1 takayama 50: for (i=0; i<MH_Mg; i++) {
1.4 ! takayama 51: sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"%lf,",MH_Beta[i]);
1.1 takayama 52: }
53: sprintf(&(MH_Dfname[strlen(MH_Dfname)]),".txt");
54: }
1.4 ! takayama 55: Df = mh_fopen(MH_Dfname,"w",MH_byFile);
1.1 takayama 56: if (Df == NULL) {
57: fprintf(stderr,"Error to open the file %s\n",MH_Dfname);
58: return(-1);
59: }
60: }
61: if (MH_Gfname != NULL) {
1.4 ! takayama 62: Gf = mh_fopen(MH_Gfname,"w",MH_byFile);
1.1 takayama 63: if (Gf == NULL) {
64: fprintf(stderr,"Error to open the file %s\n",MH_Gfname);
65: return(-1);
66: }
67: }else return(1);
1.4 ! takayama 68: if (MH_byFile) {
! 69: sprintf(fname,"%s-gp.txt",MH_Gfname);
! 70: fp = fopen(fname,"w");
! 71: fprintf(fp,"set xrange [0:20]\n");
! 72: fprintf(fp,"set yrange [0:1.2]\n");
! 73: fprintf(fp,"plot \"%s\" title \"by hg\" with lines\n",MH_Gfname);
! 74: fclose(fp);
! 75: }
1.1 takayama 76: return(0);
77: }
78: /*
79: Runge-Kutta
80: y_{n+1} = y_n + 1/6 k_1 + 1/3 k_2 + 1/3 k_3 + 1/6 k_4
81: k_1 = h f(x_n, y_n)
82: k_2 = h f(x_n + 1/2 h, y_n + 1/2 k_1)
83: k_3 = h f(x_n + 1/2 h, y_n + 1/2 k_2)
84: k_4 = h f(x_n + h, y_n + k_3)
85: */
86:
87: /* ベクトル値関数 rk is defined in tmp-code-?.c,
88: f(x, y) (x : scalar, y : vector) */
89:
1.2 takayama 90: static void show_v(double x,double *v, int n)
1.1 takayama 91: {
92: int i;
93: static int counter=0;
1.2 takayama 94: extern int MH_Dp;
1.4 ! takayama 95: char swork[MH_SSIZE];
1.1 takayama 96:
1.2 takayama 97: if ((counter % MH_Dp) != 0) { counter++; return;} else counter=1;
1.4 ! takayama 98: sprintf(swork,"%lf\n",x); mh_fputs(swork,Df);
! 99: for (i = 0; i < n; i++) {sprintf(swork," %le\n", v[i]); mh_fputs(swork,Df);}
1.1 takayama 100: }
101:
1.4 ! takayama 102: struct MH_RESULT mh_rkmain(double x0,double y0[],double xn)
1.1 takayama 103: {
104: static int initialized=0;
105: int i;
106: double h;
1.3 takayama 107: double x;
1.4 ! takayama 108: char swork[MH_SSIZE];
! 109: struct MH_RESULT result;
1.3 takayama 110: extern int MH_deallocate;
1.1 takayama 111: /*
112: double y[MH_RANK];
113: double k1[MH_RANK], k2[MH_RANK], k3[MH_RANK], k4[MH_RANK];
114: double temp[MH_RANK];
115: double ty[MH_RANK];
116: */
117: static double *y,*k1,*k2,*k3,*k4,*temp,*ty;
1.3 takayama 118: if (MH_deallocate && initialized) {
119: if (y) mh_free(y);
120: if (k1) mh_free(k1);
121: if (k2) mh_free(k2);
122: if (k3) mh_free(k3);
123: if (k4) mh_free(k4);
124: if (temp) mh_free(temp);
125: if (ty) mh_free(ty);
126: y = k1 = k2 = k3 = k4 = temp = ty = NULL;
127: initialized=0;
1.4 ! takayama 128: return(result);
1.3 takayama 129: }
1.1 takayama 130: if (!initialized) {
1.2 takayama 131: y = (double *)mh_malloc(sizeof(double)*MH_RANK);
132: k1 = (double *)mh_malloc(sizeof(double)*MH_RANK);
133: k2 = (double *)mh_malloc(sizeof(double)*MH_RANK);
134: k3 = (double *)mh_malloc(sizeof(double)*MH_RANK);
135: k4 = (double *)mh_malloc(sizeof(double)*MH_RANK);
136: temp = (double *)mh_malloc(sizeof(double)*MH_RANK);
137: ty = (double *)mh_malloc(sizeof(double)*MH_RANK);
1.1 takayama 138: initialized=1;
139: }
1.2 takayama 140: h = MH_Hg;
1.1 takayama 141: for (i = 0; i < MH_RANK; i++)
142: y[i] = y0[i];
143: for (x = x0; (h>0?(x<xn):(x>xn)); x += h) {
144: if (Df) show_v(x,y, MH_RANK);
1.4 ! takayama 145: if (Gf) {
! 146: sprintf(swork,"%lf %le\n",x,y[0]);
! 147: mh_fputs(swork,Gf);
! 148: }
1.1 takayama 149: /* Output 95% point */
1.2 takayama 150: if (MH_P95) {
151: if ((MH_P95==1) && (y[0] >= 0.95)) {
1.1 takayama 152: printf("x=%le, y[0]=%lf\n",x,y[0]);
1.2 takayama 153: MH_P95=2;
154: }else if ((MH_P95==2) && (y[0] >=0.9500000001)) {
1.1 takayama 155: printf("x=%le, y[0]=%lf\n",x,y[0]);
1.2 takayama 156: MH_P95=0;
1.1 takayama 157: }
158: }
1.2 takayama 159: mh_rf(x, y, MH_RANK, temp, MH_RANK);
1.1 takayama 160: for (i = 0; i < MH_RANK; i++)
161: k1[i] = h * temp[i];
162:
163: for (i = 0; i < MH_RANK; i++)
164: ty[i] = y[i] + 0.5 * k1[i];
1.2 takayama 165: mh_rf(x + 0.5 * h, ty, MH_RANK, temp, MH_RANK);
1.1 takayama 166: for (i = 0; i < MH_RANK; i++)
167: k2[i] = h * temp[i];
168:
169: for (i = 0; i < MH_RANK; i++)
170: ty[i] = y[i] + 0.5 * k2[i];
1.2 takayama 171: mh_rf(x + 0.5 * h, ty, MH_RANK, temp, MH_RANK);
1.1 takayama 172: for (i = 0; i < MH_RANK; i++)
173: k3[i] = h * temp[i];
174:
175: for (i = 0; i < MH_RANK; i++)
176: ty[i] = y[i] + k3[i];
1.2 takayama 177: mh_rf(x + h, ty, MH_RANK, temp, MH_RANK);
1.1 takayama 178: for (i = 0; i < MH_RANK; i++)
179: k4[i] = h * temp[i];
180:
181: for (i = 0; i < MH_RANK; i++)
182: y[i] = y[i] + 1.0/6.0 * k1[i] + 1.0/3.0 * k2[i] + 1.0/3.0 * k3[i] + 1.0/6.0 * k4[i];
183: }
184: printf("x=%lf, y[0]=%lg\n",x,y[0]);
1.4 ! takayama 185: result.x = x;
! 186: result.rank = MH_RANK;
! 187: result.y = (double *)mh_malloc(sizeof(double)*MH_RANK); /* todo, how to free it */
! 188: for (i=0; i<MH_RANK; i++) (result.y)[i] = y[i];
! 189: result.size=2;
! 190: result.sfpp = (struct SFILE **)mh_malloc(sizeof(struct SFILE *)*(result.size)); /* todo, free */
! 191: (result.sfpp)[0] = Df;
! 192: (result.sfpp)[1] = Gf;
! 193: return result;
1.1 takayama 194: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to rk.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / hgm / mh / src