/* MxN matrices and operations on them. */
/* Last edited on 2005-01-16 15:20:30 by stolfi */
/* Copyright © 2005 Jorge Stolfi, UNICAMP. See note at end of file. */

#ifndef rmxn_H
#define rmxn_H

#include <stdio.h>
#include <rn.h>

/* Matrices are assumed to be stored linearized by rows, without gaps;
  so, in a matrix {A} with {m} rows and {n} columns, the logical 
  element {A[i,j]} is actually {A[n*i + j]}.
  
  Unless said otherwise, all output matrices and vectors must be 
  disjoint from each other from all the input ones. */

void rmxn_zero(int m, int n, double *R);
  /* Stores in {R} the null matrix. */

void rmxn_ident(int m, int n, double *R);
  /* Stores in {R} an identity matrix, truncated to {m} rows
    and {n}columns. */

void rmxn_scale(int m, int n, double s, double *A, double *R);
  /* Stores in {R} the matrix {s*A}.  The matrix {R} may be the 
    same as {A}. */

void rmxn_mix (int m, int n, double s, double *A, double t, double *B, double *R);
  /* Stores in {R} the matrix {s*A + t*B}. The matrix {R} may be the 
    same as {A}.*/

void rmxn_rel_diff(int m, int n, double *A, double *B, double *R);
  /* Stores in {z = R[i,j]} the relative difference between {x = A[i,j]}
    and {y = B[i,j]}, namely {(x-y)/sqrt(x^2 + y^2)/2}.  
    (See {rel_diff} in {js.h}). The matrix {R} may be the 
    same as {A} or {B}. */

void rmxn_map_row (int m, int n, double *x, double *A, double *r);
  /* Computes {r = x * A}, where {x} is a (row) vector of size {m},
    and {A} is matrix of size {m x n}. The result is a (row) vector of
    size {n}. */

void rmxn_map_col (int m, int n, double *A, double *x, double *R);
  /* Computes {r = A * x}, where {A} is a matrix of size {m x n},
    and {x} is a (column) vector of size {n}. The result is a (column)
    vector of size {m}. */

void rmxn_mul (int m, int p, int n, double *A, double *B, double *R);
  /* Computes the matrix product {R = A * B}, where {A} has size
    {m x p} and {B} has size {p x n}. */

void rmxn_mul_tr (int m, int n, int p, double *A, double *B, double *R);
  /* Computes the matrix product {R = A * B^t}, where {A} has size
    {m x p} and {B} has size {n x p}.  (In other words, sets {R[i,j]}
    to the dot product of row {i} of {A} and row {j} of {B}.)  */

double rmxn_det (int n, double *A);
  /* Returns the determinant of the {n x n} matrix {A} */

double rmxn_cof (int n, double *A, int ix, int jx);
  /* Returns the cofactor of element {[ix,jx]} in the {n x n} matrix {A}. */

void rmxn_adj (int n, double *A, double *R);
  /* Sets {R} to the adjoint of the {n x n} matrix {A}. */

void rmxn_inv (int n, double *A, double *R);
  /* Sets {R} to the inverse of the {n x n} matrix {A}. */

/* MATRIX FACTORIZATION */

void rmxn_cholesky(int n, double *A, double *L);
  /* Factors the positive definite matrix {A}, with {n} rows and
   columns, into {L*L^t} where {L} is lower triangular and {L^t} is
   the transpose of {L}.
   
   (The name is French, thus it should be pronounced "sholeskEE"). */

/* OPERATIONS ON LOWER TRIANGULAR MATRICES */

/* The operations in this section require the matrix {L} to 
   be lower triangular, with nonzero entries in the diagonal. */

void rmxn_LT_inv_map_row(int n, double *y, double *L, double *r);
  /* Computes {r = y * (L^-1)}, where {y} is a (row) vector of size
    {n}, and {L} is a lower triangular matrix of size {n} by {n}. The
    result is a (row) vector of size {n}. In other words, finds the
    solution {r} to the linear system {r*L = y}. The vector {r} may be
    the same as {y}. */
    
void rmxn_LT_inv_map_col(int m, double *L, double *y, double *r);
  /* Computes {r = (L^-1) * y}, where {L} is a lower triangular matrix
    of size {m} by {m}, and {y} is a (column) vector of size {m}. The
    result is a (column) vector of size {m}. In other words, finds the
    solution {r} to the linear system {L*r = y}. The vector {r} may be
    the same as {y}. */

void rmxn_LT_pos_div(int m, int n, double *A, double *L, double *R);
  /* Computes {R = A * (L^-1)}, where {A} is a rectangular matrix of
    size {m} by {n}, and {L} is a lower triangular matrix of size {n}
    by {n}. In other words, finds the solution {R} to the matrix
    equation {R*L = A}. The result {R} has size {m} by {n}, and may be
    the same as {A}. */

void rmxn_LT_pre_div(int m, int n, double *L, double *A, double *R);
  /* Computes {R = (L^-1) * A}, where {L} is a lower triangular matrix
    of size {m} by {m}, and {A} is a rectangular matrix of size
    {m} by {n}. In other words, finds the solution {R} to the matrix
    equation {L*R = A}. The result has size {m} by {n}, and may be
    the same as {A}. */
    
/* MATRIX FORMATTING */

void rmxn_print (FILE *f, int m, int n, double *A);
  /* Prints the {m x n} matrix {A} to file {f}, with default format. */

void rmxn_gen_print 
  ( FILE *f, int m, int n, double *A,
    char *fmt, 
    char *olp, char *osep, char *orp, /* Outer delimiters. */
    char *ilp, char *isep, char *irp  /* Inner delimiters. */
  );
  /* Prints the {m x n} matrix {A} to file {f}, using {fmt} for each
    element. The matrix is bounded by {olp} and {orp}, and rows are
    separated by {osep}. Each row is bounded by {ilp} and {irp}, and
    elements are separated by {isep}. Defaults are provided for any of
    these strings which are NULL. */

/* HEAP ALLOCATION */

double *rmxn_alloc(int m, int n);
  /* Allocates {m*n} {double}s on the heap; bombs out if no mem. */

#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 */