Annotation of OpenXM/src/ox_cdd/ox_cdd.c, Revision 1.2
1.2 ! noro 1: /* $OpenXM: OpenXM/src/ox_cdd/ox_cdd.c,v 1.1 2005/05/25 04:42:20 noro Exp $ */
1.1 noro 2:
3: #include <stdio.h>
4: #include <stdlib.h>
5: #include <string.h>
6:
7: #if defined GMPRATIONAL
8: #include "gmp.h"
9: #define __GMP_FAKE_H__
10: #else
11: #define dd_almostzero 1.0E-7
12: #include <math.h>
13: #endif
14:
15: #include "ox_toolkit.h"
16:
17: #if defined GMPRATIONAL
18: typedef mpq_t mytype;
19: #else /* built-in C double */
20: typedef double mytype[1];
21: #endif
22:
23: typedef enum {
24: dd_LPnone=0, dd_LPmax, dd_LPmin
25: } dd_LPObjectiveType;
26:
27: OXFILE *fd_rw;
28:
29: #define INIT_S_SIZE 2048
30: #define EXT_S_SIZE 2048
31:
32: static int stack_size = 0;
33: static int stack_pointer = 0;
34: static cmo **stack = NULL;
35: extern int debug_print;
36:
37: void init_cdd(void);
38: int **redcheck(int row,int col,int **matrix,int *presult_row);
39: mytype *lpsolve(dd_LPObjectiveType type,int row,int col,int **matrix,int *obj);
1.2 ! noro 40: mytype *lpintpt(int row, int col, int **matrix);
! 41:
1.1 noro 42:
43: void dprint( int level, char *string )
44: {
45: if( debug_print >= level ){
46: printf( string );
47: fflush( stdout );
48: }
49: }
50:
51: #define LP_Q 0
52: #define LP_I 1
53: #define LP_MAX 0
54: #define LP_MIN 1
55: #define LP_MAXMIN 2
1.2 ! noro 56: #define LP_INTPT 3
1.1 noro 57:
58: int initialize_stack()
59: {
60: stack_pointer = 0;
61: stack_size = INIT_S_SIZE;
62: stack = MALLOC(stack_size*sizeof(cmo*));
63: return 0;
64: }
65:
66: static int extend_stack()
67: {
68: int size2 = stack_size + EXT_S_SIZE;
69: cmo **stack2 = MALLOC(size2*sizeof(cmo*));
70: memcpy(stack2, stack, stack_size*sizeof(cmo *));
71: free(stack);
72: stack = stack2;
73: stack_size = size2;
74: return 0;
75: }
76:
77: int push(cmo* m)
78: {
79: stack[stack_pointer] = m;
80: stack_pointer++;
81: if(stack_pointer >= stack_size) {
82: extend_stack();
83: }
84: return 0;
85: }
86:
87: cmo* pop()
88: {
89: if(stack_pointer > 0) {
90: stack_pointer--;
91: return stack[stack_pointer];
92: }
93: return new_cmo_null();
94: }
95:
96: void pops(int n)
97: {
98: stack_pointer -= n;
99: if(stack_pointer < 0) {
100: stack_pointer = 0;
101: }
102: }
103:
104: #define OX_CDD_VERSION 0x11121400
105: #define ID_STRING "2005/02/14 06:50:00"
106:
107: int sm_mathcap()
108: {
109: mathcap_init(OX_CDD_VERSION, ID_STRING, "ox_cdd", NULL, NULL);
110: push((cmo*)oxf_cmo_mathcap(fd_rw));
111: return 0;
112: }
113:
114: int sm_popCMO()
115: {
116: cmo* m = pop();
117:
118: if(m != NULL) {
119: send_ox_cmo(fd_rw, m);
120: return 0;
121: }
122: return SM_popCMO;
123: }
124:
125: cmo_error2 *make_error2(int code)
126: {
127: return (cmo_error2 *) new_cmo_int32(-1);
128: }
129:
130: int get_i()
131: {
132: cmo *c = pop();
133: if(c->tag == CMO_INT32) {
134: return ((cmo_int32 *)c)->i;
135: }else if(c->tag == CMO_ZZ) {
136: return mpz_get_si(((cmo_zz *)c)->mpz);
137: }
138: make_error2(-1);
139: return 0;
140: }
141:
142: char *get_str()
143: {
144: cmo *c = pop();
145: if(c->tag == CMO_STRING) {
146: return ((cmo_string *)c)->s;
147: }
148: make_error2(-1);
149: return "";
150: }
151:
152: int cmo2int(cmo *c)
153: {
154: if(c->tag == CMO_INT32) {
155: return ((cmo_int32 *)c)->i;
156: }else if(c->tag == CMO_ZZ) {
157: return mpz_get_si(((cmo_zz *)c)->mpz);
158: } else if(c->tag == CMO_NULL){
159: return 0;
160: }
161:
162: return 0;
163: }
164:
165: cmo *matrix2cmo( int **matrix, int row, int col )
166: {
167: cmo_list *result,*work;
168: cmo *tmp;
169: int i,j;
170:
171: result = new_cmo_list();
172: for(i=0;i<row;i++){
173:
174: work = new_cmo_list();
175: for(j=0;j<col;j++){
176: if( matrix[i][j] != 0 ){
177: tmp = (cmo*) new_cmo_zz_set_si( matrix[i][j] );
178: } else {
179: tmp = (cmo*) new_cmo_zero();
180: }
181: work = list_append( work, tmp );
182: }
183: result = list_append( result, (cmo*)work );
184:
185: }
186:
187: return (cmo *)result;
188: }
189:
1.2 ! noro 190: cmo_qq* new_cmo_qq_set_mpq(mpq_ptr q);
! 191:
! 192: cmo *vector2cmo( mytype *vector, int len)
! 193: {
! 194: cmo_list *result;
! 195: cmo *tmp;
! 196: int j;
! 197:
! 198: result = new_cmo_list();
! 199: for(j=0;j<len;j++){
! 200: if( vector[j] != 0 ){
! 201: tmp = (cmo*) new_cmo_qq_set_mpq( (mpq_ptr)vector[j] );
! 202: } else {
! 203: tmp = (cmo*) new_cmo_zero();
! 204: }
! 205: result = list_append( result, tmp );
! 206: }
! 207:
! 208: return (cmo *)result;
! 209: }
! 210:
1.1 noro 211: int mytype2int( mytype a, int i )
212: {
213: #if defined GMPRATIONAL
214: /* How to round "a" depends on "i".
215: if i < 0, then the rounding style is "floor".
216: if i = 0, then the rounding style is "truncate".
217: if i > 0, then the rounding style is "ceil". */
218: int ret;
219: mpz_t q;
220: mpz_init(q);
221:
222: if( i < 0 ){
223: mpz_fdiv_q( q, mpq_numref(a), mpq_denref(a) );
224: } else if( i > 0 ){
225: mpz_cdiv_q( q, mpq_numref(a), mpq_denref(a) );
226: } else {
227: mpz_tdiv_q( q, mpq_numref(a), mpq_denref(a) );
228: }
229:
230: ret = mpz_get_si( q );
231: mpz_clear( q );
232: return ret;
233: #else
234: /* How to round "a" depends on "i".
235: if i < 0, then the rounding style is "floor".
236: if i = 0, then the rounding style is "near".
237: if i > 0, then the rounding style is "ceil". */
238: if( i < 0 ){
239: return floor( *a + dd_almostzero );
240: } else if( i == 0 ){
241: return (int)( *a + 0.5 );
242: } else {
243: return ceil( *a - dd_almostzero );
244: }
245: #endif
246: }
247:
248: void my_redcheck(void)
249: {
250: int row,col,result_row;
251: int i,j;
252: cmo *a,*b,*c;
253: int **matrix,**result;
254:
255: pop(); /* for argc */
256: row = get_i(); /* row size */
257: col = get_i(); /* col size */
258:
259: a = pop();
260:
261: matrix = MALLOC( row * sizeof(int*) );
262: for(i=0;i<row;i++){
263: matrix[i] = MALLOC( col * sizeof(int) );
264: }
265:
266: for(i=0;i<row;i++){
267: b = list_nth( (cmo_list*)a, i );
268: for(j=0;j<col;j++){
269: c = list_nth( (cmo_list*)b, j );
270: matrix[i][j] = cmo2int( c );
271: }
272: }
273:
274: dprint(1,"redcheck...");
275: result = redcheck(row,col,matrix,&result_row);
276: dprint(1,"done\n");
277:
278: push( matrix2cmo(result,result_row,col) );
279: }
280:
281: void my_lpsolve( int resulttype, int index )
282: {
283: int row,col;
284: int i,j;
285: cmo *a,*b,*c;
286: cmo_list* ret;
287: int **matrix,*object;
288: int lpmax,lpmin;
289: mytype *tmp;
1.2 ! noro 290: cmo *cmomin,*cmomax,*cmoint;
1.1 noro 291:
292: pop(); /* for argc */
293: row = get_i(); /* row size */
294: col = get_i(); /* col size */
295:
296: matrix = MALLOC( row * sizeof(int*) );
297: for(i=0;i<row;i++){
298: matrix[i] = MALLOC( col * sizeof(int) );
299: }
300:
301: a = pop();
302:
303: /* For matrix */
304: for(i=0;i<row;i++){
305: b = list_nth( (cmo_list*)a, i );
306: for(j=0;j<col;j++){
307: c = list_nth( (cmo_list*)b, j );
308: matrix[i][j] = cmo2int( c );
309: }
310: }
311:
1.2 ! noro 312: if ( index != LP_INTPT ) {
! 313: /* For object */
! 314: object = MALLOC( col * sizeof(int) );
! 315: a = pop();
! 316:
! 317: for(i=0;i<col;i++){
! 318: c = list_nth( (cmo_list*)a, i );
! 319:
! 320: object[i] = cmo2int( c );
! 321: }
1.1 noro 322: }
323:
324: if( index == LP_MAX ){
325: dprint( 1, "lpsolve(max)...");
326: } else if( index == LP_MIN ){
327: dprint( 1, "lpsolve(min)...");
328: } else if( LP_MAXMIN ){
329: dprint( 1, "lpsolve(maxmin)...");
1.2 ! noro 330: } else if( LP_INTPT ){
! 331: dprint( 1, "lpsolve(intpt)...");
1.1 noro 332: }
333:
334: if( index == LP_MAX || index == LP_MAXMIN ){
335: tmp = lpsolve( dd_LPmax, row, col, matrix, object );
336: if( resulttype == LP_I ){
337: lpmax = mytype2int( *tmp, -1 );
338: if( lpmax == 0 ){
339: cmomax = (cmo*) new_cmo_zero();
340: } else {
341: cmomax = (cmo*) new_cmo_zz_set_si( lpmax );
342: }
343: } else if( resulttype == LP_Q ){
344: #if defined GMPRATIONAL
345: if( mpz_cmp_si( mpq_numref( *tmp ), 0 ) == 0 ){
346: cmomax = (cmo*) new_cmo_zero();
347: } else {
348: cmomax = (cmo*) new_cmo_qq_set_mpz( mpq_numref( *tmp ), mpq_denref( *tmp ) );
349: }
350: #else
351: cmomax = (cmo*) new_cmo_double( *tmp[0] );
352: #endif
353: }
354: #if defined GMPRATIONAL
355: mpq_clear( *tmp );
356: #endif
357: }
358:
359: if( index == LP_MIN || index == LP_MAXMIN ){
360: tmp = lpsolve( dd_LPmin, row, col, matrix, object );
361: if( resulttype == LP_I ){
362: lpmin = mytype2int( *tmp, 1 );
363: if( lpmin == 0 ){
364: cmomin = (cmo*) new_cmo_zero();
365: } else {
366: cmomin = (cmo*) new_cmo_zz_set_si( lpmin );
367: }
368: } else if( resulttype == LP_Q ){
369: #if defined GMPRATIONAL
370: if( mpz_cmp_si( mpq_numref( *tmp ), 0 ) == 0 ){
371: cmomin = (cmo*) new_cmo_zero();
372: } else {
373: cmomin = (cmo*) new_cmo_qq_set_mpz( mpq_numref( *tmp ), mpq_denref( *tmp ) );
374: }
375: #else
376: cmomin = (cmo*) new_cmo_double( *tmp[0] );
377: #endif
378: }
379: #if defined GMPRATIONAL
380: mpq_clear( *tmp );
381: #endif
382:
383: }
384:
1.2 ! noro 385: if( index == LP_INTPT ) {
! 386: tmp = lpintpt(row, col, matrix);
! 387: if ( tmp )
! 388: cmoint = vector2cmo(tmp,col);
! 389: else
! 390: cmoint = new_cmo_zero();
! 391: for ( i = 0; i < col; i++ )
! 392: mpq_clear( (mpq_ptr)tmp[i] );
! 393: }
! 394:
1.1 noro 395: if( index == LP_MAX ){
396: push( cmomax );
397: } else if( index == LP_MIN ){
398: push( cmomin );
399: } else if( index == LP_MAXMIN ){
400: ret = new_cmo_list();
401: ret = list_append( ret, cmomin );
402: ret = list_append( ret, cmomax );
403: push( (cmo*) ret );
1.2 ! noro 404: } else if ( index == LP_INTPT )
! 405: push( cmoint );
1.1 noro 406: }
407:
408: int sm_executeFunction()
409: {
410: cmo_string *func = (cmo_string *)pop();
411: if(func->tag != CMO_STRING) {
412: printf("sm_executeFunction : func->tag is not CMO_STRING");fflush(stdout);
413: push((cmo*)make_error2(0));
414: return -1;
415: }
416:
417: if(strcmp(func->s, "redcheck" ) == 0) {
418: my_redcheck();
419: return 0;
420: } else if( strcmp( func->s, "ilpmax" ) == 0 ){
421: my_lpsolve(LP_I,LP_MAX);
422: return 0;
423: } else if( strcmp( func->s, "ilpmin" ) == 0 ){
424: my_lpsolve(LP_I,LP_MIN);
425: return 0;
426: } else if( strcmp( func->s, "ilpmaxmin" ) == 0 ){
427: my_lpsolve(LP_I,LP_MAXMIN);
428: return 0;
429: } else if( strcmp( func->s, "lpmax" ) == 0 ){
430: my_lpsolve(LP_Q,LP_MAX);
431: return 0;
432: } else if( strcmp( func->s, "lpmin" ) == 0 ){
433: my_lpsolve(LP_Q,LP_MIN);
434: return 0;
435: } else if( strcmp( func->s, "lpmaxmin" ) == 0 ){
436: my_lpsolve(LP_Q,LP_MAXMIN);
437: return 0;
1.2 ! noro 438: #if defined GMPRATIONAL
! 439: } else if( strcmp( func->s, "intpt" ) == 0 ){
! 440: my_lpsolve(LP_Q,LP_INTPT);
! 441: return 0;
! 442: #endif
1.1 noro 443: } else if( strcmp( func->s, "debugprint" ) == 0 ){
444: pop();
445: debug_print = get_i();
446: return 0;
447: } else if( strcmp( func->s, "exit" ) == 0 ){
448: pop();
449: exit(0);
450: return 0;
451: } else {
452: push((cmo*)make_error2(0));
453: return -1;
454: }
455: }
456:
457: int receive_and_execute_sm_command()
458: {
459: int code = receive_int32(fd_rw);
460: switch(code) {
461: case SM_popCMO:
462: sm_popCMO();
463: break;
464: case SM_executeFunction:
465: sm_executeFunction();
466: break;
467: case SM_mathcap:
468: sm_mathcap();
469: break;
470: case SM_setMathCap:
471: pop();
472: break;
473: default:
474: printf("receive_and_execute_sm_command : code=%d\n",code);fflush(stdout);
475: break;
476: }
477: return 0;
478: }
479:
480: int receive()
481: {
482: int tag;
483:
484: tag = receive_ox_tag(fd_rw);
485: switch(tag) {
486: case OX_DATA:
487: push(receive_cmo(fd_rw));
488: break;
489: case OX_COMMAND:
490: receive_and_execute_sm_command();
491: break;
492: default:
493: printf("receive : tag=%d\n",tag);fflush(stdout);
494: }
495: return 0;
496: }
497:
498: int main()
499: {
500: GC_INIT();
501: ox_stderr_init(stderr);
502: initialize_stack();
503: init_cdd();
504:
505: #if defined GMPRATIONAL
506: fprintf(stderr,"ox_cddgmp\n");
507: #else
508: fprintf(stderr,"ox_cdd\n");
509: #endif
510:
511: fd_rw = oxf_open(3);
512: oxf_determine_byteorder_server(fd_rw);
513:
514: while(1){
515: receive();
516: }
517: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>