← Back to branch summary

~ubuntu-branches/ubuntu/trusty/haskell-hmatrix/trusty

« back to all changes in this revision

Viewing changes to lib/Numeric/GSL/gsl-aux.c

  • Committer: Package Import Robot
  • Author(s): Denis Laxalde
  • Date: 2013-07-06 15:37:50 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20130706153750-wxxplc788jedqvv5
Tags: 0.15.0.0-1
* New upstream release.
* Make it clear in copyright that the license is GPL-3 (as stated by the
  author in <https://github.com/albertoruiz/hmatrix/issues/45>) although
  there is no proper license file yet.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/Numeric/GSL/gsl-aux.c
32
32
#include <gsl/gsl_complex_math.h>
33
33
#include <gsl/gsl_rng.h>
34
34
#include <gsl/gsl_randist.h>
 
35
#include <gsl/gsl_roots.h>
35
36
#include <gsl/gsl_multifit_nlin.h>
36
37
#include <string.h>
37
38
#include <stdio.h>
937
938
        iter++;
938
939
        if (status) break;
939
940
        status = gsl_multimin_test_size (size, tolsize);
940
 
    } while (status == GSL_CONTINUE && iter <= maxit);
 
941
    } while (status == GSL_CONTINUE && iter < maxit);
941
942
    int i,j;
942
943
    for (i=iter; i<solr; i++) {
943
944
        solp[i*solc+0] = iter;
1033
1034
        iter++;
1034
1035
        if (status) break;
1035
1036
        status = gsl_multimin_test_gradient (s->gradient, tolgrad);
1036
 
    } while (status == GSL_CONTINUE && iter <= maxit);
 
1037
    } while (status == GSL_CONTINUE && iter < maxit);
1037
1038
    int i,j;
1038
1039
    for (i=iter; i<solr; i++) {
1039
1040
        solp[i*solc+0] = iter;
1047
1048
 
1048
1049
//---------------------------------------------------------------
1049
1050
 
 
1051
double only_f_aux_root(double x, void *pars) {
 
1052
    double (*f)(double) = (double (*)(double)) pars;
 
1053
    return f(x);
 
1054
}
 
1055
 
 
1056
int root(int method, double f(double),
 
1057
         double epsrel, int maxit,
 
1058
         double xl, double xu, RMAT(sol)) {
 
1059
    REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
 
1060
    DEBUGMSG("root_only_f");
 
1061
    gsl_function my_func;
 
1062
    // extract function from pars
 
1063
    my_func.function = only_f_aux_root;
 
1064
    my_func.params = f;
 
1065
    size_t iter = 0;
 
1066
    int status;
 
1067
    const gsl_root_fsolver_type *T;
 
1068
    gsl_root_fsolver *s;
 
1069
    // Starting point
 
1070
    switch(method) {
 
1071
        case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; }
 
1072
        case 1 : {T = gsl_root_fsolver_falsepos; break; }
 
1073
        case 2 : {T = gsl_root_fsolver_brent; break; }
 
1074
        default: ERROR(BAD_CODE);
 
1075
    }
 
1076
    s = gsl_root_fsolver_alloc (T);
 
1077
    gsl_root_fsolver_set (s, &my_func, xl, xu);
 
1078
    do {
 
1079
           double best, current_lo, current_hi;
 
1080
           status = gsl_root_fsolver_iterate (s);
 
1081
           best = gsl_root_fsolver_root (s);
 
1082
           current_lo = gsl_root_fsolver_x_lower (s);
 
1083
           current_hi = gsl_root_fsolver_x_upper (s);
 
1084
           solp[iter*solc] = iter + 1;
 
1085
           solp[iter*solc+1] = best;
 
1086
           solp[iter*solc+2] = current_lo;
 
1087
           solp[iter*solc+3] = current_hi;
 
1088
           iter++;
 
1089
           if (status)   /* check if solver is stuck */
 
1090
             break;
 
1091
 
 
1092
           status =
 
1093
               gsl_root_test_interval (current_lo, current_hi, 0, epsrel);
 
1094
        }
 
1095
        while (status == GSL_CONTINUE && iter < maxit);
 
1096
    int i;
 
1097
    for (i=iter; i<solr; i++) {
 
1098
        solp[i*solc+0] = iter;
 
1099
        solp[i*solc+1]=0.;
 
1100
        solp[i*solc+2]=0.;
 
1101
        solp[i*solc+3]=0.;
 
1102
    }
 
1103
    gsl_root_fsolver_free(s);
 
1104
    OK
 
1105
}
 
1106
 
 
1107
typedef struct {
 
1108
    double (*f)(double);
 
1109
    double (*jf)(double);
 
1110
} uniTfjf;
 
1111
 
 
1112
double f_aux_uni(double x, void *pars) {
 
1113
    uniTfjf * fjf = ((uniTfjf*) pars);
 
1114
    return (fjf->f)(x);
 
1115
}
 
1116
 
 
1117
double jf_aux_uni(double x, void * pars) {
 
1118
    uniTfjf * fjf = ((uniTfjf*) pars);
 
1119
    return (fjf->jf)(x);
 
1120
}
 
1121
 
 
1122
void fjf_aux_uni(double x, void * pars, double * f, double * g) {
 
1123
    *f = f_aux_uni(x,pars);
 
1124
    *g = jf_aux_uni(x,pars);
 
1125
}
 
