Core JavaScript

Generating a Singular Value Decomposition

The JavaScript code example demonstrates how to create the singular value decomposition of a matrix.

SingularValueDecomposition.html

<!DOCTYPE html>
<html>
	<head>
		<title>XoaX.net's JavaScript</title>
    <style>
      th, td {
        text-align: center;
        border-color:black;
        width:25px;
        padding:5px;
      }
      .cMatrixExpression {
      	float:left;
      	margin:5px;
      	border-width:1px;
      	border-radius:5px;
      	border-style:none solid;
      }
      #idA td {
        background-color:white;
      }
    </style>
		<script type="text/javascript" src="SingularValueDecomposition.js"></script>
	</head>
	<body onload="fnInitialize()">
		<button onclick="fnDecomposeMatrix()">Decompose</button>
		<button onclick="fnRandomizeMatrix()">Randomize</button>
		<label for="nRows">Rows:</label><input type="text" id="idRows" name="nRows" maxlength="2" size="1" onchange="fnResizeTheMatrix()" />
		<label for="nCols">Columns:</label><input type="text" id="idCols" name="nCols" maxlength="2" size="1" onchange="fnResizeTheMatrix()" />
		<hr style="clear:both;"/>
		<div id="idSVD">
			<table id="idA" class="cMatrixExpression">
				<tr><td contenteditable="plaintext-only">1</td><td contenteditable="plaintext-only">3</td><td contenteditable="plaintext-only">-2</td></tr>
				<tr><td contenteditable="plaintext-only">2</td><td contenteditable="plaintext-only">7</td><td contenteditable="plaintext-only">5</td></tr>
				<tr><td contenteditable="plaintext-only">-2</td><td contenteditable="plaintext-only">-3</td><td contenteditable="plaintext-only">4</td></tr>
				<tr><td contenteditable="plaintext-only">5</td><td contenteditable="plaintext-only">-3</td><td contenteditable="plaintext-only">-2</td></tr>
			</table>
		</div>
	</body>
</html>

SingularValueDecomposition.js

function fnInitialize() {
	fnFillInRowsAndColumns();
}

function fnClearTheDecomposition() {
	// The table with the matrix entries
	const qTableA = document.getElementById("idA");
	// The table with the matrix entries
	const qDivSvd = document.getElementById("idSVD");
	qDivSvd.innerHTML = "";
	qDivSvd.appendChild(qTableA);
}

function fnFillInRowsAndColumns() {
	const qRowsInput = document.getElementById("idRows");
	const qColsInput = document.getElementById("idCols");
	const qaTableA = document.getElementById("idA");
	// Get the rows of the matrix
	const kiRows = qaTableA.rows.length;
	const kiCols = qaTableA.rows[0].children.length;
	qRowsInput.value = kiRows;
	qColsInput.value = kiCols;
}

function fnResizeTheMatrix() {
	const qRowsInput = document.getElementById("idRows");
	const qColsInput = document.getElementById("idCols");
	let iNewRows = qRowsInput.value;
	let iNewCols = qColsInput.value;
	let qTableA = document.getElementById("idA");
	let qTableBody = qTableA.children[0];
	// Get the rows of the matrix
	let iRows = qTableA.rows.length;
	let iCols = qTableA.rows[0].children.length;
	
	if (iNewRows < iRows) { // Remove any extra rows
		while (iNewRows < qTableA.rows.length) {
			qTableBody.removeChild(qTableBody.children[iNewRows]);
		}
	} else if (iNewRows > iRows) { // Add any extra rows
		while (iNewRows > qTableA.rows.length) {
			var qNewTr = document.createElement("tr");
			qTableBody.appendChild(qNewTr);
			for (var j = 0; j < iCols; ++j) {
				var qNewTd = document.createElement("td");
				qNewTd.innerHTML = 0;
				qNewTr.appendChild(qNewTd);
			}
		}
	}
	if (iNewCols < iCols) { // Remove any extra columns
		for (var i = 0; i < qTableA.rows.length; ++i) {
			let qRow = qTableBody.rows[i];
			while (iNewCols < qRow.children.length) {
				qRow.removeChild(qRow.children[iNewCols]);
			}
		}
	} else if (iNewCols > iCols) { // Add any extra columns
		for (var i = 0; i < qTableA.rows.length; ++i) {
			let qRow = qTableBody.rows[i];
			while (iNewCols > qRow.children.length) {
				var qNewTd = document.createElement("td");
				qNewTd.innerHTML = 0;
				qRow.appendChild(qNewTd);
			}
		}				}
	// Clear the decomposition, if it exists.
	fnClearTheDecomposition();
}

