You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
407 lines
14 KiB
407 lines
14 KiB
/* btfint.c (interface to BT-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 "btfint.h"
|
|
|
|
BTFINT *btfint_create(void)
|
|
{ /* create interface to BT-factorization */
|
|
BTFINT *fi;
|
|
fi = talloc(1, BTFINT);
|
|
fi->n_max = 0;
|
|
fi->valid = 0;
|
|
fi->sva = NULL;
|
|
fi->btf = NULL;
|
|
fi->sgf = NULL;
|
|
fi->sva_n_max = fi->sva_size = 0;
|
|
fi->delta_n0 = fi->delta_n = 0;
|
|
fi->sgf_piv_tol = 0.10;
|
|
fi->sgf_piv_lim = 4;
|
|
fi->sgf_suhl = 1;
|
|
fi->sgf_eps_tol = DBL_EPSILON;
|
|
return fi;
|
|
}
|
|
|
|
static void factorize_triv(BTFINT *fi, int k, int (*col)(void *info,
|
|
int j, int ind[], double val[]), void *info)
|
|
{ /* compute LU-factorization of diagonal block A~[k,k] and store
|
|
* corresponding columns of matrix A except elements of A~[k,k]
|
|
* (trivial case when the block has unity size) */
|
|
SVA *sva = fi->sva;
|
|
int *sv_ind = sva->ind;
|
|
double *sv_val = sva->val;
|
|
BTF *btf = fi->btf;
|
|
int *pp_inv = btf->pp_inv;
|
|
int *qq_ind = btf->qq_ind;
|
|
int *beg = btf->beg;
|
|
int ac_ref = btf->ac_ref;
|
|
int *ac_ptr = &sva->ptr[ac_ref-1];
|
|
int *ac_len = &sva->len[ac_ref-1];
|
|
SGF *sgf = fi->sgf;
|
|
int *ind = (int *)sgf->vr_max; /* working array */
|
|
double *val = sgf->work; /* working array */
|
|
int i, j, t, len, ptr, beg_k;
|
|
/* diagonal block A~[k,k] has the only element in matrix A~,
|
|
* which is a~[beg[k],beg[k]] = a[i,j] */
|
|
beg_k = beg[k];
|
|
i = pp_inv[beg_k];
|
|
j = qq_ind[beg_k];
|
|
/* get j-th column of A */
|
|
len = col(info, j, ind, val);
|
|
/* find element a[i,j] = a~[beg[k],beg[k]] in j-th column */
|
|
for (t = 1; t <= len; t++)
|
|
{ if (ind[t] == i)
|
|
break;
|
|
}
|
|
xassert(t <= len);
|
|
/* compute LU-factorization of diagonal block A~[k,k], where
|
|
* F = (1), V = (a[i,j]), P = Q = (1) (see the module LUF) */
|
|
#if 1 /* FIXME */
|
|
xassert(val[t] != 0.0);
|
|
#endif
|
|
btf->vr_piv[beg_k] = val[t];
|
|
btf->p1_ind[beg_k] = btf->p1_inv[beg_k] = 1;
|
|
btf->q1_ind[beg_k] = btf->q1_inv[beg_k] = 1;
|
|
/* remove element a[i,j] = a~[beg[k],beg[k]] from j-th column */
|
|
memmove(&ind[t], &ind[t+1], (len-t) * sizeof(int));
|
|
memmove(&val[t], &val[t+1], (len-t) * sizeof(double));
|
|
len--;
|
|
/* and store resulting j-th column of A into BTF */
|
|
if (len > 0)
|
|
{ /* reserve locations for j-th column of A */
|
|
if (sva->r_ptr - sva->m_ptr < len)
|
|
{ sva_more_space(sva, len);
|
|
sv_ind = sva->ind;
|
|
sv_val = sva->val;
|
|
}
|
|
sva_reserve_cap(sva, ac_ref+(j-1), len);
|
|
/* store j-th column of A (except elements of A~[k,k]) */
|
|
ptr = ac_ptr[j];
|
|
memcpy(&sv_ind[ptr], &ind[1], len * sizeof(int));
|
|
memcpy(&sv_val[ptr], &val[1], len * sizeof(double));
|
|
ac_len[j] = len;
|
|
}
|
|
return;
|
|
}
|
|
|
|
static int factorize_block(BTFINT *fi, int k, int (*col)(void *info,
|
|
int j, int ind[], double val[]), void *info)
|
|
{ /* compute LU-factorization of diagonal block A~[k,k] and store
|
|
* corresponding columns of matrix A except elements of A~[k,k]
|
|
* (general case) */
|
|
SVA *sva = fi->sva;
|
|
int *sv_ind = sva->ind;
|
|
double *sv_val = sva->val;
|
|
BTF *btf = fi->btf;
|
|
int *pp_ind = btf->pp_ind;
|
|
int *qq_ind = btf->qq_ind;
|
|
int *beg = btf->beg;
|
|
int ac_ref = btf->ac_ref;
|
|
int *ac_ptr = &sva->ptr[ac_ref-1];
|
|
int *ac_len = &sva->len[ac_ref-1];
|
|
SGF *sgf = fi->sgf;
|
|
int *ind = (int *)sgf->vr_max; /* working array */
|
|
double *val = sgf->work; /* working array */
|
|
LUF luf;
|
|
int *vc_ptr, *vc_len, *vc_cap;
|
|
int i, ii, j, jj, t, len, cnt, ptr, beg_k;
|
|
/* construct fake LUF for LU-factorization of A~[k,k] */
|
|
sgf->luf = &luf;
|
|
luf.n = beg[k+1] - (beg_k = beg[k]);
|
|
luf.sva = sva;
|
|
luf.fr_ref = btf->fr_ref + (beg_k-1);
|
|
luf.fc_ref = btf->fc_ref + (beg_k-1);
|
|
luf.vr_ref = btf->vr_ref + (beg_k-1);
|
|
luf.vr_piv = btf->vr_piv + (beg_k-1);
|
|
luf.vc_ref = btf->vc_ref + (beg_k-1);
|
|
luf.pp_ind = btf->p1_ind + (beg_k-1);
|
|
luf.pp_inv = btf->p1_inv + (beg_k-1);
|
|
luf.qq_ind = btf->q1_ind + (beg_k-1);
|
|
luf.qq_inv = btf->q1_inv + (beg_k-1);
|
|
/* process columns of k-th block of matrix A~ */
|
|
vc_ptr = &sva->ptr[luf.vc_ref-1];
|
|
vc_len = &sva->len[luf.vc_ref-1];
|
|
vc_cap = &sva->cap[luf.vc_ref-1];
|
|
for (jj = 1; jj <= luf.n; jj++)
|
|
{ /* jj-th column of A~ = j-th column of A */
|
|
j = qq_ind[jj + (beg_k-1)];
|
|
/* get j-th column of A */
|
|
len = col(info, j, ind, val);
|
|
/* move elements of diagonal block A~[k,k] to the beginning of
|
|
* the column list */
|
|
cnt = 0;
|
|
for (t = 1; t <= len; t++)
|
|
{ /* i = row index of element a[i,j] */
|
|
i = ind[t];
|
|
/* i-th row of A = ii-th row of A~ */
|
|
ii = pp_ind[i];
|
|
if (ii >= beg_k)
|
|
{ /* a~[ii,jj] = a[i,j] is in diagonal block A~[k,k] */
|
|
double temp;
|
|
cnt++;
|
|
ind[t] = ind[cnt];
|
|
ind[cnt] = ii - (beg_k-1); /* local index */
|
|
temp = val[t], val[t] = val[cnt], val[cnt] = temp;
|
|
}
|
|
}
|
|
/* first cnt elements in the column list give jj-th column of
|
|
* diagonal block A~[k,k], which is initial matrix V in LUF */
|
|
/* enlarge capacity of jj-th column of V = A~[k,k] */
|
|
if (vc_cap[jj] < cnt)
|
|
{ if (sva->r_ptr - sva->m_ptr < cnt)
|
|
{ sva_more_space(sva, cnt);
|
|
sv_ind = sva->ind;
|
|
sv_val = sva->val;
|
|
}
|
|
sva_enlarge_cap(sva, luf.vc_ref+(jj-1), cnt, 0);
|
|
}
|
|
/* store jj-th column of V = A~[k,k] */
|
|
ptr = vc_ptr[jj];
|
|
memcpy(&sv_ind[ptr], &ind[1], cnt * sizeof(int));
|
|
memcpy(&sv_val[ptr], &val[1], cnt * sizeof(double));
|
|
vc_len[jj] = cnt;
|
|
/* other (len-cnt) elements in the column list are stored in
|
|
* j-th column of the original matrix A */
|
|
len -= cnt;
|
|
if (len > 0)
|
|
{ /* reserve locations for j-th column of A */
|
|
if (sva->r_ptr - sva->m_ptr < len)
|
|
{ sva_more_space(sva, len);
|
|
sv_ind = sva->ind;
|
|
sv_val = sva->val;
|
|
}
|
|
sva_reserve_cap(sva, ac_ref-1+j, len);
|
|
/* store j-th column of A (except elements of A~[k,k]) */
|
|
ptr = ac_ptr[j];
|
|
memcpy(&sv_ind[ptr], &ind[cnt+1], len * sizeof(int));
|
|
memcpy(&sv_val[ptr], &val[cnt+1], len * sizeof(double));
|
|
ac_len[j] = len;
|
|
}
|
|
}
|
|
/* compute LU-factorization of diagonal block A~[k,k]; may note
|
|
* that A~[k,k] is irreducible (strongly connected), so singleton
|
|
* phase will have no effect */
|
|
k = sgf_factorize(sgf, 0 /* disable singleton phase */);
|
|
/* now left (dynamic) part of SVA should be empty (wichtig!) */
|
|
xassert(sva->m_ptr == 1);
|
|
return k;
|
|
}
|
|
|
|
int btfint_factorize(BTFINT *fi, int n, int (*col)(void *info, int j,
|
|
int ind[], double val[]), void *info)
|
|
{ /* compute BT-factorization of specified matrix A */
|
|
SVA *sva;
|
|
BTF *btf;
|
|
SGF *sgf;
|
|
int k, rank;
|
|
xassert(n > 0);
|
|
fi->valid = 0;
|
|
/* create sparse vector area (SVA), if necessary */
|
|
sva = fi->sva;
|
|
if (sva == NULL)
|
|
{ int sva_n_max = fi->sva_n_max;
|
|
int sva_size = fi->sva_size;
|
|
if (sva_n_max == 0)
|
|
sva_n_max = 6 * n;
|
|
if (sva_size == 0)
|
|
sva_size = 10 * n;
|
|
sva = fi->sva = sva_create_area(sva_n_max, sva_size);
|
|
}
|
|
/* allocate/reallocate underlying objects, if necessary */
|
|
if (fi->n_max < n)
|
|
{ int n_max = fi->n_max;
|
|
if (n_max == 0)
|
|
n_max = fi->n_max = n + fi->delta_n0;
|
|
else
|
|
n_max = fi->n_max = n + fi->delta_n;
|
|
xassert(n_max >= n);
|
|
/* allocate/reallocate block triangular factorization (BTF) */
|
|
btf = fi->btf;
|
|
if (btf == NULL)
|
|
{ btf = fi->btf = talloc(1, BTF);
|
|
memset(btf, 0, sizeof(BTF));
|
|
btf->sva = sva;
|
|
}
|
|
else
|
|
{ tfree(btf->pp_ind);
|
|
tfree(btf->pp_inv);
|
|
tfree(btf->qq_ind);
|
|
tfree(btf->qq_inv);
|
|
tfree(btf->beg);
|
|
tfree(btf->vr_piv);
|
|
tfree(btf->p1_ind);
|
|
tfree(btf->p1_inv);
|
|
tfree(btf->q1_ind);
|
|
tfree(btf->q1_inv);
|
|
}
|
|
btf->pp_ind = talloc(1+n_max, int);
|
|
btf->pp_inv = talloc(1+n_max, int);
|
|
btf->qq_ind = talloc(1+n_max, int);
|
|
btf->qq_inv = talloc(1+n_max, int);
|
|
btf->beg = talloc(1+n_max+1, int);
|
|
btf->vr_piv = talloc(1+n_max, double);
|
|
btf->p1_ind = talloc(1+n_max, int);
|
|
btf->p1_inv = talloc(1+n_max, int);
|
|
btf->q1_ind = talloc(1+n_max, int);
|
|
btf->q1_inv = talloc(1+n_max, int);
|
|
/* allocate/reallocate factorizer workspace (SGF) */
|
|
/* (note that for SGF we could use the size of largest block
|
|
* rather than n_max) */
|
|
sgf = fi->sgf;
|
|
sgf = fi->sgf;
|
|
if (sgf == NULL)
|
|
{ sgf = fi->sgf = talloc(1, SGF);
|
|
memset(sgf, 0, sizeof(SGF));
|
|
}
|
|
else
|
|
{ tfree(sgf->rs_head);
|
|
tfree(sgf->rs_prev);
|
|
tfree(sgf->rs_next);
|
|
tfree(sgf->cs_head);
|
|
tfree(sgf->cs_prev);
|
|
tfree(sgf->cs_next);
|
|
tfree(sgf->vr_max);
|
|
tfree(sgf->flag);
|
|
tfree(sgf->work);
|
|
}
|
|
sgf->rs_head = talloc(1+n_max, int);
|
|
sgf->rs_prev = talloc(1+n_max, int);
|
|
sgf->rs_next = talloc(1+n_max, int);
|
|
sgf->cs_head = talloc(1+n_max, int);
|
|
sgf->cs_prev = talloc(1+n_max, int);
|
|
sgf->cs_next = talloc(1+n_max, int);
|
|
sgf->vr_max = talloc(1+n_max, double);
|
|
sgf->flag = talloc(1+n_max, char);
|
|
sgf->work = talloc(1+n_max, double);
|
|
}
|
|
btf = fi->btf;
|
|
btf->n = n;
|
|
sgf = fi->sgf;
|
|
#if 1 /* FIXME */
|
|
/* initialize SVA */
|
|
sva->n = 0;
|
|
sva->m_ptr = 1;
|
|
sva->r_ptr = sva->size + 1;
|
|
sva->head = sva->tail = 0;
|
|
#endif
|
|
/* store pattern of original matrix A in column-wise format */
|
|
btf->ac_ref = sva_alloc_vecs(btf->sva, btf->n);
|
|
btf_store_a_cols(btf, col, info, btf->pp_ind, btf->vr_piv);
|
|
#ifdef GLP_DEBUG
|
|
sva_check_area(sva);
|
|
#endif
|
|
/* analyze pattern of original matrix A and determine permutation
|
|
* matrices P and Q such that A = P * A~* Q, where A~ is an upper
|
|
* block triangular matrix */
|
|
rank = btf_make_blocks(btf);
|
|
if (rank != n)
|
|
{ /* original matrix A is structurally singular */
|
|
return 1;
|
|
}
|
|
#ifdef GLP_DEBUG
|
|
btf_check_blocks(btf);
|
|
#endif
|
|
#if 1 /* FIXME */
|
|
/* initialize SVA */
|
|
sva->n = 0;
|
|
sva->m_ptr = 1;
|
|
sva->r_ptr = sva->size + 1;
|
|
sva->head = sva->tail = 0;
|
|
#endif
|
|
/* allocate sparse vectors in SVA */
|
|
btf->ar_ref = sva_alloc_vecs(btf->sva, btf->n);
|
|
btf->ac_ref = sva_alloc_vecs(btf->sva, btf->n);
|
|
btf->fr_ref = sva_alloc_vecs(btf->sva, btf->n);
|
|
btf->fc_ref = sva_alloc_vecs(btf->sva, btf->n);
|
|
btf->vr_ref = sva_alloc_vecs(btf->sva, btf->n);
|
|
btf->vc_ref = sva_alloc_vecs(btf->sva, btf->n);
|
|
/* setup factorizer control parameters */
|
|
sgf->updat = 0; /* wichtig! */
|
|
sgf->piv_tol = fi->sgf_piv_tol;
|
|
sgf->piv_lim = fi->sgf_piv_lim;
|
|
sgf->suhl = fi->sgf_suhl;
|
|
sgf->eps_tol = fi->sgf_eps_tol;
|
|
/* compute LU-factorizations of diagonal blocks A~[k,k] and also
|
|
* store corresponding columns of matrix A except elements of all
|
|
* blocks A~[k,k] */
|
|
for (k = 1; k <= btf->num; k++)
|
|
{ if (btf->beg[k+1] - btf->beg[k] == 1)
|
|
{ /* trivial case (A~[k,k] has unity order) */
|
|
factorize_triv(fi, k, col, info);
|
|
}
|
|
else
|
|
{ /* general case */
|
|
if (factorize_block(fi, k, col, info) != 0)
|
|
return 2; /* factorization of A~[k,k] failed */
|
|
}
|
|
}
|
|
#ifdef GLP_DEBUG
|
|
sva_check_area(sva);
|
|
#endif
|
|
/* build row-wise representation of matrix A */
|
|
btf_build_a_rows(fi->btf, fi->sgf->rs_head);
|
|
#ifdef GLP_DEBUG
|
|
sva_check_area(sva);
|
|
#endif
|
|
/* BT-factorization has been successfully computed */
|
|
fi->valid = 1;
|
|
return 0;
|
|
}
|
|
|
|
void btfint_delete(BTFINT *fi)
|
|
{ /* delete interface to BT-factorization */
|
|
SVA *sva = fi->sva;
|
|
BTF *btf = fi->btf;
|
|
SGF *sgf = fi->sgf;
|
|
if (sva != NULL)
|
|
sva_delete_area(sva);
|
|
if (btf != NULL)
|
|
{ tfree(btf->pp_ind);
|
|
tfree(btf->pp_inv);
|
|
tfree(btf->qq_ind);
|
|
tfree(btf->qq_inv);
|
|
tfree(btf->beg);
|
|
tfree(btf->vr_piv);
|
|
tfree(btf->p1_ind);
|
|
tfree(btf->p1_inv);
|
|
tfree(btf->q1_ind);
|
|
tfree(btf->q1_inv);
|
|
tfree(btf);
|
|
}
|
|
if (sgf != NULL)
|
|
{ tfree(sgf->rs_head);
|
|
tfree(sgf->rs_prev);
|
|
tfree(sgf->rs_next);
|
|
tfree(sgf->cs_head);
|
|
tfree(sgf->cs_prev);
|
|
tfree(sgf->cs_next);
|
|
tfree(sgf->vr_max);
|
|
tfree(sgf->flag);
|
|
tfree(sgf->work);
|
|
tfree(sgf);
|
|
}
|
|
tfree(fi);
|
|
return;
|
|
}
|
|
|
|
/* eof */
|