/* scf.h (sparse updatable Schur-complement-based factorization) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2013-2014 Andrew Makhorin, Department for Applied
*  Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
*  reserved. E-mail: <mao@gnu.org>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#ifndef SCF_H
#define SCF_H

#include "btf.h"
#include "ifu.h"
#include "luf.h"

/***********************************************************************
*  The structure SCF describes sparse updatable factorization based on
*  Schur complement.
*
*  The SCF-factorization has the following format:
*
*     ( A   A1~ )     ( A0  A1 )       ( R0    ) ( S0  S )
*     (         ) = P (        ) Q = P (       ) (       ) Q,        (1)
*     ( A2~ A3~ )     ( A2  A3 )       ( R   I ) (     C )
*
*  where:
*
*  A is current (unsymmetric) square matrix (not stored);
*
*  A1~, A2~, A3~ are some additional matrices (not stored);
*
*  A0 is initial (unsymmetric) square matrix (not stored);
*
*  A1, A2, A3 are some additional matrices (not stored);
*
*  R0 and S0 are matrices that define factorization of the initial
*  matrix A0 = R0 * S0 (stored in an invertable form);
*
*  R is a matrix defined from R * S0 = A2, so R = A2 * inv(S0) (stored
*  in row-wise sparse format);
*
*  S is a matrix defined from R0 * S = A1, so S = inv(R0) * A1 (stored
*  in column-wise sparse format);
*
*  C is Schur complement (to matrix A0) defined from R * S + C = A3,
*  so C = A3 - R * S = A3 - A2 * inv(A0) * A1 (stored in an invertable
*  form).
*
*  P, Q are permutation matrices (stored in both row- and column-like
*  formats). */

typedef struct SCF SCF;

struct SCF
{     /* Schur-complement-based factorization */
      int n;
      /* order of current matrix A */
      /*--------------------------------------------------------------*/
      /* initial matrix A0 = R0 * S0 of order n0 in invertable form */
      int n0;
      /* order of matrix A0 */
      int type;
      /* type of factorization used:
       * 1 - LU-factorization (R0 = F0, S0 = V0)
       * 2 - BT-factorization (R0 = I, S0 = A0) */
      union
      {  LUF *luf; /* type = 1 */
         BTF *btf; /* type = 2 */
      }  a0;
      /* factorization of matrix A0 */
      /*--------------------------------------------------------------*/
      /* augmented matrix (A0, A1; A2, A3) of order n0+nn */
      int nn_max;
      /* maximal number of additional rows and columns in the augmented
       * matrix (this limits the number of updates) */
      int nn;
      /* current number of additional rows and columns in the augmented
       * matrix, 0 <= nn <= nn_max */
      SVA *sva;
      /* associated sparse vector area (SVA) used to store rows of
       * matrix R and columns of matrix S */
      /*--------------------------------------------------------------*/
      /* nn*n0-matrix R in row-wise format */
      int rr_ref;
      /* reference number of sparse vector in SVA, which is the first
       * row of matrix R */
#if 0 + 0
      int *rr_ptr = &sva->ptr[rr_ref-1];
      /* rr_ptr[0] is not used;
       * rr_ptr[i], 1 <= i <= nn, is pointer to i-th row in SVA;
       * rr_ptr[nn+1,...,nn_max] are reserved locations */
      int *rr_len = &sva->len[rr_ref-1];
      /* rr_len[0] is not used;
       * rr_len[i], 1 <= i <= nn, is length of i-th row;
       * rr_len[nn+1,...,nn_max] are reserved locations */
#endif
      /*--------------------------------------------------------------*/
      /* n0*nn-matrix S in column-wise format */
      int ss_ref;
      /* reference number of sparse vector in SVA, which is the first
       * column of matrix S */
#if 0 + 0
      int *ss_ptr = &sva->ptr[ss_ref-1];
      /* ss_ptr[0] is not used;
       * ss_ptr[j], 1 <= j <= nn, is pointer to j-th column in SVA;
       * ss_ptr[nn+1,...,nn_max] are reserved locations */
      int *ss_len = &sva->len[ss_ref-1];
      /* ss_len[0] is not used;
       * ss_len[j], 1 <= j <= nn, is length of j-th column;
       * ss_len[nn+1,...,nn_max] are reserved locations */
#endif
      /*--------------------------------------------------------------*/
      /* Schur complement C of order nn in invertable form */
      IFU ifu;
      /* IFU-factorization of matrix C */
      /*--------------------------------------------------------------*/
      /* permutation matrix P of order n0+nn */
      int *pp_ind; /* int pp_ind[1+n0+nn_max]; */
      /* pp_ind[i] = j means that P[i,j] = 1 */
      int *pp_inv; /* int pp_inv[1+n0+nn_max]; */
      /* pp_inv[j] = i means that P[i,j] = 1 */
      /*--------------------------------------------------------------*/
      /* permutation matrix Q of order n0+nn */
      int *qq_ind; /* int qq_ind[1+n0+nn_max]; */
      /* qq_ind[i] = j means that Q[i,j] = 1 */
      int *qq_inv; /* int qq_inv[1+n0+nn_max]; */
      /* qq_inv[j] = i means that Q[i,j] = 1 */
};