function fnRandomizeMatrix() {
	// Get the table with the matrix entries
	const qaTableA = document.getElementById("idA");
	for (var i = 0; i < qaTableA.rows.length; ++i) {
		for (var j = 0; j < qaTableA.rows[i].children.length; ++j) {
			qaTableA.rows[i].children[j].innerHTML = (Math.random()*-10.0).toFixed(4).toString();
		}
	}
}

function fnDecomposeMatrix() {
	fnClearTheDecomposition();
	// Get the rows of the matrix
	// The table with the matrix entries
	const qaTableA = document.getElementById("idA");
	// Get the rows of the matrix
	const qaTableRowsA = qaTableA.rows;
	// Allocate the row array of A
	let daaA = new Array();
	for (var i = 0; i < qaTableRowsA.length; ++i) {
		// Add a new row
		daaA.push(new Array());
		// Read the entries from the elements and put them into a 2D array
		for (var j = 0; j < qaTableRowsA[i].children.length; ++j) {
			// Convert the matrix entry into from a string to a float and add it the 2D array 
			daaA[i][j] = parseFloat(qaTableRowsA[i].children[j].innerHTML);
		}
	}
	// Send the matrix as a JSON string
	const sJsonA = JSON.stringify(daaA);	
	// The factorization occurs on the server
	var qHttpRequest = new XMLHttpRequest();
	qHttpRequest.onreadystatechange = function() {
		if (this.readyState == 4 && this.status == 200) {
			qaTableA.insertAdjacentHTML("afterend", this.responseText);
		}
	};
	qHttpRequest.open("POST", "SingularValueDecomposition.php", true);
	qHttpRequest.setRequestHeader("Content-type", "application/x-www-form-urlencoded");
	qHttpRequest.send("sJsonA="+sJsonA);
}

SingularValueDecomposition.php

<?php

$sJsonA = $_POST['sJsonA'];

// Decode the JSON into a 2D array
$daaA = json_decode($sJsonA);


$iRows = count($daaA);
$iCols = count($daaA[0]);

$daAtoU = new SplFixedArray($iRows*$iCols);

for ($i = 0; $i < $iRows; ++$i) {
	for ($j = 0; $j < $iCols; ++$j) {
		$daAtoU[$i*$iCols + $j] = $daaA[$i][$j];
	}
}

$daW	 = new SplFixedArray($iCols);
$daVt	 = new SplFixedArray($iCols*$iCols);

SVD($daAtoU, $iRows, $iCols, $daW, $daVt);
ReorderSVD($daAtoU, $iRows, $iCols, $daW, $daVt);

$daaU = array();
$daaW = array();
$daaVt = array();

// Print an equals sign at half height
$sEquals = '<table style="float:left;">';
$iHalfRows = round($iRows/2) - 1;
$iHalfLineHeight = round(imagefontheight(1));
if (2*($iHalfRows + 1) == $iRows) {
	// Insert a blank half line offset
	$sEquals .= '<tr><td style="height:'.$iHalfLineHeight.'px;"></td></tr>';
}
for ($i = 0; $i < $iRows; ++$i) {
	if ($i == $iHalfRows) {
		$sEquals .= '<tr><td>=</td></tr>';
	} else {
		$sEquals .= '<tr><td>&nbsp;</td></tr>';
	}
}
$sEquals .= '</table>';

