alpar@9: /* glplux.h (LU-factorization, bignum arithmetic) */ alpar@9: alpar@9: /*********************************************************************** alpar@9: * This code is part of GLPK (GNU Linear Programming Kit). alpar@9: * alpar@9: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@9: * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, alpar@9: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@9: * E-mail: . alpar@9: * alpar@9: * GLPK is free software: you can redistribute it and/or modify it alpar@9: * under the terms of the GNU General Public License as published by alpar@9: * the Free Software Foundation, either version 3 of the License, or alpar@9: * (at your option) any later version. alpar@9: * alpar@9: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@9: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@9: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@9: * License for more details. alpar@9: * alpar@9: * You should have received a copy of the GNU General Public License alpar@9: * along with GLPK. If not, see . alpar@9: ***********************************************************************/ alpar@9: alpar@9: #ifndef GLPLUX_H alpar@9: #define GLPLUX_H alpar@9: alpar@9: #include "glpdmp.h" alpar@9: #include "glpgmp.h" alpar@9: alpar@9: /*---------------------------------------------------------------------- alpar@9: // The structure LUX defines LU-factorization of a square matrix A, alpar@9: // which is the following quartet: alpar@9: // alpar@9: // [A] = (F, V, P, Q), (1) alpar@9: // alpar@9: // where F and V are such matrices that alpar@9: // alpar@9: // A = F * V, (2) alpar@9: // alpar@9: // and P and Q are such permutation matrices that the matrix alpar@9: // alpar@9: // L = P * F * inv(P) (3) alpar@9: // alpar@9: // is lower triangular with unity diagonal, and the matrix alpar@9: // alpar@9: // U = P * V * Q (4) alpar@9: // alpar@9: // is upper triangular. All the matrices have the order n. alpar@9: // alpar@9: // The matrices F and V are stored in row/column-wise sparse format as alpar@9: // row and column linked lists of non-zero elements. Unity elements on alpar@9: // the main diagonal of the matrix F are not stored. Pivot elements of alpar@9: // the matrix V (that correspond to diagonal elements of the matrix U) alpar@9: // are also missing from the row and column lists and stored separately alpar@9: // in an ordinary array. alpar@9: // alpar@9: // The permutation matrices P and Q are stored as ordinary arrays using alpar@9: // both row- and column-like formats. alpar@9: // alpar@9: // The matrices L and U being completely defined by the matrices F, V, alpar@9: // P, and Q are not stored explicitly. alpar@9: // alpar@9: // It is easy to show that the factorization (1)-(3) is some version of alpar@9: // LU-factorization. Indeed, from (3) and (4) it follows that: alpar@9: // alpar@9: // F = inv(P) * L * P, alpar@9: // alpar@9: // V = inv(P) * U * inv(Q), alpar@9: // alpar@9: // and substitution into (2) gives: alpar@9: // alpar@9: // A = F * V = inv(P) * L * U * inv(Q). alpar@9: // alpar@9: // For more details see the program documentation. */ alpar@9: alpar@9: typedef struct LUX LUX; alpar@9: typedef struct LUXELM LUXELM; alpar@9: typedef struct LUXWKA LUXWKA; alpar@9: alpar@9: struct LUX alpar@9: { /* LU-factorization of a square matrix */ alpar@9: int n; alpar@9: /* the order of matrices A, F, V, P, Q */ alpar@9: DMP *pool; alpar@9: /* memory pool for elements of matrices F and V */ alpar@9: LUXELM **F_row; /* LUXELM *F_row[1+n]; */ alpar@9: /* F_row[0] is not used; alpar@9: F_row[i], 1 <= i <= n, is a pointer to the list of elements in alpar@9: i-th row of matrix F (diagonal elements are not stored) */ alpar@9: LUXELM **F_col; /* LUXELM *F_col[1+n]; */ alpar@9: /* F_col[0] is not used; alpar@9: F_col[j], 1 <= j <= n, is a pointer to the list of elements in alpar@9: j-th column of matrix F (diagonal elements are not stored) */ alpar@9: mpq_t *V_piv; /* mpq_t V_piv[1+n]; */ alpar@9: /* V_piv[0] is not used; alpar@9: V_piv[p], 1 <= p <= n, is a pivot element v[p,q] corresponding alpar@9: to a diagonal element u[k,k] of matrix U = P*V*Q (used on k-th alpar@9: elimination step, k = 1, 2, ..., n) */ alpar@9: LUXELM **V_row; /* LUXELM *V_row[1+n]; */ alpar@9: /* V_row[0] is not used; alpar@9: V_row[i], 1 <= i <= n, is a pointer to the list of elements in alpar@9: i-th row of matrix V (except pivot elements) */ alpar@9: LUXELM **V_col; /* LUXELM *V_col[1+n]; */ alpar@9: /* V_col[0] is not used; alpar@9: V_col[j], 1 <= j <= n, is a pointer to the list of elements in alpar@9: j-th column of matrix V (except pivot elements) */ alpar@9: int *P_row; /* int P_row[1+n]; */ alpar@9: /* P_row[0] is not used; alpar@9: P_row[i] = j means that p[i,j] = 1, where p[i,j] is an element alpar@9: of permutation matrix P */ alpar@9: int *P_col; /* int P_col[1+n]; */ alpar@9: /* P_col[0] is not used; alpar@9: P_col[j] = i means that p[i,j] = 1, where p[i,j] is an element alpar@9: of permutation matrix P */ alpar@9: /* if i-th row or column of matrix F is i'-th row or column of alpar@9: matrix L = P*F*inv(P), or if i-th row of matrix V is i'-th row alpar@9: of matrix U = P*V*Q, then P_row[i'] = i and P_col[i] = i' */ alpar@9: int *Q_row; /* int Q_row[1+n]; */ alpar@9: /* Q_row[0] is not used; alpar@9: Q_row[i] = j means that q[i,j] = 1, where q[i,j] is an element alpar@9: of permutation matrix Q */ alpar@9: int *Q_col; /* int Q_col[1+n]; */ alpar@9: /* Q_col[0] is not used; alpar@9: Q_col[j] = i means that q[i,j] = 1, where q[i,j] is an element alpar@9: of permutation matrix Q */ alpar@9: /* if j-th column of matrix V is j'-th column of matrix U = P*V*Q, alpar@9: then Q_row[j] = j' and Q_col[j'] = j */ alpar@9: int rank; alpar@9: /* the (exact) rank of matrices A and V */ alpar@9: }; alpar@9: alpar@9: struct LUXELM alpar@9: { /* element of matrix F or V */ alpar@9: int i; alpar@9: /* row index, 1 <= i <= m */ alpar@9: int j; alpar@9: /* column index, 1 <= j <= n */ alpar@9: mpq_t val; alpar@9: /* numeric (non-zero) element value */ alpar@9: LUXELM *r_prev; alpar@9: /* pointer to previous element in the same row */ alpar@9: LUXELM *r_next; alpar@9: /* pointer to next element in the same row */ alpar@9: LUXELM *c_prev; alpar@9: /* pointer to previous element in the same column */ alpar@9: LUXELM *c_next; alpar@9: /* pointer to next element in the same column */ alpar@9: }; alpar@9: alpar@9: struct LUXWKA alpar@9: { /* working area (used only during factorization) */ alpar@9: /* in order to efficiently implement Markowitz strategy and Duff alpar@9: search technique there are two families {R[0], R[1], ..., R[n]} alpar@9: and {C[0], C[1], ..., C[n]}; member R[k] is a set of active alpar@9: rows of matrix V having k non-zeros, and member C[k] is a set alpar@9: of active columns of matrix V having k non-zeros (in the active alpar@9: submatrix); each set R[k] and C[k] is implemented as a separate alpar@9: doubly linked list */ alpar@9: int *R_len; /* int R_len[1+n]; */ alpar@9: /* R_len[0] is not used; alpar@9: R_len[i], 1 <= i <= n, is the number of non-zero elements in alpar@9: i-th row of matrix V (that is the length of i-th row) */ alpar@9: int *R_head; /* int R_head[1+n]; */ alpar@9: /* R_head[k], 0 <= k <= n, is the number of a first row, which is alpar@9: active and whose length is k */ alpar@9: int *R_prev; /* int R_prev[1+n]; */ alpar@9: /* R_prev[0] is not used; alpar@9: R_prev[i], 1 <= i <= n, is the number of a previous row, which alpar@9: is active and has the same length as i-th row */ alpar@9: int *R_next; /* int R_next[1+n]; */ alpar@9: /* R_prev[0] is not used; alpar@9: R_prev[i], 1 <= i <= n, is the number of a next row, which is alpar@9: active and has the same length as i-th row */ alpar@9: int *C_len; /* int C_len[1+n]; */ alpar@9: /* C_len[0] is not used; alpar@9: C_len[j], 1 <= j <= n, is the number of non-zero elements in alpar@9: j-th column of the active submatrix of matrix V (that is the alpar@9: length of j-th column in the active submatrix) */ alpar@9: int *C_head; /* int C_head[1+n]; */ alpar@9: /* C_head[k], 0 <= k <= n, is the number of a first column, which alpar@9: is active and whose length is k */ alpar@9: int *C_prev; /* int C_prev[1+n]; */ alpar@9: /* C_prev[0] is not used; alpar@9: C_prev[j], 1 <= j <= n, is the number of a previous column, alpar@9: which is active and has the same length as j-th column */ alpar@9: int *C_next; /* int C_next[1+n]; */ alpar@9: /* C_next[0] is not used; alpar@9: C_next[j], 1 <= j <= n, is the number of a next column, which alpar@9: is active and has the same length as j-th column */ alpar@9: }; alpar@9: alpar@9: #define lux_create _glp_lux_create alpar@9: #define lux_decomp _glp_lux_decomp alpar@9: #define lux_f_solve _glp_lux_f_solve alpar@9: #define lux_v_solve _glp_lux_v_solve alpar@9: #define lux_solve _glp_lux_solve alpar@9: #define lux_delete _glp_lux_delete alpar@9: alpar@9: LUX *lux_create(int n); alpar@9: /* create LU-factorization */ alpar@9: alpar@9: int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[], alpar@9: mpq_t val[]), void *info); alpar@9: /* compute LU-factorization */ alpar@9: alpar@9: void lux_f_solve(LUX *lux, int tr, mpq_t x[]); alpar@9: /* solve system F*x = b or F'*x = b */ alpar@9: alpar@9: void lux_v_solve(LUX *lux, int tr, mpq_t x[]); alpar@9: /* solve system V*x = b or V'*x = b */ alpar@9: alpar@9: void lux_solve(LUX *lux, int tr, mpq_t x[]); alpar@9: /* solve system A*x = b or A'*x = b */ alpar@9: alpar@9: void lux_delete(LUX *lux); alpar@9: /* delete LU-factorization */ alpar@9: alpar@9: #endif alpar@9: alpar@9: /* eof */