aboutsummaryrefslogtreecommitdiffstats
path: root/test/monniaux/glpk-4.65/src/bflib/scfint.c
diff options
context:
space:
mode:
Diffstat (limited to 'test/monniaux/glpk-4.65/src/bflib/scfint.c')
-rw-r--r--test/monniaux/glpk-4.65/src/bflib/scfint.c255
1 files changed, 255 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/bflib/scfint.c b/test/monniaux/glpk-4.65/src/bflib/scfint.c
new file mode 100644
index 00000000..06aa8f7d
--- /dev/null
+++ b/test/monniaux/glpk-4.65/src/bflib/scfint.c
@@ -0,0 +1,255 @@
+/* scfint.c (interface to Schur-complement-based factorization) */
+
+/***********************************************************************
+* This code is part of GLPK (GNU Linear Programming Kit).
+*
+* Copyright (C) 2013-2014 Andrew Makhorin, Department for Applied
+* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
+* reserved. E-mail: <mao@gnu.org>.
+*
+* GLPK is free software: you can redistribute it and/or modify it
+* under the terms of the GNU General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* GLPK is distributed in the hope that it will be useful, but WITHOUT
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+* License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+
+#include "env.h"
+#include "scfint.h"
+
+SCFINT *scfint_create(int type)
+{ /* create interface to SC-factorization */
+ SCFINT *fi;
+ fi = talloc(1, SCFINT);
+ memset(fi, 0, sizeof(SCFINT));
+ switch ((fi->scf.type = type))
+ { case 1:
+ fi->u.lufi = lufint_create();
+ break;
+ case 2:
+ fi->u.btfi = btfint_create();
+ break;
+ default:
+ xassert(type != type);
+ }
+ return fi;
+}
+
+int scfint_factorize(SCFINT *fi, int n, int (*col)(void *info, int j,
+ int ind[], double val[]), void *info)
+{ /* compute SC-factorization of specified matrix A */
+ int nn_max, old_n0_max, n0_max, k, ret;
+ xassert(n > 0);
+ fi->valid = 0;
+ /* get required value of nn_max */
+ nn_max = fi->nn_max;
+ if (nn_max == 0)
+ nn_max = 100;
+ xassert(nn_max > 0);
+ /* compute factorization of specified matrix A */
+ switch (fi->scf.type)
+ { case 1:
+ old_n0_max = fi->u.lufi->n_max;
+ fi->u.lufi->sva_n_max = 4 * n + 2 * nn_max;
+ ret = lufint_factorize(fi->u.lufi, n, col, info);
+ n0_max = fi->u.lufi->n_max;
+ fi->scf.sva = fi->u.lufi->sva;
+ fi->scf.a0.luf = fi->u.lufi->luf;
+ break;
+ case 2:
+ old_n0_max = fi->u.btfi->n_max;
+ fi->u.btfi->sva_n_max = 6 * n + 2 * nn_max;
+ ret = btfint_factorize(fi->u.btfi, n, col, info);
+ n0_max = fi->u.btfi->n_max;
+ fi->scf.sva = fi->u.btfi->sva;
+ fi->scf.a0.btf = fi->u.btfi->btf;
+ break;
+ default:
+ xassert(fi != fi);
+ }
+ /* allocate/reallocate arrays, if necessary */
+ if (old_n0_max < n0_max)
+ { if (fi->w1 != NULL)
+ tfree(fi->w1);
+ if (fi->w2 != NULL)
+ tfree(fi->w2);
+ if (fi->w3 != NULL)
+ tfree(fi->w3);
+ fi->w1 = talloc(1+n0_max, double);
+ fi->w2 = talloc(1+n0_max, double);
+ fi->w3 = talloc(1+n0_max, double);
+ }
+ if (fi->scf.nn_max != nn_max)
+ { if (fi->scf.ifu.f != NULL)
+ tfree(fi->scf.ifu.f);
+ if (fi->scf.ifu.u != NULL)
+ tfree(fi->scf.ifu.u);
+ fi->scf.ifu.f = talloc(nn_max * nn_max, double);
+ fi->scf.ifu.u = talloc(nn_max * nn_max, double);
+ }
+ if (old_n0_max < n0_max || fi->scf.nn_max != nn_max)
+ { if (fi->scf.pp_ind != NULL)
+ tfree(fi->scf.pp_ind);
+ if (fi->scf.pp_inv != NULL)
+ tfree(fi->scf.pp_inv);
+ if (fi->scf.qq_ind != NULL)
+ tfree(fi->scf.qq_ind);
+ if (fi->scf.qq_inv != NULL)
+ tfree(fi->scf.qq_inv);
+ if (fi->w4 != NULL)
+ tfree(fi->w4);
+ if (fi->w5 != NULL)
+ tfree(fi->w5);
+ fi->scf.pp_ind = talloc(1+n0_max+nn_max, int);
+ fi->scf.pp_inv = talloc(1+n0_max+nn_max, int);
+ fi->scf.qq_ind = talloc(1+n0_max+nn_max, int);
+ fi->scf.qq_inv = talloc(1+n0_max+nn_max, int);
+ fi->w4 = talloc(1+n0_max+nn_max, double);
+ fi->w5 = talloc(1+n0_max+nn_max, double);
+ }
+ /* initialize SC-factorization */
+ fi->scf.n = n;
+ fi->scf.n0 = n;
+ fi->scf.nn_max = nn_max;
+ fi->scf.nn = 0;
+ fi->scf.rr_ref = sva_alloc_vecs(fi->scf.sva, nn_max);
+ fi->scf.ss_ref = sva_alloc_vecs(fi->scf.sva, nn_max);
+ fi->scf.ifu.n_max = nn_max;
+ fi->scf.ifu.n = 0;
+ for (k = 1; k <= n; k++)
+ { fi->scf.pp_ind[k] = k;
+ fi->scf.pp_inv[k] = k;
+ fi->scf.qq_ind[k] = k;
+ fi->scf.qq_inv[k] = k;
+ }
+ /* set validation flag */
+ if (ret == 0)
+ fi->valid = 1;
+ return ret;
+}
+
+int scfint_update(SCFINT *fi, int upd, int j, int len, const int ind[],
+ const double val[])
+{ /* update SC-factorization after replacing j-th column of A */
+ int n = fi->scf.n;
+ int n0 = fi->scf.n0;
+ int nn = fi->scf.nn;
+ int *pp_ind = fi->scf.pp_ind;
+ int *qq_ind = fi->scf.qq_ind;
+ int *qq_inv = fi->scf.qq_inv;
+ double *bf = fi->w4;
+ double *dg = fi->w5;
+ int k, t, ret;
+ xassert(fi->valid);
+ xassert(0 <= n && n <= n0+nn);
+ /* (b, f) := inv(P) * (beta, 0) */
+ for (k = 1; k <= n0+nn; k++)
+ bf[k] = 0.0;
+ for (t = 1; t <= len; t++)
+ { k = ind[t];
+ xassert(1 <= k && k <= n);
+#if 1 /* FIXME: currently P = I */
+ xassert(pp_ind[k] == k);
+#endif
+ xassert(bf[k] == 0.0);
+ xassert(val[t] != 0.0);
+ bf[k] = val[t];
+ }
+ /* (d, g) := Q * (cj, 0) */
+ for (k = 1; k <= n0+nn; k++)
+ dg[k] = 0.0;
+ xassert(1 <= j && j <= n);
+ dg[fi->scf.qq_inv[j]] = 1;
+ /* update factorization of augmented matrix */
+ ret = scf_update_aug(&fi->scf, &bf[0], &dg[0], &bf[n0], &dg[n0],
+ 0.0, upd, fi->w1, fi->w2, fi->w3);
+ if (ret == 0)
+ { /* swap j-th and last columns of new matrix Q */
+ scf_swap_q_cols(j, n0+nn+1);
+ }
+ else
+ { /* updating failed */
+ fi->valid = 0;
+ }
+ return ret;
+}
+
+void scfint_ftran(SCFINT *fi, double x[])
+{ /* solve system A * x = b */
+ xassert(fi->valid);
+ scf_a_solve(&fi->scf, x, fi->w4, fi->w5, fi->w1, fi->w2);
+ return;
+}
+
+void scfint_btran(SCFINT *fi, double x[])
+{ /* solve system A'* x = b */
+ xassert(fi->valid);
+ scf_at_solve(&fi->scf, x, fi->w4, fi->w5, fi->w1, fi->w2);
+ return;
+}
+
+double scfint_estimate(SCFINT *fi)
+{ /* estimate 1-norm of inv(A) */
+ double norm;
+ xassert(fi->valid);
+ xassert(fi->scf.n == fi->scf.n0);
+ switch (fi->scf.type)
+ { case 1:
+ norm = luf_estimate_norm(fi->scf.a0.luf, fi->w1, fi->w2);
+ break;
+ case 2:
+ norm = btf_estimate_norm(fi->scf.a0.btf, fi->w1, fi->w2,
+ fi->w3, fi->w4);
+ break;
+ default:
+ xassert(fi != fi);
+ }
+ return norm;
+}
+
+void scfint_delete(SCFINT *fi)
+{ /* delete interface to SC-factorization */
+ switch (fi->scf.type)
+ { case 1:
+ lufint_delete(fi->u.lufi);
+ break;
+ case 2:
+ btfint_delete(fi->u.btfi);
+ break;
+ default:
+ xassert(fi != fi);
+ }
+ if (fi->scf.ifu.f != NULL)
+ tfree(fi->scf.ifu.f);
+ if (fi->scf.ifu.u != NULL)
+ tfree(fi->scf.ifu.u);
+ if (fi->scf.pp_ind != NULL)
+ tfree(fi->scf.pp_ind);
+ if (fi->scf.pp_inv != NULL)
+ tfree(fi->scf.pp_inv);
+ if (fi->scf.qq_ind != NULL)
+ tfree(fi->scf.qq_ind);
+ if (fi->scf.qq_inv != NULL)
+ tfree(fi->scf.qq_inv);
+ if (fi->w1 != NULL)
+ tfree(fi->w1);
+ if (fi->w2 != NULL)
+ tfree(fi->w2);
+ if (fi->w3 != NULL)
+ tfree(fi->w3);
+ if (fi->w4 != NULL)
+ tfree(fi->w4);
+ if (fi->w5 != NULL)
+ tfree(fi->w5);
+ tfree(fi);
+ return;
+}
+
+/* eof */
generated by cgit (git 2.34.1) at 2024-06-11 10:26:36 +0000