1126
 
 
1127
int rootj(int method, double f(double),
 
1128
          double df(double),
 
1129
         double epsrel, int maxit,
 
1130
         double x, RMAT(sol)) {
 
1131
    REQUIRES(solr == maxit && solc == 2,BAD_SIZE);
 
1132
    DEBUGMSG("root_fjf");
 
1133
    gsl_function_fdf my_func;
 
1134
    // extract function from pars
 
1135
    my_func.f = f_aux_uni;
 
1136
    my_func.df = jf_aux_uni;
 
1137
    my_func.fdf = fjf_aux_uni;
 
1138
    uniTfjf stfjf;
 
1139
    stfjf.f = f;
 
1140
    stfjf.jf = df;
 
1141
    my_func.params = &stfjf;
 
1142
    size_t iter = 0;
 
1143
    int status;
 
1144
    const gsl_root_fdfsolver_type *T;
 
1145
    gsl_root_fdfsolver *s;
 
1146
    // Starting point
 
1147
    switch(method) {
 
1148
        case 0 : {T = gsl_root_fdfsolver_newton;; break; }
 
1149
        case 1 : {T = gsl_root_fdfsolver_secant; break; }
 
1150
        case 2 : {T = gsl_root_fdfsolver_steffenson; break; }
 
1151
        default: ERROR(BAD_CODE);
 
1152
    }
 
1153
    s = gsl_root_fdfsolver_alloc (T);
 
1154
 
 
1155
    gsl_root_fdfsolver_set (s, &my_func, x);
 
1156
 
 
1157
    do {
 
1158
           double x0;
 
1159
           status = gsl_root_fdfsolver_iterate (s);
 
1160
           x0 = x;
 
1161
           x = gsl_root_fdfsolver_root(s);
 
1162
           solp[iter*solc+0] = iter+1;
 
1163
           solp[iter*solc+1] = x;
 
1164
 
 
1165
           iter++;
 
1166
           if (status)   /* check if solver is stuck */
 
1167
             break;
 
1168
 
 
1169
           status =
 
1170
               gsl_root_test_delta (x, x0, 0, epsrel);
 
1171
        }
 
1172
        while (status == GSL_CONTINUE && iter < maxit);
 
1173
 
 
1174
    int i;
 
1175
    for (i=iter; i<solr; i++) {
 
1176
        solp[i*solc+0] = iter;
 
1177
        solp[i*solc+1]=0.;
 
1178
    }
 
1179
    gsl_root_fdfsolver_free(s);
 
1180
    OK
 
1181
}
 
1182
 
 
1183
 
 
1184
//---------------------------------------------------------------
 
1185
 
1050
1186
typedef void TrawfunV(int, double*, int, double*);
1051
1187
 
1052
 
int only_f_aux_root(const gsl_vector*x, void *pars, gsl_vector*y) {
 
1188
int only_f_aux_multiroot(const gsl_vector*x, void *pars, gsl_vector*y) {
1053
1189
    TrawfunV * f = (TrawfunV*) pars;
1054
1190
    double* p = (double*)calloc(x->size,sizeof(double));
1055
1191
    double* q = (double*)calloc(y->size,sizeof(double));
1066
1202
    return 0; //hmmm
1067
1203
}
1068
1204
 
1069
 
int root(int method, void f(int, double*, int, double*),
 
1205
int multiroot(int method, void f(int, double*, int, double*),
1070
1206
         double epsabs, int maxit,
1071
1207
         KRVEC(xi), RMAT(sol)) {
1072
1208
    REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
1073
1209
    DEBUGMSG("root_only_f");
1074
1210
    gsl_multiroot_function my_func;
1075
1211
    // extract function from pars
1076
 
    my_func.f = only_f_aux_root;
 
1212
    my_func.f = only_f_aux_multiroot;
1077
1213
    my_func.n = xin;
1078
1214
    my_func.params = f;
1079
1215
    size_t iter = 0;
1112
1248
           status =
1113
1249
             gsl_multiroot_test_residual (s->f, epsabs);
1114
1250
        }
1115
 
        while (status == GSL_CONTINUE && iter <= maxit);
 
1251
        while (status == GSL_CONTINUE && iter < maxit);
1116
1252
 
1117
1253
    int i,j;
1118
1254
    for (i=iter; i<solr; i++) {
1175
1311
    return 0;
1176
1312
}
1177
1313
 
1178
 
int rootj(int method, int f(int, double*, int, double*),
 
1314
int multirootj(int method, int f(int, double*, int, double*),
1179
1315
                      int jac(int, double*, int, int, double*),
1180
1316
         double epsabs, int maxit,
1181
1317
         KRVEC(xi), RMAT(sol)) {
1228
1364
           status =
1229
1365
             gsl_multiroot_test_residual (s->f, epsabs);
1230
1366
        }
1231
 
        while (status == GSL_CONTINUE && iter <= maxit);
 
1367
        while (status == GSL_CONTINUE && iter < maxit);
1232
1368
 
1233
1369
    int i,j;
1234
1370
    for (i=iter; i<solr; i++) {
1293
1429
 
1294
1430
           status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel);
1295
1431
        }
1296
 
        while (status == GSL_CONTINUE && iter <= maxit);
 
1432
        while (status == GSL_CONTINUE && iter < maxit);
1297
1433
 
1298
1434
    int i,j;
1299
1435
    for (i=iter; i<solr; i++) {