// Put the entries into a 2-d array for convenience
// U should be m-by-m, but it is only m by n, because the last columns map to zero and can be anything.
for ($i = 0; $i < $iRows; ++$i) {
	array_push($daaU, array());
	for ($j = 0; $j < $iRows; ++$j) {
		// Is this from the zeroed columns?
		$dNext = (($j < $iCols) ? $daAtoU[$i*$iCols + $j] : 0);
		array_push($daaU[$i], $dNext);
	}
}

// W is diagonal and it is m-by-n, but we only have the diagonal entries
for ($i = 0; $i < $iRows; ++$i) {
	array_push($daaW, array());
	for ($j = 0; $j < $iCols; ++$j) {
		// Insert the a diagonal entries from the array.
		$dNext = (($i == $j) ? $daW[$i] : 0);
		array_push($daaW[$i], $dNext);
	}
}

// V transpose is n-by-n
for ($i = 0; $i < $iCols; ++$i) {
	array_push($daaVt, array());
	for ($j = 0; $j < $iCols; ++$j) {
		// Swap the row and column to get V from V transpose
		array_push($daaVt[$i], $daVt[$j*$iCols + $i]);
	}
}

// Print out the matrices U, W, and V as a product.
$sU = CreateMatrixTable($daaU, "idU");
$sW = CreateMatrixTable($daaW, "idW");
$sV = CreateMatrixTable($daaVt, "idVt");
echo $sEquals.$sU.$sW.$sV;
echo '<hr style="clear:both;" />';

function CreateTensorSumMatrix($daAtoU, $iRows, $iCols, $daW, $daVt) {
	$daaProduct = array();
	for ($i = 0; $i < $iRows; ++$i) {
		array_push($daaProduct, array());
		for ($j = 0; $j < $iCols; ++$j) {
			array_push($daaProduct[$i], 0.0);
		}
	}
	$kiCompCount = count($daW);
	for ($iC = 0; $iC < $kiCompCount; ++$iC) {
		for ($i = 0; $i < $iRows; ++$i) {
			for ($j = 0; $j < $iCols; ++$j) {
				// This transposes the entries of V transpose to use V instead.
				$daaProduct[$i][$j] += $daAtoU[$i*$iCols + $iC]*$daW[$iC]*$daVt[$j*$iCols + $iC];
			}
		}
	}
	return $daaProduct;
}

function CreateMatrixTable1D($daM, $iRows, $iCols) {
	$sOutput = '<table class="cMatrixExpression">';
	for ($i = 0; $i < $iRows; ++$i) {
		$sOutput .= '<tr>';
		for ($j = 0; $j < $iCols; ++$j) {
			$sOutput .= '<td>'.$daM[$i*$iCols + $j].'</td>';
		}
		$sOutput .= '</tr>';
	}
	$sOutput .= '</table>';
	return $sOutput;
}

function CreateMatrixTable($daaM, $sID) {
	$iRows = count($daaM);
	$iCols = count($daaM[0]);
	$sOutput = '<table id="{$sID}" class="cMatrixExpression">';
	for ($i = 0; $i < $iRows; ++$i) {
		$sOutput .= '<tr>';
		for ($j = 0; $j < $iCols; ++$j) {
			$sOutput .= '<td>'.round($daaM[$i][$j], 4).'</td>';
		}
		$sOutput .= '</tr>';
	}
	$sOutput .= '</table>';
	return $sOutput;
}

