lemon-project-template-glpk
diff deps/glpk/src/glpluf.c @ 9:33de93886c88
Import GLPK 4.47
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 20:59:10 +0100 |
| parents | |
| children |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/deps/glpk/src/glpluf.c Sun Nov 06 20:59:10 2011 +0100 1.3 @@ -0,0 +1,1846 @@ 1.4 +/* glpluf.c (LU-factorization) */ 1.5 + 1.6 +/*********************************************************************** 1.7 +* This code is part of GLPK (GNU Linear Programming Kit). 1.8 +* 1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 1.10 +* 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, 1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved. 1.12 +* E-mail: <mao@gnu.org>. 1.13 +* 1.14 +* GLPK is free software: you can redistribute it and/or modify it 1.15 +* under the terms of the GNU General Public License as published by 1.16 +* the Free Software Foundation, either version 3 of the License, or 1.17 +* (at your option) any later version. 1.18 +* 1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT 1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public 1.22 +* License for more details. 1.23 +* 1.24 +* You should have received a copy of the GNU General Public License 1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>. 1.26 +***********************************************************************/ 1.27 + 1.28 +#include "glpenv.h" 1.29 +#include "glpluf.h" 1.30 +#define xfault xerror 1.31 + 1.32 +/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ 1.33 + 1.34 +#define N_MAX 100000000 /* = 100*10^6 */ 1.35 +/* maximal order of the original matrix */ 1.36 + 1.37 +/*********************************************************************** 1.38 +* NAME 1.39 +* 1.40 +* luf_create_it - create LU-factorization 1.41 +* 1.42 +* SYNOPSIS 1.43 +* 1.44 +* #include "glpluf.h" 1.45 +* LUF *luf_create_it(void); 1.46 +* 1.47 +* DESCRIPTION 1.48 +* 1.49 +* The routine luf_create_it creates a program object, which represents 1.50 +* LU-factorization of a square matrix. 1.51 +* 1.52 +* RETURNS 1.53 +* 1.54 +* The routine luf_create_it returns a pointer to the object created. */ 1.55 + 1.56 +LUF *luf_create_it(void) 1.57 +{ LUF *luf; 1.58 + luf = xmalloc(sizeof(LUF)); 1.59 + luf->n_max = luf->n = 0; 1.60 + luf->valid = 0; 1.61 + luf->fr_ptr = luf->fr_len = NULL; 1.62 + luf->fc_ptr = luf->fc_len = NULL; 1.63 + luf->vr_ptr = luf->vr_len = luf->vr_cap = NULL; 1.64 + luf->vr_piv = NULL; 1.65 + luf->vc_ptr = luf->vc_len = luf->vc_cap = NULL; 1.66 + luf->pp_row = luf->pp_col = NULL; 1.67 + luf->qq_row = luf->qq_col = NULL; 1.68 + luf->sv_size = 0; 1.69 + luf->sv_beg = luf->sv_end = 0; 1.70 + luf->sv_ind = NULL; 1.71 + luf->sv_val = NULL; 1.72 + luf->sv_head = luf->sv_tail = 0; 1.73 + luf->sv_prev = luf->sv_next = NULL; 1.74 + luf->vr_max = NULL; 1.75 + luf->rs_head = luf->rs_prev = luf->rs_next = NULL; 1.76 + luf->cs_head = luf->cs_prev = luf->cs_next = NULL; 1.77 + luf->flag = NULL; 1.78 + luf->work = NULL; 1.79 + luf->new_sva = 0; 1.80 + luf->piv_tol = 0.10; 1.81 + luf->piv_lim = 4; 1.82 + luf->suhl = 1; 1.83 + luf->eps_tol = 1e-15; 1.84 + luf->max_gro = 1e+10; 1.85 + luf->nnz_a = luf->nnz_f = luf->nnz_v = 0; 1.86 + luf->max_a = luf->big_v = 0.0; 1.87 + luf->rank = 0; 1.88 + return luf; 1.89 +} 1.90 + 1.91 +/*********************************************************************** 1.92 +* NAME 1.93 +* 1.94 +* luf_defrag_sva - defragment the sparse vector area 1.95 +* 1.96 +* SYNOPSIS 1.97 +* 1.98 +* #include "glpluf.h" 1.99 +* void luf_defrag_sva(LUF *luf); 1.100 +* 1.101 +* DESCRIPTION 1.102 +* 1.103 +* The routine luf_defrag_sva defragments the sparse vector area (SVA) 1.104 +* gathering all unused locations in one continuous extent. In order to 1.105 +* do that the routine moves all unused locations from the left part of 1.106 +* SVA (which contains rows and columns of the matrix V) to the middle 1.107 +* part (which contains free locations). This is attained by relocating 1.108 +* elements of rows and columns of the matrix V toward the beginning of 1.109 +* the left part. 1.110 +* 1.111 +* NOTE that this "garbage collection" involves changing row and column 1.112 +* pointers of the matrix V. */ 1.113 + 1.114 +void luf_defrag_sva(LUF *luf) 1.115 +{ int n = luf->n; 1.116 + int *vr_ptr = luf->vr_ptr; 1.117 + int *vr_len = luf->vr_len; 1.118 + int *vr_cap = luf->vr_cap; 1.119 + int *vc_ptr = luf->vc_ptr; 1.120 + int *vc_len = luf->vc_len; 1.121 + int *vc_cap = luf->vc_cap; 1.122 + int *sv_ind = luf->sv_ind; 1.123 + double *sv_val = luf->sv_val; 1.124 + int *sv_next = luf->sv_next; 1.125 + int sv_beg = 1; 1.126 + int i, j, k; 1.127 + /* skip rows and columns, which do not need to be relocated */ 1.128 + for (k = luf->sv_head; k != 0; k = sv_next[k]) 1.129 + { if (k <= n) 1.130 + { /* i-th row of the matrix V */ 1.131 + i = k; 1.132 + if (vr_ptr[i] != sv_beg) break; 1.133 + vr_cap[i] = vr_len[i]; 1.134 + sv_beg += vr_cap[i]; 1.135 + } 1.136 + else 1.137 + { /* j-th column of the matrix V */ 1.138 + j = k - n; 1.139 + if (vc_ptr[j] != sv_beg) break; 1.140 + vc_cap[j] = vc_len[j]; 1.141 + sv_beg += vc_cap[j]; 1.142 + } 1.143 + } 1.144 + /* relocate other rows and columns in order to gather all unused 1.145 + locations in one continuous extent */ 1.146 + for (k = k; k != 0; k = sv_next[k]) 1.147 + { if (k <= n) 1.148 + { /* i-th row of the matrix V */ 1.149 + i = k; 1.150 + memmove(&sv_ind[sv_beg], &sv_ind[vr_ptr[i]], 1.151 + vr_len[i] * sizeof(int)); 1.152 + memmove(&sv_val[sv_beg], &sv_val[vr_ptr[i]], 1.153 + vr_len[i] * sizeof(double)); 1.154 + vr_ptr[i] = sv_beg; 1.155 + vr_cap[i] = vr_len[i]; 1.156 + sv_beg += vr_cap[i]; 1.157 + } 1.158 + else 1.159 + { /* j-th column of the matrix V */ 1.160 + j = k - n; 1.161 + memmove(&sv_ind[sv_beg], &sv_ind[vc_ptr[j]], 1.162 + vc_len[j] * sizeof(int)); 1.163 + memmove(&sv_val[sv_beg], &sv_val[vc_ptr[j]], 1.164 + vc_len[j] * sizeof(double)); 1.165 + vc_ptr[j] = sv_beg; 1.166 + vc_cap[j] = vc_len[j]; 1.167 + sv_beg += vc_cap[j]; 1.168 + } 1.169 + } 1.170 + /* set new pointer to the beginning of the free part */ 1.171 + luf->sv_beg = sv_beg; 1.172 + return; 1.173 +} 1.174 + 1.175 +/*********************************************************************** 1.176 +* NAME 1.177 +* 1.178 +* luf_enlarge_row - enlarge row capacity 1.179 +* 1.180 +* SYNOPSIS 1.181 +* 1.182 +* #include "glpluf.h" 1.183 +* int luf_enlarge_row(LUF *luf, int i, int cap); 1.184 +* 1.185 +* DESCRIPTION 1.186 +* 1.187 +* The routine luf_enlarge_row enlarges capacity of the i-th row of the 1.188 +* matrix V to cap locations (assuming that its current capacity is less 1.189 +* than cap). In order to do that the routine relocates elements of the 1.190 +* i-th row to the end of the left part of SVA (which contains rows and 1.191 +* columns of the matrix V) and then expands the left part by allocating 1.192 +* cap free locations from the free part. If there are less than cap 1.193 +* free locations, the routine defragments the sparse vector area. 1.194 +* 1.195 +* Due to "garbage collection" this operation may change row and column 1.196 +* pointers of the matrix V. 1.197 +* 1.198 +* RETURNS 1.199 +* 1.200 +* If no error occured, the routine returns zero. Otherwise, in case of 1.201 +* overflow of the sparse vector area, the routine returns non-zero. */ 1.202 + 1.203 +int luf_enlarge_row(LUF *luf, int i, int cap) 1.204 +{ int n = luf->n; 1.205 + int *vr_ptr = luf->vr_ptr; 1.206 + int *vr_len = luf->vr_len; 1.207 + int *vr_cap = luf->vr_cap; 1.208 + int *vc_cap = luf->vc_cap; 1.209 + int *sv_ind = luf->sv_ind; 1.210 + double *sv_val = luf->sv_val; 1.211 + int *sv_prev = luf->sv_prev; 1.212 + int *sv_next = luf->sv_next; 1.213 + int ret = 0; 1.214 + int cur, k, kk; 1.215 + xassert(1 <= i && i <= n); 1.216 + xassert(vr_cap[i] < cap); 1.217 + /* if there are less than cap free locations, defragment SVA */ 1.218 + if (luf->sv_end - luf->sv_beg < cap) 1.219 + { luf_defrag_sva(luf); 1.220 + if (luf->sv_end - luf->sv_beg < cap) 1.221 + { ret = 1; 1.222 + goto done; 1.223 + } 1.224 + } 1.225 + /* save current capacity of the i-th row */ 1.226 + cur = vr_cap[i]; 1.227 + /* copy existing elements to the beginning of the free part */ 1.228 + memmove(&sv_ind[luf->sv_beg], &sv_ind[vr_ptr[i]], 1.229 + vr_len[i] * sizeof(int)); 1.230 + memmove(&sv_val[luf->sv_beg], &sv_val[vr_ptr[i]], 1.231 + vr_len[i] * sizeof(double)); 1.232 + /* set new pointer and new capacity of the i-th row */ 1.233 + vr_ptr[i] = luf->sv_beg; 1.234 + vr_cap[i] = cap; 1.235 + /* set new pointer to the beginning of the free part */ 1.236 + luf->sv_beg += cap; 1.237 + /* now the i-th row starts in the rightmost location among other 1.238 + rows and columns of the matrix V, so its node should be moved 1.239 + to the end of the row/column linked list */ 1.240 + k = i; 1.241 + /* remove the i-th row node from the linked list */ 1.242 + if (sv_prev[k] == 0) 1.243 + luf->sv_head = sv_next[k]; 1.244 + else 1.245 + { /* capacity of the previous row/column can be increased at the 1.246 + expense of old locations of the i-th row */ 1.247 + kk = sv_prev[k]; 1.248 + if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur; 1.249 + sv_next[sv_prev[k]] = sv_next[k]; 1.250 + } 1.251 + if (sv_next[k] == 0) 1.252 + luf->sv_tail = sv_prev[k]; 1.253 + else 1.254 + sv_prev[sv_next[k]] = sv_prev[k]; 1.255 + /* insert the i-th row node to the end of the linked list */ 1.256 + sv_prev[k] = luf->sv_tail; 1.257 + sv_next[k] = 0; 1.258 + if (sv_prev[k] == 0) 1.259 + luf->sv_head = k; 1.260 + else 1.261 + sv_next[sv_prev[k]] = k; 1.262 + luf->sv_tail = k; 1.263 +done: return ret; 1.264 +} 1.265 + 1.266 +/*********************************************************************** 1.267 +* NAME 1.268 +* 1.269 +* luf_enlarge_col - enlarge column capacity 1.270 +* 1.271 +* SYNOPSIS 1.272 +* 1.273 +* #include "glpluf.h" 1.274 +* int luf_enlarge_col(LUF *luf, int j, int cap); 1.275 +* 1.276 +* DESCRIPTION 1.277 +* 1.278 +* The routine luf_enlarge_col enlarges capacity of the j-th column of 1.279 +* the matrix V to cap locations (assuming that its current capacity is 1.280 +* less than cap). In order to do that the routine relocates elements 1.281 +* of the j-th column to the end of the left part of SVA (which contains 1.282 +* rows and columns of the matrix V) and then expands the left part by 1.283 +* allocating cap free locations from the free part. If there are less 1.284 +* than cap free locations, the routine defragments the sparse vector 1.285 +* area. 1.286 +* 1.287 +* Due to "garbage collection" this operation may change row and column 1.288 +* pointers of the matrix V. 1.289 +* 1.290 +* RETURNS 1.291 +* 1.292 +* If no error occured, the routine returns zero. Otherwise, in case of 1.293 +* overflow of the sparse vector area, the routine returns non-zero. */ 1.294 + 1.295 +int luf_enlarge_col(LUF *luf, int j, int cap) 1.296 +{ int n = luf->n; 1.297 + int *vr_cap = luf->vr_cap; 1.298 + int *vc_ptr = luf->vc_ptr; 1.299 + int *vc_len = luf->vc_len; 1.300 + int *vc_cap = luf->vc_cap; 1.301 + int *sv_ind = luf->sv_ind; 1.302 + double *sv_val = luf->sv_val; 1.303 + int *sv_prev = luf->sv_prev; 1.304 + int *sv_next = luf->sv_next; 1.305 + int ret = 0; 1.306 + int cur, k, kk; 1.307 + xassert(1 <= j && j <= n); 1.308 + xassert(vc_cap[j] < cap); 1.309 + /* if there are less than cap free locations, defragment SVA */ 1.310 + if (luf->sv_end - luf->sv_beg < cap) 1.311 + { luf_defrag_sva(luf); 1.312 + if (luf->sv_end - luf->sv_beg < cap) 1.313 + { ret = 1; 1.314 + goto done; 1.315 + } 1.316 + } 1.317 + /* save current capacity of the j-th column */ 1.318 + cur = vc_cap[j]; 1.319 + /* copy existing elements to the beginning of the free part */ 1.320 + memmove(&sv_ind[luf->sv_beg], &sv_ind[vc_ptr[j]], 1.321 + vc_len[j] * sizeof(int)); 1.322 + memmove(&sv_val[luf->sv_beg], &sv_val[vc_ptr[j]], 1.323 + vc_len[j] * sizeof(double)); 1.324 + /* set new pointer and new capacity of the j-th column */ 1.325 + vc_ptr[j] = luf->sv_beg; 1.326 + vc_cap[j] = cap; 1.327 + /* set new pointer to the beginning of the free part */ 1.328 + luf->sv_beg += cap; 1.329 + /* now the j-th column starts in the rightmost location among 1.330 + other rows and columns of the matrix V, so its node should be 1.331 + moved to the end of the row/column linked list */ 1.332 + k = n + j; 1.333 + /* remove the j-th column node from the linked list */ 1.334 + if (sv_prev[k] == 0) 1.335 + luf->sv_head = sv_next[k]; 1.336 + else 1.337 + { /* capacity of the previous row/column can be increased at the 1.338 + expense of old locations of the j-th column */ 1.339 + kk = sv_prev[k]; 1.340 + if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur; 1.341 + sv_next[sv_prev[k]] = sv_next[k]; 1.342 + } 1.343 + if (sv_next[k] == 0) 1.344 + luf->sv_tail = sv_prev[k]; 1.345 + else 1.346 + sv_prev[sv_next[k]] = sv_prev[k]; 1.347 + /* insert the j-th column node to the end of the linked list */ 1.348 + sv_prev[k] = luf->sv_tail; 1.349 + sv_next[k] = 0; 1.350 + if (sv_prev[k] == 0) 1.351 + luf->sv_head = k; 1.352 + else 1.353 + sv_next[sv_prev[k]] = k; 1.354 + luf->sv_tail = k; 1.355 +done: return ret; 1.356 +} 1.357 + 1.358 +/*********************************************************************** 1.359 +* reallocate - reallocate LU-factorization arrays 1.360 +* 1.361 +* This routine reallocates arrays, whose size depends of n, the order 1.362 +* of the matrix A to be factorized. */ 1.363 + 1.364 +static void reallocate(LUF *luf, int n) 1.365 +{ int n_max = luf->n_max; 1.366 + luf->n = n; 1.367 + if (n <= n_max) goto done; 1.368 + if (luf->fr_ptr != NULL) xfree(luf->fr_ptr); 1.369 + if (luf->fr_len != NULL) xfree(luf->fr_len); 1.370 + if (luf->fc_ptr != NULL) xfree(luf->fc_ptr); 1.371 + if (luf->fc_len != NULL) xfree(luf->fc_len); 1.372 + if (luf->vr_ptr != NULL) xfree(luf->vr_ptr); 1.373 + if (luf->vr_len != NULL) xfree(luf->vr_len); 1.374 + if (luf->vr_cap != NULL) xfree(luf->vr_cap); 1.375 + if (luf->vr_piv != NULL) xfree(luf->vr_piv); 1.376 + if (luf->vc_ptr != NULL) xfree(luf->vc_ptr); 1.377 + if (luf->vc_len != NULL) xfree(luf->vc_len); 1.378 + if (luf->vc_cap != NULL) xfree(luf->vc_cap); 1.379 + if (luf->pp_row != NULL) xfree(luf->pp_row); 1.380 + if (luf->pp_col != NULL) xfree(luf->pp_col); 1.381 + if (luf->qq_row != NULL) xfree(luf->qq_row); 1.382 + if (luf->qq_col != NULL) xfree(luf->qq_col); 1.383 + if (luf->sv_prev != NULL) xfree(luf->sv_prev); 1.384 + if (luf->sv_next != NULL) xfree(luf->sv_next); 1.385 + if (luf->vr_max != NULL) xfree(luf->vr_max); 1.386 + if (luf->rs_head != NULL) xfree(luf->rs_head); 1.387 + if (luf->rs_prev != NULL) xfree(luf->rs_prev); 1.388 + if (luf->rs_next != NULL) xfree(luf->rs_next); 1.389 + if (luf->cs_head != NULL) xfree(luf->cs_head); 1.390 + if (luf->cs_prev != NULL) xfree(luf->cs_prev); 1.391 + if (luf->cs_next != NULL) xfree(luf->cs_next); 1.392 + if (luf->flag != NULL) xfree(luf->flag); 1.393 + if (luf->work != NULL) xfree(luf->work); 1.394 + luf->n_max = n_max = n + 100; 1.395 + luf->fr_ptr = xcalloc(1+n_max, sizeof(int)); 1.396 + luf->fr_len = xcalloc(1+n_max, sizeof(int)); 1.397 + luf->fc_ptr = xcalloc(1+n_max, sizeof(int)); 1.398 + luf->fc_len = xcalloc(1+n_max, sizeof(int)); 1.399 + luf->vr_ptr = xcalloc(1+n_max, sizeof(int)); 1.400 + luf->vr_len = xcalloc(1+n_max, sizeof(int)); 1.401 + luf->vr_cap = xcalloc(1+n_max, sizeof(int)); 1.402 + luf->vr_piv = xcalloc(1+n_max, sizeof(double)); 1.403 + luf->vc_ptr = xcalloc(1+n_max, sizeof(int)); 1.404 + luf->vc_len = xcalloc(1+n_max, sizeof(int)); 1.405 + luf->vc_cap = xcalloc(1+n_max, sizeof(int)); 1.406 + luf->pp_row = xcalloc(1+n_max, sizeof(int)); 1.407 + luf->pp_col = xcalloc(1+n_max, sizeof(int)); 1.408 + luf->qq_row = xcalloc(1+n_max, sizeof(int)); 1.409 + luf->qq_col = xcalloc(1+n_max, sizeof(int)); 1.410 + luf->sv_prev = xcalloc(1+n_max+n_max, sizeof(int)); 1.411 + luf->sv_next = xcalloc(1+n_max+n_max, sizeof(int)); 1.412 + luf->vr_max = xcalloc(1+n_max, sizeof(double)); 1.413 + luf->rs_head = xcalloc(1+n_max, sizeof(int)); 1.414 + luf->rs_prev = xcalloc(1+n_max, sizeof(int)); 1.415 + luf->rs_next = xcalloc(1+n_max, sizeof(int)); 1.416 + luf->cs_head = xcalloc(1+n_max, sizeof(int)); 1.417 + luf->cs_prev = xcalloc(1+n_max, sizeof(int)); 1.418 + luf->cs_next = xcalloc(1+n_max, sizeof(int)); 1.419 + luf->flag = xcalloc(1+n_max, sizeof(int)); 1.420 + luf->work = xcalloc(1+n_max, sizeof(double)); 1.421 +done: return; 1.422 +} 1.423 + 1.424 +/*********************************************************************** 1.425 +* initialize - initialize LU-factorization data structures 1.426 +* 1.427 +* This routine initializes data structures for subsequent computing 1.428 +* the LU-factorization of a given matrix A, which is specified by the 1.429 +* formal routine col. On exit V = A and F = P = Q = I, where I is the 1.430 +* unity matrix. (Row-wise representation of the matrix F is not used 1.431 +* at the factorization stage and therefore is not initialized.) 1.432 +* 1.433 +* If no error occured, the routine returns zero. Otherwise, in case of 1.434 +* overflow of the sparse vector area, the routine returns non-zero. */ 1.435 + 1.436 +static int initialize(LUF *luf, int (*col)(void *info, int j, int rn[], 1.437 + double aj[]), void *info) 1.438 +{ int n = luf->n; 1.439 + int *fc_ptr = luf->fc_ptr; 1.440 + int *fc_len = luf->fc_len; 1.441 + int *vr_ptr = luf->vr_ptr; 1.442 + int *vr_len = luf->vr_len; 1.443 + int *vr_cap = luf->vr_cap; 1.444 + int *vc_ptr = luf->vc_ptr; 1.445 + int *vc_len = luf->vc_len; 1.446 + int *vc_cap = luf->vc_cap; 1.447 + int *pp_row = luf->pp_row; 1.448 + int *pp_col = luf->pp_col; 1.449 + int *qq_row = luf->qq_row; 1.450 + int *qq_col = luf->qq_col; 1.451 + int *sv_ind = luf->sv_ind; 1.452 + double *sv_val = luf->sv_val; 1.453 + int *sv_prev = luf->sv_prev; 1.454 + int *sv_next = luf->sv_next; 1.455 + double *vr_max = luf->vr_max; 1.456 + int *rs_head = luf->rs_head; 1.457 + int *rs_prev = luf->rs_prev; 1.458 + int *rs_next = luf->rs_next; 1.459 + int *cs_head = luf->cs_head; 1.460 + int *cs_prev = luf->cs_prev; 1.461 + int *cs_next = luf->cs_next; 1.462 + int *flag = luf->flag; 1.463 + double *work = luf->work; 1.464 + int ret = 0; 1.465 + int i, i_ptr, j, j_beg, j_end, k, len, nnz, sv_beg, sv_end, ptr; 1.466 + double big, val; 1.467 + /* free all locations of the sparse vector area */ 1.468 + sv_beg = 1; 1.469 + sv_end = luf->sv_size + 1; 1.470 + /* (row-wise representation of the matrix F is not initialized, 1.471 + because it is not used at the factorization stage) */ 1.472 + /* build the matrix F in column-wise format (initially F = I) */ 1.473 + for (j = 1; j <= n; j++) 1.474 + { fc_ptr[j] = sv_end; 1.475 + fc_len[j] = 0; 1.476 + } 1.477 + /* clear rows of the matrix V; clear the flag array */ 1.478 + for (i = 1; i <= n; i++) 1.479 + vr_len[i] = vr_cap[i] = 0, flag[i] = 0; 1.480 + /* build the matrix V in column-wise format (initially V = A); 1.481 + count non-zeros in rows of this matrix; count total number of 1.482 + non-zeros; compute largest of absolute values of elements */ 1.483 + nnz = 0; 1.484 + big = 0.0; 1.485 + for (j = 1; j <= n; j++) 1.486 + { int *rn = pp_row; 1.487 + double *aj = work; 1.488 + /* obtain j-th column of the matrix A */ 1.489 + len = col(info, j, rn, aj); 1.490 + if (!(0 <= len && len <= n)) 1.491 + xfault("luf_factorize: j = %d; len = %d; invalid column len" 1.492 + "gth\n", j, len); 1.493 + /* check for free locations */ 1.494 + if (sv_end - sv_beg < len) 1.495 + { /* overflow of the sparse vector area */ 1.496 + ret = 1; 1.497 + goto done; 1.498 + } 1.499 + /* set pointer to the j-th column */ 1.500 + vc_ptr[j] = sv_beg; 1.501 + /* set length of the j-th column */ 1.502 + vc_len[j] = vc_cap[j] = len; 1.503 + /* count total number of non-zeros */ 1.504 + nnz += len; 1.505 + /* walk through elements of the j-th column */ 1.506 + for (ptr = 1; ptr <= len; ptr++) 1.507 + { /* get row index and numerical value of a[i,j] */ 1.508 + i = rn[ptr]; 1.509 + val = aj[ptr]; 1.510 + if (!(1 <= i && i <= n)) 1.511 + xfault("luf_factorize: i = %d; j = %d; invalid row index" 1.512 + "\n", i, j); 1.513 + if (flag[i]) 1.514 + xfault("luf_factorize: i = %d; j = %d; duplicate element" 1.515 + " not allowed\n", i, j); 1.516 + if (val == 0.0) 1.517 + xfault("luf_factorize: i = %d; j = %d; zero element not " 1.518 + "allowed\n", i, j); 1.519 + /* add new element v[i,j] = a[i,j] to j-th column */ 1.520 + sv_ind[sv_beg] = i; 1.521 + sv_val[sv_beg] = val; 1.522 + sv_beg++; 1.523 + /* big := max(big, |a[i,j]|) */ 1.524 + if (val < 0.0) val = - val; 1.525 + if (big < val) big = val; 1.526 + /* mark non-zero in the i-th position of the j-th column */ 1.527 + flag[i] = 1; 1.528 + /* increase length of the i-th row */ 1.529 + vr_cap[i]++; 1.530 + } 1.531 + /* reset all non-zero marks */ 1.532 + for (ptr = 1; ptr <= len; ptr++) flag[rn[ptr]] = 0; 1.533 + } 1.534 + /* allocate rows of the matrix V */ 1.535 + for (i = 1; i <= n; i++) 1.536 + { /* get length of the i-th row */ 1.537 + len = vr_cap[i]; 1.538 + /* check for free locations */ 1.539 + if (sv_end - sv_beg < len) 1.540 + { /* overflow of the sparse vector area */ 1.541 + ret = 1; 1.542 + goto done; 1.543 + } 1.544 + /* set pointer to the i-th row */ 1.545 + vr_ptr[i] = sv_beg; 1.546 + /* reserve locations for the i-th row */ 1.547 + sv_beg += len; 1.548 + } 1.549 + /* build the matrix V in row-wise format using representation of 1.550 + this matrix in column-wise format */ 1.551 + for (j = 1; j <= n; j++) 1.552 + { /* walk through elements of the j-th column */ 1.553 + j_beg = vc_ptr[j]; 1.554 + j_end = j_beg + vc_len[j] - 1; 1.555 + for (k = j_beg; k <= j_end; k++) 1.556 + { /* get row index and numerical value of v[i,j] */ 1.557 + i = sv_ind[k]; 1.558 + val = sv_val[k]; 1.559 + /* store element in the i-th row */ 1.560 + i_ptr = vr_ptr[i] + vr_len[i]; 1.561 + sv_ind[i_ptr] = j; 1.562 + sv_val[i_ptr] = val; 1.563 + /* increase count of the i-th row */ 1.564 + vr_len[i]++; 1.565 + } 1.566 + } 1.567 + /* initialize the matrices P and Q (initially P = Q = I) */ 1.568 + for (k = 1; k <= n; k++) 1.569 + pp_row[k] = pp_col[k] = qq_row[k] = qq_col[k] = k; 1.570 + /* set sva partitioning pointers */ 1.571 + luf->sv_beg = sv_beg; 1.572 + luf->sv_end = sv_end; 1.573 + /* the initial physical order of rows and columns of the matrix V 1.574 + is n+1, ..., n+n, 1, ..., n (firstly columns, then rows) */ 1.575 + luf->sv_head = n+1; 1.576 + luf->sv_tail = n; 1.577 + for (i = 1; i <= n; i++) 1.578 + { sv_prev[i] = i-1; 1.579 + sv_next[i] = i+1; 1.580 + } 1.581 + sv_prev[1] = n+n; 1.582 + sv_next[n] = 0; 1.583 + for (j = 1; j <= n; j++) 1.584 + { sv_prev[n+j] = n+j-1; 1.585 + sv_next[n+j] = n+j+1; 1.586 + } 1.587 + sv_prev[n+1] = 0; 1.588 + sv_next[n+n] = 1; 1.589 + /* clear working arrays */ 1.590 + for (k = 1; k <= n; k++) 1.591 + { flag[k] = 0; 1.592 + work[k] = 0.0; 1.593 + } 1.594 + /* initialize some statistics */ 1.595 + luf->nnz_a = nnz; 1.596 + luf->nnz_f = 0; 1.597 + luf->nnz_v = nnz; 1.598 + luf->max_a = big; 1.599 + luf->big_v = big; 1.600 + luf->rank = -1; 1.601 + /* initially the active submatrix is the entire matrix V */ 1.602 + /* largest of absolute values of elements in each active row is 1.603 + unknown yet */ 1.604 + for (i = 1; i <= n; i++) vr_max[i] = -1.0; 1.605 + /* build linked lists of active rows */ 1.606 + for (len = 0; len <= n; len++) rs_head[len] = 0; 1.607 + for (i = 1; i <= n; i++) 1.608 + { len = vr_len[i]; 1.609 + rs_prev[i] = 0; 1.610 + rs_next[i] = rs_head[len]; 1.611 + if (rs_next[i] != 0) rs_prev[rs_next[i]] = i; 1.612 + rs_head[len] = i; 1.613 + } 1.614 + /* build linked lists of active columns */ 1.615 + for (len = 0; len <= n; len++) cs_head[len] = 0; 1.616 + for (j = 1; j <= n; j++) 1.617 + { len = vc_len[j]; 1.618 + cs_prev[j] = 0; 1.619 + cs_next[j] = cs_head[len]; 1.620 + if (cs_next[j] != 0) cs_prev[cs_next[j]] = j; 1.621 + cs_head[len] = j; 1.622 + } 1.623 +done: /* return to the factorizing routine */ 1.624 + return ret; 1.625 +} 1.626 + 1.627 +/*********************************************************************** 1.628 +* find_pivot - choose a pivot element 1.629 +* 1.630 +* This routine chooses a pivot element in the active submatrix of the 1.631 +* matrix U = P*V*Q. 1.632 +* 1.633 +* It is assumed that on entry the matrix U has the following partially 1.634 +* triangularized form: 1.635 +* 1.636 +* 1 k n 1.637 +* 1 x x x x x x x x x x 1.638 +* . x x x x x x x x x 1.639 +* . . x x x x x x x x 1.640 +* . . . x x x x x x x 1.641 +* k . . . . * * * * * * 1.642 +* . . . . * * * * * * 1.643 +* . . . . * * * * * * 1.644 +* . . . . * * * * * * 1.645 +* . . . . * * * * * * 1.646 +* n . . . . * * * * * * 1.647 +* 1.648 +* where rows and columns k, k+1, ..., n belong to the active submatrix 1.649 +* (elements of the active submatrix are marked by '*'). 1.650 +* 1.651 +* Since the matrix U = P*V*Q is not stored, the routine works with the 1.652 +* matrix V. It is assumed that the row-wise representation corresponds 1.653 +* to the matrix V, but the column-wise representation corresponds to 1.654 +* the active submatrix of the matrix V, i.e. elements of the matrix V, 1.655 +* which doesn't belong to the active submatrix, are missing from the 1.656 +* column linked lists. It is also assumed that each active row of the 1.657 +* matrix V is in the set R[len], where len is number of non-zeros in 1.658 +* the row, and each active column of the matrix V is in the set C[len], 1.659 +* where len is number of non-zeros in the column (in the latter case 1.660 +* only elements of the active submatrix are counted; such elements are 1.661 +* marked by '*' on the figure above). 1.662 +* 1.663 +* For the reason of numerical stability the routine applies so called 1.664 +* threshold pivoting proposed by J.Reid. It is assumed that an element 1.665 +* v[i,j] can be selected as a pivot candidate if it is not very small 1.666 +* (in absolute value) among other elements in the same row, i.e. if it 1.667 +* satisfies to the stability condition |v[i,j]| >= tol * max|v[i,*]|, 1.668 +* where 0 < tol < 1 is a given tolerance. 1.669 +* 1.670 +* In order to keep sparsity of the matrix V the routine uses Markowitz 1.671 +* strategy, trying to choose such element v[p,q], which satisfies to 1.672 +* the stability condition (see above) and has smallest Markowitz cost 1.673 +* (nr[p]-1) * (nc[q]-1), where nr[p] and nc[q] are numbers of non-zero 1.674 +* elements, respectively, in the p-th row and in the q-th column of the 1.675 +* active submatrix. 1.676 +* 1.677 +* In order to reduce the search, i.e. not to walk through all elements 1.678 +* of the active submatrix, the routine exploits a technique proposed by 1.679 +* I.Duff. This technique is based on using the sets R[len] and C[len] 1.680 +* of active rows and columns. 1.681 +* 1.682 +* If the pivot element v[p,q] has been chosen, the routine stores its 1.683 +* indices to the locations *p and *q and returns zero. Otherwise, if 1.684 +* the active submatrix is empty and therefore the pivot element can't 1.685 +* be chosen, the routine returns non-zero. */ 1.686 + 1.687 +static int find_pivot(LUF *luf, int *_p, int *_q) 1.688 +{ int n = luf->n; 1.689 + int *vr_ptr = luf->vr_ptr; 1.690 + int *vr_len = luf->vr_len; 1.691 + int *vc_ptr = luf->vc_ptr; 1.692 + int *vc_len = luf->vc_len; 1.693 + int *sv_ind = luf->sv_ind; 1.694 + double *sv_val = luf->sv_val; 1.695 + double *vr_max = luf->vr_max; 1.696 + int *rs_head = luf->rs_head; 1.697 + int *rs_next = luf->rs_next; 1.698 + int *cs_head = luf->cs_head; 1.699 + int *cs_prev = luf->cs_prev; 1.700 + int *cs_next = luf->cs_next; 1.701 + double piv_tol = luf->piv_tol; 1.702 + int piv_lim = luf->piv_lim; 1.703 + int suhl = luf->suhl; 1.704 + int p, q, len, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr, 1.705 + ncand, next_j, min_p, min_q, min_len; 1.706 + double best, cost, big, temp; 1.707 + /* initially no pivot candidates have been found so far */ 1.708 + p = q = 0, best = DBL_MAX, ncand = 0; 1.709 + /* if in the active submatrix there is a column that has the only 1.710 + non-zero (column singleton), choose it as pivot */ 1.711 + j = cs_head[1]; 1.712 + if (j != 0) 1.713 + { xassert(vc_len[j] == 1); 1.714 + p = sv_ind[vc_ptr[j]], q = j; 1.715 + goto done; 1.716 + } 1.717 + /* if in the active submatrix there is a row that has the only 1.718 + non-zero (row singleton), choose it as pivot */ 1.719 + i = rs_head[1]; 1.720 + if (i != 0) 1.721 + { xassert(vr_len[i] == 1); 1.722 + p = i, q = sv_ind[vr_ptr[i]]; 1.723 + goto done; 1.724 + } 1.725 + /* there are no singletons in the active submatrix; walk through 1.726 + other non-empty rows and columns */ 1.727 + for (len = 2; len <= n; len++) 1.728 + { /* consider active columns that have len non-zeros */ 1.729 + for (j = cs_head[len]; j != 0; j = next_j) 1.730 + { /* the j-th column has len non-zeros */ 1.731 + j_beg = vc_ptr[j]; 1.732 + j_end = j_beg + vc_len[j] - 1; 1.733 + /* save pointer to the next column with the same length */ 1.734 + next_j = cs_next[j]; 1.735 + /* find an element in the j-th column, which is placed in a 1.736 + row with minimal number of non-zeros and satisfies to the 1.737 + stability condition (such element may not exist) */ 1.738 + min_p = min_q = 0, min_len = INT_MAX; 1.739 + for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) 1.740 + { /* get row index of v[i,j] */ 1.741 + i = sv_ind[j_ptr]; 1.742 + i_beg = vr_ptr[i]; 1.743 + i_end = i_beg + vr_len[i] - 1; 1.744 + /* if the i-th row is not shorter than that one, where 1.745 + minimal element is currently placed, skip v[i,j] */ 1.746 + if (vr_len[i] >= min_len) continue; 1.747 + /* determine the largest of absolute values of elements 1.748 + in the i-th row */ 1.749 + big = vr_max[i]; 1.750 + if (big < 0.0) 1.751 + { /* the largest value is unknown yet; compute it */ 1.752 + for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) 1.753 + { temp = sv_val[i_ptr]; 1.754 + if (temp < 0.0) temp = - temp; 1.755 + if (big < temp) big = temp; 1.756 + } 1.757 + vr_max[i] = big; 1.758 + } 1.759 + /* find v[i,j] in the i-th row */ 1.760 + for (i_ptr = vr_ptr[i]; sv_ind[i_ptr] != j; i_ptr++); 1.761 + xassert(i_ptr <= i_end); 1.762 + /* if v[i,j] doesn't satisfy to the stability condition, 1.763 + skip it */ 1.764 + temp = sv_val[i_ptr]; 1.765 + if (temp < 0.0) temp = - temp; 1.766 + if (temp < piv_tol * big) continue; 1.767 + /* v[i,j] is better than the current minimal element */ 1.768 + min_p = i, min_q = j, min_len = vr_len[i]; 1.769 + /* if Markowitz cost of the current minimal element is 1.770 + not greater than (len-1)**2, it can be chosen right 1.771 + now; this heuristic reduces the search and works well 1.772 + in many cases */ 1.773 + if (min_len <= len) 1.774 + { p = min_p, q = min_q; 1.775 + goto done; 1.776 + } 1.777 + } 1.778 + /* the j-th column has been scanned */ 1.779 + if (min_p != 0) 1.780 + { /* the minimal element is a next pivot candidate */ 1.781 + ncand++; 1.782 + /* compute its Markowitz cost */ 1.783 + cost = (double)(min_len - 1) * (double)(len - 1); 1.784 + /* choose between the minimal element and the current 1.785 + candidate */ 1.786 + if (cost < best) p = min_p, q = min_q, best = cost; 1.787 + /* if piv_lim candidates have been considered, there are 1.788 + doubts that a much better candidate exists; therefore 1.789 + it's time to terminate the search */ 1.790 + if (ncand == piv_lim) goto done; 1.791 + } 1.792 + else 1.793 + { /* the j-th column has no elements, which satisfy to the 1.794 + stability condition; Uwe Suhl suggests to exclude such 1.795 + column from the further consideration until it becomes 1.796 + a column singleton; in hard cases this significantly 1.797 + reduces a time needed for pivot searching */ 1.798 + if (suhl) 1.799 + { /* remove the j-th column from the active set */ 1.800 + if (cs_prev[j] == 0) 1.801 + cs_head[len] = cs_next[j]; 1.802 + else 1.803 + cs_next[cs_prev[j]] = cs_next[j]; 1.804 + if (cs_next[j] == 0) 1.805 + /* nop */; 1.806 + else 1.807 + cs_prev[cs_next[j]] = cs_prev[j]; 1.808 + /* the following assignment is used to avoid an error 1.809 + when the routine eliminate (see below) will try to 1.810 + remove the j-th column from the active set */ 1.811 + cs_prev[j] = cs_next[j] = j; 1.812 + } 1.813 + } 1.814 + } 1.815 + /* consider active rows that have len non-zeros */ 1.816 + for (i = rs_head[len]; i != 0; i = rs_next[i]) 1.817 + { /* the i-th row has len non-zeros */ 1.818 + i_beg = vr_ptr[i]; 1.819 + i_end = i_beg + vr_len[i] - 1; 1.820 + /* determine the largest of absolute values of elements in 1.821 + the i-th row */ 1.822 + big = vr_max[i]; 1.823 + if (big < 0.0) 1.824 + { /* the largest value is unknown yet; compute it */ 1.825 + for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) 1.826 + { temp = sv_val[i_ptr]; 1.827 + if (temp < 0.0) temp = - temp; 1.828 + if (big < temp) big = temp; 1.829 + } 1.830 + vr_max[i] = big; 1.831 + } 1.832 + /* find an element in the i-th row, which is placed in a 1.833 + column with minimal number of non-zeros and satisfies to 1.834 + the stability condition (such element always exists) */ 1.835 + min_p = min_q = 0, min_len = INT_MAX; 1.836 + for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) 1.837 + { /* get column index of v[i,j] */ 1.838 + j = sv_ind[i_ptr]; 1.839 + /* if the j-th column is not shorter than that one, where 1.840 + minimal element is currently placed, skip v[i,j] */ 1.841 + if (vc_len[j] >= min_len) continue; 1.842 + /* if v[i,j] doesn't satisfy to the stability condition, 1.843 + skip it */ 1.844 + temp = sv_val[i_ptr]; 1.845 + if (temp < 0.0) temp = - temp; 1.846 + if (temp < piv_tol * big) continue; 1.847 + /* v[i,j] is better than the current minimal element */ 1.848 + min_p = i, min_q = j, min_len = vc_len[j]; 1.849 + /* if Markowitz cost of the current minimal element is 1.850 + not greater than (len-1)**2, it can be chosen right 1.851 + now; this heuristic reduces the search and works well 1.852 + in many cases */ 1.853 + if (min_len <= len) 1.854 + { p = min_p, q = min_q; 1.855 + goto done; 1.856 + } 1.857 + } 1.858 + /* the i-th row has been scanned */ 1.859 + if (min_p != 0) 1.860 + { /* the minimal element is a next pivot candidate */ 1.861 + ncand++; 1.862 + /* compute its Markowitz cost */ 1.863 + cost = (double)(len - 1) * (double)(min_len - 1); 1.864 + /* choose between the minimal element and the current 1.865 + candidate */ 1.866 + if (cost < best) p = min_p, q = min_q, best = cost; 1.867 + /* if piv_lim candidates have been considered, there are 1.868 + doubts that a much better candidate exists; therefore 1.869 + it's time to terminate the search */ 1.870 + if (ncand == piv_lim) goto done; 1.871 + } 1.872 + else 1.873 + { /* this can't be because this can never be */ 1.874 + xassert(min_p != min_p); 1.875 + } 1.876 + } 1.877 + } 1.878 +done: /* bring the pivot to the factorizing routine */ 1.879 + *_p = p, *_q = q; 1.880 + return (p == 0); 1.881 +} 1.882 + 1.883 +/*********************************************************************** 1.884 +* eliminate - perform gaussian elimination. 1.885 +* 1.886 +* This routine performs elementary gaussian transformations in order 1.887 +* to eliminate subdiagonal elements in the k-th column of the matrix 1.888 +* U = P*V*Q using the pivot element u[k,k], where k is the number of 1.889 +* the current elimination step. 1.890 +* 1.891 +* The parameters p and q are, respectively, row and column indices of 1.892 +* the element v[p,q], which corresponds to the element u[k,k]. 1.893 +* 1.894 +* Each time when the routine applies the elementary transformation to 1.895 +* a non-pivot row of the matrix V, it stores the corresponding element 1.896 +* to the matrix F in order to keep the main equality A = F*V. 1.897 +* 1.898 +* The routine assumes that on entry the matrices L = P*F*inv(P) and 1.899 +* U = P*V*Q are the following: 1.900 +* 1.901 +* 1 k 1 k n 1.902 +* 1 1 . . . . . . . . . 1 x x x x x x x x x x 1.903 +* x 1 . . . . . . . . . x x x x x x x x x 1.904 +* x x 1 . . . . . . . . . x x x x x x x x 1.905 +* x x x 1 . . . . . . . . . x x x x x x x 1.906 +* k x x x x 1 . . . . . k . . . . * * * * * * 1.907 +* x x x x _ 1 . . . . . . . . # * * * * * 1.908 +* x x x x _ . 1 . . . . . . . # * * * * * 1.909 +* x x x x _ . . 1 . . . . . . # * * * * * 1.910 +* x x x x _ . . . 1 . . . . . # * * * * * 1.911 +* n x x x x _ . . . . 1 n . . . . # * * * * * 1.912 +* 1.913 +* matrix L matrix U 1.914 +* 1.915 +* where rows and columns of the matrix U with numbers k, k+1, ..., n 1.916 +* form the active submatrix (eliminated elements are marked by '#' and 1.917 +* other elements of the active submatrix are marked by '*'). Note that 1.918 +* each eliminated non-zero element u[i,k] of the matrix U gives the 1.919 +* corresponding element l[i,k] of the matrix L (marked by '_'). 1.920 +* 1.921 +* Actually all operations are performed on the matrix V. Should note 1.922 +* that the row-wise representation corresponds to the matrix V, but the 1.923 +* column-wise representation corresponds to the active submatrix of the 1.924 +* matrix V, i.e. elements of the matrix V, which doesn't belong to the 1.925 +* active submatrix, are missing from the column linked lists. 1.926 +* 1.927 +* Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal 1.928 +* elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies 1.929 +* the following elementary gaussian transformations: 1.930 +* 1.931 +* (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V), 1.932 +* 1.933 +* where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier. 1.934 +* 1.935 +* Additionally, in order to keep the main equality A = F*V, each time 1.936 +* when the routine applies the transformation to i-th row of the matrix 1.937 +* V, it also adds f[i,p] as a new element to the matrix F. 1.938 +* 1.939 +* IMPORTANT: On entry the working arrays flag and work should contain 1.940 +* zeros. This status is provided by the routine on exit. 1.941 +* 1.942 +* If no error occured, the routine returns zero. Otherwise, in case of 1.943 +* overflow of the sparse vector area, the routine returns non-zero. */ 1.944 + 1.945 +static int eliminate(LUF *luf, int p, int q) 1.946 +{ int n = luf->n; 1.947 + int *fc_ptr = luf->fc_ptr; 1.948 + int *fc_len = luf->fc_len; 1.949 + int *vr_ptr = luf->vr_ptr; 1.950 + int *vr_len = luf->vr_len; 1.951 + int *vr_cap = luf->vr_cap; 1.952 + double *vr_piv = luf->vr_piv; 1.953 + int *vc_ptr = luf->vc_ptr; 1.954 + int *vc_len = luf->vc_len; 1.955 + int *vc_cap = luf->vc_cap; 1.956 + int *sv_ind = luf->sv_ind; 1.957 + double *sv_val = luf->sv_val; 1.958 + int *sv_prev = luf->sv_prev; 1.959 + int *sv_next = luf->sv_next; 1.960 + double *vr_max = luf->vr_max; 1.961 + int *rs_head = luf->rs_head; 1.962 + int *rs_prev = luf->rs_prev; 1.963 + int *rs_next = luf->rs_next; 1.964 + int *cs_head = luf->cs_head; 1.965 + int *cs_prev = luf->cs_prev; 1.966 + int *cs_next = luf->cs_next; 1.967 + int *flag = luf->flag; 1.968 + double *work = luf->work; 1.969 + double eps_tol = luf->eps_tol; 1.970 + /* at this stage the row-wise representation of the matrix F is 1.971 + not used, so fr_len can be used as a working array */ 1.972 + int *ndx = luf->fr_len; 1.973 + int ret = 0; 1.974 + int len, fill, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr, k, 1.975 + p_beg, p_end, p_ptr, q_beg, q_end, q_ptr; 1.976 + double fip, val, vpq, temp; 1.977 + xassert(1 <= p && p <= n); 1.978 + xassert(1 <= q && q <= n); 1.979 + /* remove the p-th (pivot) row from the active set; this row will 1.980 + never return there */ 1.981 + if (rs_prev[p] == 0) 1.982 + rs_head[vr_len[p]] = rs_next[p]; 1.983 + else 1.984 + rs_next[rs_prev[p]] = rs_next[p]; 1.985 + if (rs_next[p] == 0) 1.986 + ; 1.987 + else 1.988 + rs_prev[rs_next[p]] = rs_prev[p]; 1.989 + /* remove the q-th (pivot) column from the active set; this column 1.990 + will never return there */ 1.991 + if (cs_prev[q] == 0) 1.992 + cs_head[vc_len[q]] = cs_next[q]; 1.993 + else 1.994 + cs_next[cs_prev[q]] = cs_next[q]; 1.995 + if (cs_next[q] == 0) 1.996 + ; 1.997 + else 1.998 + cs_prev[cs_next[q]] = cs_prev[q]; 1.999 + /* find the pivot v[p,q] = u[k,k] in the p-th row */ 1.1000 + p_beg = vr_ptr[p]; 1.1001 + p_end = p_beg + vr_len[p] - 1; 1.1002 + for (p_ptr = p_beg; sv_ind[p_ptr] != q; p_ptr++) /* nop */; 1.1003 + xassert(p_ptr <= p_end); 1.1004 + /* store value of the pivot */ 1.1005 + vpq = (vr_piv[p] = sv_val[p_ptr]); 1.1006 + /* remove the pivot from the p-th row */ 1.1007 + sv_ind[p_ptr] = sv_ind[p_end]; 1.1008 + sv_val[p_ptr] = sv_val[p_end]; 1.1009 + vr_len[p]--; 1.1010 + p_end--; 1.1011 + /* find the pivot v[p,q] = u[k,k] in the q-th column */ 1.1012 + q_beg = vc_ptr[q]; 1.1013 + q_end = q_beg + vc_len[q] - 1; 1.1014 + for (q_ptr = q_beg; sv_ind[q_ptr] != p; q_ptr++) /* nop */; 1.1015 + xassert(q_ptr <= q_end); 1.1016 + /* remove the pivot from the q-th column */ 1.1017 + sv_ind[q_ptr] = sv_ind[q_end]; 1.1018 + vc_len[q]--; 1.1019 + q_end--; 1.1020 + /* walk through the p-th (pivot) row, which doesn't contain the 1.1021 + pivot v[p,q] already, and do the following... */ 1.1022 + for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++) 1.1023 + { /* get column index of v[p,j] */ 1.1024 + j = sv_ind[p_ptr]; 1.1025 + /* store v[p,j] to the working array */ 1.1026 + flag[j] = 1; 1.1027 + work[j] = sv_val[p_ptr]; 1.1028 + /* remove the j-th column from the active set; this column will 1.1029 + return there later with new length */ 1.1030 + if (cs_prev[j] == 0) 1.1031 + cs_head[vc_len[j]] = cs_next[j]; 1.1032 + else 1.1033 + cs_next[cs_prev[j]] = cs_next[j]; 1.1034 + if (cs_next[j] == 0) 1.1035 + ; 1.1036 + else 1.1037 + cs_prev[cs_next[j]] = cs_prev[j]; 1.1038 + /* find v[p,j] in the j-th column */ 1.1039 + j_beg = vc_ptr[j]; 1.1040 + j_end = j_beg + vc_len[j] - 1; 1.1041 + for (j_ptr = j_beg; sv_ind[j_ptr] != p; j_ptr++) /* nop */; 1.1042 + xassert(j_ptr <= j_end); 1.1043 + /* since v[p,j] leaves the active submatrix, remove it from the 1.1044 + j-th column; however, v[p,j] is kept in the p-th row */ 1.1045 + sv_ind[j_ptr] = sv_ind[j_end]; 1.1046 + vc_len[j]--; 1.1047 + } 1.1048 + /* walk through the q-th (pivot) column, which doesn't contain the 1.1049 + pivot v[p,q] already, and perform gaussian elimination */ 1.1050 + while (q_beg <= q_end) 1.1051 + { /* element v[i,q] should be eliminated */ 1.1052 + /* get row index of v[i,q] */ 1.1053 + i = sv_ind[q_beg]; 1.1054 + /* remove the i-th row from the active set; later this row will 1.1055 + return there with new length */ 1.1056 + if (rs_prev[i] == 0) 1.1057 + rs_head[vr_len[i]] = rs_next[i]; 1.1058 + else 1.1059 + rs_next[rs_prev[i]] = rs_next[i]; 1.1060 + if (rs_next[i] == 0) 1.1061 + ; 1.1062 + else 1.1063 + rs_prev[rs_next[i]] = rs_prev[i]; 1.1064 + /* find v[i,q] in the i-th row */ 1.1065 + i_beg = vr_ptr[i]; 1.1066 + i_end = i_beg + vr_len[i] - 1; 1.1067 + for (i_ptr = i_beg; sv_ind[i_ptr] != q; i_ptr++) /* nop */; 1.1068 + xassert(i_ptr <= i_end); 1.1069 + /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] */ 1.1070 + fip = sv_val[i_ptr] / vpq; 1.1071 + /* since v[i,q] should be eliminated, remove it from the i-th 1.1072 + row */ 1.1073 + sv_ind[i_ptr] = sv_ind[i_end]; 1.1074 + sv_val[i_ptr] = sv_val[i_end]; 1.1075 + vr_len[i]--; 1.1076 + i_end--; 1.1077 + /* and from the q-th column */ 1.1078 + sv_ind[q_beg] = sv_ind[q_end]; 1.1079 + vc_len[q]--; 1.1080 + q_end--; 1.1081 + /* perform gaussian transformation: 1.1082 + (i-th row) := (i-th row) - f[i,p] * (p-th row) 1.1083 + note that now the p-th row, which is in the working array, 1.1084 + doesn't contain the pivot v[p,q], and the i-th row doesn't 1.1085 + contain the eliminated element v[i,q] */ 1.1086 + /* walk through the i-th row and transform existing non-zero 1.1087 + elements */ 1.1088 + fill = vr_len[p]; 1.1089 + for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) 1.1090 + { /* get column index of v[i,j] */ 1.1091 + j = sv_ind[i_ptr]; 1.1092 + /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */ 1.1093 + if (flag[j]) 1.1094 + { /* v[p,j] != 0 */ 1.1095 + temp = (sv_val[i_ptr] -= fip * work[j]); 1.1096 + if (temp < 0.0) temp = - temp; 1.1097 + flag[j] = 0; 1.1098 + fill--; /* since both v[i,j] and v[p,j] exist */ 1.1099 + if (temp == 0.0 || temp < eps_tol) 1.1100 + { /* new v[i,j] is closer to zero; replace it by exact 1.1101 + zero, i.e. remove it from the active submatrix */ 1.1102 + /* remove v[i,j] from the i-th row */ 1.1103 + sv_ind[i_ptr] = sv_ind[i_end]; 1.1104 + sv_val[i_ptr] = sv_val[i_end]; 1.1105 + vr_len[i]--; 1.1106 + i_ptr--; 1.1107 + i_end--; 1.1108 + /* find v[i,j] in the j-th column */ 1.1109 + j_beg = vc_ptr[j]; 1.1110 + j_end = j_beg + vc_len[j] - 1; 1.1111 + for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++); 1.1112 + xassert(j_ptr <= j_end); 1.1113 + /* remove v[i,j] from the j-th column */ 1.1114 + sv_ind[j_ptr] = sv_ind[j_end]; 1.1115 + vc_len[j]--; 1.1116 + } 1.1117 + else 1.1118 + { /* v_big := max(v_big, |v[i,j]|) */ 1.1119 + if (luf->big_v < temp) luf->big_v = temp; 1.1120 + } 1.1121 + } 1.1122 + } 1.1123 + /* now flag is the pattern of the set v[p,*] \ v[i,*], and fill 1.1124 + is number of non-zeros in this set; therefore up to fill new 1.1125 + non-zeros may appear in the i-th row */ 1.1126 + if (vr_len[i] + fill > vr_cap[i]) 1.1127 + { /* enlarge the i-th row */ 1.1128 + if (luf_enlarge_row(luf, i, vr_len[i] + fill)) 1.1129 + { /* overflow of the sparse vector area */ 1.1130 + ret = 1; 1.1131 + goto done; 1.1132 + } 1.1133 + /* defragmentation may change row and column pointers of the 1.1134 + matrix V */ 1.1135 + p_beg = vr_ptr[p]; 1.1136 + p_end = p_beg + vr_len[p] - 1; 1.1137 + q_beg = vc_ptr[q]; 1.1138 + q_end = q_beg + vc_len[q] - 1; 1.1139 + } 1.1140 + /* walk through the p-th (pivot) row and create new elements 1.1141 + of the i-th row that appear due to fill-in; column indices 1.1142 + of these new elements are accumulated in the array ndx */ 1.1143 + len = 0; 1.1144 + for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++) 1.1145 + { /* get column index of v[p,j], which may cause fill-in */ 1.1146 + j = sv_ind[p_ptr]; 1.1147 + if (flag[j]) 1.1148 + { /* compute new non-zero v[i,j] = 0 - f[i,p] * v[p,j] */ 1.1149 + temp = (val = - fip * work[j]); 1.1150 + if (temp < 0.0) temp = - temp; 1.1151 + if (temp == 0.0 || temp < eps_tol) 1.1152 + /* if v[i,j] is closer to zero; just ignore it */; 1.1153 + else 1.1154 + { /* add v[i,j] to the i-th row */ 1.1155 + i_ptr = vr_ptr[i] + vr_len[i]; 1.1156 + sv_ind[i_ptr] = j; 1.1157 + sv_val[i_ptr] = val; 1.1158 + vr_len[i]++; 1.1159 + /* remember column index of v[i,j] */ 1.1160 + ndx[++len] = j; 1.1161 + /* big_v := max(big_v, |v[i,j]|) */ 1.1162 + if (luf->big_v < temp) luf->big_v = temp; 1.1163 + } 1.1164 + } 1.1165 + else 1.1166 + { /* there is no fill-in, because v[i,j] already exists in 1.1167 + the i-th row; restore the flag of the element v[p,j], 1.1168 + which was reset before */ 1.1169 + flag[j] = 1; 1.1170 + } 1.1171 + } 1.1172 + /* add new non-zeros v[i,j] to the corresponding columns */ 1.1173 + for (k = 1; k <= len; k++) 1.1174 + { /* get column index of new non-zero v[i,j] */ 1.1175 + j = ndx[k]; 1.1176 + /* one free location is needed in the j-th column */ 1.1177 + if (vc_len[j] + 1 > vc_cap[j]) 1.1178 + { /* enlarge the j-th column */ 1.1179 + if (luf_enlarge_col(luf, j, vc_len[j] + 10)) 1.1180 + { /* overflow of the sparse vector area */ 1.1181 + ret = 1; 1.1182 + goto done; 1.1183 + } 1.1184 + /* defragmentation may change row and column pointers of 1.1185 + the matrix V */ 1.1186 + p_beg = vr_ptr[p]; 1.1187 + p_end = p_beg + vr_len[p] - 1; 1.1188 + q_beg = vc_ptr[q]; 1.1189 + q_end = q_beg + vc_len[q] - 1; 1.1190 + } 1.1191 + /* add new non-zero v[i,j] to the j-th column */ 1.1192 + j_ptr = vc_ptr[j] + vc_len[j]; 1.1193 + sv_ind[j_ptr] = i; 1.1194 + vc_len[j]++; 1.1195 + } 1.1196 + /* now the i-th row has been completely transformed, therefore 1.1197 + it can return to the active set with new length */ 1.1198 + rs_prev[i] = 0; 1.1199 + rs_next[i] = rs_head[vr_len[i]]; 1.1200 + if (rs_next[i] != 0) rs_prev[rs_next[i]] = i; 1.1201 + rs_head[vr_len[i]] = i; 1.1202 + /* the largest of absolute values of elements in the i-th row 1.1203 + is currently unknown */ 1.1204 + vr_max[i] = -1.0; 1.1205 + /* at least one free location is needed to store the gaussian 1.1206 + multiplier */ 1.1207 + if (luf->sv_end - luf->sv_beg < 1) 1.1208 + { /* there are no free locations at all; defragment SVA */ 1.1209 + luf_defrag_sva(luf); 1.1210 + if (luf->sv_end - luf->sv_beg < 1) 1.1211 + { /* overflow of the sparse vector area */ 1.1212 + ret = 1; 1.1213 + goto done; 1.1214 + } 1.1215 + /* defragmentation may change row and column pointers of the 1.1216 + matrix V */ 1.1217 + p_beg = vr_ptr[p]; 1.1218 + p_end = p_beg + vr_len[p] - 1; 1.1219 + q_beg = vc_ptr[q]; 1.1220 + q_end = q_beg + vc_len[q] - 1; 1.1221 + } 1.1222 + /* add the element f[i,p], which is the gaussian multiplier, 1.1223 + to the matrix F */ 1.1224 + luf->sv_end--; 1.1225 + sv_ind[luf->sv_end] = i; 1.1226 + sv_val[luf->sv_end] = fip; 1.1227 + fc_len[p]++; 1.1228 + /* end of elimination loop */ 1.1229 + } 1.1230 + /* at this point the q-th (pivot) column should be empty */ 1.1231 + xassert(vc_len[q] == 0); 1.1232 + /* reset capacity of the q-th column */ 1.1233 + vc_cap[q] = 0; 1.1234 + /* remove node of the q-th column from the addressing list */ 1.1235 + k = n + q; 1.1236 + if (sv_prev[k] == 0) 1.1237 + luf->sv_head = sv_next[k]; 1.1238 + else 1.1239 + sv_next[sv_prev[k]] = sv_next[k]; 1.1240 + if (sv_next[k] == 0) 1.1241 + luf->sv_tail = sv_prev[k]; 1.1242 + else 1.1243 + sv_prev[sv_next[k]] = sv_prev[k]; 1.1244 + /* the p-th column of the matrix F has been completely built; set 1.1245 + its pointer */ 1.1246 + fc_ptr[p] = luf->sv_end; 1.1247 + /* walk through the p-th (pivot) row and do the following... */ 1.1248 + for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++) 1.1249 + { /* get column index of v[p,j] */ 1.1250 + j = sv_ind[p_ptr]; 1.1251 + /* erase v[p,j] from the working array */ 1.1252 + flag[j] = 0; 1.1253 + work[j] = 0.0; 1.1254 + /* the j-th column has been completely transformed, therefore 1.1255 + it can return to the active set with new length; however 1.1256 + the special case c_prev[j] = c_next[j] = j means that the 1.1257 + routine find_pivot excluded the j-th column from the active 1.1258 + set due to Uwe Suhl's rule, and therefore in this case the 1.1259 + column can return to the active set only if it is a column 1.1260 + singleton */ 1.1261 + if (!(vc_len[j] != 1 && cs_prev[j] == j && cs_next[j] == j)) 1.1262 + { cs_prev[j] = 0; 1.1263 + cs_next[j] = cs_head[vc_len[j]]; 1.1264 + if (cs_next[j] != 0) cs_prev[cs_next[j]] = j; 1.1265 + cs_head[vc_len[j]] = j; 1.1266 + } 1.1267 + } 1.1268 +done: /* return to the factorizing routine */ 1.1269 + return ret; 1.1270 +} 1.1271 + 1.1272 +/*********************************************************************** 1.1273 +* build_v_cols - build the matrix V in column-wise format 1.1274 +* 1.1275 +* This routine builds the column-wise representation of the matrix V 1.1276 +* using its row-wise representation. 1.1277 +* 1.1278 +* If no error occured, the routine returns zero. Otherwise, in case of 1.1279 +* overflow of the sparse vector area, the routine returns non-zero. */ 1.1280 + 1.1281 +static int build_v_cols(LUF *luf) 1.1282 +{ int n = luf->n; 1.1283 + int *vr_ptr = luf->vr_ptr; 1.1284 + int *vr_len = luf->vr_len; 1.1285 + int *vc_ptr = luf->vc_ptr; 1.1286 + int *vc_len = luf->vc_len; 1.1287 + int *vc_cap = luf->vc_cap; 1.1288 + int *sv_ind = luf->sv_ind; 1.1289 + double *sv_val = luf->sv_val; 1.1290 + int *sv_prev = luf->sv_prev; 1.1291 + int *sv_next = luf->sv_next; 1.1292 + int ret = 0; 1.1293 + int i, i_beg, i_end, i_ptr, j, j_ptr, k, nnz; 1.1294 + /* it is assumed that on entry all columns of the matrix V are 1.1295 + empty, i.e. vc_len[j] = vc_cap[j] = 0 for all j = 1, ..., n, 1.1296 + and have been removed from the addressing list */ 1.1297 + /* count non-zeros in columns of the matrix V; count total number 1.1298 + of non-zeros in this matrix */ 1.1299 + nnz = 0; 1.1300 + for (i = 1; i <= n; i++) 1.1301 + { /* walk through elements of the i-th row and count non-zeros 1.1302 + in the corresponding columns */ 1.1303 + i_beg = vr_ptr[i]; 1.1304 + i_end = i_beg + vr_len[i] - 1; 1.1305 + for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) 1.1306 + vc_cap[sv_ind[i_ptr]]++; 1.1307 + /* count total number of non-zeros */ 1.1308 + nnz += vr_len[i]; 1.1309 + } 1.1310 + /* store total number of non-zeros */ 1.1311 + luf->nnz_v = nnz; 1.1312 + /* check for free locations */ 1.1313 + if (luf->sv_end - luf->sv_beg < nnz) 1.1314 + { /* overflow of the sparse vector area */ 1.1315 + ret = 1; 1.1316 + goto done; 1.1317 + } 1.1318 + /* allocate columns of the matrix V */ 1.1319 + for (j = 1; j <= n; j++) 1.1320 + { /* set pointer to the j-th column */ 1.1321 + vc_ptr[j] = luf->sv_beg; 1.1322 + /* reserve locations for the j-th column */ 1.1323 + luf->sv_beg += vc_cap[j]; 1.1324 + } 1.1325 + /* build the matrix V in column-wise format using this matrix in 1.1326 + row-wise format */ 1.1327 + for (i = 1; i <= n; i++) 1.1328 + { /* walk through elements of the i-th row */ 1.1329 + i_beg = vr_ptr[i]; 1.1330 + i_end = i_beg + vr_len[i] - 1; 1.1331 + for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) 1.1332 + { /* get column index */ 1.1333 + j = sv_ind[i_ptr]; 1.1334 + /* store element in the j-th column */ 1.1335 + j_ptr = vc_ptr[j] + vc_len[j]; 1.1336 + sv_ind[j_ptr] = i; 1.1337 + sv_val[j_ptr] = sv_val[i_ptr]; 1.1338 + /* increase length of the j-th column */ 1.1339 + vc_len[j]++; 1.1340 + } 1.1341 + } 1.1342 + /* now columns are placed in the sparse vector area behind rows 1.1343 + in the order n+1, n+2, ..., n+n; so insert column nodes in the 1.1344 + addressing list using this order */ 1.1345 + for (k = n+1; k <= n+n; k++) 1.1346 + { sv_prev[k] = k-1; 1.1347 + sv_next[k] = k+1; 1.1348 + } 1.1349 + sv_prev[n+1] = luf->sv_tail; 1.1350 + sv_next[luf->sv_tail] = n+1; 1.1351 + sv_next[n+n] = 0; 1.1352 + luf->sv_tail = n+n; 1.1353 +done: /* return to the factorizing routine */ 1.1354 + return ret; 1.1355 +} 1.1356 + 1.1357 +/*********************************************************************** 1.1358 +* build_f_rows - build the matrix F in row-wise format 1.1359 +* 1.1360 +* This routine builds the row-wise representation of the matrix F using 1.1361 +* its column-wise representation. 1.1362 +* 1.1363 +* If no error occured, the routine returns zero. Otherwise, in case of 1.1364 +* overflow of the sparse vector area, the routine returns non-zero. */ 1.1365 + 1.1366 +static int build_f_rows(LUF *luf) 1.1367 +{ int n = luf->n; 1.1368 + int *fr_ptr = luf->fr_ptr; 1.1369 + int *fr_len = luf->fr_len; 1.1370 + int *fc_ptr = luf->fc_ptr; 1.1371 + int *fc_len = luf->fc_len; 1.1372 + int *sv_ind = luf->sv_ind; 1.1373 + double *sv_val = luf->sv_val; 1.1374 + int ret = 0; 1.1375 + int i, j, j_beg, j_end, j_ptr, ptr, nnz; 1.1376 + /* clear rows of the matrix F */ 1.1377 + for (i = 1; i <= n; i++) fr_len[i] = 0; 1.1378 + /* count non-zeros in rows of the matrix F; count total number of 1.1379 + non-zeros in this matrix */ 1.1380 + nnz = 0; 1.1381 + for (j = 1; j <= n; j++) 1.1382 + { /* walk through elements of the j-th column and count non-zeros 1.1383 + in the corresponding rows */ 1.1384 + j_beg = fc_ptr[j]; 1.1385 + j_end = j_beg + fc_len[j] - 1; 1.1386 + for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) 1.1387 + fr_len[sv_ind[j_ptr]]++; 1.1388 + /* increase total number of non-zeros */ 1.1389 + nnz += fc_len[j]; 1.1390 + } 1.1391 + /* store total number of non-zeros */ 1.1392 + luf->nnz_f = nnz; 1.1393 + /* check for free locations */ 1.1394 + if (luf->sv_end - luf->sv_beg < nnz) 1.1395 + { /* overflow of the sparse vector area */ 1.1396 + ret = 1; 1.1397 + goto done; 1.1398 + } 1.1399 + /* allocate rows of the matrix F */ 1.1400 + for (i = 1; i <= n; i++) 1.1401 + { /* set pointer to the end of the i-th row; later this pointer 1.1402 + will be set to the beginning of the i-th row */ 1.1403 + fr_ptr[i] = luf->sv_end; 1.1404 + /* reserve locations for the i-th row */ 1.1405 + luf->sv_end -= fr_len[i]; 1.1406 + } 1.1407 + /* build the matrix F in row-wise format using this matrix in 1.1408 + column-wise format */ 1.1409 + for (j = 1; j <= n; j++) 1.1410 + { /* walk through elements of the j-th column */ 1.1411 + j_beg = fc_ptr[j]; 1.1412 + j_end = j_beg + fc_len[j] - 1; 1.1413 + for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) 1.1414 + { /* get row index */ 1.1415 + i = sv_ind[j_ptr]; 1.1416 + /* store element in the i-th row */ 1.1417 + ptr = --fr_ptr[i]; 1.1418 + sv_ind[ptr] = j; 1.1419 + sv_val[ptr] = sv_val[j_ptr]; 1.1420 + } 1.1421 + } 1.1422 +done: /* return to the factorizing routine */ 1.1423 + return ret; 1.1424 +} 1.1425 + 1.1426 +/*********************************************************************** 1.1427 +* NAME 1.1428 +* 1.1429 +* luf_factorize - compute LU-factorization 1.1430 +* 1.1431 +* SYNOPSIS 1.1432 +* 1.1433 +* #include "glpluf.h" 1.1434 +* int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j, 1.1435 +* int ind[], double val[]), void *info); 1.1436 +* 1.1437 +* DESCRIPTION 1.1438 +* 1.1439 +* The routine luf_factorize computes LU-factorization of a specified 1.1440 +* square matrix A. 1.1441 +* 1.1442 +* The parameter luf specifies LU-factorization program object created 1.1443 +* by the routine luf_create_it. 1.1444 +* 1.1445 +* The parameter n specifies the order of A, n > 0. 1.1446 +* 1.1447 +* The formal routine col specifies the matrix A to be factorized. To 1.1448 +* obtain j-th column of A the routine luf_factorize calls the routine 1.1449 +* col with the parameter j (1 <= j <= n). In response the routine col 1.1450 +* should store row indices and numerical values of non-zero elements 1.1451 +* of j-th column of A to locations ind[1,...,len] and val[1,...,len], 1.1452 +* respectively, where len is the number of non-zeros in j-th column 1.1453 +* returned on exit. Neither zero nor duplicate elements are allowed. 1.1454 +* 1.1455 +* The parameter info is a transit pointer passed to the routine col. 1.1456 +* 1.1457 +* RETURNS 1.1458 +* 1.1459 +* 0 LU-factorization has been successfully computed. 1.1460 +* 1.1461 +* LUF_ESING 1.1462 +* The specified matrix is singular within the working precision. 1.1463 +* (On some elimination step the active submatrix is exactly zero, 1.1464 +* so no pivot can be chosen.) 1.1465 +* 1.1466 +* LUF_ECOND 1.1467 +* The specified matrix is ill-conditioned. 1.1468 +* (On some elimination step too intensive growth of elements of the 1.1469 +* active submatix has been detected.) 1.1470 +* 1.1471 +* If matrix A is well scaled, the return code LUF_ECOND may also mean 1.1472 +* that the threshold pivoting tolerance piv_tol should be increased. 1.1473 +* 1.1474 +* In case of non-zero return code the factorization becomes invalid. 1.1475 +* It should not be used in other operations until the cause of failure 1.1476 +* has been eliminated and the factorization has been recomputed again 1.1477 +* with the routine luf_factorize. 1.1478 +* 1.1479 +* REPAIRING SINGULAR MATRIX 1.1480 +* 1.1481 +* If the routine luf_factorize returns non-zero code, it provides all 1.1482 +* necessary information that can be used for "repairing" the matrix A, 1.1483 +* where "repairing" means replacing linearly dependent columns of the 1.1484 +* matrix A by appropriate columns of the unity matrix. This feature is 1.1485 +* needed when this routine is used for factorizing the basis matrix 1.1486 +* within the simplex method procedure. 1.1487 +* 1.1488 +* On exit linearly dependent columns of the (partially transformed) 1.1489 +* matrix U have numbers rank+1, rank+2, ..., n, where rank is estimated 1.1490 +* rank of the matrix A stored by the routine to the member luf->rank. 1.1491 +* The correspondence between columns of A and U is the same as between 1.1492 +* columns of V and U. Thus, linearly dependent columns of the matrix A 1.1493 +* have numbers qq_col[rank+1], qq_col[rank+2], ..., qq_col[n], where 1.1494 +* qq_col is the column-like representation of the permutation matrix Q. 1.1495 +* It is understood that each j-th linearly dependent column of the 1.1496 +* matrix U should be replaced by the unity vector, where all elements 1.1497 +* are zero except the unity diagonal element u[j,j]. On the other hand 1.1498 +* j-th row of the matrix U corresponds to the row of the matrix V (and 1.1499 +* therefore of the matrix A) with the number pp_row[j], where pp_row is 1.1500 +* the row-like representation of the permutation matrix P. Thus, each 1.1501 +* j-th linearly dependent column of the matrix U should be replaced by 1.1502 +* column of the unity matrix with the number pp_row[j]. 1.1503 +* 1.1504 +* The code that repairs the matrix A may look like follows: 1.1505 +* 1.1506 +* for (j = rank+1; j <= n; j++) 1.1507 +* { replace the column qq_col[j] of the matrix A by the column 1.1508 +* pp_row[j] of the unity matrix; 1.1509 +* } 1.1510 +* 1.1511 +* where rank, pp_row, and qq_col are members of the structure LUF. */ 1.1512 + 1.1513 +int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j, 1.1514 + int ind[], double val[]), void *info) 1.1515 +{ int *pp_row, *pp_col, *qq_row, *qq_col; 1.1516 + double max_gro = luf->max_gro; 1.1517 + int i, j, k, p, q, t, ret; 1.1518 + if (n < 1) 1.1519 + xfault("luf_factorize: n = %d; invalid parameter\n", n); 1.1520 + if (n > N_MAX) 1.1521 + xfault("luf_factorize: n = %d; matrix too big\n", n); 1.1522 + /* invalidate the factorization */ 1.1523 + luf->valid = 0; 1.1524 + /* reallocate arrays, if necessary */ 1.1525 + reallocate(luf, n); 1.1526 + pp_row = luf->pp_row; 1.1527 + pp_col = luf->pp_col; 1.1528 + qq_row = luf->qq_row; 1.1529 + qq_col = luf->qq_col; 1.1530 + /* estimate initial size of the SVA, if not specified */ 1.1531 + if (luf->sv_size == 0 && luf->new_sva == 0) 1.1532 + luf->new_sva = 5 * (n + 10); 1.1533 +more: /* reallocate the sparse vector area, if required */ 1.1534 + if (luf->new_sva > 0) 1.1535 + { if (luf->sv_ind != NULL) xfree(luf->sv_ind); 1.1536 + if (luf->sv_val != NULL) xfree(luf->sv_val); 1.1537 + luf->sv_size = luf->new_sva; 1.1538 + luf->sv_ind = xcalloc(1+luf->sv_size, sizeof(int)); 1.1539 + luf->sv_val = xcalloc(1+luf->sv_size, sizeof(double)); 1.1540 + luf->new_sva = 0; 1.1541 + } 1.1542 + /* initialize LU-factorization data structures */ 1.1543 + if (initialize(luf, col, info)) 1.1544 + { /* overflow of the sparse vector area */ 1.1545 + luf->new_sva = luf->sv_size + luf->sv_size; 1.1546 + xassert(luf->new_sva > luf->sv_size); 1.1547 + goto more; 1.1548 + } 1.1549 + /* main elimination loop */ 1.1550 + for (k = 1; k <= n; k++) 1.1551 + { /* choose a pivot element v[p,q] */ 1.1552 + if (find_pivot(luf, &p, &q)) 1.1553 + { /* no pivot can be chosen, because the active submatrix is 1.1554 + exactly zero */ 1.1555 + luf->rank = k - 1; 1.1556 + ret = LUF_ESING; 1.1557 + goto done; 1.1558 + } 1.1559 + /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th 1.1560 + rows and k-th and j'-th columns of the matrix U = P*V*Q to 1.1561 + move the element u[i',j'] to the position u[k,k] */ 1.1562 + i = pp_col[p], j = qq_row[q]; 1.1563 + xassert(k <= i && i <= n && k <= j && j <= n); 1.1564 + /* permute k-th and i-th rows of the matrix U */ 1.1565 + t = pp_row[k]; 1.1566 + pp_row[i] = t, pp_col[t] = i; 1.1567 + pp_row[k] = p, pp_col[p] = k; 1.1568 + /* permute k-th and j-th columns of the matrix U */ 1.1569 + t = qq_col[k]; 1.1570 + qq_col[j] = t, qq_row[t] = j; 1.1571 + qq_col[k] = q, qq_row[q] = k; 1.1572 + /* eliminate subdiagonal elements of k-th column of the matrix 1.1573 + U = P*V*Q using the pivot element u[k,k] = v[p,q] */ 1.1574 + if (eliminate(luf, p, q)) 1.1575 + { /* overflow of the sparse vector area */ 1.1576 + luf->new_sva = luf->sv_size + luf->sv_size; 1.1577 + xassert(luf->new_sva > luf->sv_size); 1.1578 + goto more; 1.1579 + } 1.1580 + /* check relative growth of elements of the matrix V */ 1.1581 + if (luf->big_v > max_gro * luf->max_a) 1.1582 + { /* the growth is too intensive, therefore most probably the 1.1583 + matrix A is ill-conditioned */ 1.1584 + luf->rank = k - 1; 1.1585 + ret = LUF_ECOND; 1.1586 + goto done; 1.1587 + } 1.1588 + } 1.1589 + /* now the matrix U = P*V*Q is upper triangular, the matrix V has 1.1590 + been built in row-wise format, and the matrix F has been built 1.1591 + in column-wise format */ 1.1592 + /* defragment the sparse vector area in order to merge all free 1.1593 + locations in one continuous extent */ 1.1594 + luf_defrag_sva(luf); 1.1595 + /* build the matrix V in column-wise format */ 1.1596 + if (build_v_cols(luf)) 1.1597 + { /* overflow of the sparse vector area */ 1.1598 + luf->new_sva = luf->sv_size + luf->sv_size; 1.1599 + xassert(luf->new_sva > luf->sv_size); 1.1600 + goto more; 1.1601 + } 1.1602 + /* build the matrix F in row-wise format */ 1.1603 + if (build_f_rows(luf)) 1.1604 + { /* overflow of the sparse vector area */ 1.1605 + luf->new_sva = luf->sv_size + luf->sv_size; 1.1606 + xassert(luf->new_sva > luf->sv_size); 1.1607 + goto more; 1.1608 + } 1.1609 + /* the LU-factorization has been successfully computed */ 1.1610 + luf->valid = 1; 1.1611 + luf->rank = n; 1.1612 + ret = 0; 1.1613 + /* if there are few free locations in the sparse vector area, try 1.1614 + increasing its size in the future */ 1.1615 + t = 3 * (n + luf->nnz_v) + 2 * luf->nnz_f; 1.1616 + if (luf->sv_size < t) 1.1617 + { luf->new_sva = luf->sv_size; 1.1618 + while (luf->new_sva < t) 1.1619 + { k = luf->new_sva; 1.1620 + luf->new_sva = k + k; 1.1621 + xassert(luf->new_sva > k); 1.1622 + } 1.1623 + } 1.1624 +done: /* return to the calling program */ 1.1625 + return ret; 1.1626 +} 1.1627 + 1.1628 +/*********************************************************************** 1.1629 +* NAME 1.1630 +* 1.1631 +* luf_f_solve - solve system F*x = b or F'*x = b 1.1632 +* 1.1633 +* SYNOPSIS 1.1634 +* 1.1635 +* #include "glpluf.h" 1.1636 +* void luf_f_solve(LUF *luf, int tr, double x[]); 1.1637 +* 1.1638 +* DESCRIPTION 1.1639 +* 1.1640 +* The routine luf_f_solve solves either the system F*x = b (if the 1.1641 +* flag tr is zero) or the system F'*x = b (if the flag tr is non-zero), 1.1642 +* where the matrix F is a component of LU-factorization specified by 1.1643 +* the parameter luf, F' is a matrix transposed to F. 1.1644 +* 1.1645 +* On entry the array x should contain elements of the right-hand side 1.1646 +* vector b in locations x[1], ..., x[n], where n is the order of the 1.1647 +* matrix F. On exit this array will contain elements of the solution 1.1648 +* vector x in the same locations. */ 1.1649 + 1.1650 +void luf_f_solve(LUF *luf, int tr, double x[]) 1.1651 +{ int n = luf->n; 1.1652 + int *fr_ptr = luf->fr_ptr; 1.1653 + int *fr_len = luf->fr_len; 1.1654 + int *fc_ptr = luf->fc_ptr; 1.1655 + int *fc_len = luf->fc_len; 1.1656 + int *pp_row = luf->pp_row; 1.1657 + int *sv_ind = luf->sv_ind; 1.1658 + double *sv_val = luf->sv_val; 1.1659 + int i, j, k, beg, end, ptr; 1.1660 + double xk; 1.1661 + if (!luf->valid) 1.1662 + xfault("luf_f_solve: LU-factorization is not valid\n"); 1.1663 + if (!tr) 1.1664 + { /* solve the system F*x = b */ 1.1665 + for (j = 1; j <= n; j++) 1.1666 + { k = pp_row[j]; 1.1667 + xk = x[k]; 1.1668 + if (xk != 0.0) 1.1669 + { beg = fc_ptr[k]; 1.1670 + end = beg + fc_len[k] - 1; 1.1671 + for (ptr = beg; ptr <= end; ptr++) 1.1672 + x[sv_ind[ptr]] -= sv_val[ptr] * xk; 1.1673 + } 1.1674 + } 1.1675 + } 1.1676 + else 1.1677 + { /* solve the system F'*x = b */ 1.1678 + for (i = n; i >= 1; i--) 1.1679 + { k = pp_row[i]; 1.1680 + xk = x[k]; 1.1681 + if (xk != 0.0) 1.1682 + { beg = fr_ptr[k]; 1.1683 + end = beg + fr_len[k] - 1; 1.1684 + for (ptr = beg; ptr <= end; ptr++) 1.1685 + x[sv_ind[ptr]] -= sv_val[ptr] * xk; 1.1686 + } 1.1687 + } 1.1688 + } 1.1689 + return; 1.1690 +} 1.1691 + 1.1692 +/*********************************************************************** 1.1693 +* NAME 1.1694 +* 1.1695 +* luf_v_solve - solve system V*x = b or V'*x = b 1.1696 +* 1.1697 +* SYNOPSIS 1.1698 +* 1.1699 +* #include "glpluf.h" 1.1700 +* void luf_v_solve(LUF *luf, int tr, double x[]); 1.1701 +* 1.1702 +* DESCRIPTION 1.1703 +* 1.1704 +* The routine luf_v_solve solves either the system V*x = b (if the 1.1705 +* flag tr is zero) or the system V'*x = b (if the flag tr is non-zero), 1.1706 +* where the matrix V is a component of LU-factorization specified by 1.1707 +* the parameter luf, V' is a matrix transposed to V. 1.1708 +* 1.1709 +* On entry the array x should contain elements of the right-hand side 1.1710 +* vector b in locations x[1], ..., x[n], where n is the order of the 1.1711 +* matrix V. On exit this array will contain elements of the solution 1.1712 +* vector x in the same locations. */ 1.1713 + 1.1714 +void luf_v_solve(LUF *luf, int tr, double x[]) 1.1715 +{ int n = luf->n; 1.1716 + int *vr_ptr = luf->vr_ptr; 1.1717 + int *vr_len = luf->vr_len; 1.1718 + double *vr_piv = luf->vr_piv; 1.1719 + int *vc_ptr = luf->vc_ptr; 1.1720 + int *vc_len = luf->vc_len; 1.1721 + int *pp_row = luf->pp_row; 1.1722 + int *qq_col = luf->qq_col; 1.1723 + int *sv_ind = luf->sv_ind; 1.1724 + double *sv_val = luf->sv_val; 1.1725 + double *b = luf->work; 1.1726 + int i, j, k, beg, end, ptr; 1.1727 + double temp; 1.1728 + if (!luf->valid) 1.1729 + xfault("luf_v_solve: LU-factorization is not valid\n"); 1.1730 + for (k = 1; k <= n; k++) b[k] = x[k], x[k] = 0.0; 1.1731 + if (!tr) 1.1732 + { /* solve the system V*x = b */ 1.1733 + for (k = n; k >= 1; k--) 1.1734 + { i = pp_row[k], j = qq_col[k]; 1.1735 + temp = b[i]; 1.1736 + if (temp != 0.0) 1.1737 + { x[j] = (temp /= vr_piv[i]); 1.1738 + beg = vc_ptr[j]; 1.1739 + end = beg + vc_len[j] - 1; 1.1740 + for (ptr = beg; ptr <= end; ptr++) 1.1741 + b[sv_ind[ptr]] -= sv_val[ptr] * temp; 1.1742 + } 1.1743 + } 1.1744 + } 1.1745 + else 1.1746 + { /* solve the system V'*x = b */ 1.1747 + for (k = 1; k <= n; k++) 1.1748 + { i = pp_row[k], j = qq_col[k]; 1.1749 + temp = b[j]; 1.1750 + if (temp != 0.0) 1.1751 + { x[i] = (temp /= vr_piv[i]); 1.1752 + beg = vr_ptr[i]; 1.1753 + end = beg + vr_len[i] - 1; 1.1754 + for (ptr = beg; ptr <= end; ptr++) 1.1755 + b[sv_ind[ptr]] -= sv_val[ptr] * temp; 1.1756 + } 1.1757 + } 1.1758 + } 1.1759 + return; 1.1760 +} 1.1761 + 1.1762 +/*********************************************************************** 1.1763 +* NAME 1.1764 +* 1.1765 +* luf_a_solve - solve system A*x = b or A'*x = b 1.1766 +* 1.1767 +* SYNOPSIS 1.1768 +* 1.1769 +* #include "glpluf.h" 1.1770 +* void luf_a_solve(LUF *luf, int tr, double x[]); 1.1771 +* 1.1772 +* DESCRIPTION 1.1773 +* 1.1774 +* The routine luf_a_solve solves either the system A*x = b (if the 1.1775 +* flag tr is zero) or the system A'*x = b (if the flag tr is non-zero), 1.1776 +* where the parameter luf specifies LU-factorization of the matrix A, 1.1777 +* A' is a matrix transposed to A. 1.1778 +* 1.1779 +* On entry the array x should contain elements of the right-hand side 1.1780 +* vector b in locations x[1], ..., x[n], where n is the order of the 1.1781 +* matrix A. On exit this array will contain elements of the solution 1.1782 +* vector x in the same locations. */ 1.1783 + 1.1784 +void luf_a_solve(LUF *luf, int tr, double x[]) 1.1785 +{ if (!luf->valid) 1.1786 + xfault("luf_a_solve: LU-factorization is not valid\n"); 1.1787 + if (!tr) 1.1788 + { /* A = F*V, therefore inv(A) = inv(V)*inv(F) */ 1.1789 + luf_f_solve(luf, 0, x); 1.1790 + luf_v_solve(luf, 0, x); 1.1791 + } 1.1792 + else 1.1793 + { /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */ 1.1794 + luf_v_solve(luf, 1, x); 1.1795 + luf_f_solve(luf, 1, x); 1.1796 + } 1.1797 + return; 1.1798 +} 1.1799 + 1.1800 +/*********************************************************************** 1.1801 +* NAME 1.1802 +* 1.1803 +* luf_delete_it - delete LU-factorization 1.1804 +* 1.1805 +* SYNOPSIS 1.1806 +* 1.1807 +* #include "glpluf.h" 1.1808 +* void luf_delete_it(LUF *luf); 1.1809 +* 1.1810 +* DESCRIPTION 1.1811 +* 1.1812 +* The routine luf_delete deletes LU-factorization specified by the 1.1813 +* parameter luf and frees all the memory allocated to this program 1.1814 +* object. */ 1.1815 + 1.1816 +void luf_delete_it(LUF *luf) 1.1817 +{ if (luf->fr_ptr != NULL) xfree(luf->fr_ptr); 1.1818 + if (luf->fr_len != NULL) xfree(luf->fr_len); 1.1819 + if (luf->fc_ptr != NULL) xfree(luf->fc_ptr); 1.1820 + if (luf->fc_len != NULL) xfree(luf->fc_len); 1.1821 + if (luf->vr_ptr != NULL) xfree(luf->vr_ptr); 1.1822 + if (luf->vr_len != NULL) xfree(luf->vr_len); 1.1823 + if (luf->vr_cap != NULL) xfree(luf->vr_cap); 1.1824 + if (luf->vr_piv != NULL) xfree(luf->vr_piv); 1.1825 + if (luf->vc_ptr != NULL) xfree(luf->vc_ptr); 1.1826 + if (luf->vc_len != NULL) xfree(luf->vc_len); 1.1827 + if (luf->vc_cap != NULL) xfree(luf->vc_cap); 1.1828 + if (luf->pp_row != NULL) xfree(luf->pp_row); 1.1829 + if (luf->pp_col != NULL) xfree(luf->pp_col); 1.1830 + if (luf->qq_row != NULL) xfree(luf->qq_row); 1.1831 + if (luf->qq_col != NULL) xfree(luf->qq_col); 1.1832 + if (luf->sv_ind != NULL) xfree(luf->sv_ind); 1.1833 + if (luf->sv_val != NULL) xfree(luf->sv_val); 1.1834 + if (luf->sv_prev != NULL) xfree(luf->sv_prev); 1.1835 + if (luf->sv_next != NULL) xfree(luf->sv_next); 1.1836 + if (luf->vr_max != NULL) xfree(luf->vr_max); 1.1837 + if (luf->rs_head != NULL) xfree(luf->rs_head); 1.1838 + if (luf->rs_prev != NULL) xfree(luf->rs_prev); 1.1839 + if (luf->rs_next != NULL) xfree(luf->rs_next); 1.1840 + if (luf->cs_head != NULL) xfree(luf->cs_head); 1.1841 + if (luf->cs_prev != NULL) xfree(luf->cs_prev); 1.1842 + if (luf->cs_next != NULL) xfree(luf->cs_next); 1.1843 + if (luf->flag != NULL) xfree(luf->flag); 1.1844 + if (luf->work != NULL) xfree(luf->work); 1.1845 + xfree(luf); 1.1846 + return; 1.1847 +} 1.1848 + 1.1849 +/* eof */