"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/plugin/rq/rqfnb.c" (29 Aug 2019, 8170 Bytes) of package /linux/misc/gretl-2020e.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "rqfnb.c" see the Fossies "Dox" file reference documentation.

    1 /* rqfnb.f -- translated by f2c (version 20050501)
    2    and slightly cleaned up by Allin Cottrell
    3 */
    4 
    5 #include "libgretl.h"
    6 #include "gretl_f2c.h"
    7 
    8 /* Table of constant values */
    9 
   10 static double c_b4 = 1.;
   11 static integer one = 1;
   12 static double zero = 0.;
   13 static double c_b13 = -1.;
   14 
   15 /* lapack/blas functions called below */
   16 
   17 extern int dsyr_ (char *, integer *, double *, double *, 
   18           integer *, double *, integer *);
   19 
   20 extern int dposv_ (char *, integer *, integer *, double *, integer *, 
   21            double *, integer *, integer *);
   22 
   23 extern int dgemv_ (char *, integer *, integer *, double *, double *, 
   24            integer *, double *, integer *, double *, 
   25            double *, integer *);
   26 
   27 extern int dcopy_ (integer *, double *, integer *, double *, 
   28            integer *);
   29 
   30 extern int dswap_ (integer *, double *, integer *, double *, 
   31            integer *);
   32 
   33 extern int daxpy_ (integer *, double *, double *, integer *, 
   34            double *, integer *);
   35 
   36 extern int dpotrs_ (char *, integer *, integer *, double *, integer *, 
   37             double *, integer *, integer *);
   38 
   39 extern double ddot_ (integer *, double *, integer *, double *, 
   40              integer *);
   41 
   42 static int stepy_ (integer *n, integer *p, double *a, 
   43            double *d, double *b, double *ada, 
   44            integer *info)
   45 {
   46     integer i, m = *p * *p;
   47     int attempt = 0;
   48     int err = 0;
   49 
   50  try_again:
   51 
   52     for (i=0; i<m; i++) {
   53     ada[i] = 0.0;
   54     }
   55 
   56     for (i=0; i<*n; i++) {
   57     dsyr_("U", p, &d[i], &a[i * *p], &one, ada, p);
   58     }
   59 
   60     if (attempt == 0) {
   61     dposv_("U", p, &one, ada, p, b, p, info);
   62     if (*info != 0) {
   63         fprintf(stderr, "stepy: dposv gave info = %d\n", *info);
   64         attempt = 1;
   65         goto try_again;
   66     }
   67     } else {
   68     gretl_matrix A, B;
   69 
   70     gretl_matrix_init(&A);
   71     gretl_matrix_init(&B);
   72 
   73     A.rows = A.cols = *p;
   74     A.val = ada;
   75     B.rows = *p;
   76     B.cols = 1;
   77     B.val = b;
   78 
   79     err = gretl_LU_solve(&A, &B);
   80     if (err) {
   81         fprintf(stderr, "stepy: gretl_LU_solve: err = %d\n", err);
   82     }
   83     }
   84 
   85     return err;
   86 }
   87 
   88 #define ITERSTEP 5
   89 
   90 static int lpfnb_ (integer *n, integer *p, double *a, double *c__, 
   91            double *b, double *d__, double *u, double *beta,
   92            double *eps, double *x, double *s, double *y, 
   93            double *z__, double *w, double *dx, double *ds, 
   94            double *dy, double *dz, double *dw, double *dr, 
   95            double *rhs, double *ada, integer *nit, integer *info,
   96            void (*callback)(void))
   97 {
   98     integer a_dim1 = *p, ada_dim1 = *p;
   99     integer a_offset = 1 + a_dim1, ada_offset = 1 + ada_dim1;
  100     double d1, d2;
  101     static double g;
  102     static integer i;
  103     static double mu, gap;
  104     static double dsdw, dxdz;
  105     static double deltad, deltap;
  106     int main_iters = 0;
  107     int err = 0;
  108 
  109     /* Parameter adjustments */
  110     --dr;
  111     --dw;
  112     --dz;
  113     --ds;
  114     --dx;
  115     --w;
  116     --z__;
  117     --s;
  118     --x;
  119     --u;
  120     --d__;
  121     --c__;
  122     ada -= ada_offset;
  123     --rhs;
  124     --dy;
  125     --y;
  126     --b;
  127     a -= a_offset;
  128     --nit;
  129 
  130     /* Function Body */
  131     nit[1] = 0;
  132     nit[2] = 0;
  133     nit[3] = *n;
  134     dgemv_("N", p, n, &c_b4, &a[a_offset], p, &c__[1], &one, &zero, &y[1],
  135        &one);
  136     for (i = 1; i <= *n; ++i) {
  137     d__[i] = 1.;
  138     }
  139     err = stepy_(n, p, &a[a_offset], &d__[1], &y[1], &ada[ada_offset], info);
  140     if (err) {
  141     return err;
  142     }
  143     dcopy_(n, &c__[1], &one, &s[1], &one);
  144     dgemv_("T", p, n, &c_b13, &a[a_offset], p, &y[1], &one, &c_b4, &s[1],
  145        &one);
  146     for (i = 1; i <= *n; ++i) {
  147     if ((d1 = s[i], fabs(d1)) < *eps) {
  148         d1 = s[i];
  149         z__[i] = max(d1,0.) + *eps;
  150         d1 = -s[i];
  151         w[i] = max(d1,0.) + *eps;
  152     } else {
  153         d1 = s[i];
  154         z__[i] = max(d1, 0.);
  155         d1 = -s[i];
  156         w[i] = max(d1, 0.);
  157     }
  158     s[i] = u[i] - x[i];
  159     }
  160 
  161     gap = ddot_(n, &z__[1], &one, &x[1], &one) + 
  162     ddot_(n, &w[1], &one, &s[1], &one);
  163 
  164 looptop:
  165 
  166     if (callback != NULL && (main_iters++ % ITERSTEP == 0)) {
  167     callback();
  168     }
  169 
  170     if (gap > *eps && nit[1] < 50) {
  171     ++nit[1];
  172     for (i = 1; i <= *n; ++i) {
  173         d__[i] = 1. / (z__[i] / x[i] + w[i] / s[i]);
  174         ds[i] = z__[i] - w[i];
  175         dz[i] = d__[i] * ds[i];
  176     }
  177 
  178     dcopy_(p, &b[1], &one, &dy[1], &one);
  179     dgemv_("N", p, n, &c_b13, &a[a_offset], p, &x[1], &one, &c_b4, &dy[1],
  180         &one);
  181     dgemv_("N", p, n, &c_b4, &a[a_offset], p, &dz[1], &one, &c_b4, &dy[1],
  182         &one);
  183     dcopy_(p, &dy[1], &one, &rhs[1], &one);
  184     err = stepy_(n, p, &a[a_offset], &d__[1], &dy[1], &ada[ada_offset], info);
  185     if (err) {
  186         return err;
  187     }
  188 
  189     dgemv_("T", p, n, &c_b4, &a[a_offset], p, &dy[1], &one, &c_b13, 
  190            &ds[1], &one);
  191     deltap = 1e20;
  192     deltad = 1e20;
  193 
  194     for (i = 1; i <= *n; ++i) {
  195         dx[i] = d__[i] * ds[i];
  196         ds[i] = -dx[i];
  197         dz[i] = -z__[i] * (dx[i] / x[i] + 1.);
  198         dw[i] = -w[i] * (ds[i] / s[i] + 1.);
  199         if (dx[i] < 0.) {
  200         d1 = deltap, d2 = -x[i] / dx[i];
  201         deltap = min(d1,d2);
  202         }
  203         if (ds[i] < 0.) {
  204         d1 = deltap, d2 = -s[i] / ds[i];
  205         deltap = min(d1,d2);
  206         }
  207         if (dz[i] < 0.) {
  208         d1 = deltad, d2 = -z__[i] / dz[i];
  209         deltad = min(d1,d2);
  210         }
  211         if (dw[i] < 0.) {
  212         d1 = deltad, d2 = -w[i] / dw[i];
  213         deltad = min(d1,d2);
  214         }
  215     }
  216     d1 = *beta * deltap;
  217     deltap = min(d1,1.);
  218     d1 = *beta * deltad;
  219     deltad = min(d1,1.);
  220     if (min(deltap,deltad) < 1.) {
  221         ++nit[2];
  222         mu = ddot_(n, &x[1], &one, &z__[1], &one) + 
  223         ddot_(n, &s[1], &one, &w[1], &one);
  224         g = mu + deltap * ddot_(n, &dx[1], &one, &z__[1], &one) + 
  225         deltad * ddot_(n, &dz[1], &one, &x[1], &one) + 
  226         deltap * deltad * ddot_(n, &dz[1], &one, &dx[1], &one) + 
  227         deltap * ddot_(n, &ds[1], &one, &w[1], &one) + 
  228         deltad * ddot_(n, &dw[1], &one, &s[1], &one) + 
  229         deltap * deltad * ddot_(n, &ds[1], &one, &dw[1], &one);
  230         d1 = g / mu;
  231         mu = mu * (d1 * (d1 * d1)) / (double) (*n << 1);
  232         for (i = 1; i <= *n; ++i) {
  233         dr[i] = d__[i] * (mu * (1 / s[i] - 1 / x[i]) + dx[i]
  234                   * dz[i] / x[i] - ds[i] * dw[i] / s[i]);
  235         }
  236         dswap_(p, &rhs[1], &one, &dy[1], &one);
  237         dgemv_("N", p, n, &c_b4, &a[a_offset], p, &dr[1], &one, &c_b4,
  238            &dy[1], &one);
  239         dpotrs_("U", p, &one, &ada[ada_offset], p, &dy[1], p, info);
  240         if (*info != 0) {
  241         fprintf(stderr, "lpfnb: dpotrs_ gave info = %d\n", *info);
  242         }
  243         dgemv_("T", p, n, &c_b4, &a[a_offset], p, &dy[1], &one, &zero, 
  244            &u[1], &one);
  245         deltap = 1e20;
  246         deltad = 1e20;
  247         for (i = 1; i <= *n; ++i) {
  248         dxdz = dx[i] * dz[i];
  249         dsdw = ds[i] * dw[i];
  250         dx[i] = d__[i] * (u[i] - z__[i] + w[i]) - dr[i];
  251         ds[i] = -dx[i];
  252         dz[i] = -z__[i] + (mu - z__[i] * dx[i] - dxdz) / x[i];
  253         dw[i] = -w[i] + (mu - w[i] * ds[i] - dsdw) / s[i];
  254         if (dx[i] < 0.) {
  255             d1 = deltap, d2 = -x[i] / dx[i];
  256             deltap = min(d1, d2);
  257         }
  258         if (ds[i] < 0.) {
  259             d1 = deltap, d2 = -s[i] / ds[i];
  260             deltap = min(d1, d2);
  261         }
  262         if (dz[i] < 0.) {
  263             d1 = deltad, d2 = -z__[i] / dz[i];
  264             deltad = min(d1, d2);
  265         }
  266         if (dw[i] < 0.) {
  267             d1 = deltad, d2 = -w[i] / dw[i];
  268             deltad = min(d1, d2);
  269         }
  270         }
  271         d1 = *beta * deltap;
  272         deltap = min(d1,1.);
  273         d1 = *beta * deltad;
  274         deltad = min(d1,1.);
  275     }
  276 
  277     daxpy_(n, &deltap, &dx[1], &one, &x[1], &one);
  278     daxpy_(n, &deltap, &ds[1], &one, &s[1], &one);
  279     daxpy_(p, &deltad, &dy[1], &one, &y[1], &one);
  280     daxpy_(n, &deltad, &dz[1], &one, &z__[1], &one);
  281     daxpy_(n, &deltad, &dw[1], &one, &w[1], &one);
  282     gap = ddot_(n, &z__[1], &one, &x[1], &one) + 
  283         ddot_(n, &w[1], &one, &s[1], &one);
  284 
  285     goto looptop;
  286     }
  287 
  288     daxpy_(n, &c_b13, &w[1], &one, &z__[1], &one);
  289     dswap_(n, &z__[1], &one, &x[1], &one);
  290 
  291     return err;
  292 } /* end of lpfnb_ */
  293 
  294 int rqfnb_ (integer *n, integer *p, double *a, double *y, 
  295         double *rhs, double *d, double *u, double *beta,
  296         double *eps, double *wn, double *wp, integer *nit, 
  297         integer *info, void (*callback)(void))
  298 {
  299     integer a_dim = *p, wn_dim = *n, wp_dim = *p;
  300     int err;
  301 
  302     /* Parameter adjustments */
  303     wn -= 1 + wn_dim;
  304     wp -= 1 + wp_dim;
  305     a -= 1 + a_dim;
  306 
  307     err = lpfnb_(n, p, &a[a_dim + 1], y, rhs, d, u, beta, eps, 
  308          &wn[wn_dim + 1], &wn[(wn_dim << 1) + 1], &wp[wp_dim + 1], 
  309          &wn[wn_dim * 3 + 1], &wn[(wn_dim << 2) + 1], &wn[wn_dim * 5 + 1], 
  310          &wn[wn_dim * 6 + 1], &wp[(wp_dim << 1) + 1], &wn[wn_dim * 7 + 1],
  311          &wn[(wn_dim << 3) + 1], &wn[wn_dim * 9 + 1], &wp[wp_dim * 3 + 1], 
  312          &wp[(wp_dim << 2) + 1], nit, info, callback);
  313 
  314     return err;
  315 }