function MultiplyMatrices($daaA, $daaB) {
	$iRowsA = count($daaA);
	$iColsA = count($daaA[0]);
	$iRowsB = count($daaB);
	$iColsB = count($daaB[0]);
	$daaProduct = array();
	for ($i = 0; $i < $iRowsA; ++$i) {
		// Add a row to the product
		array_push($daaProduct, array());
		for ($j = 0; $j < $iColsB; ++$j) {
			$dNewEntry = 0.0;
			for ($k = 0; $k < $iRowsB; ++$k) {
				$dNewEntry += $daaA[$i][$k]*$daaB[$k][$j];
			}
			array_push($daaProduct[$i], $dNewEntry);
		}
	}
	return $daaProduct;
}
function MultiplyWithTranspose($daaA) {
	$iRowsA = count($daaA);
	$iColsA = count($daaA[0]);
	$daaProduct = array();
	for ($i = 0; $i < $iRowsA; ++$i) {
		// Add a row to the product
		array_push($daaProduct, array());
		for ($j = 0; $j < $iColsA; ++$j) {
			$dNewEntry = 0;
			for ($k = 0; $k < $iColsA; ++$k) {
				$dNewEntry += $daaA[$i][$k]*$daaA[$j][$k];
			}
			array_push($daaProduct[$i], $dNewEntry);
		}
	}
	return $daaProduct;
}


function Hypotenuse($dA, $dB) {
  $dAbsA = abs($dA);
  $dAbsB = abs($dB);
  if ($dAbsA > $dAbsB) {
    $dQ = ($dAbsB/$dAbsA);
    return $dAbsA*sqrt(1 + $dQ*$dQ);
  } else if ($dAbsB == 0.0) {
    return 0.0;
  } else {
    $dQ = ($dAbsA/$dAbsB);
    return $dAbsB*sqrt(1 + $dQ*$dQ);
  }
}

function Sign($dA, $dB) {
  return ($dB >= 0.0
    ? ($dA >= 0.0 ? $dA : -$dA)
    : ($dA >= 0.0 ? -$dA : $dA));
}

