You are here

public function SingularValueDecomposition::__construct in Loft Data Grids 7.2

Same name and namespace in other branches
  1. 6.2 vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php \SingularValueDecomposition::__construct()

* Construct the singular value decomposition * * Derived from LINPACK code. * *

Parameters

$A Rectangular matrix: * @return Structure to access U, S and V.

File

vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php, line 61

Class

SingularValueDecomposition
@package JAMA

Code

public function __construct($Arg) {

  // Initialize.
  $A = $Arg
    ->getArrayCopy();
  $this->m = $Arg
    ->getRowDimension();
  $this->n = $Arg
    ->getColumnDimension();
  $nu = min($this->m, $this->n);
  $e = array();
  $work = array();
  $wantu = true;
  $wantv = true;
  $nct = min($this->m - 1, $this->n);
  $nrt = max(0, min($this->n - 2, $this->m));

  // Reduce A to bidiagonal form, storing the diagonal elements
  // in s and the super-diagonal elements in e.
  for ($k = 0; $k < max($nct, $nrt); ++$k) {
    if ($k < $nct) {

      // Compute the transformation for the k-th column and
      // place the k-th diagonal in s[$k].
      // Compute 2-norm of k-th column without under/overflow.
      $this->s[$k] = 0;
      for ($i = $k; $i < $this->m; ++$i) {
        $this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
      }
      if ($this->s[$k] != 0.0) {
        if ($A[$k][$k] < 0.0) {
          $this->s[$k] = -$this->s[$k];
        }
        for ($i = $k; $i < $this->m; ++$i) {
          $A[$i][$k] /= $this->s[$k];
        }
        $A[$k][$k] += 1.0;
      }
      $this->s[$k] = -$this->s[$k];
    }
    for ($j = $k + 1; $j < $this->n; ++$j) {
      if ($k < $nct & $this->s[$k] != 0.0) {

        // Apply the transformation.
        $t = 0;
        for ($i = $k; $i < $this->m; ++$i) {
          $t += $A[$i][$k] * $A[$i][$j];
        }
        $t = -$t / $A[$k][$k];
        for ($i = $k; $i < $this->m; ++$i) {
          $A[$i][$j] += $t * $A[$i][$k];
        }

        // Place the k-th row of A into e for the
        // subsequent calculation of the row transformation.
        $e[$j] = $A[$k][$j];
      }
    }
    if ($wantu and $k < $nct) {

      // Place the transformation in U for subsequent back
      // multiplication.
      for ($i = $k; $i < $this->m; ++$i) {
        $this->U[$i][$k] = $A[$i][$k];
      }
    }
    if ($k < $nrt) {

      // Compute the k-th row transformation and place the
      // k-th super-diagonal in e[$k].
      // Compute 2-norm without under/overflow.
      $e[$k] = 0;
      for ($i = $k + 1; $i < $this->n; ++$i) {
        $e[$k] = hypo($e[$k], $e[$i]);
      }
      if ($e[$k] != 0.0) {
        if ($e[$k + 1] < 0.0) {
          $e[$k] = -$e[$k];
        }
        for ($i = $k + 1; $i < $this->n; ++$i) {
          $e[$i] /= $e[$k];
        }
        $e[$k + 1] += 1.0;
      }
      $e[$k] = -$e[$k];
      if ($k + 1 < $this->m and $e[$k] != 0.0) {

        // Apply the transformation.
        for ($i = $k + 1; $i < $this->m; ++$i) {
          $work[$i] = 0.0;
        }
        for ($j = $k + 1; $j < $this->n; ++$j) {
          for ($i = $k + 1; $i < $this->m; ++$i) {
            $work[$i] += $e[$j] * $A[$i][$j];
          }
        }
        for ($j = $k + 1; $j < $this->n; ++$j) {
          $t = -$e[$j] / $e[$k + 1];
          for ($i = $k + 1; $i < $this->m; ++$i) {
            $A[$i][$j] += $t * $work[$i];
          }
        }
      }
      if ($wantv) {

        // Place the transformation in V for subsequent
        // back multiplication.
        for ($i = $k + 1; $i < $this->n; ++$i) {
          $this->V[$i][$k] = $e[$i];
        }
      }
    }
  }

  // Set up the final bidiagonal matrix or order p.
  $p = min($this->n, $this->m + 1);
  if ($nct < $this->n) {
    $this->s[$nct] = $A[$nct][$nct];
  }
  if ($this->m < $p) {
    $this->s[$p - 1] = 0.0;
  }
  if ($nrt + 1 < $p) {
    $e[$nrt] = $A[$nrt][$p - 1];
  }
  $e[$p - 1] = 0.0;

  // If required, generate U.
  if ($wantu) {
    for ($j = $nct; $j < $nu; ++$j) {
      for ($i = 0; $i < $this->m; ++$i) {
        $this->U[$i][$j] = 0.0;
      }
      $this->U[$j][$j] = 1.0;
    }
    for ($k = $nct - 1; $k >= 0; --$k) {
      if ($this->s[$k] != 0.0) {
        for ($j = $k + 1; $j < $nu; ++$j) {
          $t = 0;
          for ($i = $k; $i < $this->m; ++$i) {
            $t += $this->U[$i][$k] * $this->U[$i][$j];
          }
          $t = -$t / $this->U[$k][$k];
          for ($i = $k; $i < $this->m; ++$i) {
            $this->U[$i][$j] += $t * $this->U[$i][$k];
          }
        }
        for ($i = $k; $i < $this->m; ++$i) {
          $this->U[$i][$k] = -$this->U[$i][$k];
        }
        $this->U[$k][$k] = 1.0 + $this->U[$k][$k];
        for ($i = 0; $i < $k - 1; ++$i) {
          $this->U[$i][$k] = 0.0;
        }
      }
      else {
        for ($i = 0; $i < $this->m; ++$i) {
          $this->U[$i][$k] = 0.0;
        }
        $this->U[$k][$k] = 1.0;
      }
    }
  }

  // If required, generate V.
  if ($wantv) {
    for ($k = $this->n - 1; $k >= 0; --$k) {
      if ($k < $nrt and $e[$k] != 0.0) {
        for ($j = $k + 1; $j < $nu; ++$j) {
          $t = 0;
          for ($i = $k + 1; $i < $this->n; ++$i) {
            $t += $this->V[$i][$k] * $this->V[$i][$j];
          }
          $t = -$t / $this->V[$k + 1][$k];
          for ($i = $k + 1; $i < $this->n; ++$i) {
            $this->V[$i][$j] += $t * $this->V[$i][$k];
          }
        }
      }
      for ($i = 0; $i < $this->n; ++$i) {
        $this->V[$i][$k] = 0.0;
      }
      $this->V[$k][$k] = 1.0;
    }
  }

  // Main iteration loop for the singular values.
  $pp = $p - 1;
  $iter = 0;
  $eps = pow(2.0, -52.0);
  while ($p > 0) {

    // Here is where a test for too many iterations would go.
    // This section of the program inspects for negligible
    // elements in the s and e arrays.  On completion the
    // variables kase and k are set as follows:
    // kase = 1  if s(p) and e[k-1] are negligible and k<p
    // kase = 2  if s(k) is negligible and k<p
    // kase = 3  if e[k-1] is negligible, k<p, and
    //           s(k), ..., s(p) are not negligible (qr step).
    // kase = 4  if e(p-1) is negligible (convergence).
    for ($k = $p - 2; $k >= -1; --$k) {
      if ($k == -1) {
        break;
      }
      if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k + 1]))) {
        $e[$k] = 0.0;
        break;
      }
    }
    if ($k == $p - 2) {
      $kase = 4;
    }
    else {
      for ($ks = $p - 1; $ks >= $k; --$ks) {
        if ($ks == $k) {
          break;
        }
        $t = ($ks != $p ? abs($e[$ks]) : 0.0) + ($ks != $k + 1 ? abs($e[$ks - 1]) : 0.0);
        if (abs($this->s[$ks]) <= $eps * $t) {
          $this->s[$ks] = 0.0;
          break;
        }
      }
      if ($ks == $k) {
        $kase = 3;
      }
      else {
        if ($ks == $p - 1) {
          $kase = 1;
        }
        else {
          $kase = 2;
          $k = $ks;
        }
      }
    }
    ++$k;

    // Perform the task indicated by kase.
    switch ($kase) {

      // Deflate negligible s(p).
      case 1:
        $f = $e[$p - 2];
        $e[$p - 2] = 0.0;
        for ($j = $p - 2; $j >= $k; --$j) {
          $t = hypo($this->s[$j], $f);
          $cs = $this->s[$j] / $t;
          $sn = $f / $t;
          $this->s[$j] = $t;
          if ($j != $k) {
            $f = -$sn * $e[$j - 1];
            $e[$j - 1] = $cs * $e[$j - 1];
          }
          if ($wantv) {
            for ($i = 0; $i < $this->n; ++$i) {
              $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p - 1];
              $this->V[$i][$p - 1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p - 1];
              $this->V[$i][$j] = $t;
            }
          }
        }
        break;

      // Split at negligible s(k).
      case 2:
        $f = $e[$k - 1];
        $e[$k - 1] = 0.0;
        for ($j = $k; $j < $p; ++$j) {
          $t = hypo($this->s[$j], $f);
          $cs = $this->s[$j] / $t;
          $sn = $f / $t;
          $this->s[$j] = $t;
          $f = -$sn * $e[$j];
          $e[$j] = $cs * $e[$j];
          if ($wantu) {
            for ($i = 0; $i < $this->m; ++$i) {
              $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k - 1];
              $this->U[$i][$k - 1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k - 1];
              $this->U[$i][$j] = $t;
            }
          }
        }
        break;

      // Perform one qr step.
      case 3:

        // Calculate the shift.
        $scale = max(max(max(max(abs($this->s[$p - 1]), abs($this->s[$p - 2])), abs($e[$p - 2])), abs($this->s[$k])), abs($e[$k]));
        $sp = $this->s[$p - 1] / $scale;
        $spm1 = $this->s[$p - 2] / $scale;
        $epm1 = $e[$p - 2] / $scale;
        $sk = $this->s[$k] / $scale;
        $ek = $e[$k] / $scale;
        $b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
        $c = $sp * $epm1 * ($sp * $epm1);
        $shift = 0.0;
        if ($b != 0.0 || $c != 0.0) {
          $shift = sqrt($b * $b + $c);
          if ($b < 0.0) {
            $shift = -$shift;
          }
          $shift = $c / ($b + $shift);
        }
        $f = ($sk + $sp) * ($sk - $sp) + $shift;
        $g = $sk * $ek;

        // Chase zeros.
        for ($j = $k; $j < $p - 1; ++$j) {
          $t = hypo($f, $g);
          $cs = $f / $t;
          $sn = $g / $t;
          if ($j != $k) {
            $e[$j - 1] = $t;
          }
          $f = $cs * $this->s[$j] + $sn * $e[$j];
          $e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
          $g = $sn * $this->s[$j + 1];
          $this->s[$j + 1] = $cs * $this->s[$j + 1];
          if ($wantv) {
            for ($i = 0; $i < $this->n; ++$i) {
              $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j + 1];
              $this->V[$i][$j + 1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j + 1];
              $this->V[$i][$j] = $t;
            }
          }
          $t = hypo($f, $g);
          $cs = $f / $t;
          $sn = $g / $t;
          $this->s[$j] = $t;
          $f = $cs * $e[$j] + $sn * $this->s[$j + 1];
          $this->s[$j + 1] = -$sn * $e[$j] + $cs * $this->s[$j + 1];
          $g = $sn * $e[$j + 1];
          $e[$j + 1] = $cs * $e[$j + 1];
          if ($wantu && $j < $this->m - 1) {
            for ($i = 0; $i < $this->m; ++$i) {
              $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j + 1];
              $this->U[$i][$j + 1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j + 1];
              $this->U[$i][$j] = $t;
            }
          }
        }
        $e[$p - 2] = $f;
        $iter = $iter + 1;
        break;

      // Convergence.
      case 4:

        // Make the singular values positive.
        if ($this->s[$k] <= 0.0) {
          $this->s[$k] = $this->s[$k] < 0.0 ? -$this->s[$k] : 0.0;
          if ($wantv) {
            for ($i = 0; $i <= $pp; ++$i) {
              $this->V[$i][$k] = -$this->V[$i][$k];
            }
          }
        }

        // Order the singular values.
        while ($k < $pp) {
          if ($this->s[$k] >= $this->s[$k + 1]) {
            break;
          }
          $t = $this->s[$k];
          $this->s[$k] = $this->s[$k + 1];
          $this->s[$k + 1] = $t;
          if ($wantv and $k < $this->n - 1) {
            for ($i = 0; $i < $this->n; ++$i) {
              $t = $this->V[$i][$k + 1];
              $this->V[$i][$k + 1] = $this->V[$i][$k];
              $this->V[$i][$k] = $t;
            }
          }
          if ($wantu and $k < $this->m - 1) {
            for ($i = 0; $i < $this->m; ++$i) {
              $t = $this->U[$i][$k + 1];
              $this->U[$i][$k + 1] = $this->U[$i][$k];
              $this->U[$i][$k] = $t;
            }
          }
          ++$k;
        }
        $iter = 0;
        --$p;
        break;
    }

    // end switch
  }

  // end while
}