alpar@1: /* glpluf.c (LU-factorization) */ alpar@1: alpar@1: /*********************************************************************** alpar@1: * This code is part of GLPK (GNU Linear Programming Kit). alpar@1: * alpar@1: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@1: * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, alpar@1: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@1: * E-mail: . alpar@1: * alpar@1: * GLPK is free software: you can redistribute it and/or modify it alpar@1: * under the terms of the GNU General Public License as published by alpar@1: * the Free Software Foundation, either version 3 of the License, or alpar@1: * (at your option) any later version. alpar@1: * alpar@1: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@1: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@1: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@1: * License for more details. alpar@1: * alpar@1: * You should have received a copy of the GNU General Public License alpar@1: * along with GLPK. If not, see . alpar@1: ***********************************************************************/ alpar@1: alpar@1: #include "glpenv.h" alpar@1: #include "glpluf.h" alpar@1: #define xfault xerror alpar@1: alpar@1: /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ alpar@1: alpar@1: #define N_MAX 100000000 /* = 100*10^6 */ alpar@1: /* maximal order of the original matrix */ alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * luf_create_it - create LU-factorization alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpluf.h" alpar@1: * LUF *luf_create_it(void); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine luf_create_it creates a program object, which represents alpar@1: * LU-factorization of a square matrix. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine luf_create_it returns a pointer to the object created. */ alpar@1: alpar@1: LUF *luf_create_it(void) alpar@1: { LUF *luf; alpar@1: luf = xmalloc(sizeof(LUF)); alpar@1: luf->n_max = luf->n = 0; alpar@1: luf->valid = 0; alpar@1: luf->fr_ptr = luf->fr_len = NULL; alpar@1: luf->fc_ptr = luf->fc_len = NULL; alpar@1: luf->vr_ptr = luf->vr_len = luf->vr_cap = NULL; alpar@1: luf->vr_piv = NULL; alpar@1: luf->vc_ptr = luf->vc_len = luf->vc_cap = NULL; alpar@1: luf->pp_row = luf->pp_col = NULL; alpar@1: luf->qq_row = luf->qq_col = NULL; alpar@1: luf->sv_size = 0; alpar@1: luf->sv_beg = luf->sv_end = 0; alpar@1: luf->sv_ind = NULL; alpar@1: luf->sv_val = NULL; alpar@1: luf->sv_head = luf->sv_tail = 0; alpar@1: luf->sv_prev = luf->sv_next = NULL; alpar@1: luf->vr_max = NULL; alpar@1: luf->rs_head = luf->rs_prev = luf->rs_next = NULL; alpar@1: luf->cs_head = luf->cs_prev = luf->cs_next = NULL; alpar@1: luf->flag = NULL; alpar@1: luf->work = NULL; alpar@1: luf->new_sva = 0; alpar@1: luf->piv_tol = 0.10; alpar@1: luf->piv_lim = 4; alpar@1: luf->suhl = 1; alpar@1: luf->eps_tol = 1e-15; alpar@1: luf->max_gro = 1e+10; alpar@1: luf->nnz_a = luf->nnz_f = luf->nnz_v = 0; alpar@1: luf->max_a = luf->big_v = 0.0; alpar@1: luf->rank = 0; alpar@1: return luf; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * luf_defrag_sva - defragment the sparse vector area alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpluf.h" alpar@1: * void luf_defrag_sva(LUF *luf); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine luf_defrag_sva defragments the sparse vector area (SVA) alpar@1: * gathering all unused locations in one continuous extent. In order to alpar@1: * do that the routine moves all unused locations from the left part of alpar@1: * SVA (which contains rows and columns of the matrix V) to the middle alpar@1: * part (which contains free locations). This is attained by relocating alpar@1: * elements of rows and columns of the matrix V toward the beginning of alpar@1: * the left part. alpar@1: * alpar@1: * NOTE that this "garbage collection" involves changing row and column alpar@1: * pointers of the matrix V. */ alpar@1: alpar@1: void luf_defrag_sva(LUF *luf) alpar@1: { int n = luf->n; alpar@1: int *vr_ptr = luf->vr_ptr; alpar@1: int *vr_len = luf->vr_len; alpar@1: int *vr_cap = luf->vr_cap; alpar@1: int *vc_ptr = luf->vc_ptr; alpar@1: int *vc_len = luf->vc_len; alpar@1: int *vc_cap = luf->vc_cap; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: int *sv_next = luf->sv_next; alpar@1: int sv_beg = 1; alpar@1: int i, j, k; alpar@1: /* skip rows and columns, which do not need to be relocated */ alpar@1: for (k = luf->sv_head; k != 0; k = sv_next[k]) alpar@1: { if (k <= n) alpar@1: { /* i-th row of the matrix V */ alpar@1: i = k; alpar@1: if (vr_ptr[i] != sv_beg) break; alpar@1: vr_cap[i] = vr_len[i]; alpar@1: sv_beg += vr_cap[i]; alpar@1: } alpar@1: else alpar@1: { /* j-th column of the matrix V */ alpar@1: j = k - n; alpar@1: if (vc_ptr[j] != sv_beg) break; alpar@1: vc_cap[j] = vc_len[j]; alpar@1: sv_beg += vc_cap[j]; alpar@1: } alpar@1: } alpar@1: /* relocate other rows and columns in order to gather all unused alpar@1: locations in one continuous extent */ alpar@1: for (k = k; k != 0; k = sv_next[k]) alpar@1: { if (k <= n) alpar@1: { /* i-th row of the matrix V */ alpar@1: i = k; alpar@1: memmove(&sv_ind[sv_beg], &sv_ind[vr_ptr[i]], alpar@1: vr_len[i] * sizeof(int)); alpar@1: memmove(&sv_val[sv_beg], &sv_val[vr_ptr[i]], alpar@1: vr_len[i] * sizeof(double)); alpar@1: vr_ptr[i] = sv_beg; alpar@1: vr_cap[i] = vr_len[i]; alpar@1: sv_beg += vr_cap[i]; alpar@1: } alpar@1: else alpar@1: { /* j-th column of the matrix V */ alpar@1: j = k - n; alpar@1: memmove(&sv_ind[sv_beg], &sv_ind[vc_ptr[j]], alpar@1: vc_len[j] * sizeof(int)); alpar@1: memmove(&sv_val[sv_beg], &sv_val[vc_ptr[j]], alpar@1: vc_len[j] * sizeof(double)); alpar@1: vc_ptr[j] = sv_beg; alpar@1: vc_cap[j] = vc_len[j]; alpar@1: sv_beg += vc_cap[j]; alpar@1: } alpar@1: } alpar@1: /* set new pointer to the beginning of the free part */ alpar@1: luf->sv_beg = sv_beg; alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * luf_enlarge_row - enlarge row capacity alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpluf.h" alpar@1: * int luf_enlarge_row(LUF *luf, int i, int cap); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine luf_enlarge_row enlarges capacity of the i-th row of the alpar@1: * matrix V to cap locations (assuming that its current capacity is less alpar@1: * than cap). In order to do that the routine relocates elements of the alpar@1: * i-th row to the end of the left part of SVA (which contains rows and alpar@1: * columns of the matrix V) and then expands the left part by allocating alpar@1: * cap free locations from the free part. If there are less than cap alpar@1: * free locations, the routine defragments the sparse vector area. alpar@1: * alpar@1: * Due to "garbage collection" this operation may change row and column alpar@1: * pointers of the matrix V. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * If no error occured, the routine returns zero. Otherwise, in case of alpar@1: * overflow of the sparse vector area, the routine returns non-zero. */ alpar@1: alpar@1: int luf_enlarge_row(LUF *luf, int i, int cap) alpar@1: { int n = luf->n; alpar@1: int *vr_ptr = luf->vr_ptr; alpar@1: int *vr_len = luf->vr_len; alpar@1: int *vr_cap = luf->vr_cap; alpar@1: int *vc_cap = luf->vc_cap; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: int *sv_prev = luf->sv_prev; alpar@1: int *sv_next = luf->sv_next; alpar@1: int ret = 0; alpar@1: int cur, k, kk; alpar@1: xassert(1 <= i && i <= n); alpar@1: xassert(vr_cap[i] < cap); alpar@1: /* if there are less than cap free locations, defragment SVA */ alpar@1: if (luf->sv_end - luf->sv_beg < cap) alpar@1: { luf_defrag_sva(luf); alpar@1: if (luf->sv_end - luf->sv_beg < cap) alpar@1: { ret = 1; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* save current capacity of the i-th row */ alpar@1: cur = vr_cap[i]; alpar@1: /* copy existing elements to the beginning of the free part */ alpar@1: memmove(&sv_ind[luf->sv_beg], &sv_ind[vr_ptr[i]], alpar@1: vr_len[i] * sizeof(int)); alpar@1: memmove(&sv_val[luf->sv_beg], &sv_val[vr_ptr[i]], alpar@1: vr_len[i] * sizeof(double)); alpar@1: /* set new pointer and new capacity of the i-th row */ alpar@1: vr_ptr[i] = luf->sv_beg; alpar@1: vr_cap[i] = cap; alpar@1: /* set new pointer to the beginning of the free part */ alpar@1: luf->sv_beg += cap; alpar@1: /* now the i-th row starts in the rightmost location among other alpar@1: rows and columns of the matrix V, so its node should be moved alpar@1: to the end of the row/column linked list */ alpar@1: k = i; alpar@1: /* remove the i-th row node from the linked list */ alpar@1: if (sv_prev[k] == 0) alpar@1: luf->sv_head = sv_next[k]; alpar@1: else alpar@1: { /* capacity of the previous row/column can be increased at the alpar@1: expense of old locations of the i-th row */ alpar@1: kk = sv_prev[k]; alpar@1: if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur; alpar@1: sv_next[sv_prev[k]] = sv_next[k]; alpar@1: } alpar@1: if (sv_next[k] == 0) alpar@1: luf->sv_tail = sv_prev[k]; alpar@1: else alpar@1: sv_prev[sv_next[k]] = sv_prev[k]; alpar@1: /* insert the i-th row node to the end of the linked list */ alpar@1: sv_prev[k] = luf->sv_tail; alpar@1: sv_next[k] = 0; alpar@1: if (sv_prev[k] == 0) alpar@1: luf->sv_head = k; alpar@1: else alpar@1: sv_next[sv_prev[k]] = k; alpar@1: luf->sv_tail = k; alpar@1: done: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * luf_enlarge_col - enlarge column capacity alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpluf.h" alpar@1: * int luf_enlarge_col(LUF *luf, int j, int cap); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine luf_enlarge_col enlarges capacity of the j-th column of alpar@1: * the matrix V to cap locations (assuming that its current capacity is alpar@1: * less than cap). In order to do that the routine relocates elements alpar@1: * of the j-th column to the end of the left part of SVA (which contains alpar@1: * rows and columns of the matrix V) and then expands the left part by alpar@1: * allocating cap free locations from the free part. If there are less alpar@1: * than cap free locations, the routine defragments the sparse vector alpar@1: * area. alpar@1: * alpar@1: * Due to "garbage collection" this operation may change row and column alpar@1: * pointers of the matrix V. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * If no error occured, the routine returns zero. Otherwise, in case of alpar@1: * overflow of the sparse vector area, the routine returns non-zero. */ alpar@1: alpar@1: int luf_enlarge_col(LUF *luf, int j, int cap) alpar@1: { int n = luf->n; alpar@1: int *vr_cap = luf->vr_cap; alpar@1: int *vc_ptr = luf->vc_ptr; alpar@1: int *vc_len = luf->vc_len; alpar@1: int *vc_cap = luf->vc_cap; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: int *sv_prev = luf->sv_prev; alpar@1: int *sv_next = luf->sv_next; alpar@1: int ret = 0; alpar@1: int cur, k, kk; alpar@1: xassert(1 <= j && j <= n); alpar@1: xassert(vc_cap[j] < cap); alpar@1: /* if there are less than cap free locations, defragment SVA */ alpar@1: if (luf->sv_end - luf->sv_beg < cap) alpar@1: { luf_defrag_sva(luf); alpar@1: if (luf->sv_end - luf->sv_beg < cap) alpar@1: { ret = 1; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* save current capacity of the j-th column */ alpar@1: cur = vc_cap[j]; alpar@1: /* copy existing elements to the beginning of the free part */ alpar@1: memmove(&sv_ind[luf->sv_beg], &sv_ind[vc_ptr[j]], alpar@1: vc_len[j] * sizeof(int)); alpar@1: memmove(&sv_val[luf->sv_beg], &sv_val[vc_ptr[j]], alpar@1: vc_len[j] * sizeof(double)); alpar@1: /* set new pointer and new capacity of the j-th column */ alpar@1: vc_ptr[j] = luf->sv_beg; alpar@1: vc_cap[j] = cap; alpar@1: /* set new pointer to the beginning of the free part */ alpar@1: luf->sv_beg += cap; alpar@1: /* now the j-th column starts in the rightmost location among alpar@1: other rows and columns of the matrix V, so its node should be alpar@1: moved to the end of the row/column linked list */ alpar@1: k = n + j; alpar@1: /* remove the j-th column node from the linked list */ alpar@1: if (sv_prev[k] == 0) alpar@1: luf->sv_head = sv_next[k]; alpar@1: else alpar@1: { /* capacity of the previous row/column can be increased at the alpar@1: expense of old locations of the j-th column */ alpar@1: kk = sv_prev[k]; alpar@1: if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur; alpar@1: sv_next[sv_prev[k]] = sv_next[k]; alpar@1: } alpar@1: if (sv_next[k] == 0) alpar@1: luf->sv_tail = sv_prev[k]; alpar@1: else alpar@1: sv_prev[sv_next[k]] = sv_prev[k]; alpar@1: /* insert the j-th column node to the end of the linked list */ alpar@1: sv_prev[k] = luf->sv_tail; alpar@1: sv_next[k] = 0; alpar@1: if (sv_prev[k] == 0) alpar@1: luf->sv_head = k; alpar@1: else alpar@1: sv_next[sv_prev[k]] = k; alpar@1: luf->sv_tail = k; alpar@1: done: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * reallocate - reallocate LU-factorization arrays alpar@1: * alpar@1: * This routine reallocates arrays, whose size depends of n, the order alpar@1: * of the matrix A to be factorized. */ alpar@1: alpar@1: static void reallocate(LUF *luf, int n) alpar@1: { int n_max = luf->n_max; alpar@1: luf->n = n; alpar@1: if (n <= n_max) goto done; alpar@1: if (luf->fr_ptr != NULL) xfree(luf->fr_ptr); alpar@1: if (luf->fr_len != NULL) xfree(luf->fr_len); alpar@1: if (luf->fc_ptr != NULL) xfree(luf->fc_ptr); alpar@1: if (luf->fc_len != NULL) xfree(luf->fc_len); alpar@1: if (luf->vr_ptr != NULL) xfree(luf->vr_ptr); alpar@1: if (luf->vr_len != NULL) xfree(luf->vr_len); alpar@1: if (luf->vr_cap != NULL) xfree(luf->vr_cap); alpar@1: if (luf->vr_piv != NULL) xfree(luf->vr_piv); alpar@1: if (luf->vc_ptr != NULL) xfree(luf->vc_ptr); alpar@1: if (luf->vc_len != NULL) xfree(luf->vc_len); alpar@1: if (luf->vc_cap != NULL) xfree(luf->vc_cap); alpar@1: if (luf->pp_row != NULL) xfree(luf->pp_row); alpar@1: if (luf->pp_col != NULL) xfree(luf->pp_col); alpar@1: if (luf->qq_row != NULL) xfree(luf->qq_row); alpar@1: if (luf->qq_col != NULL) xfree(luf->qq_col); alpar@1: if (luf->sv_prev != NULL) xfree(luf->sv_prev); alpar@1: if (luf->sv_next != NULL) xfree(luf->sv_next); alpar@1: if (luf->vr_max != NULL) xfree(luf->vr_max); alpar@1: if (luf->rs_head != NULL) xfree(luf->rs_head); alpar@1: if (luf->rs_prev != NULL) xfree(luf->rs_prev); alpar@1: if (luf->rs_next != NULL) xfree(luf->rs_next); alpar@1: if (luf->cs_head != NULL) xfree(luf->cs_head); alpar@1: if (luf->cs_prev != NULL) xfree(luf->cs_prev); alpar@1: if (luf->cs_next != NULL) xfree(luf->cs_next); alpar@1: if (luf->flag != NULL) xfree(luf->flag); alpar@1: if (luf->work != NULL) xfree(luf->work); alpar@1: luf->n_max = n_max = n + 100; alpar@1: luf->fr_ptr = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->fr_len = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->fc_ptr = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->fc_len = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->vr_ptr = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->vr_len = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->vr_cap = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->vr_piv = xcalloc(1+n_max, sizeof(double)); alpar@1: luf->vc_ptr = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->vc_len = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->vc_cap = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->pp_row = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->pp_col = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->qq_row = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->qq_col = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->sv_prev = xcalloc(1+n_max+n_max, sizeof(int)); alpar@1: luf->sv_next = xcalloc(1+n_max+n_max, sizeof(int)); alpar@1: luf->vr_max = xcalloc(1+n_max, sizeof(double)); alpar@1: luf->rs_head = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->rs_prev = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->rs_next = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->cs_head = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->cs_prev = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->cs_next = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->flag = xcalloc(1+n_max, sizeof(int)); alpar@1: luf->work = xcalloc(1+n_max, sizeof(double)); alpar@1: done: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * initialize - initialize LU-factorization data structures alpar@1: * alpar@1: * This routine initializes data structures for subsequent computing alpar@1: * the LU-factorization of a given matrix A, which is specified by the alpar@1: * formal routine col. On exit V = A and F = P = Q = I, where I is the alpar@1: * unity matrix. (Row-wise representation of the matrix F is not used alpar@1: * at the factorization stage and therefore is not initialized.) alpar@1: * alpar@1: * If no error occured, the routine returns zero. Otherwise, in case of alpar@1: * overflow of the sparse vector area, the routine returns non-zero. */ alpar@1: alpar@1: static int initialize(LUF *luf, int (*col)(void *info, int j, int rn[], alpar@1: double aj[]), void *info) alpar@1: { int n = luf->n; alpar@1: int *fc_ptr = luf->fc_ptr; alpar@1: int *fc_len = luf->fc_len; alpar@1: int *vr_ptr = luf->vr_ptr; alpar@1: int *vr_len = luf->vr_len; alpar@1: int *vr_cap = luf->vr_cap; alpar@1: int *vc_ptr = luf->vc_ptr; alpar@1: int *vc_len = luf->vc_len; alpar@1: int *vc_cap = luf->vc_cap; alpar@1: int *pp_row = luf->pp_row; alpar@1: int *pp_col = luf->pp_col; alpar@1: int *qq_row = luf->qq_row; alpar@1: int *qq_col = luf->qq_col; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: int *sv_prev = luf->sv_prev; alpar@1: int *sv_next = luf->sv_next; alpar@1: double *vr_max = luf->vr_max; alpar@1: int *rs_head = luf->rs_head; alpar@1: int *rs_prev = luf->rs_prev; alpar@1: int *rs_next = luf->rs_next; alpar@1: int *cs_head = luf->cs_head; alpar@1: int *cs_prev = luf->cs_prev; alpar@1: int *cs_next = luf->cs_next; alpar@1: int *flag = luf->flag; alpar@1: double *work = luf->work; alpar@1: int ret = 0; alpar@1: int i, i_ptr, j, j_beg, j_end, k, len, nnz, sv_beg, sv_end, ptr; alpar@1: double big, val; alpar@1: /* free all locations of the sparse vector area */ alpar@1: sv_beg = 1; alpar@1: sv_end = luf->sv_size + 1; alpar@1: /* (row-wise representation of the matrix F is not initialized, alpar@1: because it is not used at the factorization stage) */ alpar@1: /* build the matrix F in column-wise format (initially F = I) */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { fc_ptr[j] = sv_end; alpar@1: fc_len[j] = 0; alpar@1: } alpar@1: /* clear rows of the matrix V; clear the flag array */ alpar@1: for (i = 1; i <= n; i++) alpar@1: vr_len[i] = vr_cap[i] = 0, flag[i] = 0; alpar@1: /* build the matrix V in column-wise format (initially V = A); alpar@1: count non-zeros in rows of this matrix; count total number of alpar@1: non-zeros; compute largest of absolute values of elements */ alpar@1: nnz = 0; alpar@1: big = 0.0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { int *rn = pp_row; alpar@1: double *aj = work; alpar@1: /* obtain j-th column of the matrix A */ alpar@1: len = col(info, j, rn, aj); alpar@1: if (!(0 <= len && len <= n)) alpar@1: xfault("luf_factorize: j = %d; len = %d; invalid column len" alpar@1: "gth\n", j, len); alpar@1: /* check for free locations */ alpar@1: if (sv_end - sv_beg < len) alpar@1: { /* overflow of the sparse vector area */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* set pointer to the j-th column */ alpar@1: vc_ptr[j] = sv_beg; alpar@1: /* set length of the j-th column */ alpar@1: vc_len[j] = vc_cap[j] = len; alpar@1: /* count total number of non-zeros */ alpar@1: nnz += len; alpar@1: /* walk through elements of the j-th column */ alpar@1: for (ptr = 1; ptr <= len; ptr++) alpar@1: { /* get row index and numerical value of a[i,j] */ alpar@1: i = rn[ptr]; alpar@1: val = aj[ptr]; alpar@1: if (!(1 <= i && i <= n)) alpar@1: xfault("luf_factorize: i = %d; j = %d; invalid row index" alpar@1: "\n", i, j); alpar@1: if (flag[i]) alpar@1: xfault("luf_factorize: i = %d; j = %d; duplicate element" alpar@1: " not allowed\n", i, j); alpar@1: if (val == 0.0) alpar@1: xfault("luf_factorize: i = %d; j = %d; zero element not " alpar@1: "allowed\n", i, j); alpar@1: /* add new element v[i,j] = a[i,j] to j-th column */ alpar@1: sv_ind[sv_beg] = i; alpar@1: sv_val[sv_beg] = val; alpar@1: sv_beg++; alpar@1: /* big := max(big, |a[i,j]|) */ alpar@1: if (val < 0.0) val = - val; alpar@1: if (big < val) big = val; alpar@1: /* mark non-zero in the i-th position of the j-th column */ alpar@1: flag[i] = 1; alpar@1: /* increase length of the i-th row */ alpar@1: vr_cap[i]++; alpar@1: } alpar@1: /* reset all non-zero marks */ alpar@1: for (ptr = 1; ptr <= len; ptr++) flag[rn[ptr]] = 0; alpar@1: } alpar@1: /* allocate rows of the matrix V */ alpar@1: for (i = 1; i <= n; i++) alpar@1: { /* get length of the i-th row */ alpar@1: len = vr_cap[i]; alpar@1: /* check for free locations */ alpar@1: if (sv_end - sv_beg < len) alpar@1: { /* overflow of the sparse vector area */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* set pointer to the i-th row */ alpar@1: vr_ptr[i] = sv_beg; alpar@1: /* reserve locations for the i-th row */ alpar@1: sv_beg += len; alpar@1: } alpar@1: /* build the matrix V in row-wise format using representation of alpar@1: this matrix in column-wise format */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { /* walk through elements of the j-th column */ alpar@1: j_beg = vc_ptr[j]; alpar@1: j_end = j_beg + vc_len[j] - 1; alpar@1: for (k = j_beg; k <= j_end; k++) alpar@1: { /* get row index and numerical value of v[i,j] */ alpar@1: i = sv_ind[k]; alpar@1: val = sv_val[k]; alpar@1: /* store element in the i-th row */ alpar@1: i_ptr = vr_ptr[i] + vr_len[i]; alpar@1: sv_ind[i_ptr] = j; alpar@1: sv_val[i_ptr] = val; alpar@1: /* increase count of the i-th row */ alpar@1: vr_len[i]++; alpar@1: } alpar@1: } alpar@1: /* initialize the matrices P and Q (initially P = Q = I) */ alpar@1: for (k = 1; k <= n; k++) alpar@1: pp_row[k] = pp_col[k] = qq_row[k] = qq_col[k] = k; alpar@1: /* set sva partitioning pointers */ alpar@1: luf->sv_beg = sv_beg; alpar@1: luf->sv_end = sv_end; alpar@1: /* the initial physical order of rows and columns of the matrix V alpar@1: is n+1, ..., n+n, 1, ..., n (firstly columns, then rows) */ alpar@1: luf->sv_head = n+1; alpar@1: luf->sv_tail = n; alpar@1: for (i = 1; i <= n; i++) alpar@1: { sv_prev[i] = i-1; alpar@1: sv_next[i] = i+1; alpar@1: } alpar@1: sv_prev[1] = n+n; alpar@1: sv_next[n] = 0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { sv_prev[n+j] = n+j-1; alpar@1: sv_next[n+j] = n+j+1; alpar@1: } alpar@1: sv_prev[n+1] = 0; alpar@1: sv_next[n+n] = 1; alpar@1: /* clear working arrays */ alpar@1: for (k = 1; k <= n; k++) alpar@1: { flag[k] = 0; alpar@1: work[k] = 0.0; alpar@1: } alpar@1: /* initialize some statistics */ alpar@1: luf->nnz_a = nnz; alpar@1: luf->nnz_f = 0; alpar@1: luf->nnz_v = nnz; alpar@1: luf->max_a = big; alpar@1: luf->big_v = big; alpar@1: luf->rank = -1; alpar@1: /* initially the active submatrix is the entire matrix V */ alpar@1: /* largest of absolute values of elements in each active row is alpar@1: unknown yet */ alpar@1: for (i = 1; i <= n; i++) vr_max[i] = -1.0; alpar@1: /* build linked lists of active rows */ alpar@1: for (len = 0; len <= n; len++) rs_head[len] = 0; alpar@1: for (i = 1; i <= n; i++) alpar@1: { len = vr_len[i]; alpar@1: rs_prev[i] = 0; alpar@1: rs_next[i] = rs_head[len]; alpar@1: if (rs_next[i] != 0) rs_prev[rs_next[i]] = i; alpar@1: rs_head[len] = i; alpar@1: } alpar@1: /* build linked lists of active columns */ alpar@1: for (len = 0; len <= n; len++) cs_head[len] = 0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { len = vc_len[j]; alpar@1: cs_prev[j] = 0; alpar@1: cs_next[j] = cs_head[len]; alpar@1: if (cs_next[j] != 0) cs_prev[cs_next[j]] = j; alpar@1: cs_head[len] = j; alpar@1: } alpar@1: done: /* return to the factorizing routine */ alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * find_pivot - choose a pivot element alpar@1: * alpar@1: * This routine chooses a pivot element in the active submatrix of the alpar@1: * matrix U = P*V*Q. alpar@1: * alpar@1: * It is assumed that on entry the matrix U has the following partially alpar@1: * triangularized form: alpar@1: * alpar@1: * 1 k n alpar@1: * 1 x x x x x x x x x x alpar@1: * . x x x x x x x x x alpar@1: * . . x x x x x x x x alpar@1: * . . . x x x x x x x alpar@1: * k . . . . * * * * * * alpar@1: * . . . . * * * * * * alpar@1: * . . . . * * * * * * alpar@1: * . . . . * * * * * * alpar@1: * . . . . * * * * * * alpar@1: * n . . . . * * * * * * alpar@1: * alpar@1: * where rows and columns k, k+1, ..., n belong to the active submatrix alpar@1: * (elements of the active submatrix are marked by '*'). alpar@1: * alpar@1: * Since the matrix U = P*V*Q is not stored, the routine works with the alpar@1: * matrix V. It is assumed that the row-wise representation corresponds alpar@1: * to the matrix V, but the column-wise representation corresponds to alpar@1: * the active submatrix of the matrix V, i.e. elements of the matrix V, alpar@1: * which doesn't belong to the active submatrix, are missing from the alpar@1: * column linked lists. It is also assumed that each active row of the alpar@1: * matrix V is in the set R[len], where len is number of non-zeros in alpar@1: * the row, and each active column of the matrix V is in the set C[len], alpar@1: * where len is number of non-zeros in the column (in the latter case alpar@1: * only elements of the active submatrix are counted; such elements are alpar@1: * marked by '*' on the figure above). alpar@1: * alpar@1: * For the reason of numerical stability the routine applies so called alpar@1: * threshold pivoting proposed by J.Reid. It is assumed that an element alpar@1: * v[i,j] can be selected as a pivot candidate if it is not very small alpar@1: * (in absolute value) among other elements in the same row, i.e. if it alpar@1: * satisfies to the stability condition |v[i,j]| >= tol * max|v[i,*]|, alpar@1: * where 0 < tol < 1 is a given tolerance. alpar@1: * alpar@1: * In order to keep sparsity of the matrix V the routine uses Markowitz alpar@1: * strategy, trying to choose such element v[p,q], which satisfies to alpar@1: * the stability condition (see above) and has smallest Markowitz cost alpar@1: * (nr[p]-1) * (nc[q]-1), where nr[p] and nc[q] are numbers of non-zero alpar@1: * elements, respectively, in the p-th row and in the q-th column of the alpar@1: * active submatrix. alpar@1: * alpar@1: * In order to reduce the search, i.e. not to walk through all elements alpar@1: * of the active submatrix, the routine exploits a technique proposed by alpar@1: * I.Duff. This technique is based on using the sets R[len] and C[len] alpar@1: * of active rows and columns. alpar@1: * alpar@1: * If the pivot element v[p,q] has been chosen, the routine stores its alpar@1: * indices to the locations *p and *q and returns zero. Otherwise, if alpar@1: * the active submatrix is empty and therefore the pivot element can't alpar@1: * be chosen, the routine returns non-zero. */ alpar@1: alpar@1: static int find_pivot(LUF *luf, int *_p, int *_q) alpar@1: { int n = luf->n; alpar@1: int *vr_ptr = luf->vr_ptr; alpar@1: int *vr_len = luf->vr_len; alpar@1: int *vc_ptr = luf->vc_ptr; alpar@1: int *vc_len = luf->vc_len; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: double *vr_max = luf->vr_max; alpar@1: int *rs_head = luf->rs_head; alpar@1: int *rs_next = luf->rs_next; alpar@1: int *cs_head = luf->cs_head; alpar@1: int *cs_prev = luf->cs_prev; alpar@1: int *cs_next = luf->cs_next; alpar@1: double piv_tol = luf->piv_tol; alpar@1: int piv_lim = luf->piv_lim; alpar@1: int suhl = luf->suhl; alpar@1: int p, q, len, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr, alpar@1: ncand, next_j, min_p, min_q, min_len; alpar@1: double best, cost, big, temp; alpar@1: /* initially no pivot candidates have been found so far */ alpar@1: p = q = 0, best = DBL_MAX, ncand = 0; alpar@1: /* if in the active submatrix there is a column that has the only alpar@1: non-zero (column singleton), choose it as pivot */ alpar@1: j = cs_head[1]; alpar@1: if (j != 0) alpar@1: { xassert(vc_len[j] == 1); alpar@1: p = sv_ind[vc_ptr[j]], q = j; alpar@1: goto done; alpar@1: } alpar@1: /* if in the active submatrix there is a row that has the only alpar@1: non-zero (row singleton), choose it as pivot */ alpar@1: i = rs_head[1]; alpar@1: if (i != 0) alpar@1: { xassert(vr_len[i] == 1); alpar@1: p = i, q = sv_ind[vr_ptr[i]]; alpar@1: goto done; alpar@1: } alpar@1: /* there are no singletons in the active submatrix; walk through alpar@1: other non-empty rows and columns */ alpar@1: for (len = 2; len <= n; len++) alpar@1: { /* consider active columns that have len non-zeros */ alpar@1: for (j = cs_head[len]; j != 0; j = next_j) alpar@1: { /* the j-th column has len non-zeros */ alpar@1: j_beg = vc_ptr[j]; alpar@1: j_end = j_beg + vc_len[j] - 1; alpar@1: /* save pointer to the next column with the same length */ alpar@1: next_j = cs_next[j]; alpar@1: /* find an element in the j-th column, which is placed in a alpar@1: row with minimal number of non-zeros and satisfies to the alpar@1: stability condition (such element may not exist) */ alpar@1: min_p = min_q = 0, min_len = INT_MAX; alpar@1: for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) alpar@1: { /* get row index of v[i,j] */ alpar@1: i = sv_ind[j_ptr]; alpar@1: i_beg = vr_ptr[i]; alpar@1: i_end = i_beg + vr_len[i] - 1; alpar@1: /* if the i-th row is not shorter than that one, where alpar@1: minimal element is currently placed, skip v[i,j] */ alpar@1: if (vr_len[i] >= min_len) continue; alpar@1: /* determine the largest of absolute values of elements alpar@1: in the i-th row */ alpar@1: big = vr_max[i]; alpar@1: if (big < 0.0) alpar@1: { /* the largest value is unknown yet; compute it */ alpar@1: for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) alpar@1: { temp = sv_val[i_ptr]; alpar@1: if (temp < 0.0) temp = - temp; alpar@1: if (big < temp) big = temp; alpar@1: } alpar@1: vr_max[i] = big; alpar@1: } alpar@1: /* find v[i,j] in the i-th row */ alpar@1: for (i_ptr = vr_ptr[i]; sv_ind[i_ptr] != j; i_ptr++); alpar@1: xassert(i_ptr <= i_end); alpar@1: /* if v[i,j] doesn't satisfy to the stability condition, alpar@1: skip it */ alpar@1: temp = sv_val[i_ptr]; alpar@1: if (temp < 0.0) temp = - temp; alpar@1: if (temp < piv_tol * big) continue; alpar@1: /* v[i,j] is better than the current minimal element */ alpar@1: min_p = i, min_q = j, min_len = vr_len[i]; alpar@1: /* if Markowitz cost of the current minimal element is alpar@1: not greater than (len-1)**2, it can be chosen right alpar@1: now; this heuristic reduces the search and works well alpar@1: in many cases */ alpar@1: if (min_len <= len) alpar@1: { p = min_p, q = min_q; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* the j-th column has been scanned */ alpar@1: if (min_p != 0) alpar@1: { /* the minimal element is a next pivot candidate */ alpar@1: ncand++; alpar@1: /* compute its Markowitz cost */ alpar@1: cost = (double)(min_len - 1) * (double)(len - 1); alpar@1: /* choose between the minimal element and the current alpar@1: candidate */ alpar@1: if (cost < best) p = min_p, q = min_q, best = cost; alpar@1: /* if piv_lim candidates have been considered, there are alpar@1: doubts that a much better candidate exists; therefore alpar@1: it's time to terminate the search */ alpar@1: if (ncand == piv_lim) goto done; alpar@1: } alpar@1: else alpar@1: { /* the j-th column has no elements, which satisfy to the alpar@1: stability condition; Uwe Suhl suggests to exclude such alpar@1: column from the further consideration until it becomes alpar@1: a column singleton; in hard cases this significantly alpar@1: reduces a time needed for pivot searching */ alpar@1: if (suhl) alpar@1: { /* remove the j-th column from the active set */ alpar@1: if (cs_prev[j] == 0) alpar@1: cs_head[len] = cs_next[j]; alpar@1: else alpar@1: cs_next[cs_prev[j]] = cs_next[j]; alpar@1: if (cs_next[j] == 0) alpar@1: /* nop */; alpar@1: else alpar@1: cs_prev[cs_next[j]] = cs_prev[j]; alpar@1: /* the following assignment is used to avoid an error alpar@1: when the routine eliminate (see below) will try to alpar@1: remove the j-th column from the active set */ alpar@1: cs_prev[j] = cs_next[j] = j; alpar@1: } alpar@1: } alpar@1: } alpar@1: /* consider active rows that have len non-zeros */ alpar@1: for (i = rs_head[len]; i != 0; i = rs_next[i]) alpar@1: { /* the i-th row has len non-zeros */ alpar@1: i_beg = vr_ptr[i]; alpar@1: i_end = i_beg + vr_len[i] - 1; alpar@1: /* determine the largest of absolute values of elements in alpar@1: the i-th row */ alpar@1: big = vr_max[i]; alpar@1: if (big < 0.0) alpar@1: { /* the largest value is unknown yet; compute it */ alpar@1: for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) alpar@1: { temp = sv_val[i_ptr]; alpar@1: if (temp < 0.0) temp = - temp; alpar@1: if (big < temp) big = temp; alpar@1: } alpar@1: vr_max[i] = big; alpar@1: } alpar@1: /* find an element in the i-th row, which is placed in a alpar@1: column with minimal number of non-zeros and satisfies to alpar@1: the stability condition (such element always exists) */ alpar@1: min_p = min_q = 0, min_len = INT_MAX; alpar@1: for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) alpar@1: { /* get column index of v[i,j] */ alpar@1: j = sv_ind[i_ptr]; alpar@1: /* if the j-th column is not shorter than that one, where alpar@1: minimal element is currently placed, skip v[i,j] */ alpar@1: if (vc_len[j] >= min_len) continue; alpar@1: /* if v[i,j] doesn't satisfy to the stability condition, alpar@1: skip it */ alpar@1: temp = sv_val[i_ptr]; alpar@1: if (temp < 0.0) temp = - temp; alpar@1: if (temp < piv_tol * big) continue; alpar@1: /* v[i,j] is better than the current minimal element */ alpar@1: min_p = i, min_q = j, min_len = vc_len[j]; alpar@1: /* if Markowitz cost of the current minimal element is alpar@1: not greater than (len-1)**2, it can be chosen right alpar@1: now; this heuristic reduces the search and works well alpar@1: in many cases */ alpar@1: if (min_len <= len) alpar@1: { p = min_p, q = min_q; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* the i-th row has been scanned */ alpar@1: if (min_p != 0) alpar@1: { /* the minimal element is a next pivot candidate */ alpar@1: ncand++; alpar@1: /* compute its Markowitz cost */ alpar@1: cost = (double)(len - 1) * (double)(min_len - 1); alpar@1: /* choose between the minimal element and the current alpar@1: candidate */ alpar@1: if (cost < best) p = min_p, q = min_q, best = cost; alpar@1: /* if piv_lim candidates have been considered, there are alpar@1: doubts that a much better candidate exists; therefore alpar@1: it's time to terminate the search */ alpar@1: if (ncand == piv_lim) goto done; alpar@1: } alpar@1: else alpar@1: { /* this can't be because this can never be */ alpar@1: xassert(min_p != min_p); alpar@1: } alpar@1: } alpar@1: } alpar@1: done: /* bring the pivot to the factorizing routine */ alpar@1: *_p = p, *_q = q; alpar@1: return (p == 0); alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * eliminate - perform gaussian elimination. alpar@1: * alpar@1: * This routine performs elementary gaussian transformations in order alpar@1: * to eliminate subdiagonal elements in the k-th column of the matrix alpar@1: * U = P*V*Q using the pivot element u[k,k], where k is the number of alpar@1: * the current elimination step. alpar@1: * alpar@1: * The parameters p and q are, respectively, row and column indices of alpar@1: * the element v[p,q], which corresponds to the element u[k,k]. alpar@1: * alpar@1: * Each time when the routine applies the elementary transformation to alpar@1: * a non-pivot row of the matrix V, it stores the corresponding element alpar@1: * to the matrix F in order to keep the main equality A = F*V. alpar@1: * alpar@1: * The routine assumes that on entry the matrices L = P*F*inv(P) and alpar@1: * U = P*V*Q are the following: alpar@1: * alpar@1: * 1 k 1 k n alpar@1: * 1 1 . . . . . . . . . 1 x x x x x x x x x x alpar@1: * x 1 . . . . . . . . . x x x x x x x x x alpar@1: * x x 1 . . . . . . . . . x x x x x x x x alpar@1: * x x x 1 . . . . . . . . . x x x x x x x alpar@1: * k x x x x 1 . . . . . k . . . . * * * * * * alpar@1: * x x x x _ 1 . . . . . . . . # * * * * * alpar@1: * x x x x _ . 1 . . . . . . . # * * * * * alpar@1: * x x x x _ . . 1 . . . . . . # * * * * * alpar@1: * x x x x _ . . . 1 . . . . . # * * * * * alpar@1: * n x x x x _ . . . . 1 n . . . . # * * * * * alpar@1: * alpar@1: * matrix L matrix U alpar@1: * alpar@1: * where rows and columns of the matrix U with numbers k, k+1, ..., n alpar@1: * form the active submatrix (eliminated elements are marked by '#' and alpar@1: * other elements of the active submatrix are marked by '*'). Note that alpar@1: * each eliminated non-zero element u[i,k] of the matrix U gives the alpar@1: * corresponding element l[i,k] of the matrix L (marked by '_'). alpar@1: * alpar@1: * Actually all operations are performed on the matrix V. Should note alpar@1: * that the row-wise representation corresponds to the matrix V, but the alpar@1: * column-wise representation corresponds to the active submatrix of the alpar@1: * matrix V, i.e. elements of the matrix V, which doesn't belong to the alpar@1: * active submatrix, are missing from the column linked lists. alpar@1: * alpar@1: * Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal alpar@1: * elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies alpar@1: * the following elementary gaussian transformations: alpar@1: * alpar@1: * (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V), alpar@1: * alpar@1: * where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier. alpar@1: * alpar@1: * Additionally, in order to keep the main equality A = F*V, each time alpar@1: * when the routine applies the transformation to i-th row of the matrix alpar@1: * V, it also adds f[i,p] as a new element to the matrix F. alpar@1: * alpar@1: * IMPORTANT: On entry the working arrays flag and work should contain alpar@1: * zeros. This status is provided by the routine on exit. alpar@1: * alpar@1: * If no error occured, the routine returns zero. Otherwise, in case of alpar@1: * overflow of the sparse vector area, the routine returns non-zero. */ alpar@1: alpar@1: static int eliminate(LUF *luf, int p, int q) alpar@1: { int n = luf->n; alpar@1: int *fc_ptr = luf->fc_ptr; alpar@1: int *fc_len = luf->fc_len; alpar@1: int *vr_ptr = luf->vr_ptr; alpar@1: int *vr_len = luf->vr_len; alpar@1: int *vr_cap = luf->vr_cap; alpar@1: double *vr_piv = luf->vr_piv; alpar@1: int *vc_ptr = luf->vc_ptr; alpar@1: int *vc_len = luf->vc_len; alpar@1: int *vc_cap = luf->vc_cap; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: int *sv_prev = luf->sv_prev; alpar@1: int *sv_next = luf->sv_next; alpar@1: double *vr_max = luf->vr_max; alpar@1: int *rs_head = luf->rs_head; alpar@1: int *rs_prev = luf->rs_prev; alpar@1: int *rs_next = luf->rs_next; alpar@1: int *cs_head = luf->cs_head; alpar@1: int *cs_prev = luf->cs_prev; alpar@1: int *cs_next = luf->cs_next; alpar@1: int *flag = luf->flag; alpar@1: double *work = luf->work; alpar@1: double eps_tol = luf->eps_tol; alpar@1: /* at this stage the row-wise representation of the matrix F is alpar@1: not used, so fr_len can be used as a working array */ alpar@1: int *ndx = luf->fr_len; alpar@1: int ret = 0; alpar@1: int len, fill, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr, k, alpar@1: p_beg, p_end, p_ptr, q_beg, q_end, q_ptr; alpar@1: double fip, val, vpq, temp; alpar@1: xassert(1 <= p && p <= n); alpar@1: xassert(1 <= q && q <= n); alpar@1: /* remove the p-th (pivot) row from the active set; this row will alpar@1: never return there */ alpar@1: if (rs_prev[p] == 0) alpar@1: rs_head[vr_len[p]] = rs_next[p]; alpar@1: else alpar@1: rs_next[rs_prev[p]] = rs_next[p]; alpar@1: if (rs_next[p] == 0) alpar@1: ; alpar@1: else alpar@1: rs_prev[rs_next[p]] = rs_prev[p]; alpar@1: /* remove the q-th (pivot) column from the active set; this column alpar@1: will never return there */ alpar@1: if (cs_prev[q] == 0) alpar@1: cs_head[vc_len[q]] = cs_next[q]; alpar@1: else alpar@1: cs_next[cs_prev[q]] = cs_next[q]; alpar@1: if (cs_next[q] == 0) alpar@1: ; alpar@1: else alpar@1: cs_prev[cs_next[q]] = cs_prev[q]; alpar@1: /* find the pivot v[p,q] = u[k,k] in the p-th row */ alpar@1: p_beg = vr_ptr[p]; alpar@1: p_end = p_beg + vr_len[p] - 1; alpar@1: for (p_ptr = p_beg; sv_ind[p_ptr] != q; p_ptr++) /* nop */; alpar@1: xassert(p_ptr <= p_end); alpar@1: /* store value of the pivot */ alpar@1: vpq = (vr_piv[p] = sv_val[p_ptr]); alpar@1: /* remove the pivot from the p-th row */ alpar@1: sv_ind[p_ptr] = sv_ind[p_end]; alpar@1: sv_val[p_ptr] = sv_val[p_end]; alpar@1: vr_len[p]--; alpar@1: p_end--; alpar@1: /* find the pivot v[p,q] = u[k,k] in the q-th column */ alpar@1: q_beg = vc_ptr[q]; alpar@1: q_end = q_beg + vc_len[q] - 1; alpar@1: for (q_ptr = q_beg; sv_ind[q_ptr] != p; q_ptr++) /* nop */; alpar@1: xassert(q_ptr <= q_end); alpar@1: /* remove the pivot from the q-th column */ alpar@1: sv_ind[q_ptr] = sv_ind[q_end]; alpar@1: vc_len[q]--; alpar@1: q_end--; alpar@1: /* walk through the p-th (pivot) row, which doesn't contain the alpar@1: pivot v[p,q] already, and do the following... */ alpar@1: for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++) alpar@1: { /* get column index of v[p,j] */ alpar@1: j = sv_ind[p_ptr]; alpar@1: /* store v[p,j] to the working array */ alpar@1: flag[j] = 1; alpar@1: work[j] = sv_val[p_ptr]; alpar@1: /* remove the j-th column from the active set; this column will alpar@1: return there later with new length */ alpar@1: if (cs_prev[j] == 0) alpar@1: cs_head[vc_len[j]] = cs_next[j]; alpar@1: else alpar@1: cs_next[cs_prev[j]] = cs_next[j]; alpar@1: if (cs_next[j] == 0) alpar@1: ; alpar@1: else alpar@1: cs_prev[cs_next[j]] = cs_prev[j]; alpar@1: /* find v[p,j] in the j-th column */ alpar@1: j_beg = vc_ptr[j]; alpar@1: j_end = j_beg + vc_len[j] - 1; alpar@1: for (j_ptr = j_beg; sv_ind[j_ptr] != p; j_ptr++) /* nop */; alpar@1: xassert(j_ptr <= j_end); alpar@1: /* since v[p,j] leaves the active submatrix, remove it from the alpar@1: j-th column; however, v[p,j] is kept in the p-th row */ alpar@1: sv_ind[j_ptr] = sv_ind[j_end]; alpar@1: vc_len[j]--; alpar@1: } alpar@1: /* walk through the q-th (pivot) column, which doesn't contain the alpar@1: pivot v[p,q] already, and perform gaussian elimination */ alpar@1: while (q_beg <= q_end) alpar@1: { /* element v[i,q] should be eliminated */ alpar@1: /* get row index of v[i,q] */ alpar@1: i = sv_ind[q_beg]; alpar@1: /* remove the i-th row from the active set; later this row will alpar@1: return there with new length */ alpar@1: if (rs_prev[i] == 0) alpar@1: rs_head[vr_len[i]] = rs_next[i]; alpar@1: else alpar@1: rs_next[rs_prev[i]] = rs_next[i]; alpar@1: if (rs_next[i] == 0) alpar@1: ; alpar@1: else alpar@1: rs_prev[rs_next[i]] = rs_prev[i]; alpar@1: /* find v[i,q] in the i-th row */ alpar@1: i_beg = vr_ptr[i]; alpar@1: i_end = i_beg + vr_len[i] - 1; alpar@1: for (i_ptr = i_beg; sv_ind[i_ptr] != q; i_ptr++) /* nop */; alpar@1: xassert(i_ptr <= i_end); alpar@1: /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] */ alpar@1: fip = sv_val[i_ptr] / vpq; alpar@1: /* since v[i,q] should be eliminated, remove it from the i-th alpar@1: row */ alpar@1: sv_ind[i_ptr] = sv_ind[i_end]; alpar@1: sv_val[i_ptr] = sv_val[i_end]; alpar@1: vr_len[i]--; alpar@1: i_end--; alpar@1: /* and from the q-th column */ alpar@1: sv_ind[q_beg] = sv_ind[q_end]; alpar@1: vc_len[q]--; alpar@1: q_end--; alpar@1: /* perform gaussian transformation: alpar@1: (i-th row) := (i-th row) - f[i,p] * (p-th row) alpar@1: note that now the p-th row, which is in the working array, alpar@1: doesn't contain the pivot v[p,q], and the i-th row doesn't alpar@1: contain the eliminated element v[i,q] */ alpar@1: /* walk through the i-th row and transform existing non-zero alpar@1: elements */ alpar@1: fill = vr_len[p]; alpar@1: for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) alpar@1: { /* get column index of v[i,j] */ alpar@1: j = sv_ind[i_ptr]; alpar@1: /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */ alpar@1: if (flag[j]) alpar@1: { /* v[p,j] != 0 */ alpar@1: temp = (sv_val[i_ptr] -= fip * work[j]); alpar@1: if (temp < 0.0) temp = - temp; alpar@1: flag[j] = 0; alpar@1: fill--; /* since both v[i,j] and v[p,j] exist */ alpar@1: if (temp == 0.0 || temp < eps_tol) alpar@1: { /* new v[i,j] is closer to zero; replace it by exact alpar@1: zero, i.e. remove it from the active submatrix */ alpar@1: /* remove v[i,j] from the i-th row */ alpar@1: sv_ind[i_ptr] = sv_ind[i_end]; alpar@1: sv_val[i_ptr] = sv_val[i_end]; alpar@1: vr_len[i]--; alpar@1: i_ptr--; alpar@1: i_end--; alpar@1: /* find v[i,j] in the j-th column */ alpar@1: j_beg = vc_ptr[j]; alpar@1: j_end = j_beg + vc_len[j] - 1; alpar@1: for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++); alpar@1: xassert(j_ptr <= j_end); alpar@1: /* remove v[i,j] from the j-th column */ alpar@1: sv_ind[j_ptr] = sv_ind[j_end]; alpar@1: vc_len[j]--; alpar@1: } alpar@1: else alpar@1: { /* v_big := max(v_big, |v[i,j]|) */ alpar@1: if (luf->big_v < temp) luf->big_v = temp; alpar@1: } alpar@1: } alpar@1: } alpar@1: /* now flag is the pattern of the set v[p,*] \ v[i,*], and fill alpar@1: is number of non-zeros in this set; therefore up to fill new alpar@1: non-zeros may appear in the i-th row */ alpar@1: if (vr_len[i] + fill > vr_cap[i]) alpar@1: { /* enlarge the i-th row */ alpar@1: if (luf_enlarge_row(luf, i, vr_len[i] + fill)) alpar@1: { /* overflow of the sparse vector area */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* defragmentation may change row and column pointers of the alpar@1: matrix V */ alpar@1: p_beg = vr_ptr[p]; alpar@1: p_end = p_beg + vr_len[p] - 1; alpar@1: q_beg = vc_ptr[q]; alpar@1: q_end = q_beg + vc_len[q] - 1; alpar@1: } alpar@1: /* walk through the p-th (pivot) row and create new elements alpar@1: of the i-th row that appear due to fill-in; column indices alpar@1: of these new elements are accumulated in the array ndx */ alpar@1: len = 0; alpar@1: for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++) alpar@1: { /* get column index of v[p,j], which may cause fill-in */ alpar@1: j = sv_ind[p_ptr]; alpar@1: if (flag[j]) alpar@1: { /* compute new non-zero v[i,j] = 0 - f[i,p] * v[p,j] */ alpar@1: temp = (val = - fip * work[j]); alpar@1: if (temp < 0.0) temp = - temp; alpar@1: if (temp == 0.0 || temp < eps_tol) alpar@1: /* if v[i,j] is closer to zero; just ignore it */; alpar@1: else alpar@1: { /* add v[i,j] to the i-th row */ alpar@1: i_ptr = vr_ptr[i] + vr_len[i]; alpar@1: sv_ind[i_ptr] = j; alpar@1: sv_val[i_ptr] = val; alpar@1: vr_len[i]++; alpar@1: /* remember column index of v[i,j] */ alpar@1: ndx[++len] = j; alpar@1: /* big_v := max(big_v, |v[i,j]|) */ alpar@1: if (luf->big_v < temp) luf->big_v = temp; alpar@1: } alpar@1: } alpar@1: else alpar@1: { /* there is no fill-in, because v[i,j] already exists in alpar@1: the i-th row; restore the flag of the element v[p,j], alpar@1: which was reset before */ alpar@1: flag[j] = 1; alpar@1: } alpar@1: } alpar@1: /* add new non-zeros v[i,j] to the corresponding columns */ alpar@1: for (k = 1; k <= len; k++) alpar@1: { /* get column index of new non-zero v[i,j] */ alpar@1: j = ndx[k]; alpar@1: /* one free location is needed in the j-th column */ alpar@1: if (vc_len[j] + 1 > vc_cap[j]) alpar@1: { /* enlarge the j-th column */ alpar@1: if (luf_enlarge_col(luf, j, vc_len[j] + 10)) alpar@1: { /* overflow of the sparse vector area */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* defragmentation may change row and column pointers of alpar@1: the matrix V */ alpar@1: p_beg = vr_ptr[p]; alpar@1: p_end = p_beg + vr_len[p] - 1; alpar@1: q_beg = vc_ptr[q]; alpar@1: q_end = q_beg + vc_len[q] - 1; alpar@1: } alpar@1: /* add new non-zero v[i,j] to the j-th column */ alpar@1: j_ptr = vc_ptr[j] + vc_len[j]; alpar@1: sv_ind[j_ptr] = i; alpar@1: vc_len[j]++; alpar@1: } alpar@1: /* now the i-th row has been completely transformed, therefore alpar@1: it can return to the active set with new length */ alpar@1: rs_prev[i] = 0; alpar@1: rs_next[i] = rs_head[vr_len[i]]; alpar@1: if (rs_next[i] != 0) rs_prev[rs_next[i]] = i; alpar@1: rs_head[vr_len[i]] = i; alpar@1: /* the largest of absolute values of elements in the i-th row alpar@1: is currently unknown */ alpar@1: vr_max[i] = -1.0; alpar@1: /* at least one free location is needed to store the gaussian alpar@1: multiplier */ alpar@1: if (luf->sv_end - luf->sv_beg < 1) alpar@1: { /* there are no free locations at all; defragment SVA */ alpar@1: luf_defrag_sva(luf); alpar@1: if (luf->sv_end - luf->sv_beg < 1) alpar@1: { /* overflow of the sparse vector area */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* defragmentation may change row and column pointers of the alpar@1: matrix V */ alpar@1: p_beg = vr_ptr[p]; alpar@1: p_end = p_beg + vr_len[p] - 1; alpar@1: q_beg = vc_ptr[q]; alpar@1: q_end = q_beg + vc_len[q] - 1; alpar@1: } alpar@1: /* add the element f[i,p], which is the gaussian multiplier, alpar@1: to the matrix F */ alpar@1: luf->sv_end--; alpar@1: sv_ind[luf->sv_end] = i; alpar@1: sv_val[luf->sv_end] = fip; alpar@1: fc_len[p]++; alpar@1: /* end of elimination loop */ alpar@1: } alpar@1: /* at this point the q-th (pivot) column should be empty */ alpar@1: xassert(vc_len[q] == 0); alpar@1: /* reset capacity of the q-th column */ alpar@1: vc_cap[q] = 0; alpar@1: /* remove node of the q-th column from the addressing list */ alpar@1: k = n + q; alpar@1: if (sv_prev[k] == 0) alpar@1: luf->sv_head = sv_next[k]; alpar@1: else alpar@1: sv_next[sv_prev[k]] = sv_next[k]; alpar@1: if (sv_next[k] == 0) alpar@1: luf->sv_tail = sv_prev[k]; alpar@1: else alpar@1: sv_prev[sv_next[k]] = sv_prev[k]; alpar@1: /* the p-th column of the matrix F has been completely built; set alpar@1: its pointer */ alpar@1: fc_ptr[p] = luf->sv_end; alpar@1: /* walk through the p-th (pivot) row and do the following... */ alpar@1: for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++) alpar@1: { /* get column index of v[p,j] */ alpar@1: j = sv_ind[p_ptr]; alpar@1: /* erase v[p,j] from the working array */ alpar@1: flag[j] = 0; alpar@1: work[j] = 0.0; alpar@1: /* the j-th column has been completely transformed, therefore alpar@1: it can return to the active set with new length; however alpar@1: the special case c_prev[j] = c_next[j] = j means that the alpar@1: routine find_pivot excluded the j-th column from the active alpar@1: set due to Uwe Suhl's rule, and therefore in this case the alpar@1: column can return to the active set only if it is a column alpar@1: singleton */ alpar@1: if (!(vc_len[j] != 1 && cs_prev[j] == j && cs_next[j] == j)) alpar@1: { cs_prev[j] = 0; alpar@1: cs_next[j] = cs_head[vc_len[j]]; alpar@1: if (cs_next[j] != 0) cs_prev[cs_next[j]] = j; alpar@1: cs_head[vc_len[j]] = j; alpar@1: } alpar@1: } alpar@1: done: /* return to the factorizing routine */ alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * build_v_cols - build the matrix V in column-wise format alpar@1: * alpar@1: * This routine builds the column-wise representation of the matrix V alpar@1: * using its row-wise representation. alpar@1: * alpar@1: * If no error occured, the routine returns zero. Otherwise, in case of alpar@1: * overflow of the sparse vector area, the routine returns non-zero. */ alpar@1: alpar@1: static int build_v_cols(LUF *luf) alpar@1: { int n = luf->n; alpar@1: int *vr_ptr = luf->vr_ptr; alpar@1: int *vr_len = luf->vr_len; alpar@1: int *vc_ptr = luf->vc_ptr; alpar@1: int *vc_len = luf->vc_len; alpar@1: int *vc_cap = luf->vc_cap; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: int *sv_prev = luf->sv_prev; alpar@1: int *sv_next = luf->sv_next; alpar@1: int ret = 0; alpar@1: int i, i_beg, i_end, i_ptr, j, j_ptr, k, nnz; alpar@1: /* it is assumed that on entry all columns of the matrix V are alpar@1: empty, i.e. vc_len[j] = vc_cap[j] = 0 for all j = 1, ..., n, alpar@1: and have been removed from the addressing list */ alpar@1: /* count non-zeros in columns of the matrix V; count total number alpar@1: of non-zeros in this matrix */ alpar@1: nnz = 0; alpar@1: for (i = 1; i <= n; i++) alpar@1: { /* walk through elements of the i-th row and count non-zeros alpar@1: in the corresponding columns */ alpar@1: i_beg = vr_ptr[i]; alpar@1: i_end = i_beg + vr_len[i] - 1; alpar@1: for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) alpar@1: vc_cap[sv_ind[i_ptr]]++; alpar@1: /* count total number of non-zeros */ alpar@1: nnz += vr_len[i]; alpar@1: } alpar@1: /* store total number of non-zeros */ alpar@1: luf->nnz_v = nnz; alpar@1: /* check for free locations */ alpar@1: if (luf->sv_end - luf->sv_beg < nnz) alpar@1: { /* overflow of the sparse vector area */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* allocate columns of the matrix V */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { /* set pointer to the j-th column */ alpar@1: vc_ptr[j] = luf->sv_beg; alpar@1: /* reserve locations for the j-th column */ alpar@1: luf->sv_beg += vc_cap[j]; alpar@1: } alpar@1: /* build the matrix V in column-wise format using this matrix in alpar@1: row-wise format */ alpar@1: for (i = 1; i <= n; i++) alpar@1: { /* walk through elements of the i-th row */ alpar@1: i_beg = vr_ptr[i]; alpar@1: i_end = i_beg + vr_len[i] - 1; alpar@1: for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) alpar@1: { /* get column index */ alpar@1: j = sv_ind[i_ptr]; alpar@1: /* store element in the j-th column */ alpar@1: j_ptr = vc_ptr[j] + vc_len[j]; alpar@1: sv_ind[j_ptr] = i; alpar@1: sv_val[j_ptr] = sv_val[i_ptr]; alpar@1: /* increase length of the j-th column */ alpar@1: vc_len[j]++; alpar@1: } alpar@1: } alpar@1: /* now columns are placed in the sparse vector area behind rows alpar@1: in the order n+1, n+2, ..., n+n; so insert column nodes in the alpar@1: addressing list using this order */ alpar@1: for (k = n+1; k <= n+n; k++) alpar@1: { sv_prev[k] = k-1; alpar@1: sv_next[k] = k+1; alpar@1: } alpar@1: sv_prev[n+1] = luf->sv_tail; alpar@1: sv_next[luf->sv_tail] = n+1; alpar@1: sv_next[n+n] = 0; alpar@1: luf->sv_tail = n+n; alpar@1: done: /* return to the factorizing routine */ alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * build_f_rows - build the matrix F in row-wise format alpar@1: * alpar@1: * This routine builds the row-wise representation of the matrix F using alpar@1: * its column-wise representation. alpar@1: * alpar@1: * If no error occured, the routine returns zero. Otherwise, in case of alpar@1: * overflow of the sparse vector area, the routine returns non-zero. */ alpar@1: alpar@1: static int build_f_rows(LUF *luf) alpar@1: { int n = luf->n; alpar@1: int *fr_ptr = luf->fr_ptr; alpar@1: int *fr_len = luf->fr_len; alpar@1: int *fc_ptr = luf->fc_ptr; alpar@1: int *fc_len = luf->fc_len; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: int ret = 0; alpar@1: int i, j, j_beg, j_end, j_ptr, ptr, nnz; alpar@1: /* clear rows of the matrix F */ alpar@1: for (i = 1; i <= n; i++) fr_len[i] = 0; alpar@1: /* count non-zeros in rows of the matrix F; count total number of alpar@1: non-zeros in this matrix */ alpar@1: nnz = 0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { /* walk through elements of the j-th column and count non-zeros alpar@1: in the corresponding rows */ alpar@1: j_beg = fc_ptr[j]; alpar@1: j_end = j_beg + fc_len[j] - 1; alpar@1: for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) alpar@1: fr_len[sv_ind[j_ptr]]++; alpar@1: /* increase total number of non-zeros */ alpar@1: nnz += fc_len[j]; alpar@1: } alpar@1: /* store total number of non-zeros */ alpar@1: luf->nnz_f = nnz; alpar@1: /* check for free locations */ alpar@1: if (luf->sv_end - luf->sv_beg < nnz) alpar@1: { /* overflow of the sparse vector area */ alpar@1: ret = 1; alpar@1: goto done; alpar@1: } alpar@1: /* allocate rows of the matrix F */ alpar@1: for (i = 1; i <= n; i++) alpar@1: { /* set pointer to the end of the i-th row; later this pointer alpar@1: will be set to the beginning of the i-th row */ alpar@1: fr_ptr[i] = luf->sv_end; alpar@1: /* reserve locations for the i-th row */ alpar@1: luf->sv_end -= fr_len[i]; alpar@1: } alpar@1: /* build the matrix F in row-wise format using this matrix in alpar@1: column-wise format */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { /* walk through elements of the j-th column */ alpar@1: j_beg = fc_ptr[j]; alpar@1: j_end = j_beg + fc_len[j] - 1; alpar@1: for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) alpar@1: { /* get row index */ alpar@1: i = sv_ind[j_ptr]; alpar@1: /* store element in the i-th row */ alpar@1: ptr = --fr_ptr[i]; alpar@1: sv_ind[ptr] = j; alpar@1: sv_val[ptr] = sv_val[j_ptr]; alpar@1: } alpar@1: } alpar@1: done: /* return to the factorizing routine */ alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * luf_factorize - compute LU-factorization alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpluf.h" alpar@1: * int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j, alpar@1: * int ind[], double val[]), void *info); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine luf_factorize computes LU-factorization of a specified alpar@1: * square matrix A. alpar@1: * alpar@1: * The parameter luf specifies LU-factorization program object created alpar@1: * by the routine luf_create_it. alpar@1: * alpar@1: * The parameter n specifies the order of A, n > 0. alpar@1: * alpar@1: * The formal routine col specifies the matrix A to be factorized. To alpar@1: * obtain j-th column of A the routine luf_factorize calls the routine alpar@1: * col with the parameter j (1 <= j <= n). In response the routine col alpar@1: * should store row indices and numerical values of non-zero elements alpar@1: * of j-th column of A to locations ind[1,...,len] and val[1,...,len], alpar@1: * respectively, where len is the number of non-zeros in j-th column alpar@1: * returned on exit. Neither zero nor duplicate elements are allowed. alpar@1: * alpar@1: * The parameter info is a transit pointer passed to the routine col. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * 0 LU-factorization has been successfully computed. alpar@1: * alpar@1: * LUF_ESING alpar@1: * The specified matrix is singular within the working precision. alpar@1: * (On some elimination step the active submatrix is exactly zero, alpar@1: * so no pivot can be chosen.) alpar@1: * alpar@1: * LUF_ECOND alpar@1: * The specified matrix is ill-conditioned. alpar@1: * (On some elimination step too intensive growth of elements of the alpar@1: * active submatix has been detected.) alpar@1: * alpar@1: * If matrix A is well scaled, the return code LUF_ECOND may also mean alpar@1: * that the threshold pivoting tolerance piv_tol should be increased. alpar@1: * alpar@1: * In case of non-zero return code the factorization becomes invalid. alpar@1: * It should not be used in other operations until the cause of failure alpar@1: * has been eliminated and the factorization has been recomputed again alpar@1: * with the routine luf_factorize. alpar@1: * alpar@1: * REPAIRING SINGULAR MATRIX alpar@1: * alpar@1: * If the routine luf_factorize returns non-zero code, it provides all alpar@1: * necessary information that can be used for "repairing" the matrix A, alpar@1: * where "repairing" means replacing linearly dependent columns of the alpar@1: * matrix A by appropriate columns of the unity matrix. This feature is alpar@1: * needed when this routine is used for factorizing the basis matrix alpar@1: * within the simplex method procedure. alpar@1: * alpar@1: * On exit linearly dependent columns of the (partially transformed) alpar@1: * matrix U have numbers rank+1, rank+2, ..., n, where rank is estimated alpar@1: * rank of the matrix A stored by the routine to the member luf->rank. alpar@1: * The correspondence between columns of A and U is the same as between alpar@1: * columns of V and U. Thus, linearly dependent columns of the matrix A alpar@1: * have numbers qq_col[rank+1], qq_col[rank+2], ..., qq_col[n], where alpar@1: * qq_col is the column-like representation of the permutation matrix Q. alpar@1: * It is understood that each j-th linearly dependent column of the alpar@1: * matrix U should be replaced by the unity vector, where all elements alpar@1: * are zero except the unity diagonal element u[j,j]. On the other hand alpar@1: * j-th row of the matrix U corresponds to the row of the matrix V (and alpar@1: * therefore of the matrix A) with the number pp_row[j], where pp_row is alpar@1: * the row-like representation of the permutation matrix P. Thus, each alpar@1: * j-th linearly dependent column of the matrix U should be replaced by alpar@1: * column of the unity matrix with the number pp_row[j]. alpar@1: * alpar@1: * The code that repairs the matrix A may look like follows: alpar@1: * alpar@1: * for (j = rank+1; j <= n; j++) alpar@1: * { replace the column qq_col[j] of the matrix A by the column alpar@1: * pp_row[j] of the unity matrix; alpar@1: * } alpar@1: * alpar@1: * where rank, pp_row, and qq_col are members of the structure LUF. */ alpar@1: alpar@1: int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j, alpar@1: int ind[], double val[]), void *info) alpar@1: { int *pp_row, *pp_col, *qq_row, *qq_col; alpar@1: double max_gro = luf->max_gro; alpar@1: int i, j, k, p, q, t, ret; alpar@1: if (n < 1) alpar@1: xfault("luf_factorize: n = %d; invalid parameter\n", n); alpar@1: if (n > N_MAX) alpar@1: xfault("luf_factorize: n = %d; matrix too big\n", n); alpar@1: /* invalidate the factorization */ alpar@1: luf->valid = 0; alpar@1: /* reallocate arrays, if necessary */ alpar@1: reallocate(luf, n); alpar@1: pp_row = luf->pp_row; alpar@1: pp_col = luf->pp_col; alpar@1: qq_row = luf->qq_row; alpar@1: qq_col = luf->qq_col; alpar@1: /* estimate initial size of the SVA, if not specified */ alpar@1: if (luf->sv_size == 0 && luf->new_sva == 0) alpar@1: luf->new_sva = 5 * (n + 10); alpar@1: more: /* reallocate the sparse vector area, if required */ alpar@1: if (luf->new_sva > 0) alpar@1: { if (luf->sv_ind != NULL) xfree(luf->sv_ind); alpar@1: if (luf->sv_val != NULL) xfree(luf->sv_val); alpar@1: luf->sv_size = luf->new_sva; alpar@1: luf->sv_ind = xcalloc(1+luf->sv_size, sizeof(int)); alpar@1: luf->sv_val = xcalloc(1+luf->sv_size, sizeof(double)); alpar@1: luf->new_sva = 0; alpar@1: } alpar@1: /* initialize LU-factorization data structures */ alpar@1: if (initialize(luf, col, info)) alpar@1: { /* overflow of the sparse vector area */ alpar@1: luf->new_sva = luf->sv_size + luf->sv_size; alpar@1: xassert(luf->new_sva > luf->sv_size); alpar@1: goto more; alpar@1: } alpar@1: /* main elimination loop */ alpar@1: for (k = 1; k <= n; k++) alpar@1: { /* choose a pivot element v[p,q] */ alpar@1: if (find_pivot(luf, &p, &q)) alpar@1: { /* no pivot can be chosen, because the active submatrix is alpar@1: exactly zero */ alpar@1: luf->rank = k - 1; alpar@1: ret = LUF_ESING; alpar@1: goto done; alpar@1: } alpar@1: /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th alpar@1: rows and k-th and j'-th columns of the matrix U = P*V*Q to alpar@1: move the element u[i',j'] to the position u[k,k] */ alpar@1: i = pp_col[p], j = qq_row[q]; alpar@1: xassert(k <= i && i <= n && k <= j && j <= n); alpar@1: /* permute k-th and i-th rows of the matrix U */ alpar@1: t = pp_row[k]; alpar@1: pp_row[i] = t, pp_col[t] = i; alpar@1: pp_row[k] = p, pp_col[p] = k; alpar@1: /* permute k-th and j-th columns of the matrix U */ alpar@1: t = qq_col[k]; alpar@1: qq_col[j] = t, qq_row[t] = j; alpar@1: qq_col[k] = q, qq_row[q] = k; alpar@1: /* eliminate subdiagonal elements of k-th column of the matrix alpar@1: U = P*V*Q using the pivot element u[k,k] = v[p,q] */ alpar@1: if (eliminate(luf, p, q)) alpar@1: { /* overflow of the sparse vector area */ alpar@1: luf->new_sva = luf->sv_size + luf->sv_size; alpar@1: xassert(luf->new_sva > luf->sv_size); alpar@1: goto more; alpar@1: } alpar@1: /* check relative growth of elements of the matrix V */ alpar@1: if (luf->big_v > max_gro * luf->max_a) alpar@1: { /* the growth is too intensive, therefore most probably the alpar@1: matrix A is ill-conditioned */ alpar@1: luf->rank = k - 1; alpar@1: ret = LUF_ECOND; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* now the matrix U = P*V*Q is upper triangular, the matrix V has alpar@1: been built in row-wise format, and the matrix F has been built alpar@1: in column-wise format */ alpar@1: /* defragment the sparse vector area in order to merge all free alpar@1: locations in one continuous extent */ alpar@1: luf_defrag_sva(luf); alpar@1: /* build the matrix V in column-wise format */ alpar@1: if (build_v_cols(luf)) alpar@1: { /* overflow of the sparse vector area */ alpar@1: luf->new_sva = luf->sv_size + luf->sv_size; alpar@1: xassert(luf->new_sva > luf->sv_size); alpar@1: goto more; alpar@1: } alpar@1: /* build the matrix F in row-wise format */ alpar@1: if (build_f_rows(luf)) alpar@1: { /* overflow of the sparse vector area */ alpar@1: luf->new_sva = luf->sv_size + luf->sv_size; alpar@1: xassert(luf->new_sva > luf->sv_size); alpar@1: goto more; alpar@1: } alpar@1: /* the LU-factorization has been successfully computed */ alpar@1: luf->valid = 1; alpar@1: luf->rank = n; alpar@1: ret = 0; alpar@1: /* if there are few free locations in the sparse vector area, try alpar@1: increasing its size in the future */ alpar@1: t = 3 * (n + luf->nnz_v) + 2 * luf->nnz_f; alpar@1: if (luf->sv_size < t) alpar@1: { luf->new_sva = luf->sv_size; alpar@1: while (luf->new_sva < t) alpar@1: { k = luf->new_sva; alpar@1: luf->new_sva = k + k; alpar@1: xassert(luf->new_sva > k); alpar@1: } alpar@1: } alpar@1: done: /* return to the calling program */ alpar@1: return ret; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * luf_f_solve - solve system F*x = b or F'*x = b alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpluf.h" alpar@1: * void luf_f_solve(LUF *luf, int tr, double x[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine luf_f_solve solves either the system F*x = b (if the alpar@1: * flag tr is zero) or the system F'*x = b (if the flag tr is non-zero), alpar@1: * where the matrix F is a component of LU-factorization specified by alpar@1: * the parameter luf, F' is a matrix transposed to F. alpar@1: * alpar@1: * On entry the array x should contain elements of the right-hand side alpar@1: * vector b in locations x[1], ..., x[n], where n is the order of the alpar@1: * matrix F. On exit this array will contain elements of the solution alpar@1: * vector x in the same locations. */ alpar@1: alpar@1: void luf_f_solve(LUF *luf, int tr, double x[]) alpar@1: { int n = luf->n; alpar@1: int *fr_ptr = luf->fr_ptr; alpar@1: int *fr_len = luf->fr_len; alpar@1: int *fc_ptr = luf->fc_ptr; alpar@1: int *fc_len = luf->fc_len; alpar@1: int *pp_row = luf->pp_row; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: int i, j, k, beg, end, ptr; alpar@1: double xk; alpar@1: if (!luf->valid) alpar@1: xfault("luf_f_solve: LU-factorization is not valid\n"); alpar@1: if (!tr) alpar@1: { /* solve the system F*x = b */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { k = pp_row[j]; alpar@1: xk = x[k]; alpar@1: if (xk != 0.0) alpar@1: { beg = fc_ptr[k]; alpar@1: end = beg + fc_len[k] - 1; alpar@1: for (ptr = beg; ptr <= end; ptr++) alpar@1: x[sv_ind[ptr]] -= sv_val[ptr] * xk; alpar@1: } alpar@1: } alpar@1: } alpar@1: else alpar@1: { /* solve the system F'*x = b */ alpar@1: for (i = n; i >= 1; i--) alpar@1: { k = pp_row[i]; alpar@1: xk = x[k]; alpar@1: if (xk != 0.0) alpar@1: { beg = fr_ptr[k]; alpar@1: end = beg + fr_len[k] - 1; alpar@1: for (ptr = beg; ptr <= end; ptr++) alpar@1: x[sv_ind[ptr]] -= sv_val[ptr] * xk; alpar@1: } alpar@1: } alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * luf_v_solve - solve system V*x = b or V'*x = b alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpluf.h" alpar@1: * void luf_v_solve(LUF *luf, int tr, double x[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine luf_v_solve solves either the system V*x = b (if the alpar@1: * flag tr is zero) or the system V'*x = b (if the flag tr is non-zero), alpar@1: * where the matrix V is a component of LU-factorization specified by alpar@1: * the parameter luf, V' is a matrix transposed to V. alpar@1: * alpar@1: * On entry the array x should contain elements of the right-hand side alpar@1: * vector b in locations x[1], ..., x[n], where n is the order of the alpar@1: * matrix V. On exit this array will contain elements of the solution alpar@1: * vector x in the same locations. */ alpar@1: alpar@1: void luf_v_solve(LUF *luf, int tr, double x[]) alpar@1: { int n = luf->n; alpar@1: int *vr_ptr = luf->vr_ptr; alpar@1: int *vr_len = luf->vr_len; alpar@1: double *vr_piv = luf->vr_piv; alpar@1: int *vc_ptr = luf->vc_ptr; alpar@1: int *vc_len = luf->vc_len; alpar@1: int *pp_row = luf->pp_row; alpar@1: int *qq_col = luf->qq_col; alpar@1: int *sv_ind = luf->sv_ind; alpar@1: double *sv_val = luf->sv_val; alpar@1: double *b = luf->work; alpar@1: int i, j, k, beg, end, ptr; alpar@1: double temp; alpar@1: if (!luf->valid) alpar@1: xfault("luf_v_solve: LU-factorization is not valid\n"); alpar@1: for (k = 1; k <= n; k++) b[k] = x[k], x[k] = 0.0; alpar@1: if (!tr) alpar@1: { /* solve the system V*x = b */ alpar@1: for (k = n; k >= 1; k--) alpar@1: { i = pp_row[k], j = qq_col[k]; alpar@1: temp = b[i]; alpar@1: if (temp != 0.0) alpar@1: { x[j] = (temp /= vr_piv[i]); alpar@1: beg = vc_ptr[j]; alpar@1: end = beg + vc_len[j] - 1; alpar@1: for (ptr = beg; ptr <= end; ptr++) alpar@1: b[sv_ind[ptr]] -= sv_val[ptr] * temp; alpar@1: } alpar@1: } alpar@1: } alpar@1: else alpar@1: { /* solve the system V'*x = b */ alpar@1: for (k = 1; k <= n; k++) alpar@1: { i = pp_row[k], j = qq_col[k]; alpar@1: temp = b[j]; alpar@1: if (temp != 0.0) alpar@1: { x[i] = (temp /= vr_piv[i]); alpar@1: beg = vr_ptr[i]; alpar@1: end = beg + vr_len[i] - 1; alpar@1: for (ptr = beg; ptr <= end; ptr++) alpar@1: b[sv_ind[ptr]] -= sv_val[ptr] * temp; alpar@1: } alpar@1: } alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * luf_a_solve - solve system A*x = b or A'*x = b alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpluf.h" alpar@1: * void luf_a_solve(LUF *luf, int tr, double x[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine luf_a_solve solves either the system A*x = b (if the alpar@1: * flag tr is zero) or the system A'*x = b (if the flag tr is non-zero), alpar@1: * where the parameter luf specifies LU-factorization of the matrix A, alpar@1: * A' is a matrix transposed to A. alpar@1: * alpar@1: * On entry the array x should contain elements of the right-hand side alpar@1: * vector b in locations x[1], ..., x[n], where n is the order of the alpar@1: * matrix A. On exit this array will contain elements of the solution alpar@1: * vector x in the same locations. */ alpar@1: alpar@1: void luf_a_solve(LUF *luf, int tr, double x[]) alpar@1: { if (!luf->valid) alpar@1: xfault("luf_a_solve: LU-factorization is not valid\n"); alpar@1: if (!tr) alpar@1: { /* A = F*V, therefore inv(A) = inv(V)*inv(F) */ alpar@1: luf_f_solve(luf, 0, x); alpar@1: luf_v_solve(luf, 0, x); alpar@1: } alpar@1: else alpar@1: { /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */ alpar@1: luf_v_solve(luf, 1, x); alpar@1: luf_f_solve(luf, 1, x); alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * luf_delete_it - delete LU-factorization alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpluf.h" alpar@1: * void luf_delete_it(LUF *luf); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * The routine luf_delete deletes LU-factorization specified by the alpar@1: * parameter luf and frees all the memory allocated to this program alpar@1: * object. */ alpar@1: alpar@1: void luf_delete_it(LUF *luf) alpar@1: { if (luf->fr_ptr != NULL) xfree(luf->fr_ptr); alpar@1: if (luf->fr_len != NULL) xfree(luf->fr_len); alpar@1: if (luf->fc_ptr != NULL) xfree(luf->fc_ptr); alpar@1: if (luf->fc_len != NULL) xfree(luf->fc_len); alpar@1: if (luf->vr_ptr != NULL) xfree(luf->vr_ptr); alpar@1: if (luf->vr_len != NULL) xfree(luf->vr_len); alpar@1: if (luf->vr_cap != NULL) xfree(luf->vr_cap); alpar@1: if (luf->vr_piv != NULL) xfree(luf->vr_piv); alpar@1: if (luf->vc_ptr != NULL) xfree(luf->vc_ptr); alpar@1: if (luf->vc_len != NULL) xfree(luf->vc_len); alpar@1: if (luf->vc_cap != NULL) xfree(luf->vc_cap); alpar@1: if (luf->pp_row != NULL) xfree(luf->pp_row); alpar@1: if (luf->pp_col != NULL) xfree(luf->pp_col); alpar@1: if (luf->qq_row != NULL) xfree(luf->qq_row); alpar@1: if (luf->qq_col != NULL) xfree(luf->qq_col); alpar@1: if (luf->sv_ind != NULL) xfree(luf->sv_ind); alpar@1: if (luf->sv_val != NULL) xfree(luf->sv_val); alpar@1: if (luf->sv_prev != NULL) xfree(luf->sv_prev); alpar@1: if (luf->sv_next != NULL) xfree(luf->sv_next); alpar@1: if (luf->vr_max != NULL) xfree(luf->vr_max); alpar@1: if (luf->rs_head != NULL) xfree(luf->rs_head); alpar@1: if (luf->rs_prev != NULL) xfree(luf->rs_prev); alpar@1: if (luf->rs_next != NULL) xfree(luf->rs_next); alpar@1: if (luf->cs_head != NULL) xfree(luf->cs_head); alpar@1: if (luf->cs_prev != NULL) xfree(luf->cs_prev); alpar@1: if (luf->cs_next != NULL) xfree(luf->cs_next); alpar@1: if (luf->flag != NULL) xfree(luf->flag); alpar@1: if (luf->work != NULL) xfree(luf->work); alpar@1: xfree(luf); alpar@1: return; alpar@1: } alpar@1: alpar@1: /* eof */