// A is the m rows by n columns matrix that is to be decomposed. The matrix A = UxWx(TransposeOfV)
// U is returned in the matrix A. A and Vt are 1D arrays that are large enough to hold matrices.
function SVD($daAtoU, $iM, $iN, $daW, $daVt) {
  // Allocate an array for calculations
  $daRV1 = new SplFixedArray($iN);
  $dNorm = 0.0;
  // These three are probably not needed outside the loop
  $dScale = 0.0;
  $dG = 0.0;
  $dS = 0.0;
  $iL = 0;
  $iNM = 0;
  $dZ = 0.0;
  // Householder reflections to reduce the matrix to a bidiagonal matrix.
  for ($i = 0; $i < $iN; ++$i) {
    // It looks like $iL is always $i + 1
    $iL         = $i + 1;
    $daRV1[$i]  = ($dG*$dScale);
    $dScale     = 0.0;
    $dG         = 0.0;
    $dS         = 0.0;
    if ($i < $iM) {
      for ($k = $i; $k < $iM; ++$k) {
        $dScale += abs($daAtoU[$k*$iN + $i]);
      }
      if ($dScale != 0.0) {
        for ($k = $i; $k < $iM; ++$k) {
          $daAtoU[$k*$iN + $i] /= $dScale;
          $dS += ($daAtoU[$k*$iN + $i]*$daAtoU[$k*$iN + $i]);
        }
        $dF = $daAtoU[$i*$iN + $i];
        // This seems to always be the opposite value. Check it! Perhaps, just $dG = -$dF
        $dG = -Sign(sqrt($dS), $dF);
        $dH = ($dF*$dG - $dS);
        $daAtoU[$i*$iN + $i] = ($dF - $dG);
        for ($j = $iL; $j < $iN; ++$j) {
          $dS = 0.0;
          for ($k = $i; $k < $iM; ++$k) {
            $dS += ($daAtoU[$k*$iN + $i]*$daAtoU[$k*$iN + $j]);
          }
          $dF = ($dS/$dH);
          for ($k = $i; $k < $iM; ++$k) {
            $daAtoU[$k*$iN + $j] += ($dF*$daAtoU[$k*$iN + $i]);
          }
        }
        for ($k = $i; $k < $iM; ++$k) {
          $daAtoU[$k*$iN + $i] *= $dScale;
        }
      }
    } // if ($i < $iM)
    $daW[$i] = ($dScale*$dG);
    $dScale  = 0.0;
    $dG      = 0.0;
    $dS      = 0.0;
    if (($i < $iM) && ($i != ($iN - 1))) {
      for ($k = $iL; $k < $iN; ++$k) {
        $dScale += abs($daAtoU[$i*$iN + $k]);
      }
      if ($dScale != 0.0) {
        for ($k = $iL; $k < $iN; ++$k) {
          $daAtoU[$i*$iN + $k] /= $dScale;
          $dS += ($daAtoU[$i*$iN + $k]*$daAtoU[$i*$iN + $k]);
        }
        $dF = $daAtoU[$i*$iN + $iL];
        // This seems to always be the opposite value. Check it! Perhaps, just $dG = -$dF
        $dG = -Sign(sqrt($dS), $dF);
        $dH = ($dF*$dG - $dS);
        $daAtoU[$i*$iN + $iL] = ($dF - $dG);
        for ($k = $iL; $k < $iN; ++$k) {
          $daRV1[$k] = ($daAtoU[$i*$iN + $k]/$dH);
        }
        for ($j = $iL; $j < $iM; ++$j) {
          $dS = 0.0;
          for ($k = $iL; $k < $iN; ++$k) {
            $dS += ($daAtoU[$j*$iN + $k]*$daAtoU[$i*$iN + $k]);
          }
          for ($k = $iL; $k < $iN; ++$k) {
            $daAtoU[$j*$iN + $k] += ($dS*$daRV1[$k]);
          }
        }
        for ($k = $iL; $k < $iN; ++$k) {
          $daAtoU[$i*$iN + $k] *= $dScale;
        }
      }
    } //if (($i < $iM) && ($i != ($iN - 1))) {
    $dNorm = max($dNorm, abs($daW[$i]) + abs($daRV1[$i]));
  }
  // Accumulation of right-hand transformations
  for ($i = $iN - 1; $i >= 0; --$i) {
    if ($i < ($iN - 1)) {
      if ($dG != 0.0) {
        for ($j = $iL; $j < $iN; ++$j) {
          $daVt[$j*$iN + $i] = (($daAtoU[$i*$iN + $j]/$daAtoU[$i*$iN + $iL])/$dG);
        }
        for ($j = $iL; $j < $iN; ++$j) {
          $dS = 0.0;
          for ($k = $iL; $k < $iN; ++$k) {
            $dS += ($daAtoU[$i*$iN + $k]*$daVt[$k*$iN + $j]);
          }
          for ($k = $iL; $k < $iN; ++$k) {
            $daVt[$k*$iN + $j] += ($dS*$daVt[$k*$iN + $i]);
          }
        }
      } // if ($dG != 0.0) {
      for ($j = $iL; $j < $iN; ++$j) {
        $daVt[$i*$iN + $j] = 0.0;
        $daVt[$j*$iN + $i] = 0.0;
      }
    }
    $daVt[$i*$iN + $i] = 1.0;
    $dG = $daRV1[$i];
    $iL = $i;
  }
  // Accumulation of left-hand transformations
  for ($i = IntVal(Min($iN-1, $iM-1)); $i >= 0; --$i) {
    $iL = $i + 1;
    $dG = $daW[$i];
    for ($j = $iL; $j < $iN; ++$j) {
      $daAtoU[$i*$iN + $j] = 0.0;
    }
    if ($dG != 0.0) {
      $dG = (1.0/$dG);
      for ($j = $iL; $j < $iN; ++$j) {
        $dS = 0.0;
        for ($k = $iL; $k < $iM; ++$k) {
          $dS += ($daAtoU[$k*$iN + $i]*$daAtoU[$k*$iN + $j]);
        }
        $dF = ($dS/$daAtoU[$i*$iN + $i])*$dG;
        for ($k = $i; $k < $iM; ++$k) {
          $daAtoU[$k*$iN + $j] += ($dF*$daAtoU[$k*$iN + $i]);
        }
      }
      for ($j = $i; $j < $iM; ++$j) {
        $daAtoU[$j*$iN + $i] *= $dG;
      }
    } else {
      for ($j = $i; $j < $iM; ++$j) {
        $daAtoU[$j*$iN + $i] = 0.0;
      }
    }
    $daAtoU[$i*$iN + $i] += 1.0;
  }
  // Diagonalization of the bidiagonal form via rotations
  for ($k = $iN - 1; $k >= 0; --$k) {
    for ($iIter = 0; $iIter < 30; ++$iIter) {
      $bFlag = true;
      // Test for splitting. Note: $daRV1[$i] is always zero
      // This is the while loop equivalent of the for loop in the code.
      $iL = $k;
      $bLoop = ($iL >= 0);
      // Test for splitting. $daRV1[0] is always zero
      while ($bLoop) {
        $iNM = $iL - 1;
        if ((abs($daRV1[$iL]) + $dNorm) == $dNorm) {
          $bFlag = FALSE;
          $bLoop = FALSE;
        } else if ((abs($daW[$iNM]) + $dNorm) == $dNorm) {
          $bLoop = FALSE;
        } else {
          --$iL;
          if ($iL < 0) {
            $bLoop = FALSE;
          }
        }
      }
      if ($bFlag) {
        // Cancellation of $daRV1[$iL], if L > 0
        $dC = 0.0;
        $dS = 1.0;
        for ($i = $iL; $i <= $k; ++$i) {
          $dF = ($dS*$daRV1[$i]);
          $daRV1[$i] = ($dC*$daRV1[$i]);
          if ((abs($dF)+$dNorm) == $dNorm) {
            break;
          }
          $dG = $daW[$i];
          $dH = Hypotenuse($dF, $dG);
          $daW[$i] = $dH;
          $dH = (1.0/$dH);
          $dC = ($dG*$dH);
          $dS = (-$dF*$dH);
          // Perform rotations
          for ($j = 0; $j < $iM; ++$j) {
            $dY = $daAtoU[$j*$iN + $iNM];
            $dZ = $daAtoU[$j*$iN + $i];
            $daAtoU[$j*$iN + $iNM] = (($dY*$dC) + ($dZ*$dS));
            $daAtoU[$j*$iN + $i] = (($dZ*$dC) - ($dY*$dS));
          }
        }
      } // if ($bFlag)
      $dZ = $daW[$k];
      // Convergence
      if ($iL == $k) {
        // Singular value is made nonnegative
        if ($dZ < 0.0) {
          $daW[$k] = -$dZ;
          for ($j = 0; $j < $iN; ++$j) {
            $daVt[$j*$iN + $k] = -$daVt[$j*$iN + $k];
          }
        }
        // Break out of the outer loop
        break;
      }
      if ($iIter == 29) { // No convergence. What should we do??
      }
      $dX = $daW[$iL];
      $iNM = $k - 1;
      $dY = $daW[$iNM];
      $dG = $daRV1[$iNM];
      $dH = $daRV1[$k];
      $dF = ((($dY-$dZ)*($dY+$dZ)+($dG-$dH)*($dG+$dH))/(2.0*$dH*$dY));
      $dG = Hypotenuse($dF, 1.0);
      $dF = ( ($dX-$dZ)*($dX+$dZ)+$dH*( ($dY/($dF+Sign($dG,$dF))) -$dH) )/$dX;
      $dC = 1.0;
      $dS = 1.0;
      // Next QR Transformation
      for ($j = $iL; $j <= $iNM; ++$j) {
        $i = $j + 1;
        $dG = $daRV1[$i];
        $dY = $daW[$i];
        $dH = $dS*$dG;
        $dG = $dC*$dG;
        $dZ = Hypotenuse($dF, $dH);
        $daRV1[$j] = $dZ;
        $dC = $dF/$dZ;
        $dS = $dH/$dZ;
        $dF = ($dX*$dC + $dG*$dS);
        $dG = ($dG*$dC - $dX*$dS);
        $dH = $dY*$dS;
        $dY *= $dC;
        for ($jj = 0; $jj < $iN; ++$jj) {
          $dX = $daVt[$jj*$iN + $j];
          $dZ = $daVt[$jj*$iN + $i];
          $daVt[$jj*$iN + $j] = $dX*$dC + $dZ*$dS;
          $daVt[$jj*$iN + $i] = $dZ*$dC - $dX*$dS;
        }
        $dZ = Hypotenuse($dF, $dH);
        $daW[$j] = $dZ;
        if ($dZ != 0.0) {
          $dZ = (1.0/$dZ);
          $dC = $dF*$dZ;
          $dS = $dH*$dZ;
        }
        $dF = ($dG*$dC + $dY*$dS);
        $dX = ($dY*$dC - $dG*$dS);
        for ($jj = 0; $jj < $iM; ++$jj) {
          $dY = $daAtoU[$jj*$iN + $j];
          $dZ = $daAtoU[$jj*$iN + $i];
          $daAtoU[$jj*$iN + $j] = ($dY*$dC + $dZ*$dS);
          $daAtoU[$jj*$iN + $i] = ($dZ*$dC - $dY*$dS);
        }
      }
      $daRV1[$iL] = 0.0;
      $daRV1[$k] = $dF;
      $daW[$k] = $dX;
    } //for ($iIter = 0; $iIter < 30; ++$iIter) {
  }
} // End of function
//*/