#define scf_swap_q_cols(j1, j2) \
      do \
      {  int i1, i2; \
         i1 = qq_inv[j1], i2 = qq_inv[j2]; \
         qq_ind[i1] = j2, qq_inv[j2] = i1; \
         qq_ind[i2] = j1, qq_inv[j1] = i2; \
      }  while (0)
/* swap columns j1 and j2 of permutation matrix Q */

#define scf_r0_solve _glp_scf_r0_solve
void scf_r0_solve(SCF *scf, int tr, double x[/*1+n0*/]);
/* solve system R0 * x = b or R0'* x = b */

#define scf_s0_solve _glp_scf_s0_solve
void scf_s0_solve(SCF *scf, int tr, double x[/*1+n0*/],
      double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]);
/* solve system S0 * x = b or S0'* x = b */

#define scf_r_prod _glp_scf_r_prod
void scf_r_prod(SCF *scf, double y[/*1+nn*/], double a, const double
      x[/*1+n0*/]);
/* compute product y := y + alpha * R * x */

#define scf_rt_prod _glp_scf_rt_prod
void scf_rt_prod(SCF *scf, double y[/*1+n0*/], double a, const double
      x[/*1+nn*/]);
/* compute product y := y + alpha * R'* x */

#define scf_s_prod _glp_scf_s_prod
void scf_s_prod(SCF *scf, double y[/*1+n0*/], double a, const double
      x[/*1+nn*/]);
/* compute product y := y + alpha * S * x */

#define scf_st_prod _glp_scf_st_prod
void scf_st_prod(SCF *scf, double y[/*1+nn*/], double a, const double
      x[/*1+n0*/]);
/* compute product y := y + alpha * S'* x */

#define scf_a_solve _glp_scf_a_solve
void scf_a_solve(SCF *scf, double x[/*1+n*/],
      double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/],
      double work2[/*1+n*/], double work3[/*1+n*/]);
/* solve system A * x = b */

#define scf_at_solve _glp_scf_at_solve
void scf_at_solve(SCF *scf, double x[/*1+n*/],
      double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/],
      double work2[/*1+n*/], double work3[/*1+n*/]);
/* solve system A'* x = b */

#define scf_add_r_row _glp_scf_add_r_row
void scf_add_r_row(SCF *scf, const double w[/*1+n0*/]);
/* add new row to matrix R */

#define scf_add_s_col _glp_scf_add_s_col
void scf_add_s_col(SCF *scf, const double v[/*1+n0*/]);
/* add new column to matrix S */

#define scf_update_aug _glp_scf_update_aug
int scf_update_aug(SCF *scf, double b[/*1+n0*/], double d[/*1+n0*/],
      double f[/*1+nn*/], double g[/*1+nn*/], double h, int upd,
      double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]);
/* update factorization of augmented matrix */

#endif

/* eof */