Annotation of OpenXM/src/phc/phc6.c, Revision 1.2
1.1 maekawa 1: /* phc6.c , yama:1999/sm1-prog/phc6.c */
2: /* This is a simple C interface to the black-box solver of phc.
3: ** Requirements:
4: ** 1) executable version of phc will be searched in the following order:
5: ** OpenXM_HOME/bin/, /usr/local/bin/phc, /tmp/phc, your command search path.
6: ** Here, PHC_LIBDIR is an environment variable.
7: ** 2) user of this program has write permissions to create
8: ** the files "tmp.input.xx" and "tmp.output.xx" in the directory where
9: ** this program is executed. xx are numbers.
10: */
11:
12: #include <stdio.h>
13: #include <sys/stat.h>
14: #include <unistd.h>
15: #include <stdlib.h>
16:
17: /* Definition of class identifiers. */
18: #define Snull 0
19: #define SstringObject 5
20: #define Sarray 6
21: #define SlongdoubleComplex 601
22:
23: /* Definition of Object */
24: union cell {
25: int ival;
26: char *str;
27: struct phc_object *op;
28: long double longdouble;
29: };
30: struct phc_object{
31: int tag; /* class identifier */
32: union cell lc; /* left cell */
33: union cell rc; /* right cell */
34: };
35:
36: /* Memory allocation function.
37: Use your favorite memory allocation function.
38: I recommend not to use malloc and to use gc4.14 for large applications. */
39: #define sGC_malloc(n) malloc(n)
40:
41: /********** macros to use Sarray **************/
42: /* put to Object Array */
43: #define phc_putoa(ob,i,cc) {\
44: if ((ob).tag != Sarray) {fprintf(stderr,"Warning: PUTOA is for an array of objects\n");} else \
45: {if ((0 <= (i)) && ((i) < (ob).lc.ival)) {\
46: (ob.rc.op)[i] = cc;\
47: }else{\
48: fprintf(stderr,"Warning: PUTOA, the size is %d.\n",(ob).lc.ival);\
49: }}}
50: #define phc_getoa(ob,i) ((ob.rc.op)[i])
51: #define phc_getoaSize(ob) ((ob).lc.ival)
52:
53: /* prototypes */
54: struct phc_object phc_newObjectArray(int size);
55: void phc_printObject(FILE *fp,struct phc_object ob);
56: char *phc_generateUniqueFileName(char *s);
57: char *phc_which(char *s); /* search a path for the file s */
58: struct phc_object phc_complexTo(long double r, long double i);
59:
60:
1.2 ! takayama 61: int phc_scan_for_string(FILE *fp, char str[]);
1.1 maekawa 62: struct phc_object phc_scan_solutions(FILE *fp, int npaths, int dim );
63: struct phc_object phc_scan_output_of_phc(char *fname);
64: struct phc_object phc_call_phc(char *sys);
65:
66: int phc_verbose = 0;
67: int phc_overwrite = 1; /* Always use tmp.input.0 and tmp.output.0
68: for work files. */
69:
70: main(int argc, char *argv[]) {
71: struct phc_object ob;
72: int n,i,dim;
73: #define INPUTSIZE 4096
74: char input[INPUTSIZE];
75: #define A_SIZE 1024
76: char a[A_SIZE];
77: int message = 0;
78: for (i=1; i<argc; i++) {
79: if (strcmp(argv[i],"-v") == 0) {
1.2 ! takayama 80: phc_verbose = 1; message=1;
1.1 maekawa 81: }else if (strcmp(argv[i],"-g") == 0) {
82: phc_overwrite = 0;
83: }else if (strcmp(argv[i],"-i") == 0) {
84: ob = phc_call_phc(argv[i+1]);
85: n = phc_getoaSize(ob);
86: printf("[\n");
87: for (i=0; i<n; i++) {
88: phc_printObject(stdout,phc_getoa(ob,i));
89: if (i != n-1) printf(" ,\n"); else printf(" \n");
90: }
91: printf("]\n");
92: exit(0);
93: }
94: }
95: if (message) {
96: printf("Input example:\n 2 \n x**2 + y**2 - 1;\n x**2 + y**2 - 8*x - 3;\n");
97: printf("Note that input length is limited.\n");
98: }
99: while (1) {
100: if (message) printf("dim= ");
1.2 ! takayama 101: if (fgets(input,INPUTSIZE,stdin) <= 0) break;
! 102: sscanf(input,"%d",&dim);
1.1 maekawa 103: sprintf(input,"%d\n",dim);
104: if (message) printf("Input %d equations please.\n",dim);
105: for (i=0; i<dim; i++) {
106: if (message) {printf("eq[%d] = ",i); fflush(stdout);}
107: do {
108: fgets(a,A_SIZE-1, stdin);
109: } while (strlen(a) == 0);
110: if (strlen(a) >= A_SIZE-3) {
111: fprintf(stderr,"Too big input for the input buffer a.\n"); exit(10);
112: }
113: if (strlen(input)+strlen(a) >= INPUTSIZE) {
114: fprintf(stderr,"Too big input for the input buffer input.\n"); exit(10);
115: }
116: sprintf(input+strlen(input),"%s\n",a);
117: }
118: ob = phc_call_phc(input);
119: if (message) {
120: printf("-----------------------------------------------------------\n");
121: }
122: n = phc_getoaSize(ob);
123: printf("[\n");
124: for (i=0; i<n; i++) {
125: phc_printObject(stdout,phc_getoa(ob,i));
126: if (i != n-1) printf(" ,\n"); else printf(" \n");
127: }
128: printf("]\n");
129: }
130: }
131:
1.2 ! takayama 132: int phc_scan_for_string(FILE *fp, char str[])
1.1 maekawa 133: /*
134: ** Scans the file fp for a certain string str of length lenstr+1.
135: ** Reading stops when the string has been found, then the variable
136: ** on return equals 1, otherwise 0 is returned.
137: */
138: {
1.2 ! takayama 139: #define BUF_SIZE 1024
! 140: char buf[BUF_SIZE];
1.1 maekawa 141: char ch;
142: int index,i,compare,npaths,dim,found;
1.2 ! takayama 143: int lenstr;
! 144: lenstr = strlen(str);
! 145: if (lenstr >= BUF_SIZE-1) {
! 146: fprintf(stderr,"Too long string in phc_scan_for_string\n");
! 147: exit(-1);
! 148: }
1.1 maekawa 149: index = -1;
150: found = 0;
1.2 ! takayama 151: while (((ch = fgetc(fp))!=EOF) && found == 0)
1.1 maekawa 152: {
1.2 ! takayama 153: if (index == -1 && ch == str[0])
1.1 maekawa 154: {
1.2 ! takayama 155: index = 0;
! 156: buf[index] = ch;
1.1 maekawa 157: }
1.2 ! takayama 158: else
1.1 maekawa 159: {
1.2 ! takayama 160: if (index == lenstr)
! 161: {
! 162: compare = 0;
! 163: for (i=0; str[i] != '\0'; i++)
! 164: {
! 165: if (buf[i]!=str[i])
! 166: {
! 167: compare = compare+1;
! 168: }
! 169: }
! 170: if (compare == 0)
! 171: {
! 172: found = 1;
! 173: }
! 174: index = -1;
! 175: }
! 176: else
! 177: if (index > -1 && index < lenstr)
! 178: {
! 179: index = index+1;
! 180: buf[index] = ch;
! 181: }
1.1 maekawa 182: }
1.2 ! takayama 183: if (found == 1) break;
1.1 maekawa 184: }
185: return found;
186: }
187: struct phc_object phc_scan_solutions(FILE *fp, int npaths, int dim )
188: /*
189: ** Scans the file for the solutions, from a list of length npaths,
190: ** of complex vectors with dim entries.
191: ** The tolerance for the residual to a solution is set to 1.0E-12.
192: ** Returns solutions.
193: */
194: {
195: struct phc_object rob,sob;
196: char ch;
197: int fnd,i,j,nsols;
1.2 ! takayama 198: double res;
1.1 maekawa 199: long double realpart;
200: long double imagpart;
201: long double realparts[npaths][dim];
202: long double imagparts[npaths][dim];
203: nsols = 0;
1.2 ! takayama 204: while ((ch = fgetc(fp)) != EOF)
1.1 maekawa 205: {
1.2 ! takayama 206: fnd = phc_scan_for_string(fp,"start residual :");
1.1 maekawa 207: if (fnd==1)
208: {
209: fscanf(fp,"%E",&res);
1.2 ! takayama 210: /* printf(" residual = "); printf("%E\n",res); */
! 211: if (res < 1.0E-12) {
! 212: nsols = nsols+1;
! 213: if (nsols > npaths) {
! 214: fprintf(stderr,"Something is wrong in phc_scan_solutions\n");
! 215: fprintf(stderr,"npaths=%d, nsols=%d \n",npaths,nsols);
! 216: exit(-1);
! 217: }
! 218: }
! 219: fnd = phc_scan_for_string(fp,"the solution for t :");
1.1 maekawa 220: for (i=0;i<dim;i++)
221: {
1.2 ! takayama 222: fnd = phc_scan_for_string(fp,":");
1.1 maekawa 223: fscanf(fp,"%LE",&realpart);
224: fscanf(fp,"%LE",&imagpart);
225: if (res < 1.0E-12)
226: {
227: realparts[nsols-1][i] = realpart;
228: imagparts[nsols-1][i] = imagpart;
229: }
230: }
231: }
232: }
233: if(phc_verbose) fprintf(stderr," number of solutions = %i\n",nsols);
234: rob = phc_newObjectArray(nsols);
235: for (i=0;i<nsols;i++)
236: {
237: /* fprintf(stderr,"Solution %i :\n",i+1); */
238: sob = phc_newObjectArray(dim);
239: for (j=0;j<dim;j++)
240: {
241: /*
242: printf("%20.14LE",realparts[i][j]); printf(" ");
243: printf("%20.14LE",imagparts[i][j]); printf("\n");
244: */
245: phc_putoa(sob,j,phc_complexTo(realparts[i][j],imagparts[i][j]));
246: }
247: phc_putoa(rob,i,sob);
248: }
249: return(rob);
250: }
251: struct phc_object phc_scan_output_of_phc(char *fname)
252: /*
253: ** Scans the file "output" of phc in two stages to get
254: ** 1) the number of paths and the dimension;
255: ** 2) the solutions, vectors with residuals < 1.0E-12.
256: */
257: {
258: FILE *otp;
259: char ch;
260: int fnd,npaths,dim,i,nsols;
261: otp = fopen(fname,"r");
262: if (phc_verbose) fprintf(stderr,"Scanning the %s of phc.\n",fname);
1.2 ! takayama 263: fnd = phc_scan_for_string(otp,"THE SOLUTIONS :");
1.1 maekawa 264: fscanf(otp,"%i",&npaths);
265: if (phc_verbose) fprintf(stderr," number of paths traced = %i\n",npaths);
266: fscanf(otp,"%i",&dim);
267: if (phc_verbose) fprintf(stderr," dimension of solutions = %i\n",dim);
268: return(phc_scan_solutions(otp,npaths,dim));
269: }
270: struct phc_object phc_call_phc(char *sys) /* call phc, system as string */
271: {
272: FILE *inp;
273: char *f,*outf;
274: char cmd[1024];
275: char *w;
276: struct phc_object phc_NullObject ;
277: struct stat statbuf;
278:
279: phc_NullObject.tag = Snull;
280: f = phc_generateUniqueFileName("tmp.input");
281: if (phc_verbose) fprintf(stderr,"Creating file with name %s.\n",f);
282: outf = phc_generateUniqueFileName("tmp.output");
283: if (stat(outf,&statbuf) == 0) {
284: sprintf(cmd,"/bin/rm -f %s",outf);
285: system(cmd);
286: }
287: inp = fopen(f,"w");
288: fprintf(inp,sys);
289: fclose(inp);
290: if ((w = phc_which("phc")) != NULL) {
291: sprintf(cmd,"%s -b %s %s",w,f,outf);
292: }else{
293: sprintf(cmd,"phc -b %s %s",f,outf);
294: }
295: if (phc_verbose)fprintf(stderr,"Calling %s, black-box solver of phc.\n",cmd);
296: system(cmd);
297: if (stat(outf,&statbuf) < 0) {
298: fprintf(stderr,"execution error of phc.\n");
299: return(phc_NullObject);
300: }else{
301: if (phc_verbose) fprintf(stderr,"See the file %s for results.\n",outf);
302: return(phc_scan_output_of_phc(outf));
303: }
304: }
305:
306:
307: struct phc_object phc_newObjectArray(size)
308: int size;
309: {
310: struct phc_object rob;
311: struct phc_object *op;
312: if (size > 0) {
313: op = (struct phc_object *)sGC_malloc(size*sizeof(struct phc_object));
314: if (op == (struct phc_object *)NULL) {fprintf(stderr,"No memory\n");exit(1);}
315: }else{
316: op = (struct phc_object *)NULL;
317: }
318: rob.tag = Sarray;
319: rob.lc.ival = size;
320: rob.rc.op = op;
321: return(rob);
322: }
323:
324: void phc_printObject(FILE *fp,struct phc_object ob)
325: {
326: int n,i;
327: if (ob.tag == Snull) {
328: fprintf(fp,"null");
329: }else if (ob.tag == SstringObject) {
330: fprintf(fp,"%s",ob.lc.str);
331: }else if (ob.tag == Sarray) {
332: n = phc_getoaSize(ob);
333: fprintf(fp,"[");
334: for (i=0; i<n; i++) {
335: phc_printObject(fp,phc_getoa(ob,i));
336: if (i <n-1) fprintf(fp," , ");
337: }
338: fprintf(fp,"]");
339: }else if (ob.tag == SlongdoubleComplex) {
340: /* Try your favorite way to print complex numbers. */
341: /*fprintf(fp,"(%20.14LE)+I*(%20.14LE)",ob.lc.longdouble,ob.rc.longdouble);*/
342: fprintf(fp,"[%Lf , %Lf]",ob.lc.longdouble,ob.rc.longdouble);
343: }else{
344: fprintf(stderr,"Unknown phc_object tag %d",ob.tag);
345: }
346: }
347:
348:
349: char *phc_generateUniqueFileName(char *s)
350: {
351: char *t;
352: int i;
353: struct stat statbuf;
354: t = (char *)sGC_malloc(sizeof(char)*strlen(s)+4+2);
355: for (i=0; i<1000; i++) {
356: /* Give up if we failed for 1000 names. */
357: sprintf(t,"%s.%d",s,i);
358: if (phc_overwrite) return(t);
359: if (stat(t,&statbuf) < 0) {
360: return(t);
361: }
362: }
363: fprintf(stderr,"Could not generate a unique file name.");
364: return(NULL);
365: }
366:
367: char *phc_which(char *s) {
368: struct stat statbuf;
369: char *cmd,*a;
370: a = getenv("OpenXM_HOME");
371: if (a != NULL) {
372: cmd = (char *) sGC_malloc(sizeof(char)*(strlen(s)+strlen(a)
373: +strlen("/usr/local/bin/")+1));
374: if (cmd == NULL) {fprintf(stderr,"No more memory.\n"); exit(1);}
375: strcpy(cmd,a); strcat(cmd,"/bin/"); strcat(cmd,s);
376: if (stat(cmd,&statbuf) == 0) {
377: return(cmd);
378: }
379: }
380: cmd = (char *) sGC_malloc(sizeof(char)*(strlen(s)
381: +strlen("/usr/local/bin/")+1));
382: if (cmd == NULL) {fprintf(stderr,"No more memory.\n"); exit(1);}
383: strcpy(cmd,"/usr/local/bin/"); strcat(cmd,s);
384: if (stat(cmd,&statbuf) == 0) {
385: return(cmd);
386: }
387: strcpy(cmd,"/tmp/"); strcat(cmd,s);
388: if (stat(cmd,&statbuf) == 0) {
389: return(cmd);
390: }
391: return(NULL);
392: }
393:
394:
395: struct phc_object phc_complexTo(long double r, long double i)
396: {
397: struct phc_object rob;
398: rob.tag = SlongdoubleComplex;
399: rob.lc.longdouble = r;
400: rob.rc.longdouble = i;
401: return(rob);
402: }
403:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to phc6.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / phc