function ReorderSVD($daAtoU, $iM, $iN, $daW, $daVt) {
	$iIncr = 1;
	$daSU = new SplFixedArray($iM);
	$daSV = new SplFixedArray($iN);
	// Find the first increment larger than N
	do {
		$iIncr *= 3;
		++$iIncr;	
	} while ($iIncr <= $iN);
	do {
		$iIncr = intdiv($iIncr, 3);
		for ($i = $iIncr; $i < $iN; ++$i) {
			$dSW = $daW[$i];
			for ($k = 0; $k < $iM; ++$k) {
				$daSU[$k] = $daAtoU[$k*$iN + $i];
			}
			for ($k = 0; $k < $iN; ++$k) {
				$daSV[$k] = $daVt[$k*$iN + $i];
			}
			$j = $i;
			while ($daW[$j-$iIncr] < $dSW) {
				$daW[$j] = $daW[$j-$iIncr];
				for ($k = 0; $k < $iM; ++$k) {
					$daAtoU[$k*$iN + $j] = $daAtoU[$k*$iN + $j-$iIncr];
				}
				for ($k = 0; $k < $iN; ++$k) {
					$daVt[$k*$iN + $j] = $daVt[$k*$iN + $j-$iIncr];
				}
				$j -= $iIncr;
				if ($j < $iIncr) {
					break;
				}
			}
			$daW[$j] = $dSW;
			for ($k = 0; $k < $iM; ++$k) {
				$daAtoU[$k*$iN + $j] = $daSU[$k];
			}
			for ($k = 0; $k < $iN; ++$k) {
				$daVt[$k*$iN + $j] = $daSV[$k];
			}
		}
	} while ($iIncr > 1);
	for ($k = 0; $k < $iN; ++$k) {
		$s = 0;
		for ($i = 0; $i < $iM; ++$i) {
			if ($daAtoU[$i*$iN + $k] < 0.0) {
				++$s;
			}
		}
		for ($j = 0; $j < $iN; ++$j) {
			if ($daVt[$j*$iN + $k] < 0.0) {
				++$s;
			}
		}
		if ($s > ($iM + $iN)/2.0) {
			for	($i	=	0; $i	<	$iM; ++$i) {
				$daAtoU[$i*$iN + $k] = -$daAtoU[$i*$iN + $k];
			}
			for	($j	=	0; $j	<	$iN; ++$j) {
				$daVt[$j*$iN + $k] = -$daVt[$j*$iN + $k];
			}
		}
	}
}

?>
 
 

Output

 
 

© 2007–2025 XoaX.net LLC. All rights reserved.