<?php

/**
* @package NLA
*
* NLA is short for NUMAL Linear Algebra.  The concept behind this package is 
* still in flux.  I am just playing around with the higher level NUMAL routines.  
* I ported it because the inv method is something I want to play with and it gets 
* its pivots from the dec method.  Also, the dec and sol approach is a neat unix 
* like toolkit approach that potentially offers and easy to use and learn API.  
*
* @author Paul Meagher 
* @author Hang T. Lau 
* @license PHP License v3.0
* @version 0.1
*
* The dec method was ported from:
*
* Hang T. Lau (2003) A Numerical Library in Java for Scientists & Engineers. 
* CRC Press. [http://www.crcpress.com]
*/

include "../Basic.php";

$basic = new NUMAL_Basic;

/**
* 3.1.1.
*
* A. dec
*
* Decomposes the nxn matrix A in the form LU=PA, where L is lower triangular, U is unit 
* upper triangular, and P is a permutation matrix.
*
* The decomposition only uses partial pivoting [Dec68, Wi63, Wi65].  Since in exceptional
* cases, partial pivoting may yield useless results, even for well-conditioned matrices, 
* the user is advised to use the procedure gsselm. However, if the number of variables
* is small relative to the number of binary digits in the mantissa of the number 
* representation then the procedure dec may also be used.  Refer to the procedure gsselm
* for more details.

* @param $a matrix to be decomposed
* @param $n order of the matrix
* @param $aux[1] whether determinant is positive (1) or negative (-1)
* @param $aux[2] a relative tolerance
* @param $aux[3] number of elimination steps performed
* @param $p the pivot indices; row i and row p[i] are interchanged in the i-th iteration. 
* @uses matmat, mattam, ichrow
*/
function dec(&$a, $n, &$aux, &$p)
{
  global $basic;
      
  $pk=0;
  $r = -1.0;
  for ($i=1; $i<=$n; $i++) {
    $s=sqrt($basic->mattam(1,$n,$i,$i,$a,$a));
    if ($s > $r) $r=$s;
    $v[$i]=1.0/$s;
  }
  $eps=$aux[2]*$r;
  $d=1;
  for ($k=1; $k<=$n; $k++) {
    $r = -1.0;
    $k1=$k-1;
    for ($i=$k; $i<=$n; $i++) {
      $a[$i][$k] -= $basic->matmat(1,$k1,$i,$k,$a,$a);
      $s=abs($a[$i][$k])*$v[$i];
      if ($s > $r) {
        $r=$s;
        $pk=$i;
      }
    }
    $p[$k]=$pk;
    $v[$pk]=$v[$k];
    $s=$a[$pk][$k];
    if (abs($s) < $eps) break;
    if ($s < 0.0) $d = -$d;
    if ($pk != $k) {
      $d = -$d;
      $basic->ichrow(1,$n,$k,$pk,$a);
    }
    for ($i=$k+$1; $i<=$n; $i++) 
      $a[$k][$i]=($a[$k][$i]-$basic->matmat(1,$k1,$k,$i,$a,$a))/$s;
  }
  $aux[1]=$d;
  $aux[3]=$k-1;
}

?>