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