summaryrefslogtreecommitdiff
path: root/packages/glpk/lib/Numeric/LinearProgramming/glpk.c
blob: 0f4ad79626d4fab00aa7fce53e5af305ef260b97 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
typedef struct { double r, i; } doublecomplex;

#define DVEC(A) int A##n, double*A##p
#define CVEC(A) int A##n, doublecomplex*A##p
#define DMAT(A) int A##r, int A##c, double*A##p
#define CMAT(A) int A##r, int A##c, doublecomplex*A##p

#define AT(M,r,co) (M##p[(r)*M##c+(co)])

/*-----------------------------------------------------*/

#include <stdlib.h>
#include <stdio.h>
#include <glpk.h>
#include <math.h>

int c_simplex_dense(DMAT(c), DMAT(b), DVEC(s)) {
    glp_prob *lp;
    lp = glp_create_prob();
    glp_set_obj_dir(lp, GLP_MAX);

    int m = cr-1;
    int n = cc;
    glp_add_rows(lp, m);
    glp_add_cols(lp, n);
    int * ia = malloc((1+m*n)*sizeof(int));
    int * ja = malloc((1+m*n)*sizeof(int));
    double * ar = malloc((1+m*n)*sizeof(double));
    int i,j,k;
    k=0;
    for (i=1;i<=m;i++) {
        for(j=1;j<=n;j++) {
            k++;
            ia[k] = i;
            ja[k] = j;
            ar[k] = AT(c,i,j-1);
            //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
        }
    }
    glp_load_matrix(lp, k, ia, ja, ar);
    for (j=1;j<=n;j++) {
        glp_set_obj_coef(lp, j, AT(c,0,j-1));
    }
    int t;
    for (i=1;i<=m;i++) {
    switch((int)rint(AT(b,i-1,0))) {
        case 0: { t = GLP_FR; break; } 
        case 1: { t = GLP_LO; break; }
        case 2: { t = GLP_UP; break; }
        case 3: { t = GLP_DB; break; }
       default: { t = GLP_FX; break; }  
    }
    glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
    }
    for (j=1;j<=n;j++) {
    switch((int)rint(AT(b,m+j-1,0))) {
        case 0: { t = GLP_FR; break; } 
        case 1: { t = GLP_LO; break; }
        case 2: { t = GLP_UP; break; }
        case 3: { t = GLP_DB; break; }
       default: { t = GLP_FX; break; }  
    }
    glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
    }
    glp_term_out(0);
    glp_simplex(lp, NULL);
    sp[0] = glp_get_status(lp);
    sp[1] = glp_get_obj_val(lp);
    for (k=1; k<=n; k++) {
        sp[k+1] = glp_get_col_prim(lp, k);
    }
    
    glp_delete_prob(lp);
    free(ia);
    free(ja);
    free(ar);

    return 0;
}

/* ---------------------------------------------------- */

int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {
    glp_prob *lp;
    lp = glp_create_prob();
    glp_set_obj_dir(lp, GLP_MAX);
    int i,j,k;
    int tot = cr - n;
    glp_add_rows(lp, m);
    glp_add_cols(lp, n);

    //printf("%d %d\n",m,n);

    // the first n values
    for (k=1;k<=n;k++) {
        glp_set_obj_coef(lp, k, AT(c, k-1, 2));
        //printf("%d %f\n",k,AT(c, k-1, 2));
    }

    int * ia = malloc((1+tot)*sizeof(int));
    int * ja = malloc((1+tot)*sizeof(int));
    double * ar = malloc((1+tot)*sizeof(double));

    for (k=1; k<= tot; k++) {
        ia[k] = rint(AT(c,k-1+n,0));
        ja[k] = rint(AT(c,k-1+n,1));
        ar[k] =      AT(c,k-1+n,2);
        //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
    }
    glp_load_matrix(lp, tot, ia, ja, ar);
    
    int t;
    for (i=1;i<=m;i++) {
    switch((int)rint(AT(b,i-1,0))) {
        case 0: { t = GLP_FR; break; } 
        case 1: { t = GLP_LO; break; }
        case 2: { t = GLP_UP; break; }
        case 3: { t = GLP_DB; break; }
       default: { t = GLP_FX; break; }  
    }
    glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
    }
    for (j=1;j<=n;j++) {
    switch((int)rint(AT(b,m+j-1,0))) {
        case 0: { t = GLP_FR; break; } 
        case 1: { t = GLP_LO; break; }
        case 2: { t = GLP_UP; break; }
        case 3: { t = GLP_DB; break; }
       default: { t = GLP_FX; break; }  
    }
    glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
    }
    glp_term_out(0);
    glp_simplex(lp, NULL);
    sp[0] = glp_get_status(lp);
    sp[1] = glp_get_obj_val(lp);
    for (k=1; k<=n; k++) {
        sp[k+1] = glp_get_col_prim(lp, k);
    }
    glp_delete_prob(lp);
    free(ia);
    free(ja);
    free(ar);

    return 0;
}