Annotation of OpenXM/src/ox_cdd/ox_cdd.c, Revision 1.3
1.3 ! noro 1: /* $OpenXM: OpenXM/src/ox_cdd/ox_cdd.c,v 1.2 2007/09/12 07:07:36 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:
1.3 ! noro 281: cmo_qq* new_cmo_qq_set_mpq(mpq_ptr q);
! 282:
1.1 noro 283: void my_lpsolve( int resulttype, int index )
284: {
285: int row,col;
286: int i,j;
287: cmo *a,*b,*c;
288: cmo_list* ret;
289: int **matrix,*object;
290: int lpmax,lpmin;
291: mytype *tmp;
1.2 noro 292: cmo *cmomin,*cmomax,*cmoint;
1.1 noro 293:
294: pop(); /* for argc */
295: row = get_i(); /* row size */
296: col = get_i(); /* col size */
297:
298: matrix = MALLOC( row * sizeof(int*) );
299: for(i=0;i<row;i++){
300: matrix[i] = MALLOC( col * sizeof(int) );
301: }
302:
303: a = pop();
304:
305: /* For matrix */
306: for(i=0;i<row;i++){
307: b = list_nth( (cmo_list*)a, i );
308: for(j=0;j<col;j++){
309: c = list_nth( (cmo_list*)b, j );
310: matrix[i][j] = cmo2int( c );
311: }
312: }
313:
1.2 noro 314: if ( index != LP_INTPT ) {
315: /* For object */
316: object = MALLOC( col * sizeof(int) );
317: a = pop();
318:
319: for(i=0;i<col;i++){
320: c = list_nth( (cmo_list*)a, i );
321:
322: object[i] = cmo2int( c );
323: }
1.1 noro 324: }
325:
326: if( index == LP_MAX ){
327: dprint( 1, "lpsolve(max)...");
328: } else if( index == LP_MIN ){
329: dprint( 1, "lpsolve(min)...");
330: } else if( LP_MAXMIN ){
331: dprint( 1, "lpsolve(maxmin)...");
1.2 noro 332: } else if( LP_INTPT ){
333: dprint( 1, "lpsolve(intpt)...");
1.1 noro 334: }
335:
336: if( index == LP_MAX || index == LP_MAXMIN ){
337: tmp = lpsolve( dd_LPmax, row, col, matrix, object );
338: if( resulttype == LP_I ){
339: lpmax = mytype2int( *tmp, -1 );
340: if( lpmax == 0 ){
341: cmomax = (cmo*) new_cmo_zero();
342: } else {
343: cmomax = (cmo*) new_cmo_zz_set_si( lpmax );
344: }
345: } else if( resulttype == LP_Q ){
346: #if defined GMPRATIONAL
347: if( mpz_cmp_si( mpq_numref( *tmp ), 0 ) == 0 ){
348: cmomax = (cmo*) new_cmo_zero();
349: } else {
350: cmomax = (cmo*) new_cmo_qq_set_mpz( mpq_numref( *tmp ), mpq_denref( *tmp ) );
351: }
352: #else
353: cmomax = (cmo*) new_cmo_double( *tmp[0] );
354: #endif
355: }
356: #if defined GMPRATIONAL
357: mpq_clear( *tmp );
358: #endif
359: }
360:
361: if( index == LP_MIN || index == LP_MAXMIN ){
362: tmp = lpsolve( dd_LPmin, row, col, matrix, object );
363: if( resulttype == LP_I ){
364: lpmin = mytype2int( *tmp, 1 );
365: if( lpmin == 0 ){
366: cmomin = (cmo*) new_cmo_zero();
367: } else {
368: cmomin = (cmo*) new_cmo_zz_set_si( lpmin );
369: }
370: } else if( resulttype == LP_Q ){
371: #if defined GMPRATIONAL
372: if( mpz_cmp_si( mpq_numref( *tmp ), 0 ) == 0 ){
373: cmomin = (cmo*) new_cmo_zero();
374: } else {
375: cmomin = (cmo*) new_cmo_qq_set_mpz( mpq_numref( *tmp ), mpq_denref( *tmp ) );
376: }
377: #else
378: cmomin = (cmo*) new_cmo_double( *tmp[0] );
379: #endif
380: }
381: #if defined GMPRATIONAL
382: mpq_clear( *tmp );
383: #endif
384:
385: }
386:
1.2 noro 387: if( index == LP_INTPT ) {
388: tmp = lpintpt(row, col, matrix);
1.3 ! noro 389: if ( tmp ) {
! 390: cmoint = (cmo *)new_cmo_list();
! 391: for(j=0;j < col;j++){
! 392: if( tmp[j] != 0 ){
! 393: #if defined GMPRATIONAL
! 394: a = (cmo*) new_cmo_qq_set_mpq( (mpq_ptr)tmp[j] );
! 395: #else
! 396: a = (cmo*) new_cmo_double( *tmp[j] );
! 397: #endif
! 398: } else {
! 399: a = (cmo*) new_cmo_zero();
! 400: }
! 401: cmoint = (cmo *)list_append( (cmo_list *) cmoint, a );
! 402: }
! 403: } else
1.2 noro 404: cmoint = new_cmo_zero();
1.3 ! noro 405: #if defined GMPRATIONAL
1.2 noro 406: for ( i = 0; i < col; i++ )
407: mpq_clear( (mpq_ptr)tmp[i] );
1.3 ! noro 408: #endif
1.2 noro 409: }
410:
1.1 noro 411: if( index == LP_MAX ){
412: push( cmomax );
413: } else if( index == LP_MIN ){
414: push( cmomin );
415: } else if( index == LP_MAXMIN ){
416: ret = new_cmo_list();
417: ret = list_append( ret, cmomin );
418: ret = list_append( ret, cmomax );
419: push( (cmo*) ret );
1.2 noro 420: } else if ( index == LP_INTPT )
421: push( cmoint );
1.1 noro 422: }
423:
424: int sm_executeFunction()
425: {
426: cmo_string *func = (cmo_string *)pop();
427: if(func->tag != CMO_STRING) {
428: printf("sm_executeFunction : func->tag is not CMO_STRING");fflush(stdout);
429: push((cmo*)make_error2(0));
430: return -1;
431: }
432:
433: if(strcmp(func->s, "redcheck" ) == 0) {
434: my_redcheck();
435: return 0;
436: } else if( strcmp( func->s, "ilpmax" ) == 0 ){
437: my_lpsolve(LP_I,LP_MAX);
438: return 0;
439: } else if( strcmp( func->s, "ilpmin" ) == 0 ){
440: my_lpsolve(LP_I,LP_MIN);
441: return 0;
442: } else if( strcmp( func->s, "ilpmaxmin" ) == 0 ){
443: my_lpsolve(LP_I,LP_MAXMIN);
444: return 0;
445: } else if( strcmp( func->s, "lpmax" ) == 0 ){
446: my_lpsolve(LP_Q,LP_MAX);
447: return 0;
448: } else if( strcmp( func->s, "lpmin" ) == 0 ){
449: my_lpsolve(LP_Q,LP_MIN);
450: return 0;
451: } else if( strcmp( func->s, "lpmaxmin" ) == 0 ){
452: my_lpsolve(LP_Q,LP_MAXMIN);
453: return 0;
1.2 noro 454: } else if( strcmp( func->s, "intpt" ) == 0 ){
455: my_lpsolve(LP_Q,LP_INTPT);
456: return 0;
1.1 noro 457: } else if( strcmp( func->s, "debugprint" ) == 0 ){
458: pop();
459: debug_print = get_i();
460: return 0;
461: } else if( strcmp( func->s, "exit" ) == 0 ){
462: pop();
463: exit(0);
464: return 0;
465: } else {
466: push((cmo*)make_error2(0));
467: return -1;
468: }
469: }
470:
471: int receive_and_execute_sm_command()
472: {
473: int code = receive_int32(fd_rw);
474: switch(code) {
475: case SM_popCMO:
476: sm_popCMO();
477: break;
478: case SM_executeFunction:
479: sm_executeFunction();
480: break;
481: case SM_mathcap:
482: sm_mathcap();
483: break;
484: case SM_setMathCap:
485: pop();
486: break;
487: default:
488: printf("receive_and_execute_sm_command : code=%d\n",code);fflush(stdout);
489: break;
490: }
491: return 0;
492: }
493:
494: int receive()
495: {
496: int tag;
497:
498: tag = receive_ox_tag(fd_rw);
499: switch(tag) {
500: case OX_DATA:
501: push(receive_cmo(fd_rw));
502: break;
503: case OX_COMMAND:
504: receive_and_execute_sm_command();
505: break;
506: default:
507: printf("receive : tag=%d\n",tag);fflush(stdout);
508: }
509: return 0;
510: }
511:
512: int main()
513: {
514: GC_INIT();
515: ox_stderr_init(stderr);
516: initialize_stack();
517: init_cdd();
518:
519: #if defined GMPRATIONAL
520: fprintf(stderr,"ox_cddgmp\n");
521: #else
522: fprintf(stderr,"ox_cdd\n");
523: #endif
524:
525: fd_rw = oxf_open(3);
526: oxf_determine_byteorder_server(fd_rw);
527:
528: while(1){
529: receive();
530: }
531: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>