/* Gaussian triangulation and elimination. */
/* Last edited on 2005-01-16 15:20:11 by stolfi */
/* Copyright © 2005 Jorge Stolfi, UNICAMP. See note at end of file. */

#ifndef gauss_elim_H
#define gauss_elim_H

#include <js.h>

/* In all the procedures below, tw-dimensional matrices are stored
  into one-dimensional vectors, in row-by-row order. That is, an {m×n}
  matrix {A[0..m-1,0..n-1]} is stored as a vector {A[0..m*n-1]}, 
  with entry {A[i,j]} of the matrix in element {A[n*i+j]} of the
  vector. */

void gsel_triangularize(int m, int n, double *A, int total);
  /* Applies the Gaussian elimination method to the {m×n} matrix
    {A}, leaving it upper triangular. 
    
    Specifically, the matrix {A} is modified by row operations that do
    not change its determinant. Let {lead(A,i)} be the column index of
    the first non-zero element in row {i} of matrix {A}, or {+oo} if
    that row is all zeros. Upon exit, in any case, we will have
    {lead(A,i) >= i} for every {i}. Thus, in particular, if {m = n},
    the determinant of {A} will be the product of its diagonal
    elements.
    
    If {total == TRUE}, the output matrix satisfies a stronger
    condition: for all {i < k}, {lead(A,i) < lead(A,k)}, or {lead(A,i)
    == lead(A,k) = +oo}. In this case, the non-zero rows of {A} are a
    basis for the space spanned by the original rows. */

void gsel_diagonalize(int m, int n, double *A, int total);
  /* Assumes that the {m×n} matrix {A} has been triangularized,
    namely that {lead(A,i) >= i} for every {i}.  Applies 
    row operations to {A} so that it becomes diagonal-like,
    without change to its determinant.
    
    If {total == 0}, then, whenever {A[i,i] != 0}, all other elements
    on column {i} are set to zero. 
    
    If {total == 1}, then, whenever {j := lead(A,i)} is finite,
    all elements {A[k,j]} in rows {k != i} are set to zero. */
    
double gsel_det(int m, int n, double *A);
  /* Assumes that the {m×n} matrix {A} has been triangularized.
    Returns the product of the elements on the main diagonal.
    Note that when {m == n} the result is the determinant of {A}. */

void gsel_normalize(int m, int n, double *A, int total);
  /* Assumes that the {m×n} matrix {A} has been 
    triangularized (and possibly diagonalized), namely that {lead(A,i) >= i}
    for every {i}.  
    
    If {total == 0}, scales each row so that the diagonal element is 1.0.
    This will fail if any diagonal element is zero.
    
    If {total == 1}, then, whenever {j := lead(A,i)} is finite,
    scales the row so that {A[i,j] == 1.0}.  This operation will
    always succeed. */

void gsel_extract_solution(int m, int n, double *A, int p, double *X, int total);
  /* Assumes that {A} is an {m×n} matrix that has been
    triangularized, diagonalized, and normalized with the 
    specified {total} flag.
    
    The procedure also assumes that {n >= m+p}, and interprets {A} as
    a composite of two blocks {A = (M U)} where {M} is {m×n-p} and
    {U} is {m×p}. The procedure tries to solve the {p} linear
    systems {M x = u} where {x} is an unknown vector of size {n-p}
    and {u} is replaced in turn by each of the columns of {U}. 
    That is, it tries to compute an {(n-p)×p} matrix {X} such
    that {M X + U = 0}.  
    
    However, this goal may fail if the matrix {M} is not diagonal.
    Also, the solution may not be unique, in which case some elements
    of {X} are arbitrarily set to 0. Specifically, let {Y := M X + U}.
    If {total = FALSE}:
    
      * If {lead(A,i) > i} for any {i}, then {X[i,0..p-1]} will be
      set to 0, and {Y[i,0..p-1]} may be nonzero.

    If {total = TRUE}:
    
      * If {j} in {0..n-p-1} is not {lead(A,i)} for any
      {i}, then {X[j,0..p-1]} will be set to 0.
      
      * If {i} in {0..m} is such that {lead(A,i) >= n-p},
      then {Y[i,0..p-1]} may be nonzero.
    
    The matrix {X[i,k]} is stored in the given vector {X},
    linearized by rows. */

#endif
/* COPYRIGHT AND AUTHORSHIP NOTICE

  Copyright © 2005 Jorge Stolfi, Universidade Estadual de Campinas (UNICAMP).
  Created by Jorge Stolfi in 1992--2005.
  
  This source file can be freely distributed, used, and modified,
  provided that this copyright and authorship notice is preserved in
  all copies, and that any modified versions of this file are clearly
  marked as such.
  
  This software has NO WARRANTY of correctness or applicability for
  any purpose. Neither the author nor his employers shall be held
  responsible for any losses or damages that may result from its use.

  END OF NOTICE */